A Boundary conditions

Cavitation inception of a van der Waals fluid at a sack-wall obstacle


Cavitation in a liquid moving past a constraint is numerically investigated by means of a free-energy lattice Boltzmann simulation based on the van der Waals equation of state. The fluid is streamed past an obstacle and, depending on the pressure drop between inlet and outlet, vapor formation underneath the corner of the sack-wall is observed. The circumstances of cavitation formation are investigated and it is found that the local bulk pressure and mean stress are insufficient to explain the phenomenon. Results obtained in this study strongly suggest that the viscous stress, interfacial contributions to the local pressure, and the Laplace pressure are relevant to the opening of a vapor cavity. This can be described by a generalization of Joseph’s criterion that includes these contributions. A macroscopic investigation measuring mass flow rate behavior and discharge coefficient was also performed. As theoretically predicted, mass flow rate increases linearly with the square root of the pressure drop. However, when cavitation occurs, the mass flow growth rate is reduced and eventually it collapses into a choked flow state. In the cavitating regime, as theoretically predicted and experimentally verified, the discharge coefficient grows with the Nurick cavitation number.


I Introduction

The analysis of cavitation is a significant concern for both fundamental fluid-dynamics and engineering applications such as high-pressure atomizers, spray generators, pumps, propellers, hydraulic turbines, etc. Cavitation is defined as the rupture of a liquid due to a pressure drop, which falls below a certain critical value, at approximately constant liquid temperature, with the subsequent formation of vapor bubbles  Brennen (1995); Dabiri, Sirignano, and Joseph (2007). These bubbles collapse suddenly when they encounter a region with higher pressure  Koivula (), thus releasing a large amount of energy in the form of shock and possibly light waves Koivula (); Brennen (1995); Lohse (2005). This phenomenon causes noise, vibrations and above all it can significantly affect the performance and damage the solid structures (cavitation erosion  Koivula ()) of engineering devices. Despite the passing of more than one century since the first studies of Reynolds Reynolds (1873) and Parsons Parsons (1906) on the effects of cavitation on ship propellers, a comprehensive understanding of the physics of the phenomenon is still lacking as its analysis is challenging from both an experimental and a numerical point of view.

From an experimental point of view, cavitation is a relevant subject addressed by several researchers for both external Kawanami et al. (1997); Weiland and Vlachos (2012); Huang, Ducoin, and Young (2013) and internal flows Winklhofer et al. (); Payri et al. (2013); Sou, Hosokawa, and Tomiyama (2007); Sou et al. (2006); Mishra and Peles (2005); Tamaki et al. (1998); Hiroyasu (2000); Payri et al. (2004); De Giorgi, Ficarella, and Tarantino (2013); Gavaises et al. (2009); Giannadakis, Gavaises, and Arcoumanis (2008). In the field of high-pressure injection processes an inclusive study of this phenomenon is fundamental, since the disturbance induced by cavity bubbles promotes the jet atomization Tamaki et al. (1998); Hiroyasu (2000), thus affecting the mixing and combustion process. In this context a very accurate and complete documentation of the flow phenomena, which occur in a transparent quasi two-dimensional real size nozzle, was provided by Winklhofer et al. Winklhofer et al. (). In their work, the flow properties are described with both the conventional hydraulic characterisation, i.e., with data for mass flow as a function of the pressure drop, and with measurements of the vapor field, of the pressure field, and of the velocity profiles. The results show as the onset of cavitation is located near the wall at the entrance of the hole, specifically, in the shear layer between the incoming fluid, which separates from the entrance corner, and the recirculating fluid, which is near the nozzle wall. Velocity peaks at the boundary of the recirculation zone were found. These peaks are very large near the liquid-gas interface for choked flow conditions. Payri et al. Payri et al. (2013) also studied the behavior of cavitation in a transparent real size cylindrical nozzle by using four different fluids. They performed both hydraulic characterization and flow visualization and showed that cavitation begins before the mass flow collapse. Additionally their experiments suggest that fluids with smaller viscosity tend to cavitate sooner. Flow visualization and laser Doppler velocity measurements were carried out by Suo et al. Sou, Hosokawa, and Tomiyama (2007) at various conditions of cavitation and Reynolds numbers. They assert that both the cavitation regime, inside the nozzle, and the liquid jet, near the nozzle exit, are strongly affected by the cavitation number and less by the Reynolds number. The role of surface tension effects on cavitation in microdevices was pointed out by Mishra and Peles Mishra and Peles (2005). They examined cavitating flows through a two-dimensional micro-orifice, wide, and found that with respect to macroscale counterparts a much larger effort is needed to promote cavitation.

From a numerical point of view, cavitation is part of the wider family of multiphase flows Kuo and Acharya (2012). In this context, several efforts have been addressed in order to develop models able to capture the complex nature of such flows Chung, Pak, and Chang (2004); Niu and Lin (2006); Darbandi and Sadeghi (2010); Giannadakis, Gavaises, and Arcoumanis (2008); Yeom and Chang (2006). The main difficulties, dealing with multiphase flows, are related to the presence of an interface between the different phases, which quickly changes its shape and where huge variations of the fluid properties are located. Traditional approaches are the Eulerian-Eulerian and Eulerian-Lagrangian modeling, with the vapor phase treated as a continuum or a discrete phase, respectively. The Eulerian-Eulerian approach requires ‘jump conditions’, in order to take into account the sharp changes across the interface, and techniques, such as the volume of fluid (VOF) Hirt and Nichols (1981), to track the interface. Another widely used approach for the simulation of cavitating flows, is the homogeneous equilibrium flow model Wang et al. (2001); Utturkar et al. (2005); Tseng et al. (2010). This model considers the two phases perfectly mixed, thus solving the conservation equations for the whole mixture. The liquid (or vapor) mass (or volume) fraction is introduced to model the density change Utturkar et al. (2005). The problem can be closed by using some specific equation that couples density with pressure or a transport equation for the liquid or gas phase  Wang et al. (2001); Utturkar et al. (2005); Tseng et al. (2010). The latter approach requires the definition of source and sink terms which are modelled by using either empirical correlation, interfacial dynamics or bubble dynamics Tseng et al. (2010).

In the last decades several studies Chung, Pak, and Chang (2004); Niu and Lin (2006); Darbandi and Sadeghi (2010); Giannadakis, Gavaises, and Arcoumanis (2008); Salvador et al. (2013); Goncalvès and Charrière (2014) made use of these models in order to simulate cavitating flows. Niu and Lin Niu and Lin (2006) developed a multiscale model in order to predict slurry erosion in ducts subjected to cavitating flows. In their work an Eulerian two-phase mixture model is used to simulate the cavitating flow field. The model employs three equations of state for water, vapor, and saturation states in order to take into account the phase-change process. Darbandi and Sadeghi Darbandi and Sadeghi (2010) examined the fluid dynamic behavior of a cavitating flow through a sharp-edged circular orifice, by using a homogeneous approach, in order to untangle the question whether it is necessary to consider the effect of non-condensable gas when dealing with cavitation. Another important study, based on the Eulerian-Lagrangian approach, was performed by Giannadakis et al. Giannadakis, Gavaises, and Arcoumanis (2008). Here, cavitation is triggered by pre-existing nuclei and a full Rayleigh-Plesset equation is used to simulate the bubble growth and collapse dynamics. The authors claim that for the first time a model which takes into account the intrinsic stochastic nature of cavitation and also physical phenomena such as bubble breakup, coalescence and turbulent dispersion, is presented. The predictive capability was assessed through comparisons with experimental data obtained in real-size and enlarged models of diesel nozzle holes. For the sake of brevity and considering that the body of literature in this field is quite large, we prefer to refer the reader to other publications (Refs.  Brennen, 1995; Knapp, Daily, and Hammitt, 1970; Lecoffre, 1999; Wang et al., 2001; Utturkar et al., 2005; Tseng et al., 2010 ), where an extensive review of many aspects of cavitation can be found from both physical and numerical point of view.

The above mentioned studies have clarified many features of cavitation processes from the phenomenological point of view. However, fundamental aspects of cavitation, such as the physics that underlies cavitation inception, remain not satisfactorily understood. Important questions concern the role played by viscous stress in the formation of the first cavitating bubbles. Indeed, while in the literature it is typically assumed that the bulk value of the pressure determines the creation of a cavitating bubble, Joseph Joseph. (1998) argued that, in fact, the significant parameter is the total stress where both bulk pressure and viscous stress are taken into account. Experiments are not able to analyse this phenomenon, since, in most of the experimental tests, nucleation has a stochastic behavior due to its heterogeneous nature. Numerical simulations therefore appear to be a more promising approach for the purpose of comparing the traditional pressure criterion with Joseph’s expression Joseph. (1998). In this context, it would be also important to clarify the role of interfacial energies, contributing to the total stress balance, and expected to be relevant in any process of bubble formationBray (1994); Onuki (2002).

