Low-order modelling of three-dimensional surface waves in liquid film flow on a rotating disk

Abstract

Using low-dimensional numerical simulations, we investigate the characteristics of complex and three-dimensional surface waves in a liquid film flowing over a rotating disk, focusing on large flow rates from a nozzle. Existing integral boundary layer (IBL) models, which are based on spatially averaged variables along the direction normal to the disk surface, have primarily focused on the formation of axisymmetric waves under relatively small flow rates. In this study, an extended IBL model that accounts for both laminar and turbulent regimes is developed by considering the non-uniformity of the local flow rate in the spreading film flow and incorporating closure models dependent on the local Reynolds number. Our numerical results successfully capture the generation of concentric waves by an impinging circular liquid jet and their transition into three-dimensional solitary waves. These findings are in good agreement with visualization images and time-series data of free-surface fluctuations from a displacement sensor. The backscattering of small-scale three-dimensional turbulence into large-scale horizontal turbulence inside the film plays a critical role in determining the transition of wave modes and the nonlinear dynamics of the waves in the turbulent regime. Furthermore, the behaviour of three-dimensional waves in the downstream region, including frequent wave coalescence in the transition region and the breakup of small-scale solitons, is distinct from that of gravity-driven falling film flows. The amplitude of the three-dimensional waves is inversely related to the generalized Reynolds number defined for rotating films.

1. Introduction

Viscous liquid films spreading over a horizontally rotating disk are widely encountered in engineering applications such as spin coating, chemical reactions, and heat and mass transfer enhancement because they rapidly reach a fully wetted state. However, in typical scenarios, the flow of these films is accompanied by the formation of large-amplitude surface waves near the inlet region. The surface waves undergo a series of transitions as they propagate outwards along the film. The modelling and control of these waves are crucial, significantly affecting the performance and effectiveness of applications based on liquid film flow over a rotating disk (Aoune & Ramshaw 1999).

Low-order modelling of surface waves in film flow typically relies on integral boundary layer (IBL) models. These models exclude the vertical components of the velocity and pressure fields from the full Navier–Stokes equation, introduce the assumption of a polynomial velocity profile and describe the dynamics of a liquid film as a function of thickness and depth-averaged velocity. Integral boundary layer models are often used because they require lower computational cost than full three-dimensional simulations. Since Shkadov (1968) derived pioneering two-dimensional equations, the application of the integral formulation has been popular in modelling gravity-driven liquid film flows. Ruyer-Quil & Manneville (2000) presented an improved weighted residual integral boundary layer (WRIBL) model that implements the Galerkin method with polynomial basis functions in assuming the local velocity profile. The WRIBL model is a valuable tool for modelling the dynamics of solitary waves in gravity-driven films (Chang & Demekhin 2002; Craster & Matar 2009; Kalliadasis et al. 2012; Ruyer-Quil et al. 2014). Integral boundary layer equations for resolving three-dimensional waves have also been developed. Since Demekhin & Shkadov (1985) suggested the IBL approach for both the spanwise and streamwise directions in computing three-dimensional surface waves, the dynamics of three-dimensional solitary waves have been studied within the IBL framework (Scheid, Ruyer-Quil & Manneville 2006; Demekhin et al. 2007). Recently, the IBL models have been applied for the investigation of three-dimensional waves in more practical conditions such as film flows on moving (Ivanova et al. 2023) or cylindrical (Ding & Wong 2017) surfaces and film flows with surfactant (Batchvarov et al. 2021).

As for the modelling of surface waves in film flow over a rotating disk, integral formulations based on two-dimensional equations have been used to obtain steady axisymmetric solutions (Miyasaka 1974a; Sisoev, Tal'drik & Shkadov 1986) and explore the characteristics of linear stability (Charwat, Kelly & Gazley 1972; Sisoev & Shkadov 1987, 1990). Furthermore, Sisoev, Matar & Lawrence (2003) introduced a set of two-dimensional IBL equations under axisymmetric assumptions. Numerical solutions of this axisymmetric IBL model revealed the emergence of a group of surface waves known as axisymmetric waves, which arise due to convective instability in the downstream (outward) region. Subsequent investigations have focused on the spatiotemporal evolution of these axisymmetric waves and their impacts on heat and mass transfer (Matar, Sisoev & Lawrence 2004; Sisoev, Matar & Lawrence 2005; Kim & Kim 2009; Prieling & Steiner 2013). The axisymmetric waves predominantly occur under relatively small flow rates and high viscosities, where the dynamics of surface waves are strongly influenced by centrifugal and Coriolis forces. However, the structures of surface waves are more intricate and typically three dimensional in larger inlet-flow-rate conditions. That is, existing IBL models have only been successful in capturing axisymmetric wave patterns.

Dynamics of surface waves under large inlet-flow-rate conditions are much different from those of axisymmetric waves. The spatial distribution of the local flow rate is non-uniform due to the radial spreading of the liquid from the inlet, and so various wave regimes emerge locally along the flow path. These are clearly different from the aforementioned axisymmetric waves. The regimes observed sequentially from upstream to downstream in the radial direction include input, first laminar-wave, turbulent, and second laminar-wave regimes (Butuzov & Pukhovoi 1976). The wave regimes observed under large-flow-rate conditions can be simplified as large-scale laminar waves that form upstream and subsequently undergo breakup, leading to the generation of three-dimensional waves. The breakup of concentric laminar waves, known as wave turbulization, gives rise to a chain of three-dimensional solitary pulses, often referred to as $\varLambda$ solitons (Charwat et al. 1972; Miyasaka 1974b; Butuzov & Pukhovoi 1976; Thomas, Faghri & Hankey 1991; Leneweit, Roesner & Koehler 1999; Li et al. 2019). The primary objective of the present work is to establish a novel IBL model capable of capturing the formation of the $\varLambda$-soliton regime, thereby providing valuable physical insights into the dynamics and characteristics of these three-dimensional waves.

In this study the classical IBL equations are extended to encompass the modelling of wave turbulization in the film flow spreading over a rotating disk. It is widely acknowledged that the behaviours of $\varLambda$ solitons in gravity-driven falling films are governed by the Reynolds number, which is defined as the ratio of the volumetric flow rate per unit transverse width to the kinematic viscosity (Adomeit & Renz 2000; Demekhin et al. 2007, 2010). Given that the flow rate in a spreading film is not uniformly distributed (unlike in falling films), we employ the local Reynolds number, which is proportional to the non-uniform local flow rate. In the region with low local Reynolds number, the flow is considered laminar, and the assumptions of existing IBL models are valid. However, in the region near the inlet with high local Reynolds number, the flow can be approximated as turbulent shallow water with strong surface tension. A similar approach was suggested by Mendez et al. (2021) for high-Reynolds-number film flow in jet wiping processes. In line with this approach, the present work suggests a shallow-water model specifically tailored for turbulent films spreading over a rotating disk. Furthermore, our modelling considers, for the first time in the modelling of film flows, the effect of internal turbulent structures, referred to as sub-depth scale turbulent structures. A depth-averaged model proposed for turbulent shallow water (Hinterberger, Fröhlich & Rodi 2007) is modified and integrated into the governing equations. The IBL formulation proposed in this paper is a comprehensive model that covers both laminar and turbulent regimes, with the aim of describing the transition from two- to three-dimensional surface waves. By analysing the wave structures obtained from low-order numerical simulations, the distinct features of the three-dimensional waves in rotating films are elucidated.

In § 2, our problem and experimental set-up are described, followed by an explanation of the proposed modelling strategy. A detailed derivation of the IBL equations is provided in § 3. The numerical methods adopted to solve the IBL equations are presented in § 4. In § 5, numerical results are validated against experimentally obtained visualization images and film thickness data, and the effects of the backscatter model on the dynamics of the film flow are examined. The interaction modes between three-dimensional waves and variations in their scales under different flow conditions are also discussed. Finally, concluding remarks are presented in § 6.

2. Problem description and modelling strategy

2.1. Problem description

The flow configuration considered throughout the present study is illustrated in figure 1(a). Liquid is supplied by a circular jet discharged vertically from a nozzle. The circular liquid jet is chosen as the inflow configuration because it is more extensively used than collars and slot jets. Hereafter, dimensional parameters are denoted with a hat symbol ($\hat {\cdot }$). The liquid jet from the nozzle, with inner diameter $\hat {d}$ and volume flow rate $\hat {Q}$, impinges on the centre of a disk rotating with angular velocity $\hat {\varOmega }$. The distance from the nozzle exit to the surface of the disk is $\hat {H}$. The origin of the coordinate system is fixed at the centre of the rotating disk, and the coordinate system rotates together with the disk. Here $\hat {u}_{r}$, $\hat {u}_{\theta }$ and $\hat {u}_{z}$ represent the radial, angular and vertical components of velocity, respectively. The ranges of the input parameters are presented in table 1. The selected ranges of the nozzle flow rate $\hat {Q}$ and disk angular velocity $\hat {\varOmega }$ encompass conditions that induce the generation of three-dimensional waves for the working fluid. For smaller values of $\hat {Q}$ and $\hat {\varOmega }$, other regimes of waves, including axisymmetric waves and gravity waves, may occur downstream. Exact ranges of $\hat {Q}$ and $\hat {\varOmega }$ for each wave regime are presented in Appendix A. When a circular liquid jet impinges on a flat surface, a thin viscous boundary layer is formed in the vicinity of the stagnation zone (Wang & Khayat 2018). The thickness of the boundary layer increases radially until it merges with the liquid–gas interface at $\hat {r} = \hat {r}_{0}$ (figure 1a). Upstream of the merger point, the long-wave approximation required for integral modelling is no longer valid. Therefore, modelling the region of $\hat {r}<\hat {r}_{0}$ relies on the theory of impinging liquid jets.

Figure 1. (a) Schematic and (b) visualization image for flow configuration and surface wave regimes. In panel (b) the angular velocity of the disk is $\hat {\varOmega } = 52.4\,{\rm rad}\,{\rm s}^{-1}$, and the volume flow rate of the liquid jet at the nozzle exit is $\hat {Q} = 12.5\,{\rm mL}\,{\rm s}^{-1}$.

Table 1. Input parameters (liquid: water at $20\,^{\circ }{\rm C}$).

2.2. Experimental set-up

To validate the IBL model, visualization and sensor experiments are conducted. The experimental apparatus (figure 2) provided by SEMES Co., Ltd.consists of a rotating-disk module and an impinging-jet module, along with their respective control units. In the rotating-disk module a silicon disk of diameter 0.3 m is mounted on a vacuum chuck. The centre of the disk is securely fixed to the vacuum chuck by imposing a negative pressure of $-95\,{\rm Pa}$ through a vacuum pump. A servo motor rotates the chuck and disk at the designated angular velocity $\hat {\varOmega }$. Deionized water at a temperature of $20\,^{\circ }{\rm C}$ is discharged from the nozzle. The control units for the motor motion, nozzle position and mass flow rate are all integrated into the apparatus.

Figure 2. Experimental set-up.

The film flow is visualized using a high-speed camera (FASTCAM MINI-UX 50, Photron, Inc.) at a sampling rate of 2000 frames per second. The camera is positioned at a vertical distance of 50 cm above the disk, and an LED light source is used to illuminate the disk surface. To measure the displacement of the water–air interface, a confocal displacement sensor (CL-P015N, Keyence Co., Ltd.) with a resolution of $0.25\,\mathrm {\mu }{\rm m}$ is installed above the rotating disk. The vertical displacement of the interface is measured at multiple points with a sampling frequency of $10^{4}\,{\rm Hz}$. The measurement points are located at $\hat {r} = 43$ and 109 mm from the disk centre. For detailed information on the measurement method and data processing of the confocal chromatic sensors in film flows, see Hu et al. (2021) and Ubara, Sugimoto & Asano (2022). The time series of interface displacement include the mechanical vibration of the motor–disk system, which has a much slower time scale than the dynamics of surface waves. To remove the signal corresponding to mechanical vibration, a baseline filtering method based on sparsity (Ning, Selesnick & Duval 2014) is employed. A cutoff frequency of 10 times the rate of disk rotation ($10\hat {\varOmega }/2{\rm \pi}$) is applied. As a result, the data obtained from the sensor experiment capture the vertical fluctuations of film thickness.

2.3. Modelling strategy

Figure 1(b) presents a visualization of the typical transitional surface waves formed on a rotating film flow. In the upstream region, concentric waves are generated from the impingement zone of the liquid jet and propagate downstream (outwards). These two-dimensional waves are accompanied by azimuthal instabilities, which become more pronounced as $\hat {\varOmega }$ increases. Once the concentric waves reach a critical radius, three-dimensional solitary waves begin to develop. In the downstream region, coherent wave structures known as $\varLambda$ solitons emerge, characterized by their resemblance to the Greek letter lambda (Demekhin et al. 2007). As the film spreads outwards, these solitons undergo horizontal expansion while their peak thickness decreases. Additionally, smaller-scale solitons are formed between the larger ones.

Aforementioned observations are consistent with the experimental study of Butuzov & Pukhovoi (1976), where the transition from the concentric-wave regime to the $\varLambda$-soliton regime is referred to as wave turbulization. The criteria for this transition have been described in terms of the local Reynolds numbers, which are defined based on the radial flow rate ($Re_{L,r}=\hat {Q}/2{\rm \pi} \hat {r}\hat {\nu }$) and tangential velocity ($Re_{L,t}=\hat {\varOmega }\hat {r}^{2}/\hat {\nu }$) at a given location. In this study we consider the local flow rate per unit width, denoted as $\hat {q}=\rvert \int _{0}^{\hat {h}}(\hat {u}_{r}, \hat {u}_{\theta })\,{\rm d}{z} \rvert$, as a pivotal variable governing the wave behaviours. The time-averaged horizontal flow rate, $\langle \hat {q} \rangle$, is the root mean square of the radial and tangential flow rates, which are roughly proportional to $Re_{L,r}$ and $Re_{L,t}$, respectively (Kim & Kim 2009). Since the local flow rate $\hat {q}$ is determined by disk angular velocity $\hat {\varOmega }$ as well as nozzle flow rate $\hat {Q}$, the choice of $\hat {q}$ as a key variable allows us to explore the influence of both the flow conditions on the dynamics of three-dimensional waves, as presented in § 5.3.

According to Demekhin et al. (2007), in gravity-driven falling film flows, the transition from two-dimensional waves to three-dimensional waves and roll waves occurs as the global Reynolds number ($Re_{G}=\hat {Q}/\hat {\nu }\hat {L}$) increases. The global Reynolds number is calculated based on the total inlet discharge rate divided by the transverse width $\hat {L}$ of the domain. Flows with $Re_{G}$ above a certain threshold are considered to be turbulent or within the transition regime, and are characterized by non-polynomial velocity profiles along the direction normal to the disk surface and vortical structures in regions of high flow rate. Regarding the threshold value for this transition, it is generally accepted that laminar film flow occurs for Reynolds numbers below $Re_{G} \approx 75$ (Ishigai et al. 1972; Karimi & Kawaji 1999; Demekhin et al. 2007). At higher flow rates, occasional turbulent spots are observed when wave solitons interact. A fully turbulent film flow is typically considered for Reynolds numbers above $Re_{G} \approx 200$ (Adomeit & Renz 2000).

The investigation of a rotating film flow in this work accounts for the transition of wave regimes by considering the local flow rate $\hat {q}$, which significantly decreases as the film spreads outwards; see figure 3. Hence, we adopt the local Reynolds number $Re_{L}$ proportional to $\hat {q}$ to determine whether a region is subjected to laminar or turbulent conditions: $Re_{L} = \hat {q}/\hat {\nu }$. Specifically, $Re_{L}$ in proximity to the jet impingement point exceeds the criterion of transition from the laminar to turbulent regime given in previous studies on falling films, and so the region near the jet impingement is turbulent. By contrast, the downstream region with lower local flow rates (and local Reynolds numbers) is predominantly laminar. Consequently, to develop accurate IBL equations for a rotating film with three-dimensional solitary waves, it is essential to incorporate both laminar and turbulent models, taking into account the local Reynolds number.

Figure 3. Time-averaged local Reynolds number $\langle Re_{L} \rangle$ with respect to radial position $\hat {r}$.

The concept of the turbulent regime in the film flow discussed in this study should not be confused with the wave turbulization proposed by Butuzov & Pukhovoi (1976). The term interfacial turbulence is commonly employed to describe the complex dynamics of solitary waves. As the system of solitary waves exhibits robust coherent structures, continuously interacting with each other as quasiparticles, it is often referred to as weak and dissipative turbulence (Denner et al. 2018). Therefore, wave turbulization indicates the transformation of concentric waves into a group of three-dimensional waves exhibiting spatiotemporal chaos. The system of irregular solitary waves is typically associated with low-flow-rate conditions. By contrast, the turbulence discussed in this study emerges when the local flow rate is beyond the range in which solitary waves are observed. In this context, turbulence refers to a hydrodynamic state characterized by the presence of internal vortex structures (Hwang & Chang 1987; Adomeit & Renz 2000) and the dominance of roll waves, along with a power-law velocity distribution in the vertical direction (Brauner 1989). This study emphasizes the presence of turbulent film flow in the upstream region where the local flow rate is large. The turbulent structures in the upstream region are correlated with the generation of chaotic solitary waves, known as wave turbulization, in the downstream region.

Before proceeding with integral modelling, it is crucial to understand the source of concentric waves in the upstream region. Although the exact mechanism remains unclear, it is widely accepted that various interfacial flow phenomena arising from different nozzle configurations affect the formation of concentric waves (Charwat et al. 1972). In this study the concentric waves in the upstream region are caused by shear-induced disturbances on the surface of the circular liquid jet prior to impingement. The inception of disturbances on the jet interface is related with turbulent jets that are characterized by the jet Reynolds number $4\hat {Q}/(\hat {\nu } {\rm \pi}\hat {d})$ exceeding approximately 4000 for a smooth circular nozzle (Lienhard 2006). In the current study the liquid jets under consideration are turbulent, falling within the range of $4\hat {Q}/(\hat {\nu } {\rm \pi}\hat {d})=4200\unicode{x2013}21\,000$. For scenarios with smaller $\hat {Q}$ or $\hat {H}$, concentric waves can be suppressed according to Wang et al. (2023). The comprehensive analysis about the effect of disk angular velocity $\hat {\varOmega }$ on the initiation of concentric waves has not been conducted. Nevertheless, we assume that the effect of $\hat {\varOmega }$ is negligible for the following reasons. Generation of concentric waves in the impinging zone occurs before the interface is influenced by disk rotation, and the tangential velocity by disk rotation is small near the disk centre; the effect of $\hat {\varOmega }$ becomes significant only after the boundary layer merger ($\hat {r} > \hat {r}_{0}$). Moreover, concentric waves are observed on a liquid film even without disk rotation ($\hat {\varOmega }=0$) when a turbulent liquid jet impinges on a solid surface (Lienhard 2006).

To investigate the effects of liquid jet disturbance, concentric waves are visualized in the upstream region under various nozzle heights (figure 4). The volume flow rate $\hat {Q}$ and disk angular velocity $\hat {\varOmega }$ are the same in all cases. As the vertical distance between the nozzle and the disk increases, concentric waves with larger amplitudes and wavelengths appear. Furthermore, the concentric waves exhibit greater azimuthal disturbances when the nozzle is placed closer to the disk. Previous studies reported that surface disturbances on circular liquid jets are significantly influenced by the distance from the nozzle (Lienhard, Liu & Gabour 1992; Bhunia & Lienhard 1994). Thus, our experimental findings indicate that the concentric waves in the upstream region are induced by the temporal fluctuations of the impinging jet due to disturbances on the jet interface. The concentric waves are modelled with an oscillatory input, having statistical characteristics that match the results from the sensor experiment. The modelling of the oscillatory input is discussed in § 4.3.

Figure 4. Formation of concentric surface waves in the upstream region under different nozzle heights $\hat {H}$; for all cases, $\hat {Q}= 16.7\,{\rm mL}\,{\rm s}^{-1}$ and $\hat {\varOmega }=52.4\,{\rm rad}\,{\rm s}^{-1}$. Results are shown for (a) $\hat {H}=40\,{\rm mm}$, (b) $\hat {H}=20\,{\rm mm}$, (c) $\hat {H}=15\,{\rm mm}$.

3. Integral modelling of film flow

The derivation of low-dimensional models for both laminar and turbulent regimes follows a similar process involving depth integration and long-wave assumptions. The governing equations integrated along the vertical direction include closure terms that should be determined by different assumptions regarding the velocity profiles in the vertical direction. The closure terms depend on whether the flow is in the laminar or turbulent regime. In this section we first present the depth integration and long-wave formulation, and then consider the application of closure models based on specific velocity profiles for laminar and turbulent conditions.

3.1. Integral formulation of governing equations

The film flow is described by mass and momentum conservation equations, with boundary conditions applied on the liquid–gas interface ($\hat {z}=\hat {h}$) and the disk surface ($\hat {z}=0$). The governing equations and boundary conditions are transformed into a set of depth-integrated equations, employing the methodology presented by Kim & Kim (2009) to derive IBL equations. As mentioned in § 2.1, the IBL equations presented below are applicable for $\hat {r}>\hat {r}_{0}$, irrespective of the local Reynolds number. To make the governing equations and boundary conditions dimensionless, the following normalization is adopted:

(3.1)\begin{equation} \left. \begin{aligned} r & = \frac{\hat{r}}{\hat{l}}, \quad z = \frac{\hat{z}}{\hat{\delta}}, \quad t = \frac{1}{\hat{\varOmega}},\\ u_{r} & = \frac{\hat{u}_{r}}{\hat{u}_{0}}, \quad u_{\theta} = \frac{\hat{u}_{\theta}}{\hat{u}_{0}}, \quad h = \frac{\hat{h}}{\hat{\delta}}, \quad u_{z} = \frac{\hat{u}_{z}}{\hat{\delta}\hat{\varOmega}}, \quad p = \frac{\hat{p}}{\hat{\rho}\hat{l}^{2}\hat{\varOmega}^{2}}. \end{aligned} \right\} \end{equation}