We explore in this work the potential of the lattice Boltzmann method (LBM)Succi (2001, 2015); Benzi, Succi, and Vergassola (1992); Gonnella and Yeomans (2009) in simulating cavitating flows. We analyze the flow and thermodynamic behavior underlying the inception of cavitation in a constrained geometry and we will examine the validity of Joseph’s criterion. LBM is an algorithm that solves the Boltzmann equation in a discretized phase space. The method allows to recover the continuity and the Navier-Stokes equations in the continuum limit at the second order in the Knudsen number. This method has been widely applied in several contexts to study simple and complex fluids Dünweg and Ladd (2009). Here the main benefit of the LBM approach is that it does not need an ad hoc model to capture the distinctive features of cavitating flows, such as the phase transition process, with its multiscale nature, and the dynamically varying liquid-vapor interface. Indeed, by using an appropriate non-ideal equation of state Sofonea et al. (2004); Cristea et al. (2006); Gonnella, Lamura, and Sofonea (2007, 2009); Cristea et al. (2010), which takes into account the real fluid behavior, cavitation emerges spontaneously from the dynamical approach Falcucci et al. (2013a). In our simulations we considered homogeneous cavitation induced by a fast flow past a sack-wall obstacle. This standard geometry is widely used in numerical studies of cavitation Dabiri, Sirignano, and Joseph (2007); Darbandi and Sadeghi (2010); Giannadakis, Gavaises, and Arcoumanis (2008); Falcucci et al. (2013a, b). It can be considered a benchmark problem in order to have a comprehensive understanding of flow-induced cavitation as well as a simplified geometry useful for analyzing high pressure injectors.

LBM was first used for studying cavitation by Sukop and Or Sukop and Or (2005). In that paper the cavitation of a single bubble was induced on a square geometry by decreasing the external pressure on top and bottom walls in absence of any imposed flow. A similar study was later addressed in Refs. Chen, Zhong, and Yuan, 2011; Zhong, Zhong, and Bai, 2012. Only recently, the much more complex problem of cavitation in constrained flow environments has been successfully investigated Falcucci et al. (2013b, a). Falcucci et al. Falcucci et al. (2013a) applied a lattice Boltzmann model with a two-belt implementation Falcucci et al. (2007) of a Shan-Chen pseudo-potential Shan and Chen (1993) to a geometry similar to the one considered in this work and found cavitation underneath the obstacle. They relate their results to the stress criterion put forward by Joseph Joseph. (1998) once the cavitation threshold is properly modified by considering surface-tension contribution. In the model of Falcucci et al. Falcucci et al. (2013a) the chosen pseudo-potential mimics a non-ideal equation of state. Here we use a lattice Boltzmann equation (LBE) model Coclite, Gonnella, and Lamura (2014) which is derived by using a Gauss-Hermite projection of the continuum equation Shan, Yuan, and Chen (2006). A body force in the LBE Guo, Zheng, and Shi (2002); Tiribocchi et al. (2009); Gonnella et al. (2010); Gonnella, Lamura, and Tiribocchi (2011), based on the van der Waals (vdW) free-energy functionalRowlinson and Widom (1982), guarantees that the fluid locally satisfies the van der Waals equation of state. The present approach allows to rigorously obtain in the continuum limit the continuity and the generalized Navier-Stokes equations without spurious terms. However, the presence of unphysical spurious velocities at interfaces cannot be avoided and a general 9-point second-order finite difference scheme to compute the spatial derivatives was used. This approach is known to reduce these unwanted effects Shan (2006); Sbragaglia et al. (2007). Another advantage of this model is that the phase diagram does not depend on the relaxation time of the LBE which controls fluid viscosity. Moreover, the surface tension is properly included in the pressure tensor so that it would be interesting to analyse the inception of cavitation and the role of stress and interfacial contributions by our method.

Additionally, for the first time, we employ LBM for direct hydraulic characterization of the constrained duct and compare the results with those obtained in experimental studies in similar geometries.

The present work is organized as follows: In section II the mathematical model is introduced, together with some remarks on the implementation. In section III a local and a global criterion for predicting and measuring cavitation are introduced. Numerical results for the former are discussed in section IV and for the latter in section V. Finally, the conclusions are summarized.

Ii Model and Implementation

We briefly describe here the two-dimensional isothermal lattice Boltzmann model for non-ideal fluids introduced in Ref. Coclite, Gonnella, and Lamura, 2014. More details on the method and on the algorithm implementation are given in Ref. Coclite, Gonnella, and Lamura, 2014 and in Appendix A. The ability of the model in reproducing the correct phase diagram of a van der Waals fluid and the surface tension behavior are also described in Ref. Coclite, Gonnella, and Lamura, 2014.

ii.1 Lattice Boltzmann Equation for a van der Waals Fluid

The time evolution of the system is defined by means of the dimensionless lattice Boltzmann equation (LBE)


where the distribution functions are associated to the lattice velocities at lattice site and discrete time . Here the square lattice is considered, where 9 distribution functions are defined at each node with links to nearest and next-to-nearest neighbors. The forcing terms depend on the the force density acting on the fluid and have to be conveniently determined, is the relaxation time, and is the time step. The equilibrium distribution functions are given by a second-order Hermite expansion of the Maxwell-Boltzmann distribution function Shan, Yuan, and Chen (2006)


where is the fluid density, is the fluid momentum, and is the temperature, which is not a dynamical variable in this model. The Hermite expansion fixes the values with being the lattice unit (lu) for horizontal and vertical links with , for diagonal links with , for the rest velocity, and the weights for , for , and .

The forcing term in Eq. (1) is given by Guo, Zheng, and Shi (2002); Coclite, Gonnella, and Lamura (2014)


In order to describe a van der Waals fluid the force has to be


where is the ideal gas equation of state (EOS),


is the van der Waals EOS with critical point at , and


is the pressure tensor. It can be computed from the free-energy functional Rowlinson and Widom (1982)


with bulk free-energy density given by


The term proportional to takes into account the energy cost for the formation of interfaces and allows to change the surface tension independently of the temperature.

Chapman-Enskog expansion at second order shows that the continuity equation


and the Navier-Stokes equation


can be properly recovered with being the fluid viscosity.

ii.2 Geometry and Boundary Conditions

The model is used to simulate a fluid flowing in a channel from the left to the right past a sack-wall obstacle. This geometry is frequently adopted in numerical simulations of cavitation Darbandi and Sadeghi (2010); Falcucci et al. (2013b, a).

The simulated system consists of a rectangular simulation box where the ratio of length and height is given by . The fluid moves into the system from the left hand side and, after travelling for a distance , finds an obstacle of height which is connected to the top of the channel. The channel is thus constrained to a height , half of the initial one, for the remaining length .

In the following we will use two different sets of inflow/outflow boundary conditions (BC). We will refer to the first choice as fixed-density BC and to the second one as fixed-pressure BC. In the first case a fixed density and a fixed velocity are used at the inflow layer. On this layer the distribution functions are given by the values of the equilibrium distribution functions . At the outflow a zero-gradient velocity condition is enforced with fixed density . In general we use where is the initial density condition applied to all of the bulk fluid domain.

To compare with experimental conditions Payri et al. (2013); Nurick (1976), also fixed-pressure BC at inlet and outlet are considered. At the inflow the total upstream pressure , which is the sum of the kinetic contribution and of the van der Waals pressure, is fixed and the -component of the velocity is calculated enforcing a zero-gradient velocity condition. At the outflow a similar approach is used but the fixed-pressure used at the boundary is just the van der Waals one. (More details on the implementation of the inlet and outlet boundary conditions are given in Appendix A.1)

At walls and corners no-slip BC with neutral wetting and local density conservation are enforced Lamura and Gonnella (2001); Tiribocchi et al. (2011) (see details in Appendix A.2).

ii.3 Thermodynamic Consistency of the Model

In all simulation runs the lattice unit was fixed to 1 and the value was used unless specified otherwise. Several quantities, which will be later useful, were measured at different values of the temperature and of the interface term and we will give their values in the following. They will be expressed in units of , , , and .

In table I the equilibrium vapor and liquid densities measured in simulations are compared to the corresponding analytic values and , respectively. Runs were performed on a lattice of size where a vapor droplet of radius at density was immersed in liquid at density and let to equilibrate. Both equilibrium vapor and liquid densities are closer to the analytic values at larger and .

0.9 0.0 1.65727 1.59266 0.425742 0.39654
0.1 1.60999 0.40368
0.2 1.62043 0.40830
0.3 1.62731 0.41143
0.85 0.0 1.80714 1.7067 0.31973 0.28552
0.1 1.7331 0.29401
0.2 1.7492 0.29967
0.3 1.7599 0.30353
0.8 0.0 1.93270 1.81931 0.23966 0.20405
0.1 1.82595 0.21424
0.2 1.85206 0.22049
0.3 1.86695 0.22476
Table 1: Values of the equilibrium vapor and liquid densities, and , measured from simulations, are compared to the analytical values and for different values of the temperature and of the interface term .

Table II reports the maximum density from and simulations at which spinodal decomposition on the liquid branch was observed. For this purpose a very slightly perturbed system is initialized at a density close to the analytic spinodal value and it is observed whether the system phase separates. was measured in a lattice while in a system. The analytic value of the liquid branch spinodal density is well reproduced in both the one- and the two-dimensional cases. In Appendix B we investigated the impact of shear on the spinodal density. The results suggest that shear does not significantly impact the following results on cavitation inception of section IV.

0.9 0.0 1.39160 1.3911 1.3895
0.1 1.3907 1.3881
0.2 1.3904 1.3868
0.3 1.3901 1.3856
0.85 0.0 1.48880 1.4884 1.4873
0.1 1.4881 1.4862
0.2 1.4879 1.4852
0.3 1.4876 1.4842
0.8 0.0 1.57428 1.5739 1.5731
0.1 1.5737 1.5722
0.2 1.5735 1.5713
0.3 1.5733 1.5705
Table 2: Values of liquid-branch spinodal density from () and () simulations are compared to the analytic values for different values of the temperature and of the interface term .

Surface tension was measured directly through the free energy in a one-dimensional system. For this purpose two quantities were computed:

  • The bulk free energy


    calculated with half of the system in the pure liquid phase and with the other half in the pure vapor phase. The densities used are the numerical equilibrium values reported in table I.

  • The total free energy including gradient contributions after the system had fully equilibrated


The difference would then yield the total free energy change due to the interface including numerical contributions that cannot be attributed to the -term, but that effectively act as an interface energy contribution. The surface tension is then calculated as


Additionally the surface tension was also calculated in a 1D system according to its definition


by using the numerical values of in the equilibrated system () and also the approximate equilibrium profile Wagner and Pooley (2007) ()


The results are given in table III. For all three temperatures, the values of and coincide fairly well and the agreement improves with growing . The numerical contribution added to any value of or for non-zero is expected to give the equivalent surface tension for non-zero . This appears to be reasonably consistent, moreso for higher temperature and lower surface tension .

0.9 0.0 0.0418 0 0
0.1 0.0634 0.0361 0.0272
0.2 0.0814 0.0516 0.0493
0.3 0.0971 0.0636 0.0681
0.85 0.0 0.0738 0 0
0.1 0.1124 0.0649 0.0472
0.2 0.1449 0.0929 0.0868
0.3 0.1733 0.1150 0.1208
0.8 0.0 0.1555 0 0
0.1 0.1665 0.0968 0.0687
0.2 0.2157 0.1403 0.1277
0.3 0.2589 0.1739 0.1791
Table 3: Measured values of the surface tension in a one-dimensional system for different values of the temperature and of the interface term (see the text for details).

Iii Cavitation Prediction Criteria

Cavitation criteria are quantitative relations used to predict the occurrence of cavitation. Macroscopic criteria are based on the hydraulic characterization of the flow and on the definition of the cavitation number (see below). They are frequently used in engineering applications. Microscopic criteria are based on the behavior of local variables. They relate the local stress to a maximum pressure at which the cavity can be formed. We describe now two cavitation criteria and some relations which were investigated numerically and that will be useful in the following.

Hydraulic characterization originates in the experimental field. This is in part because when dealing with engineering devices, direct local detection of cavitation is difficult due to lack of visual access to the relevant area of the experiment Koivula (). Then cavitation can be detected by looking at the stationary mass flow rate as a function of the pressure drop. Indeed, for an incompressible flow, the theoretical mass flow rate obtained from Bernoulli’s equation is directly proportional to . However, the actual mass flow rate is smaller than the theoretical one because of losses due to boundary layer, vena contracta, turbulence and cavitation. Therefore, it is common to define a discharge coefficient as the ratio between the actual flow rate and the theoretical one


where is the actual mass flow rate and the equilibrium liquid density. Experimental works Winklhofer et al. (); Payri et al. (2013) showed that by fixing and decreasing mass flow rate grows with until it reaches the so-called critical cavitation point after which, though decreasing , mass flow remains constant and the flow is called choked. Therefore, by looking at the behavior of the discharge coefficient, it is possible to distinguish a region where cavitation is absent or negligible and another region where cavitation significantly reduces or completely removes a further increase in mass flow rate.

When the flow is not cavitating, is primarily a function of the Reynolds number, which is defined as


where is the theoretical Bernoulli velocity measured at the end of the channel and is the kinematic viscosity. In contrast, for cavitating flow, the cavitation number we will define is the main influence on Payri et al. (2013); Nurick (1976).

The definition of the cavitation number is based on Bernoulli’s equation. It relates the difference between the pressure at a certain location and the maximum pressure at which a cavity can be formed, often the vapor pressure , to a dynamic pressure. However, the specific location at which the relevant observables are measured, and in particular, the geometry are not uniquely defined. Different approaches have been proposed Brennen (1995); Nurick (1976); Payri et al. (2013). Following Nurick Nurick (1976) we choose here a definition appropriate for a constrained orifice where the upstream pressure is considered and the dynamic pressure is expressed as :


This definition originates from the one-dimensional model discussed in Ref. Nurick, 1976 and allows the calculation of for a cavitating flow as


Here is the contraction coefficient given by the ratio , where is the vena contracta flow section. The above is valid only in the presence of cavitation and therefore represents an extremely useful tool for indirect identification of cavitation in practical applications as also shown in the experimental work of Payri et al Payri et al. (2013).

In contrast to the macroscopic considerations one can also establish local criteria based on the behavior of the total stress of the fluid, with respect to some critical pressure, at a given position where cavitation might occur. Joseph et al. Dabiri, Sirignano, and Joseph (2007) argue that it is the largest directional stress overcoming a pressure threshold that gives rise to cavity formation. For a simple fluid the total or Cauchy stress tensor can be expressed by the following constitutive equation


where is the bulk pressure, is the viscous stress tensor, and is the unit tensor. The original version of Joseph’s maximum tension criterion states that the quantity to be compared to a pressure threshold is the largest eigenvalue of . Cavitation would occur at points where


Taking into account that is a negative number one may interpret this as follows: If the magnitude of the maximum principal stress drops below the critical pressure threshold for vapor formation, the formation of a cavity is expected. This generalizes the usual assumption Brennen (1995) that it is the pressure or the local mean stress that determines the formation of a cavity in a fluid. The mean stress is defined in terms of the trace of as


where is the spatial dimension. For a Newtonian and incompressible fluid is the stress deviator of the total stress tensor, since , and therefore the mean stress is equal to .

In the microscopic view there is one caveat here. The bulk pressure is not sufficient to express all pressure contributions, which also include interfacial terms. For this reason the full pressure tensor Eq. (6) must be considered, so that we introduce as the total stress the expression


with the corresponding mean stress denoted as . In this context it is natural to assume that also the pressure threshold is affected by interfacial contributions. Falcucci et al. Falcucci et al. (2013a) include the Laplace pressure in the pressure threshold so that


Here is the surface tension and the minimal radius for the formation of a vapor bubble. In experimental studies and numerical engineering works is usually identified with the vapor pressure of the fluid. While this is a reasonable assumption in experimental investigations due to the almost inevitable presence of nucleation kernels, in this particular study the equilibrium vapor pressure is likely not a very good estimate. In the case of a fluid without any kind of nucleation kernels or other sources of local symmetry breaking such as the system used here, one expects to be close to the spinodal pressure rather than the equilibrium vapor value.

Iv Cavitation Inception

iv.1 Overview

Figure 1: Contour plots of the total pressure at different times in simulation on a lattice with fixed-density BC together with the flow field. The parameters() are chosen such that the typical phases of the cavitation can be observed. In a) the transient initial wave front propagating to the right is visible. b) shows the formation of the stable flow profile after the wave has made contact with the outflow boundary layer. In c) the flow profile has matured and the pressure depletion underneath the obstacle corner becomes visible in dark blue color. In d) the initial cavity has formed which grows in e) to finally, after a long time, extend to the end of the channel while covering a significant section of the lower obstacle boundary in f). The vapor fraction in d), e), and f) has been colored in cyan. The pressure of the vapor area is about and is between the equilibrium pressures calculated from the numerical equilibrium densities of table 1 where and .
Figure 2: Total, static and dynamic pressures as function of lattice position at , i.e. half the constrained channel length, in the channel of Fig. 1 at . Here, is the lower boundary of the channel whereas is at the wall of the obstacle.
Figure 3: Contour plot of the density and macroscopic velocity shown at four different times shortly after the pressure reached the liquid-branch spinodal value. Simulation parameters: mesh with lattice spacing of and time unit adjusted accordingly, fixed-density BC, , and inflow velocity .

As previously mentioned, the simulations are performed in a two-dimensional channel constrained by a sack-wall. The simulation box has the size of lattice units (lu). Resolutions of times and twice the number of lattice nodes with correspondingly reduced and have been tested and, where stability allows, the results are consistent with the findings at the standard resolution shown here.