The dimensional pressure and time are denoted by $\hat {p}$ and $\hat {t}$, respectively. The characteristic length scales for the horizontal direction, $\hat {l}=(9\hat {Q}^{2}/ 4{\rm \pi} ^{2}\hat {\nu }\hat {\varOmega })^{1/4}$, and the vertical direction, $\hat {\delta }=(\hat {\nu }/\hat {\varOmega })^{1/2}$, are chosen following Rauscher, Kelly & Cole (1973). The reference horizontal velocity $\hat {u}_{0}$ is defined as $\hat {l}\hat {\varOmega }$. For a detailed discussion on scaling parameters and normalization, see Kim & Kim (2009). The scaling proposed by Kim & Kim (2009), which uses a universal length scale for both the radial and azimuthal directions, is more suitable for modelling three-dimensional waves than scaling based on the Ekman number (Sisoev et al. 2003). This is because the velocity fluctuations resulting from continuous interactions among solitary waves in the downstream region means that the velocity scale in the azimuthal direction is no longer distinct from that in the radial direction. The same scaling holds upstream because the influence of the inertia force surpasses that of the rotational body forces.

The governing equations and boundary conditions are formulated using the long-wave approach, which assumes that the ratio between the horizontal length scale $\hat {l}$ and the depth scale $\hat {\delta }$ is much smaller than unity. In this study a dimensionless long-wave parameter $\epsilon = \hat {\delta }/\hat {l} \ll 1$ is introduced. The product of $\epsilon$ and the global Reynolds number, defined as $Re_{G} = \hat {u}_{0}\hat {\delta }/\hat {\nu }$, is of order unity. The global Reynolds number, therefore, does not appear explicitly in the normalized governing equations and boundary conditions. The dimensionless governing equations in the reference frame rotating with the disk are then presented as

(3.2a)$$\begin{gather} \frac{1}{r}\frac{\partial}{\partial r}\left(ru_{r}\right)+\frac{1}{r}\frac{\partial u_{\theta}}{\partial \theta}+\frac{\partial u_{z}}{\partial z} =0, \end{gather}$$

(3.2b)$$\begin{gather} \frac{\partial u_{r}}{\partial t} + u_{r}\frac{\partial u_{r}}{\partial r} + \frac{u_{\theta}}{r}\frac{\partial u_{r}}{\partial \theta} + u_{z}\frac{\partial u_{r}}{\partial z} - \frac{u_{\theta}^{2}}{r} ={-} \frac{\partial p}{\partial r} + r + 2u_{\theta} + \frac{\partial^{2} u_{r}}{\partial z^{2}} \nonumber\\ +\ \epsilon^{2}\left(\frac{\partial^{2} u_{r}}{\partial r^{2}} + \frac{1}{r}\frac{\partial u_{r}}{\partial r} + \frac{1}{r^{2}}\frac{\partial^{2} u_{r}}{\partial \theta^{2}} - \frac{2}{r^{2}}\frac{\partial u_{\theta}}{\partial \theta} - \frac{u_{r}}{r^{2}}\right), \end{gather}$$

(3.2c)$$\begin{gather} \frac{\partial u_{\theta}}{\partial t} + u_{r}\frac{\partial u_{\theta}}{\partial r} + \frac{u_{\theta}}{r}\frac{\partial u_{\theta}}{\partial \theta} + u_{z}\frac{\partial u_{\theta}}{\partial z} + \frac{u_{r}u_{\theta}}{r} ={-} \frac{1}{r}\frac{\partial p}{\partial \theta} - 2u_{r} + \frac{\partial^{2} u_{\theta}}{\partial z^{2}} \nonumber\\ +\ \epsilon^{2}\left(\frac{\partial^{2} u_{\theta}}{\partial r^{2}} + \frac{1}{r}\frac{\partial u_{\theta}}{\partial r} + \frac{1}{r^{2}}\frac{\partial^{2} u_{\theta}}{\partial \theta^{2}} + \frac{2}{r^{2}}\frac{\partial u_{r}}{\partial \theta} - \frac{u_{\theta}}{r^{2}}\right), \end{gather}$$

(3.2d)$$\begin{gather} \epsilon^{2}\left(\frac{\partial u_{z}}{\partial t} + u_{r}\frac{\partial u_{z}}{\partial r} + \frac{u_{\theta}}{r}\frac{\partial u_{z}}{\partial \theta} + u_{z}\frac{\partial u_{z}}{\partial z}\right) ={-} \frac{\partial p}{\partial z} - \epsilon Fr^{{-}1} + \epsilon^{2}\frac{\partial^{2} u_{z}}{\partial z^{2}} \nonumber\\ +\ \epsilon^{4}\left(\frac{\partial^{2} u_{z}}{\partial r^{2}} + \frac{1}{r}\frac{\partial u_{z}}{\partial r} + \frac{1}{r^{2}}\frac{\partial^{2} u_{z}}{\partial \theta^{2}}\right), \end{gather}$$

where $Fr=\hat {\varOmega }^{2}\hat {l}/\hat {g}$ is the Froude number. Equation (3.2a) is the dimensionless continuity equation, whereas (3.2b–d) are the dimensionless momentum equations in the $r$, $\theta$ and $z$ directions, respectively.

The boundary condition on the disk surface ($z = 0$) is the no-slip condition, and the kinematic boundary condition is imposed on the film surface ($z = h$). The dimensionless forms of these boundary conditions are

(3.3a)$$\begin{gather} u_{r}=u_{\theta}=u_{z}=0,\quad \text{at} \ z=0, \end{gather}$$

(3.3b)$$\begin{gather}\frac{\partial h}{\partial t} + u_{r}\frac{\partial h}{\partial r} + \frac{u_{\theta}}{r}\frac{\partial h}{\partial \theta} = u_{z} ,\quad \text{at}\ z=h. \end{gather}$$

Tangential and normal stress balance conditions should also be imposed on the film surface ($z = h$). Here, the shear stress induced by the gas flow above the film surface is neglected (stress-free condition). These boundary conditions are

(3.4a)\begin{align} &\left(\epsilon^{2} \frac{\partial u_{z}}{\partial r} + \frac{\partial u_{r}}{\partial z}\right)\left[1 - \epsilon^{2}\left(\frac{\partial h}{\partial r}\right)^{2}\right] + 2\epsilon^{2}\left(\frac{\partial u_{z}}{\partial z} - \frac{\partial u_{r}}{\partial r}\right)\frac{\partial h}{\partial r} \nonumber\\ &\qquad-\epsilon^{2} \frac{1}{r}\frac{\partial h}{\partial \theta}\left[\frac{1}{r}\frac{\partial u_{r}}{\partial \theta} + \frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r} + \left(\frac{\partial u_{\theta}}{\partial z} + \epsilon^{2}\frac{1}{r}\frac{\partial u_{z}}{\partial \theta}\right)\frac{\partial h}{\partial r}\right] = 0,\quad \text{at} \ z = h, \end{align}

(3.4b)\begin{align} &\left(\epsilon^{2} \frac{1}{r}\frac{\partial u_{z}}{\partial \theta} + \frac{\partial u_{\theta}}{\partial z}\right)\left[1 - \epsilon^{2}\left(\frac{1}{r}\frac{\partial h}{\partial \theta}\right)^{2}\right] + 2\epsilon^{2}\left(\frac{\partial u_{z}}{\partial z} - \frac{1}{r}\frac{\partial u_{\theta}}{\partial \theta} - \frac{u_{r}}{r}\right)\frac{1}{r}\frac{\partial h}{\partial \theta} \nonumber\\ &\qquad-\epsilon^{2}\frac{\partial h}{\partial r}\left[\frac{1}{r}\frac{\partial u_{r}}{\partial \theta} + \frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r} + \left(\frac{\partial u_{r}}{\partial z} + \epsilon^{2}\frac{\partial u_{z}}{\partial r}\right)\frac{1}{r}\frac{\partial h}{\partial \theta}\right] = 0,\quad \text{at} \ z = h, \end{align}

(3.4c) \begin{align} &{-} p - \epsilon^{3}\kappa We + 2\epsilon^{2}\left[1 + \epsilon^{2}\left(\frac{\partial h}{\partial r}\right)^{2} + \frac{\epsilon^{2}}{r^{2}}\left(\frac{\partial h}{\partial \theta}\right)^{2}\right]^{{-}1} \nonumber\\ &\qquad \times \left[\frac{\partial u_{z}}{\partial z} - \left(\frac{\partial u_{r}}{\partial z} + \epsilon^{2}\frac{\partial u_{z}}{\partial r}\right)\frac{\partial h}{\partial r} + \epsilon^{2}\frac{\partial u_{r}}{\partial r}\left(\frac{\partial h}{\partial r}\right)^{2} - \left(\frac{\partial u_{\theta}}{\partial z} + \frac{\epsilon^{2}}{r}\frac{\partial u_{z}}{\partial \theta}\right)\frac{1}{r}\frac{\partial h}{\partial \theta}\right. \nonumber\\ &\qquad +\left. \epsilon^{2}\left(\frac{1}{r}\frac{\partial u_{r}}{\partial \theta} + \frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r}\right)\frac{1}{r}\frac{\partial h}{\partial \theta}\frac{\partial h}{\partial r} + \epsilon^{2}\left(\frac{1}{r}\frac{\partial u_{\theta}}{\partial \theta} + \frac{\partial u_{r}}{\partial r}\right)\frac{1}{r^{2}}\left(\frac{\partial h}{\partial \theta}\right)^{2} \right] = 0,\nonumber\\&\quad \text{at} \ z = h, \end{align}

where $\kappa$ denotes the mean curvature of the free surface and $We = \hat {\sigma }/\hat {\rho }\hat {\varOmega }^{2}\hat {l}\hat {\delta }^{2}$ is the Weber number. Boundary conditions (3.4a,b) concern the stress balance in the tangential direction, whereas (3.4c) is the stress balance in the normal direction. The unknown film thickness and small $\epsilon$ make numerical solutions from governing equations (3.2) subject to boundary conditions (3.3) and (3.4) challenging. Hence, the boundary layer approximation is introduced along with depth integration to simplify the equations.

Before the integral approach is employed to relate the film thickness $h$ with the horizontal velocity components $u_r$ and $u_\theta$, terms of $O(\epsilon ^{2})$ and higher orders are neglected in accordance with the first-order long-wave approximation. The stress-free boundary condition (3.4) then yields

(3.5a–c)\begin{equation} \left.\frac{\partial u_{r}}{\partial z}\right\rvert_{z = h} = 0, \quad \left.\frac{\partial u_{\theta}}{\partial z}\right\rvert_{z = h} = 0, \quad p\rvert_{z=h} ={-}\epsilon \kappa \mathbb{W}, \end{equation}

where $\mathbb {W} = \epsilon ^{2}We (= \hat {\sigma }/\hat {\rho }\hat {\varOmega }^{2}\hat {l}^{3})$. By applying the pressure on the free surface (3.5c) and the long-wave approximation ($\epsilon \ll 1$), the following pressure distribution is obtained by integrating the vertical component of the momentum balance equation (3.2d) from an arbitrary location $z$ to $z=h$ along the $z$ direction:

(3.6)\begin{equation} p ={-}\epsilon \kappa \mathbb{W} + \epsilon Fr^{{-}1}(h-z). \end{equation}

Incorporating the pressure equation (3.6), simplified stress-free boundary conditions (3.5a,b) and other boundary conditions (3.3), the continuity equation (3.2a) and momentum conservation equations in the $r$ and $\theta$ directions (3.2b,c) are integrated along the $z$ direction from $z=0$ to $z=h$ using the Leibniz integral rule. To facilitate the presentation of the integrated equations, we introduce the depth-averaged velocity components $\bar {u}_{r}(=({1}/{h})\int _0^h u_r\,\mathrm {d}z)$ and $\bar {u}_{\theta }(=({1}/{h})\int _0^h u_\theta \,\mathrm {d}z)$. Consequently, the depth-integrated equations, which establish the relationship between $h$, $\bar {u}_{r}$ and $\bar {u}_{\theta }$, are expressed as

(3.7a)$$\begin{gather} \frac{\partial h}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}(r\bar{u}_{r}h) + \frac{1}{r}\frac{\partial}{\partial \theta}(\bar{u}_{\theta}h) = 0, \end{gather}$$

(3.7b)$$\begin{gather} \frac{\partial \bar{u}_{r}h}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}(r\overline{u_{r}^{2}}h) + \frac{1}{r}\frac{\partial}{\partial \theta}(\overline{u_{r}u_{\theta}}h) - \frac{\overline{u_{\theta}^{2}}h}{r} ={-} \frac{\partial h\bar{p}}{\partial r} & + p\rvert_{z=h}\frac{\partial h}{\partial r} + rh \nonumber\\ & +\ 2 \bar{u}_{\theta}h -\left. \frac{\partial u_{r}}{\partial z}\right\rvert_{z=0}, \end{gather}$$

(3.7c)$$\begin{gather} \frac{\partial \bar{u}_{\theta}h}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}(r\overline{u_{r}u_{\theta}}h) + \frac{1}{r}\frac{\partial}{\partial \theta}(\overline{u_{\theta}^{2}}h) + \frac{\overline{u_{r}u_{\theta}}h}{r} & ={-} \frac{1}{r}\frac{\partial h\bar{p}}{\partial \theta} + p\rvert_{z=h}\frac{1}{r}\frac{\partial h}{\partial \theta} \nonumber\\ & -\ 2 \bar{u}_{r}h - \left.\frac{\partial u_{\theta}}{\partial z}\right\rvert_{z=0}, \end{gather}$$

(3.7d)$$\begin{gather}\bar{p} ={-}\epsilon\kappa\mathbb{W} + \frac{1}{2}\epsilon Fr^{{-}1} h, \end{gather}$$

where the terms $rh$, $2\bar {u}_{r}h$ and $2\bar {u}_{\theta }h$ on the right-hand side come from the centrifugal and Coriolis forces, respectively. The gravitational body force is incorporated into the pressure equation (3.7d) in terms of $Fr^{-1}$, which comes from (3.6). The mean curvature $\kappa$ of the free surface is expressed as $\kappa =\boldsymbol {\nabla } \boldsymbol {\cdot } [\boldsymbol {\nabla } h/(|\boldsymbol {\nabla }h|^{2}+1)^{1/2}]$. Here, the gradient operator $\boldsymbol {\nabla }$ is defined in the horizontal coordinates as $\boldsymbol {\nabla } = ({\partial }/{\partial r})\boldsymbol {e}_{r} + ({1}/{r})({\partial }/{\partial \theta })\boldsymbol {e}_{\theta }$.

The variables of the governing equations (3.7) are $h$, $\bar {u}_{r}$ and $\bar {u}_{\theta }$; the depth-averaged pressure $\bar {p}$ in (3.7d) is a function of $h$. Thus, the equations are not closed because of the convective terms ($\overline {u_{r}^{2}}$, $\overline {u_{r} u_{\theta }}$, $\overline {u_{\theta }^{2}}$) and bottom shear terms ($\partial u_{r}/\partial z\rvert _{z=0}$, $\partial u_{\theta }/\partial z\rvert _{z=0}$). In the Kapitza–Shkadov IBL equations, the convective and bottom shear terms are closed by assuming a self-similar velocity profile (Kalliadasis et al. 2012). A specific form of the velocity profile enables the closure terms to be defined in terms of the film thickness and depth-averaged velocity components:

(3.8a–c)$$\begin{gather} \overline{u_{r}^{2}} = k_{A}\bar{u}_{r}^{2}, \quad \overline{u_{r}u_{\theta}} = k_{B}\bar{u}_{r}\,\bar{u}_{\theta}, \quad \overline{u_{\theta}^{2}} = k_{C}\bar{u}_{\theta}^{2}, \end{gather}$$

(3.8d,e)$$\begin{gather}\left.\frac{\partial u_{r}}{\partial z}\right\rvert_{z=0} = f_{r}\frac{\bar{u}_{r}}{h}, \quad \left.\frac{\partial u_{\theta}}{\partial z}\right\rvert_{z=0} = f_{\theta}\frac{\bar{u}_{\theta}}{h}. \end{gather}$$

Here $k_{A}$, $k_{B}$, $k_{C}$, $f_{r}$ and $f_{\theta }$ are constants determined by the velocity profile along the vertical direction. These closure models, based on a single self-similar polynomial velocity profile, have successfully predicted the laminar film flow up to the global Reynolds number $Re_{G} \approx 75$ in falling films (Demekhin et al. 2007; Denner et al. 2018). However, the velocity profile deviates from the polynomial distribution for $Re_{G} > O(10^{2})$, as mentioned in § 2.3. Given the high local flow rate in the region near the jet impingement zone, it is essential to encompass the vertical velocity profile under high-Reynolds-number conditions in the closure models.

In this study we propose that the velocity profile along the $z$ direction should be chosen locally at any instant, in terms of the local Reynolds number $Re_{L} = \hat {q}/\hat {\nu }$ based on the local horizontal flow rate. We can express $Re_{L}$, using the IBL equation variables $\hat {h}$ and $\bar {\hat {\boldsymbol {u}}}$ as

(3.9)\begin{equation} Re_{L}=\frac{\hat{h}\rvert\bar{\hat{\boldsymbol{u}}}\rvert}{\hat{\nu}}, \end{equation}

where $\bar {\hat {\boldsymbol {u}}}$ denotes the depth-averaged horizontal velocity vector in dimensional form: $({1}/{\hat {h}})\int _{0}^{\hat {h}} (\hat {u}_{r},\hat {u}_{\theta }) \,\mathrm {d}\hat {z}$. As mentioned in § 2.3, the flow field is divided into laminar and turbulent film flow based on the local Reynolds number $Re_{L}$. The critical Reynolds number $Re_{*}$, which acts as a threshold between these two regimes, is set to 100, following the criterion proposed by Mendez et al. (2021). Figure 5 conceptually illustrates velocity profiles in the laminar and turbulent film flow regimes. For regions with $Re_{L} \leq 100$, the closure models for the laminar film flow are based on a polynomial velocity profile. By contrast, a power-law velocity profile is used to describe turbulent film flow for regions with $Re_{L} > 100$, where sub-depth scale vortical structures emerge near the bottom.

Figure 5. Schematics of velocity profiles along the vertical direction for two flow regimes based on the local Reynolds number $Re_{L}$. The solid line denotes the instantaneous velocity profile and the dashed line is the time-averaged velocity profile. Results are shown for (a) $Re_{L}\leq 100$ and (b) $Re_{L}>100$.

3.2. Closure model for laminar regime

In the laminar regime the closure terms are determined based on a self-similar polynomial velocity profile in the vertical direction. The shape of the velocity profile is considered to be independent of $Re_{L}$ in this regime, as experimentally verified under a low flow rate (Dietze, Al-sibai & Kneer 2009). The following velocity profile assumption introduced by Sisoev et al. (2003) in the Kapitza–Shkadov method for rotating films has been widely adopted: $u_{r} = 3\bar {u}_{r}(\eta - \tfrac {1}{2}\eta ^{2})$ and $u_{\theta } = 5\bar {u}_{\theta }\tfrac {1}{4}(2\eta - \eta ^{3} + \tfrac {1}{4}\eta ^{4})$, where $\eta = z/h$ varies from 0 to 1. The parabolic and quartic velocity profiles are obtained as asymptotic solutions of the governing equations in the limit of a large Ekman number (Shkadov 1973). The primary assumption underlying the derivation of separate profiles in the $r$ and $\theta$ directions is the dominance of viscous and body forces over inertia forces. However, the present study is characterized by relatively large flow rates, so this assumption is no longer applicable. Thus, we suggest a single quartic velocity profile for both radial and tangential directions.

Both the radial and tangential velocity profiles are assumed to follow the quartic formulations

(3.10a)$$\begin{gather} u_{r}(r,\theta,z,t)=\frac{5}{4}\bar{u}_{r}(r,\theta,t)\left(2\eta-\eta^{3}+\frac{\eta^{4}}{4}\right), \end{gather}$$

(3.10b)$$\begin{gather}u_{\theta}(r,\theta,z,t)=\frac{5}{4}\bar{u}_{\theta}(r,\theta,t)\left(2\eta-\eta^{3}+\frac{\eta^{4}}{4}\right). \end{gather}$$

These velocity profiles satisfy the boundary conditions (3.3a) and (3.5a,b). Quartic profiles were partially adopted by Kim & Kim (2009), where the time-averaged film thickness acquired from the numerical results of IBL equations was found to be in good agreement with experimental results. A rationale behind adopting the quartic profile is extensively discussed in Appendix B. The closure terms in (3.8) for the laminar regime ($Re_{L} \leq Re_{*}$, where $Re_{*} = 100$) are calculated from (3.10) as

(3.11a,b)\begin{equation} k_{A}=k_{B}=k_{C}=k_{l}=\frac{155}{126}, \quad f_{r}=f_{\theta}=\frac{5}{2}, \end{equation}

where $k_{l}$ denotes the advection closure constant in the laminar regime.

Despite numerous benefits of using more advanced and complicated mathematical models, including the wall-resolved inner boundary layer approach that accurately predicts the threshold of linear stability, the Kármán–Pohlhausen-type integral approach employed in the present work is sufficient to capture the dynamics of three-dimensional solitary waves (Chang, Demekhin & Kopelevitch 2006), which still remain elusive for rotating films. Moreover, it excels in integrating the power-law profiles in the high-$Re_{L}$ regime into the IBL framework.

3.3. Closure model for turbulent regime

Next, we present closure models for the turbulent regime ($Re_{L} > Re_{*}$). Analysis of the film flows in this regime often involves mixing length theory (King 1966; Geshev 2014) and shallow-water assumptions (Mukhopadhyay, Chhay & Ruyer-Quil 2017; James et al. 2019; Mendez et al. 2021), as well as investigation of the characteristics of roll waves (Nakoryakov, Ostapenko & Bartashevich 2012; Yu & Chu 2022). Among these approaches, the one particularly suitable to IBL formulations is the shallow-water dynamics (James et al. 2019). The shallow-water equations, which are derived through the long-wave approximation and depth averaging, exhibit compatibility with the framework of the present study. In turbulent shallow water, the power-law velocity profile is assumed instead of the polynomial profile. We refer to the shallow-water model proposed for liquid films undergoing jet wiping (Mendez et al. 2021) and develop a turbulent closure model that is specifically tailored for rotating films. In addition to the velocity profile, the influence of sub-depth scale turbulent structures on the horizontal momentum balance, as suggested in studies on shallow-water flows, is also incorporated. Regarding film flows over a rotating disk, this is the first application of the shallow-water model to the IBL equations, to the best of our knowledge.