An example of the behavior of the system from initialization at rest to the formation of a stable velocity profile and eventually the formation of a cavity is given in Fig. 1. The total pressure and the flow lines are given at different times. In this example, as well as in most measurements in this work, the temperature is , the surface tension parameter , and the viscosity . The inflow and outflow boundaries are computed according to the fixed-density BC at an upstream velocity of and , which is close to the equilibrium liquid density. The static pressure corresponding to this density is . Non-boundary nodes are initialized at rest with the same density . The fluid begins to stream in from the left hand side. The initial pressure front first collides with the obstacle wall and is, as seen in Fig. 1 (a), partially reflected. The attenuation of the inflow velocity near the upper and lower channel walls, according to Eq. (26), reduces the pressure near the corners. Instabilities due to the rebounding pressure wave that are otherwise prominent are thus eliminated.

Once the pressure wave approaches the outflow boundary the zero velocity gradient boundary at the outflow layer allows material to leave the system to the right. As the channel walls enforce zero velocity the fluid near the center of the channel accelerates and the static pressure is reduced. This reduction first occurs near the outflow boundary and then propagates towards the left as the material in the channel begins to move. The outflow BC fix the density at the outflow layer and, as , no vapor can originate from the outflow layer in this example. The stationary flow pattern of later stages of the simulation begins to form, see Fig. 1 (b).

The configuration after another 1500 time steps is shown in Fig. 1 (c). Pressure and the flow pattern have become almost stationary and, as expected, the fluid velocity profile reaches a maximum at the center of the constrained channel. The static pressure is constant perpendicular to the flow direction and any change in total pressure is due to the increased dynamic pressure contribution near the center of the channel. Fig. 2 shows these pressure profiles in the -direction at half the constrained channel length. At this time it is also visible that the lowest static pressure in the system is just underneath the obstacle. The pressure distribution and flow pattern observed here are rather typical for the sack-wall geometry (see, e.g., Ref. Falcucci et al., 2013a).

With the static pressure falling below a threshold value, a nucleation kernel is formed which then evolves into a vapor cavity. This happens first underneath the obstacle. This small initial cavity is first visible in Fig. 1 (d). In the last three images of Fig. 1 the developed vapor is colored in cyan. The formation of this initial cavity is investigated in detail in section IV.2 and a higher resolution view of the density of the initial bubble formation is shown in Fig. 3 albeit at lower inlet velocity than that of Fig. 1. In this example the area near the wall exhibits a pressure well below the equilibrium vapor pressure and thus the cavity immediately begins to grow. The exact shape and direction of this growth depend on the pressure or density distribution in the channel. The cavity typically grows towards the center of the channel first. Then, due to the increase of the velocity towards the center of the channel, it is stretched along the flow lines of the fluid as seen in Fig. 1 (e). The bubble continues to grow until it hits the outflow boundary layer. The de-wetting of the gas bubble at the obstacle is a slower process but in almost all observed cases, a sufficiently long running simulation will exhibit a fully formed gas bubble covering most of the lower obstacle surface. An intermediate state of this bubble growth is shown in Fig. 1 (f). In this example, at the point of smallest diameter for fluid flow, the vena contracta, less than of the channel width is available for fluid transport resulting in a reduced mass transport. This effect and its consequences are discussed in section V. We finally mention that, unlike in the computation results shown in Fig. 1, if the downstream pressure remains well above the equilibrium vapor pressure the de-wetting of the wall remains partial in our numerical experiments.

iv.2 Local Cavitation Inception

In order to understand the local mechanism of vapor formation and analyze the cavitation criteria discussed in section III, the conditions under which the cavity forms were carefully investigated. Fig. 3 mentioned in the previous section gives in double standard resolution the density and fluid velocity underneath the obstacle during initial bubble formation. This higher resolution computation does not show significant quantitative or qualitative differences compared with the standard resolution case. In Fig. 3 (a) the spinodal density is reached for the first time. The formation of the interface is initiated in Fig. 3 (b) which then quickly progresses to the nucleation of a vapor cavity (Fig. 3 (c) and (d)). It is worth noting that at this level of observation configurations that will not cavitate but have otherwise very similar parameters will be visually indistinguishable from that of Fig. 3 (a). They will approach a state similar to Fig. 3 (b) and will recede to a configuration similar to Fig. 3 (a) where the density remains close to but vapor formation does not occur.

We observe that prior to interface formation a single lattice site underneath the obstacle assumes the lowest density in the channel. This site is situated two lattice spacings below the obstacle and four lattice spacings to the right of the vertical wall of the obstacle for the resolution of sites. In the following, individual measurements of pressure and density with regards to cavitation inception are taken at this lattice position. Cavitation is observed at all temperatures considered in table 1. Due to stability concerns, however, only cases of are accessible at lower temperatures. For the best compromise between correct reproduction of the equilibrium phase diagram and the best stability it is chosen to limit the analysis of cavitation inception to . The macroscopic parameters available are the mean density and the velocity imposed on the inflow boundary . For observing cavitation we choose an initial density which is, as before, slightly below the equilibrium liquid density and in the metastable regime of the phase diagram, and vary the inflow velocity. If not otherwise specified, the first set of BC with fixed-density at inlet and outlet is used.

Figure 4: Density at the lowest density site during the “transient cavitation measurement” as a function of simulation time for different initial velocities for a) and b) . Cavitation only occurs for the highest shown (red continuous line). Other parameters: , .

In Fig. 4 the density at the lowest density site underneath the obstacle is given as a function of time for different inflow velocities. In both cases, and , for sufficiently high inflow velocity, ( at , respectively), the density, independently of subsequent actual cavity formation, decreases below the spinodal value. As shown in Fig. 4, if is not large enough this is a transient drop, with the density relaxing to a value close to the spindodal one at larger times. A good estimate for the highest upstream velocity for which cavitation does not occur was found with an iterative approximation. The minimal value of density obtained in this way is significantly below the liquid branch spinodal value of as seen in Fig. 5 for the case with . It appears that cases can be constructed where the density temporarily drops very significantly below the liquid branch spinodal, but, due to the transient nature of this approach, no cavity is created. Likely the time for which this low density environment exists, is insufficient to allow for nucleation. We call this protocol for observing cavitation “transient cavitation measurement”.

Figure 5: Density at the lowest density site during the “transient cavitation measurement” as a function of simulation time for different initial velocities at . Results here were obtained using an iterative approach minimizing the inflow velocity. A similar behaviour has been found at . Other parameters: , .
Figure 6: Density measured at the lowest density site in simulation with fixed-density BC with velocity increments of every time updates. The dashed lines indicate the lowest stable density observed and corresponds to the configuration shown in Fig. 8. Other parameters: , . A similar measurement was conducted but with fixed-pressure BC although at much larger pressure steps. The inset plots do show that observed thresholds appear at similar values. Differences between inset and main figure are likely due to the much longer equilibration time of the fixed-pressure BC.

To avoid the difficulties in quantifying the density threshold by the above approach, the protocol for observing incipient cavitation was altered. The system is initialized at an inflow velocity small enough that the system is not cavitating in the initial transient approach. Then, to establish a sequence of stationary states not affected by transient behavior, the simulation is run for a time of order , see Fig. 6, significantly exceeding the one needed for transient cavitation formation, and the inflow velocity is then increased by a small . This process, that we call ”‘steady-state cavitation measurement”’, is repeated until cavity formation is observed. The lowest value for which the density still converges in the liquid regime is considered the closest pre-cavitation value we can attain and is used as a threshold. These threshold measurements are given in the main plot of Fig. 6 for fixed-density BC.

Similar measurements were performed using the alternative fixed-pressure inlet and outlet boundary conditions (discussed in Appendix A), varying the upstream pressure instead of the upstream velocity. These results are shown in the insets of Fig. 6. For both types of boundary conditions we observed that the threshold densities are well below the liquid branch spinodal value.

The above results indicate that the simple assumption that cavitation occurs when the fluid density or pressure reach the liquid branch spinodal values is insufficient. Taking the liquid branch coexistence pressure or density as critical threshold would further increase the difference between the measured value and this critical value. One might wonder if the spinodal values listed in Table 2, referring to a system at rest, would be changed by flow. It is known that multiphase fluids under applied flow, shear for example, exhibit a critical point depending on the flow rate Gonnella and Pellicoro (2000); Corberi, Gonnella, and Lamura (1998); Onuki (2002). In our model, the jump in the horizontal velocity component underneath the obstacle could be roughly approximated by a linear shear profile (see Fig. 14 in the following). Thus we decided to check how much the spinodal values of our system are affected by a shear flow. The results of these simulations are summarized in Appendix IB. We find that the spinodal values measured at rest are changed very little by shear or, in any case, in a way that cannot explain the cavitation threshold measured above.

One could also consider, instead of the pressure, the behavior of the mean stress defined by (22) and compare this with the spinodal critical pressure value. From Fig. 7 one sees that the behavior of the mean stress is very close to that of the static pressure so that also this quantity cannot be useful for explaining the formation of cavitation. We observe that the van der Waals pressure remains above the spinodal value as predicted by the analytical expression of van der Waals isotherms.

The Joseph criterion or its generalizations discussed in Sect. III, based on the largest stress eigenvalue and also including interface contributions, offer a natural way to further analyze the formation of cavitation in terms of local variables. Before doing this, since interfacial contributions also affect the critical threshold (24) through the Laplace term, we will evaluate the typical size of the low density region that, under proper conditions, will evolve into a vapor bubble.