The turbulent flow requires the modelling of length scales to be divided into resolved and unresolved components. The resolved scales encompass greater horizontal flow structures, whereas the unresolved scales represent smaller internal structures. To address both scales, the velocity components are initially decomposed into filtered (resolved) and sub-grid scale (unresolved) components. This decomposition allows for the derivation of filtered versions of the depth-averaged equations, drawing inspiration from the literature on depth-averaged large-eddy simulations (DA LES) for turbulent shallow waters (Hinterberger et al. 2007; van Prooijen & Ujittewaal 2009). In these simulations, the grid size is typically comparable to the water depth, leading to the term sub-depth scale instead of sub-grid scale. The filtered velocity is then assumed to have a power-law profile along the vertical direction, which is used to derive the closure equations (3.8). As for the unfiltered components, sub-depth scale stress is interpreted as a random force vector field applied to the governing equations of the resolved scale, which is referred to as a stochastic backscatter model.

3.3.1. Filtered formulation

To capture the influence of turbulent structures, the concept of filtering, commonly used in DA LES, is introduced to the IBL formulation. Let us decompose the velocity vector $\boldsymbol {u} = (u_{r},u_{\theta },u_{z})$ into $\boldsymbol {u}=\tilde {\boldsymbol {u}}+\boldsymbol {u}'$, where $\tilde {\boldsymbol {u}}$ denotes the thickness-based Favre-filtered velocity $\tilde {\boldsymbol {u}} = \underline {h\boldsymbol {u}}/\underline {h}$; the underline marker ($\underline {\cdot }$) denotes a spatial filter $\underline {\varphi }(\boldsymbol {x})= \int _D \varphi (\boldsymbol {x}')G(\boldsymbol {x}'-\boldsymbol {x}) \,\textrm {d}\,\boldsymbol {x}'$ (Hinterberger et al. 2007). The pressure is similarly decomposed. Filtering the continuity and horizontal momentum equations (3.2a–c) yields

(3.12a)$$\begin{gather} \frac{1}{r}\frac{\partial}{\partial r}(r\tilde{u}_{r})+\frac{1}{r}\frac{\partial \tilde{u}_{\theta}}{\partial \theta}+\frac{\partial \tilde{u}_{z}}{\partial z}=0, \end{gather}$$

(3.12b)$$\begin{gather}\frac{\partial \tilde{u}_{r}}{\partial t} + \frac{\partial \widetilde{u_{r}^{2}}}{\partial r} + \frac{1}{r}\frac{\partial \widetilde{u_{r}u_{\theta}}}{\partial \theta} + \frac{\partial \widetilde{u_{r}u_{z}}}{\partial z} - \frac{\widetilde{u_{r}^{2}}-\widetilde{u_{\theta}^{2}}}{r} ={-} \frac{\partial \tilde{p}}{\partial r} + r + 2\tilde{u}_{\theta} + \frac{\partial^{2} \tilde{u}_{r}}{\partial z^{2}}, \end{gather}$$

(3.12c)$$\begin{gather}\frac{\partial \tilde{u}_{\theta}}{\partial t} + \frac{\partial \widetilde{u_{\theta}u_{r}}}{\partial r} + \frac{1}{r}\frac{\partial \widetilde{u_{\theta}^{2}}}{\partial \theta} + \frac{\partial \widetilde{u_{\theta}u_{z}}}{\partial z} + \frac{2\widetilde{u_{r}u_{\theta}}}{r} ={-} \frac{1}{r}\frac{\partial \tilde{p}}{\partial \theta} - 2\tilde{u}_{r} + \frac{\partial^{2} \tilde{u}_{\theta}}{\partial z^{2}}, \end{gather}$$

where the long-wave approximation still holds for the filtered velocity scales. By applying the Favre-filtering technique to the governing momentum equation in the $z$ direction (3.2d) and the corresponding boundary conditions (3.3) and (3.4), the velocity and pressure terms of $O(\epsilon ^{2})$ and lower are replaced with their respective Favre-filtered variables. Consequently, the filtered vertical momentum equation and boundary conditions yield the filtered pressure, which is equivalent to (3.6) with $h$ replaced by $\underline{h}$:

(3.13)\begin{equation} \tilde{p} ={-} \epsilon \tilde{\kappa} \mathbb{W} + \epsilon Fr^{{-}1} (\underline{h}-z). \end{equation}

Here $\tilde {\kappa } = \boldsymbol {\nabla } \boldsymbol {\cdot } [\boldsymbol {\nabla } \underline {h}/ (\rvert \boldsymbol {\nabla } \underline {h} \rvert ^{2}+1)^{1/2}]$.

In the filtered equations (3.12), all terms except advection terms are expressed with $\tilde {u}_{r}, \tilde {u}_{\theta }, \tilde {u}_{z}$ and $\tilde {p}$. The distinction between the filtered advection and the advection of the filtered velocity, denoted as $\partial (\widetilde {u_{i}u_{j}} - \tilde {u}_{i}\tilde {u}_{j})/\partial x_{j}$ in Cartesian tensor notation, plays a crucial role in LES, giving rise to the concepts of sub-grid scale stress and eddy viscosity. In the present work, the filtered velocity products $\widetilde {u_{r}^{2}}, \widetilde {u_{\theta }^{2}}, \widetilde {u_{r}u_{\theta }}, \widetilde {u_{r}u_{z}}, \widetilde {u_{\theta }u_{z}}$ are likewise decomposed into $\tilde {u}_{r}^{2}, \tilde {u}_{\theta }^{2}, \tilde {u}_{r}\tilde {u}_{\theta }, \tilde {u}_{r}\tilde {u}_{z}, \tilde {u}_{\theta }\tilde {u}_{z}$ and the sub-grid scale stress terms $(\widetilde {u_{i}u_{j}} - \tilde {u}_{i}\,\tilde {u}_{j})$.

Following previous studies on shallow-water LES (Nadaoka & Yagi 1998; Hinterberger et al. 2007), the Favre-filtered velocity components are interpreted as being grid resolved, and are considered identical to the velocity components without filtering in the laminar regime. The scaling parameters in (3.1) and the local Reynolds number in (3.9) are defined based on the grid-resolved values. Assuming identical vertical velocity profiles for $\tilde {u}_{r}$ and $\tilde {u}_{\theta }$, the depth integration of (3.12) using the pressure equation (3.13) and boundary conditions leads to a set of first-order long-wave boundary layer equations. The governing equations are expressed in Cartesian tensor notation to maintain consistency with existing LES formulations. Additionally, describing small-scale turbulence near the disk surface is challenging within the framework of cylindrical coordinates. The depth-integrated governing equations are

(3.14a)$$\begin{gather} \frac{\partial \underline{h}}{\partial t}+\frac{\partial }{\partial x_{i}} (\underline{h} \bar{\tilde{u}}_{i})=0, \end{gather}$$

(3.14b)$$\begin{gather} \frac{\partial (\underline{h} \bar{\tilde{u}}_{i})}{\partial t} + k_{t}\frac{\partial}{\partial x_{j}} (\underline{h} \bar{\tilde{u}}_{i}\,\bar{\tilde{u}}_{j}) & ={-} \frac{\partial \underline{h} \bar{\tilde{p}}}{\partial x_{i}} - \epsilon \tilde{\kappa} \mathbb{W}\frac{\partial \underline{h}}{\partial x_{i}} - \frac{\partial}{\partial x_{j}}(\underline{h} G_{ij}) \nonumber\\ & -\ \frac{1}{2}C_{f}\bar{\tilde{u}}_{i} \left(\bar{\tilde{u}}_{i}^{2} + \bar{\tilde{u}}_{j}^{2}\right)^{1/2} + B_{i}, \end{gather}$$

(3.14c)$$\begin{gather} \bar{\tilde{p}} ={-}\epsilon \tilde{\kappa} \mathbb{W} + \frac{1}{2} \epsilon Fr^{{-}1}\underline{h}, \end{gather}$$

where $(i,j = x,y)$ and $B_{i}$ indicates filtered and depth-integrated body forces (centrifugal force and Coriolis force terms). The depth-integrated pressure $\bar {\tilde {p}}$ in (3.14c) is obtained from (3.6), with unfiltered values replaced by filtered ones.

The constant $k_{t}$ in the advection term of (3.14b) is related to the assumption of a specific profile for the grid-resolved velocity in the vertical direction, which corresponds to (3.8a–c) and (3.11a) for the laminar regime. Here $k_{t}$ is equivalent to the constants $k_{A}$, $k_{B}$ and $k_{C}$; i.e. $k_{t}\bar {\tilde {u}}_{i} \,\bar {\tilde {u}}_{j} = \overline {\tilde {u}_{i} \tilde {u}_{j}}$. The fourth term on the right-hand side of (3.14b) represents bottom friction, which corresponds to the friction closure terms in (3.8d,e) and (3.11b) for the laminar regime. The determination of the closure models for $k_{t}$ and $C_{f}$ relies on the selection of the resolved velocity profile in the vertical direction, which is depicted in § 3.3.2. Additionally, $G_{ij}$ on the right-hand side, defined as $G_{ij} = \overline {\widetilde {u_{i}u_{j}}} - \overline {\tilde {u}_{i}\tilde {u}_{j}}$, represents the sub-grid scale stress arising from the filtering operator. This closure term for the sub-grid scale is modelled as a random force field using the stochastic backscatter model. The procedure for deriving the random force term is described in § 3.3.3.

3.3.2. Closure model for resolved velocity

The closure terms $k_{t}$ and $C_{f}$ in (3.14) are formulated based on assumptions regarding the resolved (filtered) velocity profile. The friction coefficient $C_f$ is inversely proportional to $Re_{L}$ in the laminar regime ($Re_L \leq Re_{*}$); $C_{f} = 5/Re_{L}$. For the turbulent regime ($Re_L > Re_{*}$), the friction coefficient is also regarded as a function of $Re_{L}$, exhibiting continuity with the laminar model at $Re_{L} = Re_{*}$, as suggested by Mendez et al. (2021). Specifically, we adopt an empirical relation of the form $C_{f} \approx a Re_{L}^{b}$, where $b \approx -1/4$ has been documented in studies of turbulent lubrication films (Elrod & Ng 1967; Hirs 1973). To ensure continuity at the transition point $Re_{*}$, the following expressions are proposed for the friction coefficients of the laminar and turbulent regimes:

(3.15)\begin{equation} C_{f} = \begin{cases} 5/Re_{L}, & Re_{L} \leq Re_{*}, \\ 5/(Re_{L}^{1/4}Re_*^{3/4}), & Re_{L} > Re_{*}. \end{cases} \end{equation}

The trend of $C_{f}$ in the vicinity of the transition regime is depicted in figure 6(a).

Figure 6. (a) Trends of friction coefficient $C_f$ with respect to local Reynolds number $Re_L$ ($Re_{*}=100$). (b) Vertical velocity profiles for different local Reynolds numbers ($n_{T}=21$). The subscript $i (= x,y)$ indicates each horizontal direction.

To determine the closure terms for advection, it is crucial to establish a self-similar velocity profile, as outlined in § 3.2. To address transition from laminar to turbulent flow, the combination of the polynomial profile and the power-law profile, which represent the laminar and turbulent regimes, respectively, can be considered (Mendez et al. 2021). As the local Reynolds number increases, the turbulent component becomes more dominant. The following mathematical expressions give the velocity profiles for the turbulent regime:

(3.16a) \begin{align} \tilde{u}_{x}(x,y,z,t)&=a_{L}\bar{\tilde{u}}_{x}(x,y,t)\left[\left(\frac{z}{\underline{h}}\right)-\frac{1}{2}\left(\frac{z}{\underline{h}}\right)^{3}+\frac{1}{8}\left(\frac{z}{\underline{h}}\right)^{4}\right] \nonumber \\ &\quad+ a_{T}\bar{\tilde{u}}_{x}(x,y,t)\left[\vphantom{\left(\frac{z}{\underline{h}}\right)-\frac{1}{2}\left(\frac{z}{\underline{h}}\right)^{3}+\frac{1}{8}\left(\frac{z}{\underline{h}}\right)^{4}}\left(\frac{z}{\underline{h}}-1\right)^{n_{T}}+1\vphantom{\left(\frac{z}{\underline{h}}\right)-\frac{1}{2}\left(\frac{z}{\underline{h}}\right)^{3}+\frac{1}{8}\left(\frac{z}{\underline{h}}\right)^{4}}\right], \end{align}

(3.16b) \begin{align} \tilde{u}_{y}(x,y,z,t)&=a_{L}\bar{\tilde{u}}_{y}(x,y,t)\left[\left(\frac{z}{\underline{h}}\right)-\frac{1}{2}\left(\frac{z}{\underline{h}}\right)^{3}+\frac{1}{8}\left(\frac{z}{\underline{h}}\right)^{4}\right] \nonumber\\ &\quad+ a_{T}\bar{\tilde{u}}_{y}(x,y,t)\left[\vphantom{\left(\frac{z}{\underline{h}}\right)-\frac{1}{2}\left(\frac{z}{\underline{h}}\right)^{3}+\frac{1}{8}\left(\frac{z}{\underline{h}}\right)^{4}}\left(\frac{z}{\underline{h}}- 1\right)^{n_{T}}+1\vphantom{\left(\frac{z}{\underline{h}}\right)-\frac{1}{2}\left(\frac{z}{\underline{h}}\right)^{3}+\frac{1}{8}\left(\frac{z}{\underline{h}}\right)^{4}}\right]. \end{align}

Here $a_{L}$ and $a_{T}$ are the superposition constants determined by the local flow rate per unit width and the friction coefficient. The polynomial component is identical to (3.10), so the profile is continuous at $Re_{L}=Re_{*}$ or, equivalently, $a_{L}=5/2, a_{T}=0$.

The velocity profile suggested in (3.16) should satisfy the following two conditions. First, the integration of the horizontal velocity from the bottom to the free surface $(({1}/{\underline {h}})\int _{0}^{\underline {h}}\tilde {u}_{i}\,\textrm {d}z)$ should yield the depth-averaged velocity $\bar {\tilde {u}}_{i}$; second, the bottom slope $\partial \tilde {u}_{i}/\partial z\rvert _{z=0}$ should match the friction coefficient $C_{f}$ given in (3.15). These constraints lead to the following expressions for the constants $a_{L}$ and $a_{T}$:

(3.17a)$$\begin{gather} a_{L}=\frac{5\left(n_{T}+1-\dfrac{1}{2}Re_{L}C_{f}\right)}{2n_{T}-3}, \end{gather}$$

(3.17b)$$\begin{gather}a_{T}=\frac{(n_{T}+1)\left(\dfrac{1}{2}Re_{L}C_{f}-5\right)}{n_{T}(2n_{T}-3)}. \end{gather}$$

As the local Reynolds number $Re_{L}$ continues to increase beyond $Re_{*}$, the profile gradually approaches the power-law profile. The velocity should increase monotonically along the $z$ direction, and the maximum velocity should be attained at $z = \underline {h}$. These conditions impose a physical constraint that determines the valid ranges of $a_{T}$ and $a_{L}$, thereby bounding the admissible limit of $Re_{L}$ based on (3.15) and (3.17). Moreover, choosing a larger value for $n_{T}$ tends to increase the maximum admissible value of $Re_{L}$. Thus, $n_{T} = 21$ is adopted because the range of $Re_{L}$ spans $O(10^{3})$ in the present work. More details on the choice of $n_{T}$ and the valid range of $Re_{L}$ have been described by Mendez et al. (2021). The resulting velocity profile applied in this work is presented in figure 6(b). From $k_{t}\bar {\tilde {u}}_{i} \,\bar {\tilde {u}}_{j} = \overline {\tilde {u}_{i} \tilde {u}_{j}}$, the advection constant $k_{t}$ in (3.14b) is expressed as

(3.18)\begin{equation} k_{t}=\frac{62}{315}a_{L}^{2}+\frac{9107}{22880}a_{L}a_{T}+\frac{441}{473}a_{T}^{2}. \end{equation}

3.3.3. Sub-depth scale turbulence model

In the field of shallow-water turbulence, the eddy scale of turbulent structures is typically divided into two ranges: three-dimensional structures of sub-depth scale and horizontal structures of large scale (Nadaoka & Yagi 1998). One intriguing aspect of the coexistence of two-range length scales is the phenomenon of an inverse energy cascade in the spectral domain, also known as backscatter. This involves the transfer of energy from small three-dimensional eddies to large two-dimensional eddies. This backscattering process triggers the development of small-scale instabilities into large-scale horizontal motions (Hinterberger et al. 2007). Given the feature of long waves in film flows, it is reasonable to expect these two-range eddies to play a dominant role in the film flows. The backscatter effect is modelled as a random force field, i.e. stochastic backscatter model, as suggested by Hinterberger et al. (2007). Therefore, this study incorporates a random force vector term into the closure process for the sub-depth scale stress $G_{ij}$ in (3.14b), which accounts for the influence of the inverse cascade.

In the context of DA LES models used in shallow-water turbulence, backscatter modelling is typically conducted for the total stress term (Nadaoka & Yagi 1998; Hinterberger et al. 2007), rather than $G_{ij}$. The total stress $T_{ij} = \overline {\widetilde {u_{i}u_{j}}} - \bar {\tilde {u}}_{i}\, \bar {\tilde {u}}_{j}$ represents the unresolved stress resulting from depth-averaging ($\overline {\tilde {u}_{i}\tilde {u}_{j}} - \bar {\tilde {u}}_{i}\, \bar {\tilde {u}}_{j}$) added to $G_{ij} = \overline {\widetilde {u_{i}u_{j}}} - \overline {\tilde {u}_{i}\tilde {u}_{j}}$; note that $G_{ij}$ arises solely from the filtering operator. In many turbulent shallow-water models, where the time-averaged velocity is approximately constant along the vertical direction, except near solid surfaces, the effect of the depth-averaging operator ($\overline {\tilde {u}_{i}\tilde {u}_{j}} - \bar {\tilde {u}}_{i}\, \bar {\tilde {u}}_{j}$) is often neglected. However, in the present study we consider its influence through the advection closure constant $(k_{t}-1)\bar {\tilde {u}}_{i}\, \bar {\tilde {u}}_{j}(=\overline {\tilde {u}_{i}\tilde {u}_{j}} - \bar {\tilde {u}}_{i}\, \bar {\tilde {u}}_{j})$, as the effect of the velocity profile is deemed significant in rotating films. Therefore, the term $G_{ij} = T_{ij} - (\overline {\tilde {u}_{i}\tilde {u}_{j}} - \bar {\tilde {u}}_{i}\, \bar {\tilde {u}}_{j})$ arises in (3.14b) instead of $T_{ij}$.

First, we describe the decomposition of the total stress $T_{ij}$. The Favre filter is assumed to be separable to enable $\bar {\tilde {\phi }} \approx \tilde {\bar {\phi }}$. Then, the approximation

(3.19)\begin{equation} T_{ij}(=G_{ij} + \overline{\tilde{u}_{i}\tilde{u}_{j}} - \bar{\tilde{u}}_{i}\,\bar{\tilde{u}}_{j})=\overline{\widetilde{u_{i}u_{j}}} - \bar{\tilde{u}}_{i}\,\bar{\tilde{u}}_{j}\approx \widetilde{\overline{u_{i}u_{j}}} - \tilde{\bar{u}}_{i}\,\tilde{\bar{u}}_{j}\qquad (i,j=x,y) \end{equation}

holds, which is exact on the bottom of the disk and the free surface. Equation (3.19) is then decomposed into

(3.20)\begin{equation} (\widetilde{\overline{u_{i}u_{j}}} - \widetilde{\bar{u}_{i}\,\bar{u}_{j}})+(\widetilde{\bar{u}_{i}\,\bar{u}_{j}}-\tilde{\bar{u}}_{i}\,\tilde{\bar{u}}_{j}) =\tilde{D}_{ij}(\boldsymbol{u})+(\widetilde{\bar{u}_{i}\,\bar{u}_{j}}- \tilde{\bar{u}}_{i}\,\tilde{\bar{u}}_{j}), \end{equation}

where the second term $(\widetilde {\bar {u}_{i}\,\bar {u}_{j}}- \tilde {\bar {u}}_{i}\,\tilde {\bar {u}}_{j})$ corresponds to the sub-grid scale stress induced by filtering depth-averaged velocity components, which is typically modelled using eddy viscosity in studies on conventional two-dimensional turbulence (Nadaoka & Yagi 1998; Hinterberger et al. 2007). The effect of three-dimensional sub-grid scale turbulence is included in the first term $\tilde {D}_{ij}(\boldsymbol {u})$, where $D_{ij}(\boldsymbol {u})$ is defined as

(3.21)\begin{equation} D_{ij}(\boldsymbol{u})=\frac{1}{\underline{h}}\int_{0}^{\underline{h}}(u_{i}-\bar{u}_{i})(u_{j}-\bar{u}_{j})\,{\rm d}z. \end{equation}

To further decompose the sub-grid scale stress $\tilde {D}_{ij}(\boldsymbol {u})$ in (3.20), we consider splitting the instantaneous velocity $\boldsymbol {u}$ into the time-averaged component $\langle \boldsymbol {u} \rangle$ and the fluctuating component $\boldsymbol {u}'$, where the notation $\langle \boldsymbol {\cdot } \rangle$ represents the time-averaging operator. Then, from (3.21), the decomposition is

(3.22)\begin{equation} \langle D_{ij}(\boldsymbol{u}) \rangle=D_{ij}(\langle\boldsymbol{u}\rangle)+\langle\overline{u_{i}'u_{j}'}\rangle-\langle\overline{u_{i}'}\,\overline{u_{j}'}\rangle=D_{ij}(\langle\boldsymbol{u}\rangle)+\langle D_{ij}(\boldsymbol{u}')\rangle. \end{equation}

Here, $\langle \overline {u'_{i}u'_{j}}\rangle$ represents the depth-averaged Reynolds stresses of all turbulent structures, and $\langle \overline {u_{i}'}\,\overline {u_{j}'}\rangle$ is the amount of stress caused by the horizontal flow structures of all scales. Thus, the difference between them, namely $\langle D_{ij}(\boldsymbol {u}')\rangle$, indicates the Reynolds stresses induced by three-dimensional fluctuations only. Given that Favre filtering separates the flow structures into fluctuations and large-scale resolved structures, we find the decomposition presented in (3.22) to be approximately valid for $\tilde {D}_{ij}(\boldsymbol {u})$:

(3.23)\begin{equation} \tilde{D}_{ij}(\boldsymbol{u})\approx D_{ij}(\tilde{\boldsymbol{u}}) + \tilde{D}_{ij}(\boldsymbol{u}'). \end{equation}

The first term $D_{ij}(\tilde {\boldsymbol {u}})$ on the right-hand side of (3.23) is equivalent to the stress induced by depth-averaging ($\overline {\tilde {u}_{i}\tilde {u}_{j}} - \bar {\tilde {u}}_{i}\,\bar {\tilde {u}}_{j}$) on the left-hand side of (3.19). Thus, from (3.20) and (3.23), $G_{ij}$ can be approximated as

(3.24)\begin{align} G_{ij} &\approx \tilde{D}_{ij}(\boldsymbol{u}') + (\widetilde{\bar{u}_{i}\,\bar{u}_{j}}- \tilde{\bar{u}}_{i}\,\tilde{\bar{u}}_{j})\nonumber\\ &\approx \tilde{D}_{ij}(\boldsymbol{u}') - \nu_{t}2\bar{\tilde{S}}_{ij}. \end{align}

The first term on the right-hand side of (3.24), $\tilde {D}_{ij}(\boldsymbol {u}')$, represents the sub-grid scale stress generated by three-dimensional turbulent structures near the disk surface (see figure 5b). The loss of information during depth-averaging and Favre filtering means that we must construct $\tilde {D}_{ij}(\boldsymbol {u}')$. In this study a stochastic backscatter model is applied to incorporate a random two-dimensional force vector field $F_{i,bsm}$, which is based on DA LES for shallow-water turbulence (Hinterberger et al. 2007). The second term, $(\widetilde {\bar {u}_{i}\,\bar {u}_{j}} - \tilde {\bar {u}}_{i}\,\tilde {\bar {u}}_{j})$, represents the unresolved stress typically used in LES, where the velocity components are replaced with the depth-averaged velocity and modelled using eddy viscosity, i.e. $-\nu _{t}2\bar{\tilde{S}}_{ij}$; the strain rate $S_{ij}$ is $\tfrac {1}{2}(\partial u_{i} / \partial x_{j} + \partial u_{j}/\partial x_{i})$. The eddy viscosity is commonly represented as $\nu _{t}=C\underline {h}\tilde {u}_{\tau }$, where $u_{\tau }$ denotes the friction velocity and $C$ is an empirical constant. The multiplication of the eddy viscosity and the strain rate, $\nu _{t}\bar {\tilde {S}}_{ij}$, is of $O(\epsilon ^{2})$ when differentiated with respect to the horizontal length scale, as discussed earlier for the long-wave formulation. Thus, the sub-grid scale stress term in (3.14) can be simplified to a random backscatter as follows:

(3.25)\begin{equation} \frac{\partial}{\partial x_{j}}(G_{ij}) ={-} F_{i,bsm} \quad (i,j=x,y). \end{equation}

The backscatter forcing term is linked to local turbulence production, which characterizes the dynamics of internal turbulent structures, and it is challenging to describe the backscattering of three-dimensional turbulence in terms of the long-wave scaling approach. Therefore, modelling the backscatter term requires expressions using dimensional parameters. We consider the dimensional $\hat {F}_{i,bsm} = \hat {l}\hat {\varOmega }^{2} F_{i,bsm}$, where $\hat {l}$ represents the horizontal length scale ($\hat {l}=(9\hat {Q}^{2}/ 4{\rm \pi} ^{2}\hat {\nu }\hat {\varOmega })^{1/4}$). The random force vector defined over the two-dimensional Cartesian coordinate system is

(3.26)\begin{equation} \hat{\boldsymbol{F}}_{bsm} = [\hat{\boldsymbol{\nabla}} \times (\hat{F}_{rms}Z \boldsymbol{e}_{z})] \hat{\varDelta} , \end{equation}

where the scalar value $\hat {F}_{rms}$ represents the root mean square of the backscatter force. To generate the backscatter forcing term, we introduce a filtered random scalar field $Z$ with zero mean and a normal distribution. A random force vector on the two-dimensional domain is constructed by the operator $\hat {\boldsymbol {\nabla }} = (\partial / \partial \hat {x}, \partial / \partial \hat {y})$. To ensure a divergence-free random vector field, the curl operator is applied to a vertically defined scalar field $\hat {F}_{rms}Z \boldsymbol {e}_{z}$. The solenoidal field introduced in (3.26), which pertains to the two-dimensional domain, is equivalent to the random vector potential field method used in three-dimensional backscatter models (Leith 1990; Schumann 1995). The reference grid length scale $\hat {\varDelta }$ in (3.26) is obtained by averaging the grid sizes of the entire domain.

The scaling factor $\hat {F}_{rms}$ in (3.26) is derived from the two-dimensional production of turbulent kinetic energy, $\hat {P}_{2D}$: $\hat {P}_{2D} \sim \hat {F}_{rms}^{2} {\cdot } \Delta \hat {t}$, where $\Delta \hat {t}$ is the applied time step (Alvelius 1999). Here $\hat {P}_{2D}$ scales as

(3.27)\begin{equation} \hat{P}_{2D} \sim \frac{\hat{P}_{3D}}{Re_{\tau}}=\frac{\rvert \bar{\tilde{\hat{\boldsymbol{u}}}} \rvert^{2}\hat{\nu}{C_{f}}^{1/2}}{\hat{\underline{h}}^{2}}, \end{equation}

where $\hat {P}_{3D}$ represents the three-dimensional turbulence production induced from bottom friction, and $Re_{\tau }=\hat {u}_{\tau }\hat {\underline {h}}/\hat {\nu }$ is the Reynolds number defined in terms of the friction velocity $\hat {u}_{\tau } = \rvert \bar {\tilde {\hat {\boldsymbol{u}}}} \rvert (C_{f}/2 )^{1/2}$. Hence, $\hat {F}_{rms}$ is updated every time step as

(3.28)\begin{equation} \hat{F}_{rms}=c_{B}\left(\frac{\hat{P}_{2D}}{\Delta \hat{t}}\right)^{1/2}=c_{B}\frac{\rvert\bar{\tilde{\hat{\boldsymbol{u}}}}\rvert C_{f}^{1/4}}{\hat{\underline{h}}}\left(\frac{\hat{\nu}}{\Delta \hat{t}}\right)^{1/2}, \end{equation}

where the empirical backscatter constant $c_{B} = 2.0$ is determined by comparing numerical results with experimental results obtained from sensor measurements. The effects of $c_{B}$ on the wave dynamics and local thickness of the liquid film are reported in § 5.2, along with a justification for the choice of $c_{B} = 2.0$.

4. Numerical methods

This section presents numerical methods to solve the IBL equations established in § 3. The implementation of the two-regime governing equations in a two-dimensional domain, as well as the modelling of the fluctuations induced by the vertically impinging jet, is explained. The length and time scales of fluctuations induced by the jet impingement and the backscatter force $\hat {F}_{i,bsm}$ are not in accordance with the IBL formulations based on the long-wave approximation. Therefore, we use dimensional variables to integrate the governing equations and the impinging-jet model into a single numerical framework. The numerical set-up is carefully designed to capture wave propagation and regime transition accurately, ensuring that the obtained film thickness and depth-averaged velocity fields produce reliable quantitative analyses.

4.1. Solver for IBL equations

Equations (3.7) and (3.14) for the laminar and turbulent regimes, respectively, are combined into a single set of equations. The combined equations are mathematically expressed in Cartesian tensor notation as follows:

(4.1a)$$\begin{gather} \frac{\partial \hat{h}}{\partial \hat{t}}+\frac{\partial}{\partial \hat{x}_{i}}(\hat{h} \bar{\hat{u}}_{i}) =0, \end{gather}$$

(4.1b)$$\begin{gather} \frac{\partial (\hat{h}\bar{\hat{u}}_{i})}{\partial \hat{t}} + \mathcal{K}\frac{\partial (\hat{h}\bar{\hat{u}}_{i}\,\bar{\hat{u}}_{j})}{\partial \hat{x}_{j}} & ={-}\frac{1}{\hat{\rho}}\frac{\partial(\hat{h} \bar{\hat{p}})}{\partial \hat{x}_{i}} - \frac{\hat{\sigma}\hat{\kappa}}{\hat{\rho}}\frac{\partial \hat{h}}{\partial \hat{x}_{i}} + \mathcal{L}\hat{h}\hat{F}_{i,bsm} \nonumber\\ & -\ \frac{1}{2}C_{f}\bar{\hat{u}}_{i}\left( \bar{\hat{u}}_{i}^{2} + \bar{\hat{u}}_{j}^{2} \right)^{1/2} + B_{i}, \end{gather}$$

(4.1c)$$\begin{gather} \bar{\hat{p}} ={-}\hat{\sigma}\hat{\kappa} + \frac{\hat{\rho}\hat{g}\hat{h}}{2} \quad (i,j = x,y). \end{gather}$$

In (4.1) the depth-averaged velocity $\bar {\hat {u}}_{i}$ and film thickness $\hat {h}$ are considered as resolved (filtered) values in the turbulent regime ($\bar {\hat {u}}_{i} \rightarrow \bar {\tilde {\hat {u}}}_{i}, \hat {h} \rightarrow \hat {\underline {h}}$), whereas they are equivalent to the original (unfiltered) values in the laminar regime. The backscatter force term $\hat {F}_{i,bsm}$ solely represents the sub-grid scale stress.

The closure models discussed in §§ 3.2 and 3.3 are combined into the constants $\mathcal {K}$, $\mathcal {L}$ and $C_{f}$. The values of these constants are determined based on the local Reynolds number $Re_{L}$ in each cell at each time step:

(4.2)\begin{equation} \left. \begin{array}{@{}ll@{}} \mathcal{L}=0,\mathcal{K} = k_{l},C_{f} = 5/Re_{L}, & Re_{L} \leq Re_{*}, \\ \mathcal{L}=1,\mathcal{K} = k_{t},C_{f} = 5/(Re_{L}^{1/4}Re_{*}^{3/4}), & Re_{L} > Re_{*}. \end{array} \right\} \end{equation}

By substituting the constants for the laminar regime ($Re_{L} \leq Re_{*}$) into (4.1), we obtain the Kapitza–Shkadov equations featuring the quartic velocity profiles, as suggested in § 3.2. When the constants for $Re_{L} > Re_{*}$ are used, the governing equations are transformed into the shallow-water DA LES equations discussed in § 3.3.

4.2. Computational domain and discretization

To solve (4.1), the implicit and sequential finite area solver proposed by Rauter & Tuković (2018) is employed. This finite area solver is known for accurately capturing the depth-averaged velocity and film thickness of two-dimensional flows. The fluid domain has a two-dimensional disk shape with a radius of $\hat {r} = 0.15$ m, which is the same as in the experiments. The domain is discretized using structured grids, referred to as the butterfly geometry (figure 7a). The length scale of the grids, denoted as $\hat {\varDelta }$, is determined based on the local maximum of the time-averaged film thickness, $\hat {h}_{hump}$ (figure 9). Following grid convergence tests, the grid length is set as $\hat {\varDelta }\approx 0.55\hat {h}_{hump}$. The grid convergence tests were performed by comparing the standard deviation of the time series of film thickness at a single point $\hat {r}=43$ mm, which was extracted from each simulation case without inlet fluctuations over the period $\hat {t}=0.05{-}0.25\,\textrm {s}$. The convergence with various grid resolutions is presented in figure 7(b).

Figure 7. (a) Grid layout of a two-dimensional fluid domain. Grey cells in the lower inset indicate the region subjected to the jet impingement model. (b) Result of grid convergence test: standard deviation of film thickness from $\hat {t}=0.05{-}0.25\,\textrm {s}$ versus normalized grid size $\hat {\varDelta }/\hat {h}_{hump}$. Simulations were conducted without inlet fluctuations for $\hat {\varOmega }=52.4\,\textrm {rad}\,\textrm {s}^{-1}$ and $\hat {Q}=12.5\,\textrm {mL}\,\textrm {s}^{-1}$.

In the finite area method, the discretization of each term in (4.1) is accomplished by integrating them over a control area. This integration process allows us to express each term with regard to the values at the cell centres, values on the cell edges and vectors normal to the edges, using the second-order midpoint rule and Gauss’ theorem. As for the spatial interpolation to represent values on the edges of each cell, a local blend between linear and upwind interpolations is performed. For scalar fields, a normalized variable diagram scheme called Gamma interpolation is used. The interpolation constant relevant to stability and diffusivity is chosen to be $\beta = 1/5$. The interpolation of vector fields is performed in the local edge-based coordinate system (Tuković & Jasak 2012). The interpolation factor is likewise calculated using the Gamma scheme.

Equation (4.1) is solved in a time-marching manner. Discretization of the temporal derivative at the $n$th time step is conducted with the implicit second-order scheme according to

(4.3)\begin{equation} \left(\frac{\partial \phi}{\partial t}\right)^{n} \approx \frac{3 \phi^{n}-4\phi^{n-1}+\phi^{n-2}}{2 \Delta t}, \end{equation}

where the superscripts $n-1$ and $n-2$ denote the values of the previous two time steps, and the superscript $n$ indicates the implicit value of the present time step. Time step control is implemented based on the Courant–Friedrichs–Lewy (CFL) number $C = c_e \Delta t / \varDelta _e$, where $c_e = \max (\boldsymbol {m}_e \boldsymbol {\cdot } \bar {\boldsymbol {u}}_e)$ represents the maximum depth-averaged velocity along the edges of the cell and $\varDelta _e$ corresponds to the square root of the cell area. The subscript $e$ denotes edge values and $\boldsymbol {m}_{e}$ is the unit normal vector outward from the edge. The time step is adjusted to ensure that the maximum CFL number satisfies the criterion $C < 0.4$.

The set of coupled nonlinear differential equations is solved numerically using a sequential approach. Following the methodology of Rauter & Tuković (2018), the depth-averaged pressure $\bar {\hat {p}}$ is first updated with the film thickness of the previous time step using (4.1c). The momentum conservation equation (4.1b) is then solved with the updated pressure and the old values of the film thickness to obtain $\bar {\hat {u}}_{i}$. Finally, the continuity equation (4.1a) is solved for $\hat {h}$. This iterative algorithm is repeated for a given time step until the initial residual of the depth-averaged velocity reaches $10^{-6}$. The initial conditions are $\hat {h}=20\,\mathrm {\mu }\textrm {m}$ and $\bar {\hat {\boldsymbol {u}}}=0$ for every cell in the simulation domain. A von Neumann boundary condition is employed at the outlet boundary, which corresponds to the disk edge. The normal gradient of each variable $\phi$ is zero at the boundary, i.e. $\boldsymbol {\nabla } \phi (\boldsymbol {x})=0$ for all $\boldsymbol {x}$ belonging to the boundary of the domain.

4.3. Analytical modelling of inflow: impingement of fluctuating jet

The complex nature of the film flow formed by a vertically impinging jet makes it challenging to treat the impinging jet as a fixed inlet boundary. This is particularly true in the region where the boundary layer merges with the free surface ($\hat {r} < \hat {r}_{0}$ in figure 1a) and in the presence of concentric waves near the impinging jet. These waves originate from disturbances in the impinging jet prior to impingement, as discussed in § 2. To account for the effects of the impinging jet in the simulations, the local film thickness and depth-averaged velocity obtained from an analytical model are assigned to every cell in $\hat {r} < \hat {r}_{0}$ (grey region in figure 7a) at every time step.

The impingement of a circular liquid jet on a rotating disk has been extensively studied in the context of hydraulic jumps (Wang & Khayat 2018; Wang et al. 2020; Ipatova, Smirnov & Mogilevskiy 2021) and heat transfer applications (Rice, Faghri & Cetegen 2005; Lienhard 2006). In most of these studies, a low-flow-rate condition is assumed, considering a steady and axisymmetric jet. However, the present study focuses on modelling the film flow resulting from the impingement of a high-flow-rate jet, whereby large surface fluctuations are induced by shear instability. To capture the effects of temporal fluctuations of the impinging jet on the formation of concentric waves, the existing theory for the steady impinging jet is modified by introducing temporal fluctuations. The influence of the fluctuations on the time-averaged flow properties is considered to be minimal, as demonstrated in the work of Bhagat & Wilson (2016).

The first step in modelling the liquid film thickness produced by jet impingement is to introduce an analytical model under steady flow assumptions. For axisymmetric and steady flow conditions, the film thickness and the size of the boundary merging zone ($r < r_{0}$) were formulated by Wang & Khayat (2018) as

(4.4a)$$\begin{gather} \hat{h}(\hat{r}<\hat{r}_{0}) = \frac{\hat{d}}{8}\left[\frac{\hat{d}}{\hat{r}} + \left(\frac{210{\rm \pi}\hat{r}\hat{\nu}}{13\hat{Q}}\right)^{1/2} \right], \end{gather}$$

(4.4b)$$\begin{gather}\hat{r}_{0}=\frac{\hat{d}}{2}\left(\frac{156\hat{Q}}{875{\rm \pi}\hat{\nu}\hat{d}}\right)^{1/3}. \end{gather}$$

Although the analytic solutions in (4.4) were derived under the assumption of $\hat {\varOmega }=0$, they remain valid when the ratio of disk angular velocity to nozzle flow rate is small. The difference of the formula (4.4a) from the numerical solution in non-zero $\hat {\varOmega }$ condition presented by Wang & Khayat (2018) is less than 5 % at $\hat {r}=\hat {r}_{0}$ for $\hat {Q} = 16.7\,\textrm {mL}\,\textrm {s}^{-1}$ and $\hat {\varOmega } = 52.4\,\textrm {rad}\,\textrm {s}^{-1}$.

Previous studies on shear-induced fluctuations on the surface of circular liquid jets have shown that axisymmetric modes dominate over asymmetric modes when the wavenumber of the fluctuations is small ($\hat {k}\hat {d}<2$) (Yang 1992; Shi et al. 1999). This dominance of axisymmetric modes is also evident in our visualizations of the film flow formed in the jet impingement zone (figure 4). Thus, the azimuthal fluctuations on the impinging jet are neglected here. Additionally, the boundary layer near the disk surface is known to maintain a cubic velocity profile, as in laminar jet impingement, even in the presence of a turbulent and fluctuating free-stream flow (Lienhard et al. 1992). Accordingly, we introduce a temporal axisymmetric fluctuation $\mathcal {G}(\hat {t})$ multiplied by the nozzle diameter $\hat {d}$ to the steady theoretical solutions (4.4) to model the axial fluctuations of the circular liquid jet in the impingement process. The flow conditions for the region $\hat {r} < \hat {r}_{0}$ are then given by

(4.5a)$$\begin{gather} \hat{h}(\hat{r}<\hat{r}_{0}) = \frac{\mathcal{G}(\hat{t})\hat{d}}{8}\left[\frac{\mathcal{G}(\hat{t})\hat{d}}{\hat{r}} + \left(\frac{210{\rm \pi}\hat{r}\hat{\nu}}{13\hat{Q}}\right)^{1/2}\right], \end{gather}$$

(4.5b)$$\begin{gather}\hat{r}_{0} = \frac{\mathcal{G}(\hat{t})\hat{d}}{2}\left(\frac{156\hat{Q}}{875{\rm \pi}\hat{\nu}\mathcal{G}(\hat{t})\hat{d}}\right)^{1/3}, \end{gather}$$

(4.5c)$$\begin{gather}\bar{\hat{\boldsymbol{u}}}=\frac{\hat{Q}}{2{\rm \pi}\hat{r}\hat{h}}\boldsymbol{e}_{r}-\frac{14}{17}\hat{r}\hat{\varOmega}\boldsymbol{e}_{\theta}. \end{gather}$$

The temporal fluctuations $\mathcal {G}(\hat {t})=1 + A\sum _{n=0}^{100}\cos (n \hat {\omega } \hat {t} + \mathcal {Z}_{n})$ follow a widely applied method of adding noise to the film flow inlet (Chang, Demekhin & Saprikin 2002). In this expression, $\hat {\omega }$ represents the frequency unit, typically chosen to be 1/50 of the neutral frequency of the main instability. The term $\mathcal {Z}_{n}$ denotes a random phase difference that is uniformly distributed in the range $[0, 2{\rm \pi} ]$. The constant $A$ represents the magnitude of the random fluctuations. Equation (4.5) provides quasi-steady solutions of film flow for an impinging jet with an oscillating diameter.

The specific values of the forcing frequency $\hat {\omega }$ and the fluctuation magnitude $A$ are chosen based on the fluctuations in film thickness obtained through sensor experiments; see § 2.2. For $\hat {\omega }$, the power spectra of the time series in the experimental data do not exhibit a specific peak corresponding to the main instability. Therefore, $\hat {\omega }$ is set to $60\,\textrm {s}^{-1}$, which allows the wavenumber range to cover most of the large-amplitude spectra observed in the experimental results. The value of $A$ significantly affects the downstream behaviour of the liquid film thickness. In this study, $A = 5.0 \times 10^{-3}$ is selected to ensure a match between the standard deviations of the time series obtained from numerical simulations and experimental measurements. Both sets of data are measured at $\hat {r} = 43\,\textrm {mm}$ for a time span of $\hat {t} = 0.05{-}0.25\,\textrm {s}$. Finally, the expression of velocity in (4.5c) for $\hat {r} < \hat {r}_{0}$ is derived following the approach of Kim & Kim (2009). This expression represents an asymptotic solution of the laminar IBL equation in the region $\hat {r} \ll 1$ for given values of $\hat {h}$, $\hat {Q}$ and $\hat {\varOmega }$.

5. Results and discussion

5.1. Distribution of surface waves

The film flow computed by our IBL model is shown in figure 8 for disk angular velocity $\hat {\varOmega }=52.36\,\textrm {rad}\,\textrm {s}^{-1}$ and nozzle flow rate $\hat {Q}=12.5\,\textrm {mL}\,\textrm {s}^{-1}$. The waves in the numerical results are represented by the contours of dimensionless film thickness $h$ (figure 8b), and they are qualitatively compared with those of the visualization image (figure 8a); see supplementary movie available at https://doi.org/10.1017/jfm.2024.274. One of the most pronounced features is the generation of concentric waves upstream and their transition into three-dimensional solitary waves. During the initial stage of this transition, three-dimensional wave structures of various scales emerge simultaneously. These structures later converge into a group of large-scale $\varLambda$ solitons as they propagate outwards. As the three-dimensional waves propagate towards the edge of the disk, their peak amplitude decreases and they become horizontally elongated. Another intriguing feature is the generation of small-scale $\varLambda$ solitons amidst the elongated three-dimensional waves near the edge of the disk. These observations from the numerical simulations align well with the visualization results in figure 8(a), demonstrating the successful capture of wave behaviours by the present IBL model.

Figure 8. (a) Visualization of free-surface waves and (b–d) numerical results at a specific instant: (b) contours of dimensionless film thickness $h$, (c) distribution of film thickness $h$ along a straight line from disk centre ($r=0{-}4.0$) and (d) depth-averaged horizontal velocities. In panel (d) the red and blue lines denote the radial $\bar {u}_{r}$ and tangential $\bar {u}_{\theta }$ components, respectively, and the black dashed line is the magnitude of horizontal velocity, $\rvert \bar {\boldsymbol {u}}\rvert$. All data were acquired with $\hat {\varOmega }=52.4\,\textrm {rad}\,\textrm {s}^{-1}$ and $\hat {Q}=12.5\,\textrm {mL}\,\textrm {s}^{-1}$. See supplementary movie for panels (a,b).

The regional mode transition of surface waves results in variations in the local height of the waves (figure 8c). In the concentric-wave regime ($r < 1.0$), the amplitude of the fluctuations in surface elevation is smaller than that in the downstream three-dimensional waves. Additionally, the phases of the instantaneous depth-averaged velocities and film thickness fluctuation are not strongly coupled (figure 8d), indicating that the flow in the concentric-wave regime is primarily influenced by random motions induced by the backscattering of sub-depth scale turbulence. As the three-dimensional waves develop in $r = 1.5{-}3.0$, the height of the wave crest increases, and the fluctuation phases observed in the velocity components and film thickness coincide in this region because the flow is no longer dominated by backscattering as the local flow rate enters the laminar condition. Finally, a reduction in the peak height is evident in the downstream region ($r = 3.0{-}4.0$), as shown in figure 8(c). The trends observed under different flow conditions ($\hat {\varOmega }$ and $\hat {Q}$) generally match those presented in figure 8.

To further validate our numerical simulations in a quantitative manner, the film thickness is compared with previous experimental results and the data acquired from our confocal chromatic sensor. As the confocal chromatic sensor used in this study can only measure the vertical displacement of the liquid–gas interface, it provides a time series of the free-surface elevation at a measurement point. Hence, the validation strategy consists of two parts. First, the time-averaged film thickness is compared with data from previous investigations. Second, the characteristics of the instantaneous film thickness fluctuations are used to evaluate the reliability of the numerical simulations in capturing the dynamics of surface waves.

For the time-averaged film thickness $\langle h \rangle$, the instantaneous thickness is averaged over a time span of $t=15.7{-}23.5$ at a given point and plotted along a radial line from the disk centre, $r = 0{-}4.0$. In previous IBL models on axisymmetric waves (Sisoev et al. 2003; Kim & Kim 2009), the radial distribution of $\langle h \rangle$ was assumed to have a universal dimensionless profile, regardless of the disk angular velocity $\hat {\varOmega }$ and nozzle flow rate $\hat {Q}$. However, our numerical results show that the profiles no longer fully collapse onto a single line due to the inclusion of the impinging-jet model, which does not adhere to the long-wave approximation (figure 9a,b). The film thickness $\langle h \rangle$ in the upstream region $r < 1.0$ exhibits variations as $\hat {\varOmega }$ and $\hat {Q}$ increase. The local maximum of $\langle h \rangle$, known as the hump thickness $h_{hump}$, tends to increase with increasing $\hat {\varOmega }$, and its location approaches $r = 0$. The variations in $\langle h \rangle$ with respect to $\hat {Q}$ are relatively minor. The hump thickness displays negligible changes with increasing $\hat {Q}$, while the film thickness between $r = 0$ and $h_{hump}$ decreases. In contrast to the upstream tendencies, the film thickness profiles collapse onto a single line in the downstream region ($r > 1.0$), indicating a reduction in the influence of jet impingement downstream (in the time-averaged sense).

Figure 9. Time-averaged film thickness $\langle h \rangle$ for (a) several disk angular velocities $\hat {\varOmega }$ and (b) nozzle flow rates $\hat {Q}$. Dimensional film thickness and radial position are presented in panels (c,d), which correspond to panels (a,b), respectively.

To clarify the effects of $\hat {Q}$ and $\hat {\varOmega }$ on the hump formation, the variations in the dimensional location and thickness of the humps are examined (figure 9c,d). The dimensional location $\hat {r}_{hump}$ of the hump increases as $\hat {\varOmega }$ decreases and $\hat {Q}$ increases. Simultaneously, the dimensional thickness $\hat {h}_{hump}$ of the hump tends to increase as $\hat {\varOmega }$ reduces, which is similar to the result of Wang et al. (2020). Interestingly, with increasing $\hat {Q}$, $\hat {h}_{hump}$ becomes greater slightly. This result is contrary to the correlation proposed by Wang et al. (2020) who used the jet diameter $\hat {d}$ as the only length scale. Our findings suggest that the long-wave parameters ($\hat {l}$ and $\hat {\delta }$) also influence $\hat {h}_{hump}$, particularly when the jet diameter is small.

In the vicinity of the jet impingement point ($r=0{-}0.3$), $\langle h \rangle$ from our numerical simulation is consistent with the experimental data obtained from set-ups with impinging jets, although $\langle h \rangle$ is sensitive to the inflow configuration (i.e. $\hat {Q}$, $\hat {H}$ and $\hat {d}$). Furthermore, the present numerical results are in excellent agreement with the theoretical model proposed by Bhagat & Wilson (2016) for film flows generated by a circular high-flow-rate liquid jet impinging on a stationary wall (dashed line in figure 10a) in the region $r=0{-}1.0$. As the influence of the impinging jet diminishes downstream ($r > 1.0$), the IBL results are in better agreement with the previous experimental data (figure 10a); the data are also compared in a narrower $y$-axis range in figure 10(b). Notably, the present numerical results match better with the experimental data obtained under large-flow-rate conditions. The present simulation produces lower values than the solutions given by the laminar IBL model (Kim & Kim 2009). This discrepancy arises from the assumption of a power-law velocity profile in the turbulent regime, leading to greater bottom friction.

Figure 10. Time-averaged film thickness $\langle h \rangle$ along radial coordinate: (a) upstream region and (b) downstream region. The red line is our numerical result. The black dashed line is the theoretical model for film flow formed by the impinging circular liquid jet (Bhagat & Wilson 2016), and the black dotted line is the theoretical model for laminar film flow (Kim & Kim 2009). The symbols denote experimental data from previous studies (Miyasaka 1974b; Muzhilko, Rifert & Barabesh 1983; Leneweit et al. 1999; Wang et al. 2020). Our numerical results were obtained with $\hat {\varOmega }=52.4\,\textrm {rad}\,\textrm {s}^{-1}$ and $\hat {Q}=12.5\,\textrm {mL}\,\textrm {s}^{-1}$.

The wave characteristics are now assessed by analysing the temporal fluctuations in the film thickness. The time series of film thickness obtained from the numerical simulations is compared with that of the sensor experiments in figure 11(a,b). As the sensor was fixed above the rotating disk, the location of the sampling point from which the film thickness is extracted in the numerical domain is rotated in the direction opposite to the direction of the disk rotation; note that the reference frame in our numerical simulations rotates together with the disk. In the sensor experiments, the vertical displacement of the free surface was measured at two positions ($\hat {r} = 43$ and 109 mm) instead of film thickness itself. The temporal average was then subtracted from the time-series data so that the experimental values fluctuate around zero. For both upstream ($\hat {r} = 43\, \textrm {mm}$) and downstream ($\hat {r} = 109\, \textrm {mm}$) positions, the numerical results are in good agreement with the experimental data in terms of both the local maxima $h_{max}$ of the instantaneous film thickness and the waveform.

Figure 11. (a,b) Time series of film thickness $h$, (c,d) histograms of the local maxima $h_{max}$ of film thickness from numerical simulations and (e, f) $h_{max}$ histograms from experimental measurements. Panels (a,c,e) are for the upstream point ($\hat {r} = 43\,\textrm {mm}$) and panels (b,d, f) are for the downstream point ($\hat {r} = 109\,\textrm {mm}$).

Regarding the histograms of $h_{max}$ extracted from the time series (figure 11c–f), both the numerical and experimental results exhibit several common features. In the upstream position (figure 11c,e), the histogram shows a single Weibull distribution for both the numerical and experimental data. In the downstream position (figure 11d, f), $h_{max}$ displays a prominent peak in the small-thickness range and a more uniform distribution in the large-thickness range. However, there are also discrepancies between the two sets of data. These deviations are partly attributable to the limited number of numerical data and the effect of baseline filtering on the experimental data (see § 2.2). The histogram of the experimental data includes a larger number of $h_{max}$ data, as real-time sensor experiments allow for significantly longer data acquisition ($\hat {t}=0.2{-}5.0\,\textrm {s}$), whereas data acquisition from the numerical simulations is limited to $\hat {t}=0.2{-}0.7\,\textrm {s}$ by the computational cost. To reduce the discrepancy between the numerical and experimental results, first it is necessary to incorporate more rigorous backscatter models. The use of high-fidelity direct numerical simulations makes it possible to formulate the energy balance that encompasses the backscattering effect in wavy film flows. Thus, more accurate backscatter models similar to those in shallow waters (van Prooijen & Ujittewaal 2009; Klöwer et al. 2018) can be established. In addition, a more sophisticated experimental set-up capable of directly acquiring real-time film thickness data, instead of measuring the displacement of the liquid–gas interface, should be implemented.

In summary, the numerical model proposed in this study produces reliable results at large nozzle flow rates (and high global Reynolds numbers $Re_{G}$). Despite some discrepancies from the experimental data, the overall shape and scale of the surface waves is well aligned with those of the experimental data. Furthermore, the analysis of the time-averaged film thickness provides new insights into the influence of the fluctuating impinging jet and turbulent flow structures under large flow rates.

5.2. Effects of upstream backscatter

The generation and propagation of three-dimensional waves identified in the numerical simulations exhibit unique characteristics in terms of the wave initiation process. In gravity-driven film flows, three-dimensional waves are generated as two-dimensional waves enter a high-velocity regime dominated by transverse instability (Chang et al. 2002). In contrast, the generation of three-dimensional waves in the film flow spreading over a rotating disk is primarily triggered by extensive fluctuations in the concentric-wave regime under a large local flow rate near the impinging zone. These fluctuations result from sub-depth scale turbulence. As a result, the patterns of three-dimensional waves in the onset region exhibit distinct features from those reported for falling films.

The response of surface waves to the backscatter phenomenon is of particular interest because the backscatter model has not previously been included in film flow studies. The backscatter model described in (3.28) includes an empirical constant $c_{B}$ that represents the rate of energy transfer from internal sub-depth scale turbulence to fluctuations of horizontally resolved scale. Although a specific value of $c_{B}$ is suitable for each shallow-water system, we conduct a model test by comparing the results obtained with different values of $c_{B}$, which allows us to investigate the impact of backscattering on the characteristics of three-dimensional waves.

First, the numerical results with backscatter constants of $c_{B} = 2.0$ and $c_{B} = 0$ (i.e. no backscatter) are compared. For the case of $c_B = 0$ (figure 12a), three-dimensional waves are not clearly present in a substantial portion of the region where the wave regime originally transitions from concentric waves ($1.5 < r < 2.0$). Furthermore, as $\varLambda$-soliton waves travel downstream, they do not break up or coalesce with each other until reaching the far downstream region ($r > 3.0$). In the presence of backscatter effects ($c_{B} = 2.0$, figure 12b), concentric waves are azimuthally disturbed much earlier, entering the onset region of three-dimensional waves near $r = 1.5$. This demonstrates that, under a large flow rate, the backscattering of sub-depth scale turbulence in the upstream concentric-wave regime with a high local Reynolds number is responsible for the fluctuations in horizontal flow structures, which in turn induce the generation of unique three-dimensional surface waves downstream.

Figure 12. Contours of film thickness $h$ for (a) backscatter constant $c_{B} = 0$ and (b) $c_{B} = 2.0$. (c) Distribution of $h$ with respect to radial coordinate $r$ for the two cases in panels (a,b).

Compared with $c_B = 2.0$, the amplitude of the three-dimensional waves for $c_B = 0$ decays slowly in the downstream region ($r=3.0{-}4.0$), and the streamwise spacing between two large-scale solitary waves becomes narrower (figure 12c). This result is correlated with the observation that, without backscatter effects, the interaction of waves, including their coalescence, barely occurs in the region $r=1.5{-}3.0$. For $c_B = 0$, there is less likelihood of a loss of wave energy via interaction, leading to the large amplitude and narrow spacing of the waves downstream. The influence of backscatter on horizontal flow structures in long-wave systems has been examined in prior studies on turbulent shallow-water flows. For instance, Hinterberger et al. (2007) showed that vortical coherent structures induced by horizontal shear do not appear without a backscatter model in a mixing layer of two adjacent streams of shallow water with different velocities. Taking $\varLambda$ solitons to be coherent structures in the context of weak and dissipative turbulence (Kalliadasis et al. 2012), our study provides another example of the generation of horizontal coherent structures induced by backscattering.

According to figure 13(a), as the extent of bottom-induced turbulent backscattering increases (i.e. $c_{B}$ increases), the transition from concentric waves to three-dimensional waves occurs closer to the disk centre. Moreover, the horizontal size of the three-dimensional waves in their onset region ($1.5< r<2.0$) becomes smaller, while the downstream waves become larger after undergoing multiple coalescence processes. In the onset region of the three-dimensional waves, the local Reynolds number $Re_{L}$ exceeds the threshold established for gravity-driven films (known to lie in the range $Re_{G,c}=40{-}70$), where a two-dimensional solitary wave becomes unstable to transverse perturbations (Demekhin et al. 2007). Consequently, the dynamics of the waves downstream of this region are highly sensitive to the azimuthal perturbations of the concentric waves entering the transition zone ($1.5< r<2.0$). The backscatter of small-scale turbulent structures near the disk surface is responsible for such azimuthal perturbations, leading to the early generation of three-dimensional waves; a quantitative analysis of the backscatter impact on wave perturbations is beyond the scope of the present work.

Figure 13. Effects of backscatter constants $c_{B}$ on wave dynamics: (a) contours of film thickness $h$; time series of $h$ at (b) upstream ($r = 1.25$) and (c) downstream ($r = 3.18$) locations.

In the time-series data of film thickness numerically obtained at upstream ($r = 1.25$) and downstream ($r = 3.18$) locations (figures 13b and 13c, respectively), the sampling point is fixed on the bottom of the disk surface, which is equivalent to the measurement point rotating with the disk in the experimental domain (unlike the time-series data presented in figure 11). The selection of a fixed sampling point on the disk surface offers better insights into the physical characteristics of turbulent film flows (Sivashinsky & Michelson 1980). At the upstream location, a greater backscatter constant $c_B$ results in an increase in the amplitude of thickness fluctuations (figure 13b). This finding, coupled with the larger horizontal fluctuations depicted in figure 13(a), implies that turbulent kinetic energy from sub-depth scale turbulence is transformed into wave energy in the concentric-wave regime. By contrast, in the downstream region, there is no significant difference in the height of the wave peaks as $c_B$ increases (figure 13c). However, the peak values become more regular with greater values of $c_B$. In the case of $c_B = 0.5$, the amplitudes of the peaks exhibit a high level of randomness, while for $c_B = 4.0$, the wave pattern is closer to that of periodic pulsations.

The high level of spatiotemporal chaos in the film thickness downstream prevents traditional frequency-domain analyses from providing meaningful information. The time-series data in figure 13(c) are similar to those obtained by solving the generalized Kuramoto–Sivashinsky equation with the dispersion constant in the range 0.2–0.4, where no distinguishable peaks are present in the power spectral density. In such cases, the Lorenz return map is a more informative tool (Gotoda, Pradas & Kalliadasis 2015). In the present work, the Lorenz return map is constructed by plotting pairs of successive local maxima $(h_{m,k}, h_{m,k+1})$ in the time series of film thickness at $r = 3.18$ in figure 13(c). The points plotted in figure 14 form four distinct groups. There are two distinct groups with small and large film thickness on the line $h_{m,k} = h_{m,k+1}$. The other two groups consist of pairs of small and large successive peaks and their counterparts. In chaos theory, a signal is regarded as more predictable in the long term when the points are concentrated in these four groups, indicating that a signal is closer to periodic pulsations with stochastic behaviour than a chaotic one. Therefore, figure 14(a–c), obtained from the numerical results over a time span of $\hat {t}=0.2{-}0.7\,\textrm {s}$, suggests that the three-dimensional waves in the downstream region exhibit greater long-term predictability when the backscatter constant $c_B$ is larger.

Figure 14. Lorenz return map obtained from numerical simulations for $c_{B}$ = 1.0, 2.0 and 4.0: $\hat {\varOmega }=31.4\,\textrm {rad}\,\textrm {s}^{-1}$, $\hat {Q}=16.7\,\textrm {mL}\,\textrm {s}^{-1}$.

As mentioned in § 3.3.3, $c_{B} = 2.0$ was chosen for our turbulent film flow because wave patterns for $c_B = 2.0$ are closest to those of the experimental visualization. The transition from concentric to three-dimensional waves, which occurs near $r = 1.5{-}2.0$ in the experiment (figure 8a), matches with the numerical result obtained using $c_{B} = 2.0$. Moreover, the wavelength of azimuthal perturbation in the concentric waves and the horizontal amplitude of the three-dimensional surface waves in this region exhibit good agreement between the experiment and the simulation with $c_{B} = 2.0$.

The value of $c_{B} = 2.0$ adopted in this study is significantly smaller than the values used in studies on shallow-water turbulence where $c_{B}$ is often greater than 50 (Hinterberger et al. 2007; van Prooijen & Ujittewaal 2009). This lower backscatter constant can be attributed to the presence of the free surface, which absorbs turbulent energy during wave propagation. For our disk with a much smaller scale, surface tension plays a crucial role in extracting energy from internal turbulence when making an additional free-surface area (forming surface waves). A deforming liquid–gas interface is known to reduce the intensity of turbulent structures with eddy size comparable to the wavelength of the deforming interface (Chiodi, McCaslin & Desjardins 2019). Because the free surface absorbs a significant portion of turbulent energy through surface deformation, the amount of fluctuation in the horizontal flow structures induced by the turbulent backscattering is smaller than in the shallow-water systems where surface tension is negligible, which contributes to the smaller value of $c_B$.

5.3. Interaction of three-dimensional waves

The interactions among three-dimensional surface waves after the transition from concentric waves are now examined in detail. The propagation of surface waves is first analysed using the spatiotemporal evolution of film thickness $h$ along a straight line from disk centre to disk edge over a time span of $t-t_0=0{-}4.0$; $t_{0}=10.5$ is a time at which the flow is fully developed. In figure 15 the phase velocity of the waves is represented by the slope of the bright ridges, corresponding to the crests of each wave. In the concentric-wave regime ($r < 1.0$), the phase velocity of the waves appears relatively uniform. However, in the onset region of the three-dimensional waves ($r = 1.5{-}2.0$), a wide range of phase velocities and wave amplitudes are observed. This non-uniform distribution differs from the onset of the three-dimensional regime in falling films, where the dominant wavelength of transverse disturbances in the two-dimensional waves determines the scale of relatively uniform $\varLambda$ solitons (Kalliadasis et al. 2012).

Figure 15. Spatiotemporal evolution of surface waves along a straight line from disk centre to disk edge ($r=0{-}5.2$).

Solitary waves in film flows tend to travel faster when they have a larger scale (Malamataris, Vlachogiannis & Bontozoglou 2002). As a result, two waves with distinct differences in scale can merge to form a single large wave. The three-dimensional waves in the onset region, which initially possess a wide range of scales, undergo frequent coalescence with each other as they travel downstream. This process is illustrated in figure 16(a–c), where a smaller-scale soliton ‘A’ merges with the larger wave ‘B’. After multiple coalescence events in the region $r = 1.5{-}2.0$, the scales and phase velocities of the waves converge in the region $r = 2.0{-}3.5$ (figure 15). Despite the deceleration and shallowing of the base film flow, the phase velocity and peak amplitude of waves in this downstream region remain significantly greater than in the upstream region. Farther downstream ($r > 3.5$), the amplitudes of the three-dimensional waves decrease while their phase velocities remain in a similar range (figure 15). The wave interactions in this region are unique in that the large waves no longer directly interact with each other, unlike for $r < 3.5$. Instead, small-scale $\varLambda$ solitons are constantly separated from the preceding large waves and coalesce into the following ones (figure 16d–f); wave front ‘B’ separates from ‘A’ and subsequently merges with ‘C’.

Figure 16. Sequential snapshots to illustrate interactions between three-dimensional waves: (a–c) close-up of the region $r \approx 1.0{-}2.0$ and (d–f) close-up of the region $r \approx 4.0{-}5.0$ at identical instants. The time step for the three panels in the upper and lower rows is $\Delta t = 0.42$.

Because the local Reynolds number strongly influences the generation and interaction mechanisms of three-dimensional waves, correlating the amplitude of the waves with the long-wave parameter $\hat {\delta } = (\hat {\nu }/\hat {\varOmega })^{1/2}$, which represents the characteristic scale for film thickness in (3.1), is not preferable. To characterize the variations in wave amplitude, a parameter other than $\hat {\delta }$ should be considered. In the study on three-dimensional waves in falling films (Demekhin et al. 2007), the amplitude and phase velocity of a single $\varLambda$ soliton were found to be governed by a dimensionless parameter called the generalized Reynolds number $R_{g} = Re_{G}^{11/9}/(3^{7/9}5\gamma ^{1/3})$, where $\gamma = \hat {\sigma }/(\hat {\rho }\hat {\nu }^{4/3}\hat {g}^{1/3})$ is the Kapitza number for gravity-driven films. This parameter includes the flow rate of the liquid film on which the $\varLambda$ soliton lies as well as the body force acting on it. The concept of the generalized Reynolds number is introduced to combine the Kapitza number $\gamma$, which represents the ratio of surface tension force to inertial force, with the Reynolds number $Re_{G}$, which characterizes the flow rate per unit width, in a single parameter. As $\varLambda$ solitons are regarded as the fundamental coherent structures of three-dimensional waves in film flows, we propose a generalized Reynolds number for rotating films to characterize the amplitude of three-dimensional waves.

By fixing the point of observation on each three-dimensional wave following the approach of Demekhin et al. (2010), the difference between rotating films and gravity-driven films is reduced to the influence of the body force, where gravitational acceleration $\hat {g}$ in the original definition of the generalized Reynolds number $R_g$ is replaced by centrifugal acceleration $\hat {r}\hat {\varOmega }^{2}$. Regarding the Reynolds number, $Re_{G}$ in the definition of $R_g$ is replaced by the time-averaged local Reynolds number $\langle Re_{L}\rangle$. The adoption of the local Reynolds number is justified, as the point of observation is fixed on a specific three-dimensional wave travelling in a particular region on the rotating film. The generalized Reynolds number for rotating films is then formulated as follows:

(5.1)\begin{equation} R_{g} = \frac{\langle Re_{L}\rangle^{11/9}}{3^{7/9} 5 \gamma^{1/3}}, \quad \text{where} \ \gamma = \frac{\hat{\sigma}}{\hat{\rho}\hat{\nu}^{4/3}(\hat{r}\hat{\varOmega}^{2})^{1/3}}. \end{equation}

The Kapitza number $\gamma$ dependent on $\hat {r}$ is specifically defined to represent the local wave amplitude under the influence of a spatially non-uniform body force, and it differs from that used by Sisoev et al. (2003) for rotating films, which was based on liquid properties ($\hat {\rho }, \hat {\nu }$ and $\hat {\sigma }$) and $\hat {\varOmega }$.

In the time series of film thickness, the amplitude of three-dimensional waves can be represented by the standard deviation $\sigma _{h}$. The values of $\sigma _{h}$ are obtained from sensor measurements at the downstream position ($\hat {r} = 109\,\textrm {mm}$) for various input conditions, $\hat {\varOmega } = 31.4{-}314.2\,\textrm {rad}\,\textrm {s}^{-1}$ and $\hat {Q}=1.3{-}2.1\,\textrm {m}^{3}\,\textrm {s}^{-1}$ (figure 17). The data of $\sigma _{h}$ between the numerical simulations and experimental measurements are not entirely consistent, particularly for large $\hat {\varOmega }$. Thus, the experimental data are used to obtain $\sigma _{h}$ whereas the values used to calculate the generalized Reynolds number $R_g$ in (5.1) are from the numerical simulations. The standard deviation $\sigma _h$ monotonically decreases in the low-$R_{g}$ regime, where the data are fitted to $\sigma _{h}=0.35R_{g}^{-3.5}$. This result seems to contradict the trend of $\varLambda$ solitons in falling films (Demekhin et al. 2007), where the volume of the soliton increases with $R_{g}$ and saturates as the flow becomes convectively stable. However, the range of $R_{g}$ considered in this study is more than 10 times that of the previous work ($R_{g}=0{-}0.2$) because the centrifugal acceleration $\hat {r}\hat {\varOmega }^{2}$ in the definition of $\gamma$ (5.1) greatly exceeds the gravitational acceleration $\hat {g}$ used for falling films. Importantly, the magnitude of the waves, $\sigma _{h}$, stops decreasing as the local flow enters the turbulent regime, $\langle Re_{L} \rangle > 100$ (figure 17). This suggests that the horizontal fluctuation induced by the backscattering of three-dimensional turbulence inside the film flow limits the reduction in $\sigma _h$ as $R_{g}$ increases.

Figure 17. Standard deviation $\sigma _h$ of the time series for film thickness at $\hat {r} = 109\, \textrm {mm}$ with respect to the generalized Reynolds number $R_{g}$. Insets A and B are the film thickness contours for $\hat {\varOmega }=31.4\,\textrm {rad}\,\textrm {s}^{-1}$, $\hat {Q}=12.5\,\textrm {mL}\,\textrm {s}^{-1}$ and $\hat {\varOmega }=78.6\,\textrm {rad}\,\textrm {s}^{-1}$, $\hat {Q}=16.7\,\textrm {mL}\,\textrm {s}^{-1}$, respectively. The dashed line is an exponential fitting line obtained from the points of $\langle Re_{L} \rangle \leq 100$.

6. Concluding remarks

We have developed a numerical model for simulating the formation and evolution of three-dimensional surface waves in a liquid film spreading over a rotating disk. Our IBL model improves existing IBL models for rotating films to include the turbulent regime by incorporating a power-law velocity profile and the backscatter effect. The Reynolds number based on the local flow rate is introduced as a criterion for determining whether the turbulent or laminar formulation should be employed locally at a given instant. The modified IBL equations successfully capture the transition from upstream concentric waves to downstream three-dimensional waves and the behaviours of three-dimensional waves under high flow rates, as confirmed by experimental measurements. The inclusion of the backscatter model in the turbulent regime leads to a greater disturbance of the concentric waves and the earlier onset of three-dimensional waves. A greater backscatter constant produces more predictable and periodic patterns for the crests of the three-dimensional waves. The amplitude of these waves tends to decrease monotonically with increasing generalized Reynolds number based on centrifugal acceleration and the local Reynolds number. In terms of separation and coalescence, the unique interaction process among the waves is also observed downstream.

This study has elaborated the complicated behaviours of three-dimensional surface waves, emphasizing the role of backscatter. The proposed IBL model has broad applicability to other film flows with large flow rates, such as rotating films with hydraulic jumps as well as gravity- and shear-driven films on different non-rotating solid geometries. Future work should focus on enhancing the accuracy of the IBL model by exploring alternative backscatter models and higher-order numerical schemes. It is also important to extend this IBL model to investigate the effects of these surface waves on key features relevant to industrial applications, including particle transport, heat and mass transfer, and contact line dynamics. Furthermore, a quantitative analysis of the backscatter process inside the film flow, using high-fidelity three-dimensional simulations, would be valuable.