A magnified view of the density underneath the obstacle at the time of the threshold of Fig. 6 (b) is given in Fig. 8. This configuration is the lowest stable density prior to cavitation formation and can be used to estimate the minimum droplet radius for evaluating the Laplace pressure contributions. The area inside the contour contains all lattice nodes for which . This half circular shape has a diameter of so for this study it is considered .

Figure 7: Evaluation of the Joseph criterion measured in a lu simulation at the cavity inception site underneath the obstacle. Each data point corresponds to one stationary state of the type shown in Fig. 6 after increasing the upstream velocity. The last stable threshold is found at which corresponds to an inflow velocity of . Other parameters: , , .
Figure 8: Contour plot of the density observed directly underneath the obstacle at cavitation inception in a lu simulation at , and kinematic viscosity , i.e. . The black line indicates the threshold of the spinodal density .

The maximum tension criterion of Joseph of Eq. (21) compares the magnitude of the maximum stress with a critical pressure. The maximum tension for both stress definitions (Eq. (20) and Eq. (23)), and also the behavior of the mean stress are shown in Fig. 7. Moreover, two different critical thresholds corresponding to different Laplace contributions in Eq. (24) are plotted. The Laplace pressure contribution is calculated by using the surface tension values , evaluated as described by Eq. (13), and the critical radii and .

According to the criterion one would expect cavitation when the magnitude of the local stress is smaller than the pressure threshold. In Fig. 7 one can observe that the magnitude of the maximum stress eigenvalue, for both stress definitions, is significantly below the spinodal pressure also before cavitation occurs. Therefore the Joseph criterion with the critical threshold given by the spinodal pressure has not any particular meaning in our context. On the other hand, the spinodal pressure adjusted by the Laplace pressure with a minimum droplet radius of appears in close vicinity to the maximum stress magnitude . So it appears that both the maximum stress eigenvalue and the critical pressure with the Laplace contribution are relevant and in this particular experiment they are found in close vicinity to each other. Furthermore there is an appreciable difference between the two definitions of Eq. (20) and Eq. (23) indicating that the surface tension contributions due to gradient terms in the pressure tensor should be considered relevant as well. 2 We also remark that the behavior of the mean stress does not appear relevant in this context.

Simulations at higher (doubled) resolution confirm the above picture with the stress behaving as in Fig. 7. For example, at , the magnitude of the maximum stress eigenvalue almost reaches the pressure threshold before increasing due to cavitation with a plateau at and (at the differences for the values of the pressure thresholds with respect to the case at higher resolution are negligible and the stress plateau is at ).

Similar measurements were performed for the fixed-pressure inflow BC at three different viscosities: , , and and confirm the above picture with a critical threshold diminishing with increased viscosity.

V Macroscopic Cavitation Prediction

In this section cavitation is examined by using the standard approach of hydraulic characterization described in section III. This technique is widely used in experimental works, such as by Payri et al. Payri et al. (2013) and by Winklhofer et al. Winklhofer et al. (), in order to analyse the behavior of high pressure systems. Both, Payri et al. Payri et al. (2013), and Winklhofer et al. Winklhofer et al. (), used an axial-symmetric and a quasi two-dimensional nozzle, respectively, monitored the mass flow rate under steady flow conditions by keeping constant the upstream total pressure and varying the downstream static pressure. Thus they obtained data sets for mass flow rate as a function of the pressure drop independent of the total upstream conditions. The results show that mass flow rate increases by increasing the pressure drop up to a maximum value after which even further decreasing the downstream pressure mass flow remains constant. At this point the flow is said to be choked. Mass flow saturation is mainly caused by cavitation Winklhofer et al. (); Payri et al. (2013) and it accompanied by a reduction of the discharge coefficient Payri et al. (2013). Specifically, Payri et al. Payri et al. (2013), who performed a comprehensive analysis by considering four different fluids with various viscosities in order to cover a wide range in terms of Reynolds number, showed that under non-cavitating conditions mass flow rate grows asymptotically with the Reynolds number whereas under cavitating conditions,  Eq. (19) is met. Therefore, both mass flow choking and discharge coefficient reduction are indirect evidence of cavitation. Moreover, they showed that less viscous fluids tend to cavitate sooner.

The contribution of this work is to test, for the first time, the applicability of the LBM for the analysis of the cavitation phenomenon by using a macroscopic approach closer to engineering practice which relies on observing the results in terms of mass flow rate versus pressure drop and discharge coefficient versus both cavitation and Reynolds number.

Measurements reported in this section are obtained in simulations with the fixed-pressure BC discussed in Appendix A. We used the base resolution of lu, at and , to perform several simulations by fixing the total inlet pressure , and varying for kinematic viscosities . The results in terms of stationary mass flow rate as a function of the square root of the pressure drop are given in Fig.  9. At the mass flow rate increases almost linearly up to (last right filled circle on the dash-dotted line in Fig.  9). Subsequently, further reducing the growth rate is reduced. A similar behavior is observed with kinematic viscosity equal to . Here the mass flow rate increases almost linearly up to (last right filled circle on the dashed line in Fig.  9), then the growth rate is reduced and, finally, the system collapses into a choked flow state. Mass flow choking is well visible for the case with kinematic viscosity equal to when the critical cavitation point is reached at (last right filled circle on the continuous line in Fig. 9). Indeed, further reducing the downstream pressure the mass flow rate shows a small reduction at and then settles to values comparable with those achieved at .

Figure 9: Stationary mass flow rate as a function of the square root of the pressure drop for different values of the kinematic viscosity. Simulation parameters: lattice sites lu with fixed-pressure BC at , and total inlet pressure .

The mass flow choking is clearly due to cavitation. This can be seen in Fig. 10, where the density contour plots for the cases with and are given. There is no vapor with , whereas, when vapor extends attached to the wall, along the entire length of the obstacle. A similar observation is also valid for the cases with kinematic viscosities equal to and . Moreover, the flow lines, reported in the same figure, show that in both cases there is a significant recirculation zone near the front side of the obstacle.

Figure 10: Density contour plots with flow lines. Simulation parameters: lattice of lu with fixed-pressure BC at , , and (a), (b) with kinematic viscosity equal to .

Finally, by considering again Fig.  9, when the kinematic viscosity is increased to , the mass flow rate always increases, almost linearly, with the square root of the pressure drop and it does not show any reduction of the growth rate up to the smallest considered, which for this kinematic viscosity is . A further reduction of would compromise the numerical stability. The absence of a reduction in the mass flow growth rate is due to the fact that up to cavitation does not occur.

Figure 11: Discharge coefficient as a function of the square root of the Nurick Nurick (1976) cavitation number for different values of the kinematic viscosity. Simulation parameters: lattice sites lu with fixed-pressure BC at , and total inlet pressure .

Obviously, the reduction of the mass flow growth rate causes a decrease of as shown in Fig. 11, where the discharge coefficient as a function of the square root of the Nurick cavitation number of Eq. (18) is given (here is determined considering the liquid branch spinodal pressure instead of the vapor pressure in Eq. (18)). Indeed, except for the case with kinematic viscosity equal to , we can identify a critical cavitation number below which there is an abrupt reduction of . In the non-cavitating regime, for cavitation numbers higher than the critical value (filled circles in Fig. 11), slightly increases by decreasing while it reaches a plateau at large . In the cavitating regime, moving towards the critical value of the cavitation number, for kinematic viscosity equal to and , shows a small undershoot before reaching the critical point. The undershoot is not present when the kinematic viscosity is . These findings are qualitatively in good agreement with the experimental ones Payri et al. (2013) since it is possible to identify a critical point below which the discharge coefficient abruptly decreases. In this region is a linear function of , as suggested by Eq. (19) only for the case with , whereas  Eq. (19) is not met when the kinematic viscosity is equal to and .

The partial disagreement with the Nurick law of Eq. (19) can be due to the effects of the Reynolds number which are not taken into account in the derivation of Eq. (19). We show the dependence of and on . In Fig. 12, where the discharge coefficient is given as a function of , one sees that reaches a maximum and then it abruptly decreases due to cavitation effects. This result is qualitatively in good agreement with the experimental one by Payri et. al Payri et al. (2013). Fig. 13 shows that decreases with in both cavitating and non-cavitating cases. Given the dependence of and on the Reynolds number, it is reasonable to expect, coming back to discuss the cavitation regime in Fig. 11, that Reynolds number effects come into play also where cavitation is occurring so that the linear behavior in terms of of Eq. (19) is not completely fulfilled. On the other hand, in the non-cavitating regime of Fig. 11, the reduction of with is mainly triggered by the reduction of , since decreases and consistently increases with increasing .

We have seen that cavitation and Reynolds number affect the mass flow rate and therefore . It could be instructive to consider the above results also looking at the behavior of the velocity and density profiles and analyzing how they are affected by cavitation and . This will also help to understand the role of other factors like compressibility on the behavior of the mass flow rate.

Figure 12: Discharge coefficient as a function of the Reynolds number for different values of the kinematic viscosity. Simulation parameters: lattice of lu with fixed-pressure BC at , and total inlet pressure .
Figure 13: Nurick cavitation number as a function of the Reynolds number for different values of the kinematic viscosity. Simulation parameters: lattice of lu with fixed-pressure BC at , and total inlet pressure .
Figure 14: Normalized longitudinal velocity profiles along with the density profiles at the position across the channel. Simulation parameters: lattice of lu with fixed-pressure BC at , , total inlet pressure , kinematic viscosity equal to and downstream pressure , , , and .

Fig. 14 shows the normalized longitudinal velocity profiles along with the density profiles across the constrained channel at half of its length () for the cases with kinematic viscosity and downstream pressure . The figure shows that density is constant across the channel for the non-cavitating cases with downstream pressure equal to and , whereas in the event of cavitation there is a significant reduction of the effective liquid flow cross section and therefore a reduction of . The smaller effective liquid flow cross section is caused by a sharp density discontinuity which separates the liquid from the vapor. It is also visible that, by going from to , compressibility is not negligible since there is an appreciable reduction of the density which obviously affects mass flow rate and therefore .

At the considered longitudinal distance the transversal velocity is negligible and it is at most of the longitudinal one. For the non cavitating cases at , see the behaviour of the densities, the shape of the velocity profiles is very close to that of a parabola since the flow is laminar, even if they are not perfectly symmetric due to the geometry used. A larger due to smaller causes a wider velocity profile which in turn increases the effective flow cross section and consequently . On the other hand, when cavitation occurs, i.e. the cases with equal to and , the velocity shows a typical linear shear flow profile in the vapor region which persists, with a different slope, moving towards the inner of the system. For both the cavitating cases considered here the ratio between the vapor and the liquid velocity slopes, across the interface, is roughly and it can be compared with the theoretically expected value Patanastasiou, Georgiou, and Alexandrou (1999) in a system with two phases of different density separated by a horizontal interface under shear flow. It is given by the ratio between liquid and vapor density which here is approximately .

A way to evaluate effects is also looking at the boundary layer thickness which is significant to have an idea about the reduction of the effective flow cross section due to viscous forces. Smaller boundary layers means larger effective flow cross section and therefore larger . Before the onset of cavitation the presence of the boundary layer in addition to compressibility effects are the main contributions to the discharge coefficient reduction. For the non-cavitating cases we can quantify the boundary layer by using the notion of displacement thickness Schlichting (1979), here defined as


where is the maximum value of the longitudinal velocity profile . These are equal to  lu and  lu when the downstream pressure is and respectively, thus showing that by increasing the boundary layer decreases and therefore increases. This analysis indicates that the mass flow rate and therefore the discharge coefficient are affected by different factors. Cavitation causes the most significant reduction of the effective flow cross section. It is further reduced by the boundary layer, which decreases with the Reynolds number. Finally, the compressibility causes a density reduction which in turn reduces the discharge coefficient. This compressibility effect is, however, quite small compared to the previous two.

Vi Conclusions

We conducted an in-depth numerical investigation of the formation of vapor in a constrained liquid flow using an enhanced forcing-term based lattice Boltzmann model. To model the liquid-vapor transition a free-energy approach based on the van der Waals equation of state was employed. The geometry of the system is a simple two-dimensional channel with a sack-wall obstacle positioned at the top of the channel and well behind the inlet layer inside the system. The channel is thus constrained to half its original width. Two different sets of inflow and outflow boundary conditions were implemented. One where density and velocity at the inlet layer and only the density at the outlet are fixed, and another fixing total pressure at the inflow layer and the static pressure at the outflow layer. The former approach allows for easier control of the velocity profile within the system and was used to investigate the circumstances of local cavitation formation. The latter approach on the other hand better reproduces the situation of experimental work focused on macroscopic hydraulic characterization of cavitation.

Liquid was forced into the channel from the left to flow past the obstacle. At a sufficient inlet velocity or pressure difference between inlet and outlet boundaries vapor formation underneath the obstacle was observed.

Careful investigation of the local circumstances of the formation of the vapor cavity indicates that it is necessary to reach the liquid branch spinodal pressure. The obtained data imply that this is a necessary but not sufficient condition. The density goes well below the liquid branch spinodal value without inducing the formation of a vapor bubble. Our observations then indicate that additional factors, such as viscous stress, interfacial contributions to the local pressure, and the Laplace pressure, are relevant to the opening of a vapor cavity. This is well described by an appropriate generalization of Joseph’s Joseph. (1998) criterion that includes these contributions. Cavitation occurs when the magnitude of the maximum stress eigenvalue reaches a threshold given by the spinodal pressure corrected by the Laplace contribution. The results are in agreement with those of Falcucci et al. Falcucci et al. (2013a) but we stress they are obtained from a LBM implementation based on a different thermodynamic model.

This investigation is also the first application of the LBM to the engineering practice of characterizing cavitation through a macroscopic approach by analyzing the behavior of stationary mass flow rate versus pressure drop and discharge coefficient versus both cavitation and Reynolds number. The outcomes found here are qualitatively in good agreement with the experimental ones of Ref. Payri et al., 2013 and indicate that the mass flow rate and therefore the discharge coefficient are affected by cavitation but also by other factors such as boundary layer and compressibility. As theory predicts, a linear increase of the mass flow rate with the square root of the pressure drop was found. However, the occurrence of cavitation causes a reduction in the mass flow growth rate and eventually the system collapses into a choked flow state. Reduction of the mass flow growth rate coincides with a smaller discharge coefficient which shows two different behaviors in the non-cavitating and cavitating regimes. In the non-cavitating regime, the discharge coefficient grows with the Reynolds number since the reduction of the boundary layer with causes larger effective flow cross section. In the case of cavitation, however, vapor causes a significant reduction of the effective liquid flow cross section and therefore a reduction of the discharge coefficient. In this regime, as predicted by NurickNurick (1976), the discharge coefficient grows with the cavitation number . However, agreement with Nurick law is, at this point, qualitative, since we found the predicted square-root dependence of the discharge coefficient with only for a specific value of viscosity. We suspect this is due to Reynolds number effects not taken into account in Nurick’s one-dimensional model.

For future research we will investigate the impact of wetting properties of the walls on the cavity formation process and the morphology of cavities.

We warmly thank Giacomo Falcucci for detailed discussion of his work.

Appendix A Boundary conditions

We outline here some details of the implementation of the numerical model with respect to boundary conditions.

a.1 Inlet and Outlet Boundary Conditions

Two different sets of inflow boundary conditions (BC) are employed in this study. The first one is called fixed-density BC and the second one is referred to as fixed-pressure BC. For the fixed-density BC a constant density at the inflow and a fixed inflow velocity are used. The distribution functions at the inflow layer are then calculated as equilibrium distribution functions according to Eq. (2) such that . Additionally, in order to avoid instabilities at the corners the horizontal inflow velocity is attenuated near the corners according to


Calculation of the density at the outflow layer follows a similar approach with a zero gradient on the velocity with a fixed outflow density . Copying the populations from the layer on the left of the outflow layer and rescaling them by the density, is a simple way of solving this problem:


The fixed-pressure BC are introduced to ensure compatibility with experimental work Payri et al. (2013); Nurick (1976); Winklhofer et al. () and to obtain the best possible estimate for the cavitation number and its various definitions. For this approach the pressure is fixed at both the inlet and the outlet boundaries. At the inlet, the total upstream pressure is fixed. The velocity at the inlet in the -direction is calculated from a zero-gradient velocity, i.e. velocities are taken from the layer adjacent to the right of the inlet layer, . The -velocity is fixed to vanish . At initialization of the simulation the fluid is at rest and the density is typically chosen to be slightly below the equilibrium liquid density.

At the inlet layer we thus identify the definition of the total pressure with the inlet pressure, i.e. . Replacing with the equation of state of Eq. (5) and rearranging the terms we obtain a third order polynomial in


The largest real solution of this expression represents the liquid branch density that solves the equation of state for an effective van der Waals pressure of . This density value is then used to calculate the distribution functions at the boundary layer. This density is streamed to the adjacent layer to the right causing a development of a velocity profile along the x-direction assuming that .

For the outflow a similar approach is used. The exception is that now only the pressure due to the equation of state enters the boundary condition at the outlet and we have Eq. (28) but with . Solving again for the largest density fixes at the outflow layer. This density, together with the zero-velocity gradient, is used in Eq. (27) to determine the distribution functions at the outflow layer. Including the dynamic pressure contribution in the inlet but not in the outlet pressure may seem counter-intuitive but this is done specifically to reproduce the upstream and downstream pressure definitions of the experimental work in Ref. Nurick, 1976. Thus the corresponding subsequent definition of the cavitation number in Eq. (18) also taken fromNurick (1976) is maintained.

The effect of the fixed-pressure BC at the outlet is not very significant in the current configuration. However, at the inlet, unlike the fixed profile of Eq. (26) in the fixed-density case, the velocity is no longer fixed. Instead the velocity may vary along the inflow layer and even backflow can occur. This has three practical effects for the numerical study: The simulation tends to be less stable during initial formation of the flow profile, to obtain cavitation mean pressure at the inlet layer needs to be chosen higher than in the fixed-density case, and flow profiles and vapor patterns look slightly different.

a.2 Wall Boundary Conditions

For the DQ simulation we give here explicit examples for the wall boundary condition implementation for all cases: flat wall, convex corner, and concave corner. In this implementation the base velocity set is given by for , respectively.

Let be the outgoing distribution functions in a wall lattice site at time and be those ones streamed from neighboring sites at time . For the bottom wall the local fluid density is then given by


Mass and momentum conservation require


which can be solved for the missing populations if one assumes the bounce-back prescription for the population normal to the wall. It then results in


The wall velocity is included here to illustrate how to introduce shear.

In the case of the obstacle-wall concave angle, here on the example of a corner with walls at the top and to the right, the local fluid density is given by


Mass and momentum conservation laws and bounce back rules ( and ) give the unknown distribution functions


In the case of the obstacle convex angle the local fluid density is given by


Mass and momentum conservation gives the unknown distribution functions


In order to ensure a higher isotropy and reduce spurious velocities at interfaces, which are unavoidable also in our model (see Ref. Coclite, Gonnella, and Lamura, 2014 for details), we calculated the spatial derivatives in (4) by using a general 9-point second-order finite difference scheme (see Refs. Shan, 2006; Sbragaglia et al., 2007 for details).

Appendix B Spinodal decomposition under shear

velocity field
Table 4: Tests of spinodal decomposition under shear in a channel of width at , , and shear rates between and . Indicated are also the type of shear implementation in the velocity field colum, the highest values of the mean initial density at which phase separation was observed, and the corresponding time after which liquid and vapor were noted.

The nucleation process responsible for the formation of vapor shown in this paper is confined to a relatively small region at the inlet of the constrained channel. This region is subjected to shear. It is known that phase transition behavior may change if the material is shearedWagner and Yeomans (1999); Onuki (2002); Corberi, Gonnella, and Lamura (1998, 1999). Thus it is prudent to investigate whether the shear is relevant in the ranges of the parameters where cavitation was observed in measurements presented in this paper. We evaluate if and how shear impacts the spinodal density of the fluid. In the simulation discussed in section IV the shear rates underneath the obstacle near the cavity inception point are found to be . To emulate comparable circumstances we set and to see the effect of increased shear . Two simulation scenarios were prepared. The first one, abbreviated as ’free’ in table 4, is a channel of width with periodic boundaries in the x-direction and on-grid bounceback on the upper and lower channel walls. These wall boundaries are constructed according to Eqs. (33) - (35). The lower wall then is assigned a wall velocity , and the upper wall thus ensuring in stationary conditions an average shear rate of . The second type of preparation, labeled ’fixed’, refers to a configuration where in addition to the boundary walls the fluid velocities are constrainted according to the imposed linear shear profile at the time of each collision. These velocities enter the equilibrium distribution functions in the collision term Eq. (1). Material can still move through advected differences between the local and equilibrium distribution functions but the phase separation dynamics is observed to be slower. The simulation is then initialized with a mean density with a very small random offset per lattice node of amplitude of where is a random number chosen from a flat distribution with zero mean and variance of . The initial density is very slightly below the spinodal density measured in simulations of table 2. Then, if no phase separation is observed after simulation updates, the experiment is repeated at a different initial density which is lower than the previous one. In order to reproduce the conditions at which cavitation formation is observed underneath the obstacle, we choose a channel width of . This is large enough to avoid interfacial effects (the interface width is , see, for example, Ref. Wagner and Li, 2006) and it is comparable to the width of the boundary layer which is in the right end of the channel (see Figs. 2 and 14) and is less wide near the inlet of the constrained channel. The results of these measurements are given in table 4, where the highest initial densities for which phase separation is observed, are reported. The observed difference between the measured spinodal density of table 2 and the largest effective density where the system phase separates at is . This difference is significantly smaller than the one between analytic spinodal density and the density at cavitation inception which is .


  1. preprint: AIP/123-QED
  2. The final increase in the stress seen before cavitation is due to the contribution of the van der Waals pressure that, following the van der Waals isotherm, increases when the fluid density decrease below the liquid spinodal value.


  1. C. E. Brennen, Cavitation and Bubble Dynamics (Oxford University Press, 1995).
  2. S. Dabiri, W. A. Sirignano,  and D. D. Joseph, “Cavitation in an orifice flow,” Physics of Fluids 19, 072112 (2007).
  3. T. Koivula, “On cavitation in fluid power,” in Proceedings of the 1st FPNI PhD Symposium, Hamburg, Germany, September 20–22, 2000, edited by M. Ivantysynova (Technical University of Hamburg-Harburg) pp. 371–382.
  4. D. Lohse, “Cavitation hots up,” Nature 434, 33 (2005).
  5. O. Reynolds, “The causes of the racing of the engines of screw steamers investigated theoretically and by experiment,” Transactions of the Royal Institution of Naval Architects 14, 56 (1873).
  6. C. A. Parsons, “The steam turbine on land and at sea,” Lecture to the Royal Institution, London  (1906).
  7. Y. Kawanami, H. Kato, H. Yamaguchi, M. Tanimura,  and Y. Tagaya, “Mechanism and Control of Cloud Cavitation,” Journal of Fluids Engineering 119, 788 (1997).
  8. C. Weiland and P. P. Vlachos, “Time-scale for critical growth of partial and supercavitation development over impulsively translating projectiles,” International Journal of Multiphase Flow 38, 73 (2012).
  9. B. Huang, A. Ducoin,  and Y. L. Young, “Physical and numerical investigation of cavitating flows around a pitching hydrofoil,” Physics of Fluids 25, 102109 (2013).
  10. E. Winklhofer, E. Kull, E. Kelz,  and A. Morozov, “Comprehensive hydraulic and flow field documentation in model throttle experiments under cavitation conditions,” in Proceedings 17th ANNUAL CONFERENCE ON LIQUID ATOMIZATION & SPRAY SYSTEMS, Zurich, Switzerland, September 2–6, 2001, edited by B. Ineichen (I.C. Engines and Combustion Laboratory at ETH, Zurich).
  11. R. Payri, F. J. Salvador, J. Gimeno,  and O. Venegas, “Study of cavitation phenomenon using different fuels in a transparent nozzle by hydraulic characterization and visualization,” Experimental Thermal and Fluid Science 44, 235 (2013).
  12. A. Sou, S. Hosokawa,  and A. Tomiyama, “Effects of cavitation in a nozzle on liquid jet atomization,” International Journal of Heat and Mass Transfer 50, 3575 (2007).
  13. A. Sou, A. Tomiyama, S. Hosokawa, S. Nigorikawa,  and T. Maeda, “Cavitation in a two-dimensional nozzle and liquid jet atomization (LDV measurement of liquid velocity in a nozzle),” JSME International Journal Series B 49, 1253 (2006).
  14. C. Mishra and Y. Peles, “Cavitation in flow through a micro-orifice inside a silicon microchannel,” Physics of Fluids 17, 013601 (2005).
  15. N. Tamaki, M. Shimizu, K. Nishida,  and H. Hiroyasu, “Effects of cavitation and internal flow on atomization of a liquid jet,” Atomization Sprays 8, 179 (1998).
  16. H. Hiroyasu, “Spray breakup mechanism from the hole-type nozzle and its applications,” Atomization Sprays 10, 511 (2000).
  17. F. Payri, V. Bermúdez, R. Payri,  and F. J. Salvador, “The influence of cavitation on the internal flow and the spray characteristics in diesel injection nozzles,” Fuel 83, 419 (2004).
  18. M. G. De Giorgi, A. Ficarella,  and M. Tarantino, “Evaluating cavitation regimes in an internal orifice at different temperatures using frequency analysis and visualization,” International Journal of Heat and Fluid Flow 39, 160 (2013).
  19. M. Gavaises, A. Andriotis, D. Papoulias, N. Mitroglou,  and A. Theodorakakos, “Characterization of string cavitation in large-scale Diesel nozzles with tapered holes,” Physics of Fluids 21, 052107 (2009).
  20. E. Giannadakis, M. Gavaises,  and C. Arcoumanis, “Modelling of cavitation in diesel injector nozzles,” Journal of Fluid Mechanics 616, 153 (2008).
  21. K. K. Kuo and R. Acharya, Fundamentals of turbulent and multiphase combustion (John Wiley & Sons, Inc., 2012).
  22. M. S. Chung, S. K. Pak,  and K. S. Chang, “A numerical study of two-phase flow using a two-dimensional two-fluid model,” Numerical Heat Transfer, Part A 45, 1049 (2004).
  23. Y. Y. Niu and Y. M. Lin, “A multiscale model to predict slurry erosion in cavitated duct flows,” Numerical Heat Transfer, Part A 50, 545 (2006).
  24. M. Darbandi and H. Sadeghi, “Numerical simulation of orifice cavitating flows using two-fluid and three-fluid cavitation models,” Numerical Heat Transfer, Part A 58, 505 (2010).
  25. G. S. Yeom and K. S. Chang, “Numerical simulation of two-fluid two-phase flows by hll scheme using an approximate jacobian matrix,” Numerical Heat Transfer, Part B 49, 155–177 (2006).
  26. C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics 39, 201 (1981).
  27. G. Wang, I. Senocak, W. Shyy, T. Ikohagi,  and S. Cao, “Dynamics of attached turbulent cavitating flows,” Progress in Aerospace Sciences 37, 551 (2001).
  28. Y. Utturkar, J. Wu, G. Wang,  and W. Shyy, “Recent progress in modeling of cryogenic cavitation for liquid rocket propulsion,” Progress in Aerospace Sciences 41, 558 (2005).
  29. C. C. Tseng, Y. Wei, G. Wang,  and W. Shyy, “Modeling of turbulent, isothermal and cryogenic cavitation under attached conditions,” Acta Mechanica Sinica 26, 325 (2010).
  30. F. J. Salvador, J. Martìnez-López, J. V. Romero,  and M. D. Roselló, “Computational study of the cavitation phenomenon and its interaction with the turbulence developed in diesel injector nozzles by Large Eddy Simulation (LES),” Mathematical and Computer Modelling 57, 1656 (2013).
  31. E. Goncalvès and B. Charrière, “Modelling for isothermal cavitation with a four-equation model,” International Journal of Multiphase Flow 59, 54 (2014).
  32. R. T. Knapp, J. W. Daily,  and F. G. Hammitt, Cavitation (McGraw-Hill, 1970).
  33. Y. Lecoffre, Cavitation bubble trackers (CRC Press, 1999).
  34. D. D. Joseph., “Cavitation and the state of stress in a flowing liquid,” Journal of Fluid Mechanics 366, 367 (1998).
  35. A. Bray, “Theory of Phase-Ordering Kinetics,” Adv. Phys. 43, 357–459 (1994).
  36. A. Onuki, Phase Transition Dynamics, 1st ed. (Cambridge University Press, 2002).
  37. S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, 1st ed. (Clarendon, Oxford, 2001).
  38. S. Succi, “Lattice boltzmann 2038,” EPL (Europhysics Letters) 109, 50001 (2015).
  39. R. Benzi, S. Succi,  and M. Vergassola, “The lattice Boltzmann equation: theory and applications,” Physics Reports 222, 145–197 (1992).
  40. G. Gonnella and J. M. Yeomans, “Using the lattice Boltzmann algorithm to explore phase ordering in fluids.” in Kinetics of Phase Transitions, edited by S. Puri (CRC, Boca Raton, 2009) pp. 89–166.
  41. B. Dünweg and A. Ladd, “Lattice boltzmann simulations of soft matter systems,” in Advanced Computer Simulation Approaches for Soft Matter Sciences III, Advances in Polymer Science, Vol. 221, edited by C. Holm and K. Kremer (Springer Berlin / Heidelberg, 2009) pp. 89–166.
  42. V. Sofonea, A. Lamura, G. Gonnella,  and A. Cristea, “Finite-difference Lattice Boltzmann Model with Flux Limiters for Liquid-vapor Systems,” Phys. Rev. E 70, 046702 (2004).
  43. A. Cristea, G. Gonnella, A. Lamura,  and V. Sofonea, “Finite-difference Lattice Boltzmann Model for Liquid-vapor Systems,” Math. Comput. Simul. 72, 113 (2006).
  44. G. Gonnella, A. Lamura,  and V. Sofonea, “Lattice Boltzmann simulation of thermal nonideal fluids,” Phys. Rev. E 76, 036703 (2007).
  45. G. Gonnella, A. Lamura,  and V. Sofonea, “A lattice Boltzmann method for thermal nonideal fluids,” Eur. Phys. J.-Spec. Top. 171, 181 (2009).
  46. A. Cristea, G. Gonnella, A. Lamura,  and V. Sofonea, “A lattice Boltzmann study of phase separation in liquid-vapor systems with gravity,” Commun. Comput. Phys. 7, 350 (2010).
  47. G. Falcucci, E. Jannelli, S. Ubertini,  and S. Succi, “Direct numerical evidence of stress-induced cavitation,” Journal of Fluid Mechanics 728, 362 (2013a).
  48. G. Falcucci, S. Ubertini, G. Bella,  and S. Succi, “Lattice Boltzmann simulation of cavitating flows,” Commun. Comput. Phys. 13, 685–695 (2013b).
  49. M. C. Sukop and D. Or, “Lattice boltzmann method for homogeneous and heterogeneous cavitation,” Phys. Rev. E 71, 046703 (2005).
  50. X.-P. Chen, C.-W. Zhong,  and X.-L. Yuan, “Lattice boltzmann simulation of cavitating bubble growth with large density ratio,” Computers and Mathematics with Applications 61, 3577 – 3584 (2011).
  51. M. Zhong, C. Zhong,  and C. Bai, “A high-order discrete scheme of lattice boltzmann method for cavitation simulations.” Adv. Comput. Sci. Appls 1, 73–77 (2012).
  52. G. Falcucci, G. Bella, G. Chiatti, S. Chibbaro, M. Sbragaglia,  and S. Succi, “Lattice Boltzmann models with mid-range interactions,” Commun. Comput. Phys. 2, 1071–1084 (2007).
  53. X. Shan and H. Chen, “Lattice Boltzmann model for simulating flows with multiple phases and components,” Phys. Rev. E 47, 1815–1819 (1993).
  54. A. Coclite, G. Gonnella,  and A. Lamura, “Pattern formation in liquid-vapor systems under periodic potential and shear,” Physical Review E 89, 063303 (2014).
  55. X. Shan, X. Yuan,  and H. Chen, “Kinetic theory representation of hydrodynamics: a way beyond the Navier Stokes equation,” J. Fluid Mech. 550, 413–441 (2006).
  56. Z. Guo, C. Zheng,  and B. Shi, “Discrete lattice effects on the forcing term in the lattice boltzmann method,” Physical Review E 65, 46308 (2002).
  57. A. Tiribocchi, N. Stella, G. Gonnella,  and A. Lamura, “Hybrid lattice Boltzmann model for binary fluid mixtures,” Phys. Rev. E 80, 026701 (2009).
  58. G. Gonnella, A. Lamura, A. Piscitelli,  and A. Tiribocchi, “Phase separation of binary fluids with dynamic temperature,” Phys. Rev. E 82, 046302 (2010).
  59. G. Gonnella, A. Lamura,  and A. Tiribocchi, “Thermal and hydrodynamic effects in the ordering of lamellar fluids,” Philos. Trans. R. Soc. A 369, 2592 (2011).
  60. J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon, Oxford, 1982).
  61. X. Shan, “Analysis and reduction of the spurious current in a class of multiphase lattice Boltzmann models,” Phys. Rev. E 73, 047701 (2006).
  62. M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama,  and F. Toschi, “Generalized lattice Boltzmann method with multirange pseudopotential,” Phys. Rev. E 75, 026702 (2007).
  63. W. H. Nurick, “Orifice cavitation and its effect on spray mixing,” J Fluid Eng 98, 681–687 (1976).
  64. A. Lamura and G. Gonnella, “Lattice boltzmann simulations of segregating binary fluid mixtures in shear flow,” Physica A 294, 295 (2001).
  65. A. Tiribocchi, A. Piscitelli, G. Gonnella,  and A. Lamura, “Pattern study of thermal phase separation for binary fluid mixtures,” Int. J. Numer. Methods Heat Fluid Flow 21, 572 (2011).
  66. A. J. Wagner and C. M. Pooley, “Interface width and bulk stability: Requirements for the simulation of deeply quenched liquid-gas systems,” Phys. Rev. E 76, 045702(R) (2007).
  67. G. Gonnella and M. Pellicoro, “Critical temperatures in driven binary mixtures with conserved and non-conserved dynamics,” Journal of Physics. A 33, 7043–7051 (2000).
  68. F. Corberi, G. Gonnella,  and A. Lamura, “Spinodal decomposition of binary mixtures in uniform shear flow,” Phys. Rev. Lett. 81, 3852 (1998).
  69. The final increase in the stress seen before cavitation is due to the contribution of the van der Waals pressure that, following the van der Waals isotherm, increases when the fluid density decrease below the liquid spinodal value.
  70. T. Patanastasiou, G. Georgiou,  and A. Alexandrou, Viscous Fluid Flow (CRC, Boca Raton, 1999).
  71. H. Schlichting, Boundary layer theory, 7th ed. (McGraw-Hill Book Company, New York, 1979).
  72. A. J. Wagner and J. M. Yeomans, “Phase separation under shear in two-dimensional binary fluids,” Phys. Rev. E 59, 4366–4373 (1999).
  73. F. Corberi, G. Gonnella,  and A. Lamura, “Two-scale competition in phase separation with shear,” Phys. Rev. Lett. 83, 4057 (1999).
  74. A. J. Wagner and Q. Li, “Investigation of galilean invariance of multi-phase lattice boltzmann methods,” Physica A 362, 105–110 (2006).
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.