JET  43
PREDICTION OF CAVITATION AND 
PARTICLE EROSION IN A RADIAL 
DIVERGENT TEST SECTION
NAPOVED KAVITACIJSKE EROZIJE 
IN EROZIJE DELCEV V RADIALNO 
DIVERGENTNI TESTNI SEKCIJI
Luka Kevorkijan
1
, Luka Lešnik
1
, Ignacijo Biluš
1ℜ
Keywords: Cavitation, Particles ANSYS Fluent, Erosion, CFD, modelling, DPM
Abstract
The 3D unsteady, cavitating, particle-laden flow through a radial divergent test section was 
simulated with the homogeneous mixture model and Discrete Phase Model (DPM) within the 
commercial CFD code ANSYS Fluent. For turbulence, a RANS approach was adopted with the 
Reboud’s correction of turbulent viscosity in the k-ω SST model. Cavitation erosion was predicted 
with the Schenke-Melissaris-Terwisga (SMT) model, while particle erosion was predicted with 
the Det Norske Veritas (DNV) model. Two distinct erosion zones were identified, one for pure 
cavitation erosion and one for pure particle erosion. The occurrence of the pure particle erosion 
zone downstream of the cavitation erosion zone was analysed. By observing the streamlines 
downstream of the cavitation structures, it was found that vortices form in the flow and redirect 
the particles towards the wall, causing a pure particle erosion zone on the wall. Particles under 
consideration in this study were not found to alter the flow to the extent that the cavitation ero-
sion zone would be significantly altered compared with the results without solid particles which 
are reported in the literature.
JET Volume 15 (2022) p.p. 43-64
Issue 2, 2022
Type of article: 1.01 
www.fe.um.si/si/jet.html
 
ℜ Corresponding author: Associate Professor, Ignacijo Biluš, University of Maribor, Faculty of Mechanical Engineer-
ing, Smetanova ulica 17, Maribor, Slovenia, Tel.: +386 (2) 220 7742, E-mail address: ignacijo.bilus@um.si
1
 University of Maribor, Faculty of Mechanical Engineering, Smetanova ulica 17, Maribor, Slovenia
Prediction of cavitation and particle erosion in a radial divergent test section
Napoved kavitacijske erozije in erozije delcev v radialno divergentni testni sekciji
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš

44 JET 44 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
Povzetek
3D nestacionaren, kavitirajoč tok z delci skozi radialno divergentno testno sekcijo je bil simuliran 
z modelom homogene zmesi in modelom diskretne faze (DPM) s komercialno RDT kodo ANSYS 
Fluent. Turbulenca je modelirana po pristopu RANS z Reboudovim popravkom turbulentne visko-
znosti v modelu k-ω SST. Izvedeni sta bili napoved kavitacijske erozije z modelom Schenke-Me-
lissaris-Terwisga (SMT) in napoved erozije zaradi delcev z modelom Det Norske Veritas (DNV). 
Prepoznani sta bili dve različni erozijski coni, ena za zgolj kavitacijsko erozijo in ena za erozijo zgolj 
zaradi delcev. Analiziran je bil pojav cone čiste erozije zaradi delcev dolvodno od cone kavitacijske 
erozije. Z opazovanjem tokovnic dolvodno od kavitacijskih struktur je bilo ugotovljeno, da se v 
toku oblikujejo vrtinci in preusmerjajo delce proti steni, kar povzroča na steni cono erozije zgolj 
zaradi delcev. Ugotovljeno je bilo, da delci, obravnavani v tej študiji, ne spreminjajo toka do te 
mere, da bi se območje kavitacijske erozije znatno spremenilo v primerjavi z rezultati brez trdnih 
delcev, o katerih poroča literatura.
1 INTRODUCTION
Erosion of solid surfaces can occur in hydraulic systems as a result of multiple causes, including 
cavitation and solid particles.
Cavitation is a phase change phenomenon that can arise when the pressure of a liquid drops 
to vaporisation pressure. As a result of this, vapour cavities (bubbles, or more generally vapour 
structures) form inside of the liquid. When vapour cavities are subjected to a pressure that ex-
ceeds the vaporisation pressure, condensation takes place and they collapse.
After the collapse of a vapour cavity, an instantaneous pressure peak occurs in the liquid. For 
spherical bubble collapse, Lord Rayleigh [1] derived the expression for the time evolution of the 
pressure in a liquid during the collapse. The pressure wave emitted after the collapse can cause 
the erosion of a solid surface as it impacts that surface. Hammit [2] introduced the hypothesis, 
based on experimental observations of cavitation damage in a Venturi channel, that erosion is 
related to the initial potential energy spectrum of vapour cavities. Vogel et al. [3] later derived 
the equation for potential energy of cavitation bubbles based on the observations of a single 
bubble collapse near a solid wall.
However, when a bubble collapses in the vicinity of a solid wall, the collapse becomes asymmet-
ric due to the presence of the wall on one side of the bubble. As a result of an asymmetric col-
lapse, a liquid microjet develops and impacts the solid surface. Based on this observation, Plesset 
and Chapman [4] calculated the velocity of the impacting microjet and concluded that this can be 
the cause of solid surface erosion. Whether the pressure wave, microjet, or some combination of 
both causes cavitation erosion in larger flow scales (complex macroscopic cavitation structures as 
opposed to a single bubble) remains an open area of research.
In engineering applications, erosion is undesired as it can reduce efficiency and ultimately cause 
failure of a given hydraulic system, for example hydraulic turbines, pumps, fuel injectors in in-
ternal combustion engines, pipe elbows and valves. It is therefore important to incorporate the 
study of erosion risk in the design process of a hydraulic system. In the past, most of the insight 
into the cavitation erosion process was obtained by experiments. Recently, however, approaches 
based on computational fluid dynamics (CFD) have emerged.

JET  45
Prediction of cavitation and particle erosion in a radial divergent test section
From the observations made by Hammit [2] and Vogel et al. [3], an energy cascade mechanism 
has been proposed by both Pereira et al. [5] and Fortes-Patella et al. [6] to explain the transfer 
of potential energy from large cavitation structures to the solid surface. Fortes-Patella et al. [6] 
introduced a cavitation erosion model based on an assumption that the erosive aggressiveness 
of large cavitation structures is proportional to their distance from the solid surface. Melissaris 
et al. [7] cavitation erosion prediction becomes more and more critical, as the requirements for 
more efficient propellers increase. Model testing is yet the most typical way a propeller designer 
can, nowadays, get an estimation of the erosion risk on the propeller blades. However, cavita-
tion erosion prediction using computational fluid dynamics (CFD analysed the applicability of 
this cavitation erosion model, with a simplification that potential energy is calculated only in 
the first cell layer next to a solid surface. This finding is supported by previous observations by 
Philipp and Lauterborn [8], who found that cavitation structures in contact with the solid surface 
are the most aggressive. Leclercq et al. [9] which can be described as the fluid mechanical loa-
ding leading to cavitation damage. The estimation of this quantity is a challenging problem both 
in terms of modeling the cavitating flow and predicting the erosion due to cavitation. For this 
purpose, a numerical methodology was proposed to estimate cavitation intensity from 3D un-
steady cavitating flow simulations. CFD calculations were carried out using Code_Saturne, which 
enables U-RANS equations resolution for a homogeneous fluid mixture using the Merkle’s mo-
del, coupled to a k-ε turbulence model with the Reboud’s correction. A post-process cavitation 
intensity prediction model was developed based on pressure and void fraction derivatives. This 
model is applied on a flow around a hydrofoil using different physical (inlet velocities included all 
cavitation structures in cavitation erosion prediction by projecting the release of potential energy 
from the cells that are not in contact with a solid surface via the solid angle. While Leclercq et 
al. calculated solid angles for individual cells using a discrete form written by Van Osteroom and 
Strackee [10], Schenke and van Terwisga [11] instead introduced a continuous form of the solid 
angle. Based on the description of a global energy balance of a bubble in an infinite liquid [12], 
and the observed focusing of potential energy into the centre of a bubble cloud collapse [13], 
Schenke et al. [14] have introduced a model to describe the conversion of cavitation potential 
energy into kinetic energy, which is then released as a pressure wave. Melissaris et al. [15] have 
analysed three mathematical forms of a vapour volume fraction total derivative used in the cal-
culation of cavitation potential energy, concluding that the least numerical error is introduced 
when using the form with mass source term. Melissaris et al. [15] and Pezdevšek et al. [16] also 
compared the two solid angle projection approaches.
The above-described cavitation erosion models have been developed with an important limita-
tion that both the vapour and the liquid phase are considered as a homogeneous mixture. This 
approach is widely used, and several cavitation models [17]–[19] exist and have been analysed 
[20]–[25]. However, an alternative approach is to track some if not all of the vapour phase as 
individual bubbles in a Lagrangian frame [26]–[32]. Hybrid (multiscale) approaches, where only 
sufficiently small cavitation structures are converted into bubbles and tracked in a Lagrangian 
frame, are particularly promising for cavitation erosion prediction.
When the solid particles present in the flow impact a solid surface, they too can cause erosion 
of that surface. Finnie [33] identified key parameters that determine the mass loss of a solid sur-
face that is impacted by solid particles. An expression that relates the solid particle impact angle 
to the erosion of the solid surface has been written. Similarly, Bitter [34], [35] later proposed a 
different impact angle dependency. There are many empirical expressions in the literature, some 
of the more often used are by Oka et al. [36], [37], Ahlert [38] and Det Norske Veritas (DNV) [39].

46 JET 46 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
Empirical particle erosion models can be used in CFD simulations in combination with Lagrangian 
particle tracking, where particles are assumed to be point-masses. When particles hit a solid 
surface, they rebound from that surface. The rebound behaviour is described using empirical 
expressions that relate the particle velocity after the rebound with the particle impact angle 
[40], [41].
In the past, researches usually focused their attention on one phenomenon only, either the cav-
itation erosion or the particle erosion. Few studies have been conducted, where both phenom-
ena are present at the same time. These are mostly experimental studies [42]–[44], but some 
studies apply CFD simulations as well [45]–[47].
In this work a CFD simulation of cavitation in a particle-laden flow has been conducted on a geom-
etry commonly used in cavitation erosion studies using a commercial CFD software ANSYS Fluent 
2021R2. Both the cavitation and the particle erosion are predicted using the Schenke-Melissar-
is-van Terwisga cavitation erosion model and the DNV particle erosion model respectively. The 
Schenke-Melissaris-van Terwisga cavitation erosion model was previously implemented within 
the ANSYS Fluent using User Defined Function (UDF), written in a C programming language [48]. 
The mesh was generated with a Fluent meshing module version 2020 R2.
2 SIMULATION OF CAVITATION IN A PARTICLE-LADEN FLOW
2.1 Geometry and mesh
Radial divergent test section [28] geometry (Fig. 1) was chosen, as it is often the choice of numer-
ical cavitation erosion studies [31], [32], [48]–[51] and because the present authors previously 
analysed the implemented cavitation erosion model on this geometry [48].
Figure 1: Radial divergent test section geometry and dimensions
The mesh was generated using ANSYS Fluent meshing module. To reduce the computational 
effort, only a 45° sector was meshed. The resulting poly-hexcore mesh is shown in Fig. 2 and 

JET  47
Prediction of cavitation and particle erosion in a radial divergent test section
the mesh statistics are presented in Table 1. The mesh nodalization study was conducted by the 
authors in the previous work [48].
Figure 2: Mesh of a 45° sector: a) whole domain, b) zoomed in on the refined part of the 
domain, c) zoomed in on the radius
Table 1: Mesh statistics
Mesh statistic Value
Number of cells 1,396,703
Minimum orthogonal quality 0.28
Maximum aspect ratio 65
Minimum y
+
 dimensionless wall distance 0.174
Maximum y
+ 
dimensionless wall distance 51.9
2.2 Governing equations
2.2.1 Homogeneous mixture
Vapour and liquid phases are treated as a homogeneous mixture, the mixture density (ρ) and the 
mixture viscosity (µ) are determined via the mixing rule
 Prediction of cavitation and particle erosion in a radial divergent test section 5 
   
---------- 
 
Figure 2: Mesh of a 45° sector: a) whole domain, b) zoomed in on the refined part of the domain, 
c) zoomed in on the radius. 
 
Table 1: Mesh statistics. 
Mesh statistic Value 
Number of cells 1,396,703 
Minimum orthogonal quality 0.28 
Maximum aspect ratio 65 
Minimum 𝑦𝑦 �
 dimensionless wall distance 0.174 
Maximum 𝑦𝑦 �
 dimensionless wall distance 51.9 
 
2.2 Governing equations 
2.2.1 Homogeneous mixture 
Vapour and liquid phases are treated as a homogeneous mixture, the mixture density (𝜌𝜌 ) and the 
mixture viscosity (𝜇𝜇 ) are determined via the mixing rule 
𝜌𝜌 = 𝛼𝛼𝜌𝜌
�
+ (1−𝛼𝛼 )𝜌𝜌 �
               (2.1) 
𝜇𝜇 = 𝛼𝛼𝜇𝜇
�
+ (1−𝛼𝛼 )𝜇𝜇 �
              (2.2) 
where 𝛼𝛼 is the vapour volume fraction, the liquid and vapour densities of water are 𝜌𝜌 �
=
998.85 kg m
3
⁄ and 𝜌𝜌 �
= 0.01389 kg m
3
⁄ . The liquid and vapour viscosities of water are 𝜇𝜇 �
=
0.0011 Pa·s and 𝜇𝜇 �
= 9.63 ∙ 10
��
 Pa ∙ s. 
(2.1)

48 JET 48 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
 Prediction of cavitation and particle erosion in a radial divergent test section 5 
   
---------- 
 
Figure 2: Mesh of a 45° sector: a) whole domain, b) zoomed in on the refined part of the domain, 
c) zoomed in on the radius. 
 
Table 1: Mesh statistics. 
Mesh statistic Value 
Number of cells 1,396,703 
Minimum orthogonal quality 0.28 
Maximum aspect ratio 65 
Minimum 𝑦𝑦 �
 dimensionless wall distance 0.174 
Maximum 𝑦𝑦 �
 dimensionless wall distance 51.9 
 
2.2 Governing equations 
2.2.1 Homogeneous mixture 
Vapour and liquid phases are treated as a homogeneous mixture, the mixture density (𝜌𝜌 ) and the 
mixture viscosity (𝜇𝜇 ) are determined via the mixing rule 
𝜌𝜌 = 𝛼𝛼𝜌𝜌
�
+ (1−𝛼𝛼 )𝜌𝜌 �
               (2.1) 
𝜇𝜇 = 𝛼𝛼𝜇𝜇
�
+ (1−𝛼𝛼 )𝜇𝜇 �
              (2.2) 
where 𝛼𝛼 is the vapour volume fraction, the liquid and vapour densities of water are 𝜌𝜌 �
=
998.85 kg m
3
⁄ and 𝜌𝜌 �
= 0.01389 kg m
3
⁄ . The liquid and vapour viscosities of water are 𝜇𝜇 �
=
0.0011 Pa·s and 𝜇𝜇 �
= 9.63 ∙ 10
��
 Pa ∙ s. 
(2.2)
where α is the vapour volume fraction, the liquid and vapour densities of water are ρ
l
 = 998.85 kg/m³ 
and ρ
v
 = 0.01389 kg/m³. The liquid and vapour viscosities of water are μ
l
 = 0.0011 Pa∙s and 
μ
v
 = 9.63 ∙ 10
−6
 Pa∙s.
The flow of mixture is governed by the mixture mass conservation equation
6 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
where α is the vapour volume fraction, the liquid and vapour densities of water are
ρ
l
= 998.85 kg / m
3
 and ρ
v
= 0.01389 kg / m
3
. The liquid and vapour viscosities of
water are μ
l
= 0.0011 P a∙ s and μ
v
= 9.63 ∙ 10
− 6
Pa ∙ s .
The flow of mixture is governed by the mixture mass conservation equation
∂ ρ
∂ t
+ ∇ ∙ ( ρ u ) = 0 , (2.3)
where u is the mixture velocity field, and by the mixture momentum conservation equation
∂ ( ρ u )
∂ t
+ ∇ ∙ ( ρ uu ) = − ∇ p + ∇ ∙ [ μ ( ∇ u + ∇ u
T
) ] + ρ g + S
M
, (2.4)
where p is the pressure, g is the gravitational acceleration and S
M
 is the
momentum source term due to the presence of particles in the flow. 
2.2.2 Turbulence modelling
The RANS approach was adopted and the SST k-w model was chosen. In order to reproduce the
transient cavitation effects that occur due to the compressibility of a liquid-vapour mixture, the
turbulent viscosity μ
t
 needs to be modified according to Reboud et al. [52] as follows:
μ
t
=
f ( ρ ) k
ω
1
m ax
[
1
α
❑
,
S F
2
a
1
ω ]
,
(2.5)
where k is the turbulence kinetic energy, ω is the specific rate of dissipation of the
turbulence kinetic energy, α
❑
 is the damping coefficient, S is the strain rate magnitude,
F
2
 is the second blending function and a
1
 is the model constant equaling 0.31.
Correction of turbulent viscosity is achieved by modifying the mixture density function f ( ρ )
such that a lower turbulent viscosity is achieved in the presence of both the liquid and vapour
phases
f ( ρ ) = ρ
v
+
(
ρ − ρ
v )
n
(
ρ
l
− ρ
v )
n − 1
, (2.6)
where n is an arbitrary exponent with a recommended value of 10.
This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF (User
Defined Function) [48].
2.2.3 Cavitation modelling
(2.3)
where u is the mixture velocity field, and by the mixture momentum conservation equation
6 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
where α is the vapour volume fraction, the liquid and vapour densities of water are
ρ
l
= 998.85 kg / m
3
 and ρ
v
= 0.01389 kg / m
3
. The liquid and vapour viscosities of
water are μ
l
= 0.0011 P a∙ s and μ
v
= 9.63 ∙ 10
− 6
Pa ∙ s .
The flow of mixture is governed by the mixture mass conservation equation
∂ ρ
∂ t
+ ∇ ∙ ( ρ u ) = 0 , (2.3)
where u is the mixture velocity field, and by the mixture momentum conservation equation
∂ ( ρ u )
∂ t
+ ∇ ∙ ( ρ uu ) = − ∇ p + ∇ ∙ [ μ ( ∇ u + ∇ u
T
) ] + ρ g + S
M
, (2.4)
where p is the pressure, g is the gravitational acceleration and S
M
 is the
momentum source term due to the presence of particles in the flow. 
2.2.2 Turbulence modelling
The RANS approach was adopted and the SST k-w model was chosen. In order to reproduce the
transient cavitation effects that occur due to the compressibility of a liquid-vapour mixture, the
turbulent viscosity μ
t
 needs to be modified according to Reboud et al. [52] as follows:
μ
t
=
f ( ρ ) k
ω
1
m ax
[
1
α
❑
,
S F
2
a
1
ω ]
,
(2.5)
where k is the turbulence kinetic energy, ω is the specific rate of dissipation of the
turbulence kinetic energy, α
❑
 is the damping coefficient, S is the strain rate magnitude,
F
2
 is the second blending function and a
1
 is the model constant equaling 0.31.
Correction of turbulent viscosity is achieved by modifying the mixture density function f ( ρ )
such that a lower turbulent viscosity is achieved in the presence of both the liquid and vapour
phases
f ( ρ ) = ρ
v
+
(
ρ − ρ
v )
n
(
ρ
l
− ρ
v )
n − 1
, (2.6)
where n is an arbitrary exponent with a recommended value of 10.
This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF (User
Defined Function) [48].
2.2.3 Cavitation modelling
(2.4)
where p is the pressure, g is the gravitational acceleration and S
M
 is the momentum source term 
due to the presence of particles in the flow.
2.2.2 Turbulence modelling
The RANS approach was adopted and the SST k-ω model was chosen. In order to reproduce the 
transient cavitation effects that occur due to the compressibility of a liquid-vapour mixture, the 
turbulent viscosity μ
t
 needs to be modified according to Reboud et al. [52] as follows:
6 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022) 
  Issue 2 
---------- 
The flow of mixture is governed by the mixture mass conservation equation 
𝜕𝜕 𝜌𝜌 𝜕𝜕𝜕𝜕
+ ∇ ∙ (𝜌𝜌 𝒖𝒖 ) = 0, 
(2.3)
where 𝒖𝒖 is the mixture velocity field, and by the mixture momentum conservation equation 
𝜕𝜕 (𝜌𝜌 𝒖𝒖 )
𝜕𝜕𝜕𝜕
+ ∇ ∙ (𝜌𝜌 𝒖𝒖𝒖𝒖 ) = −∇𝑝𝑝 + ∇ ∙ [𝜇𝜇 (∇𝒖𝒖 + ∇𝒖𝒖 �
)] + 𝜌𝜌 𝒈𝒈 + 𝑺𝑺 𝑴𝑴 , 
(2.4)
where 𝑝𝑝 is the pressure, 𝒈𝒈 is the gravitational acceleration and 𝑺𝑺 𝑴𝑴 is the momentum source term 
due to the presence of particles in the flow.  
 
2.2.2 Turbulence modelling 
The RANS approach was adopted and the SST k-  model was chosen. In order to reproduce the 
transient cavitation effects that occur due to the compressibility of a liquid-vapour mixture, the 
turbulent viscosity 𝜇𝜇 �
 needs to be modified according to Reboud et al. [52] as follows: 
𝜇𝜇 �
=
𝑓𝑓 (𝜌𝜌 )𝑘𝑘 𝜔𝜔 1
𝑚𝑚𝑚𝑚𝑚𝑚 �
1
𝛼𝛼 ∗
,
𝑆𝑆 𝐹𝐹 �
𝑚𝑚 �
𝜔𝜔 �
, 
(2.5)
where 𝑘𝑘 is the turbulence kinetic energy, 𝜔𝜔 is the specific rate of dissipation of the turbulence 
kinetic energy, 𝛼𝛼 ∗
 is the damping coefficient, 𝑆𝑆 is the strain rate magnitude, 𝐹𝐹 �
 is the second 
blending function and 𝑚𝑚 �
 is the model constant equaling 0.31. Correction of turbulent viscosity is 
achieved by modifying the mixture density function 𝑓𝑓 (𝜌𝜌 ) such that a lower turbulent viscosity is 
achieved in the presence of both the liquid and vapour phases 
𝑓𝑓 (𝜌𝜌 ) = 𝜌𝜌 �
+
(𝜌𝜌 − 𝜌𝜌 �
)
�
(𝜌𝜌 �
− 𝜌𝜌 �
)
�� �
, (2.6)
where 𝑛𝑛 is an arbitrary exponent with a recommended value of 10. 
This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF (User Defined 
Function) [48]. 
 
2.2.3 Cavitation modelling 
To model the mass transfer between the liquid and the vapour phase, an additional equation for 
conservation of vapour mass is solved  
𝜕𝜕 (𝛼𝛼 𝜌𝜌 �
)
𝜕𝜕𝜕𝜕
= ∇ ∙ (𝛼𝛼 𝜌𝜌 �
𝒖𝒖 ) = 𝜌𝜌 �
S
�
�
, 
(2.7)
where S
�
�
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used, which 
states that mass transfer source term S
�
�
 is 
(2.5)
where k is the turbulence kinetic energy, ω is the specific rate of dissipation of the turbulence 
kinetic energy, α* is the damping coefficient, S is the strain rate magnitude, F
2
 is the second 
blending function and a
1
 is the model constant equaling 0.31. Correction of turbulent viscosity is 
achieved by modifying the mixture density function f(ρ) such that a lower turbulent viscosity is 
achieved in the presence of both the liquid and vapour phases
6 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
where α is the vapour volume fraction, the liquid and vapour densities of water are
ρ
l
= 998.85 kg / m
3
 and ρ
v
= 0.01389 kg / m
3
. The liquid and vapour viscosities of
water are μ
l
= 0.0011 P a∙ s and μ
v
= 9.63 ∙ 10
− 6
Pa ∙ s .
The flow of mixture is governed by the mixture mass conservation equation
∂ ρ
∂ t
+ ∇ ∙ ( ρ u ) = 0 , (2.3)
where u is the mixture velocity field, and by the mixture momentum conservation equation
∂ ( ρ u )
∂ t
+ ∇ ∙ ( ρ uu ) = − ∇ p + ∇ ∙ [ μ ( ∇ u + ∇ u
T
) ] + ρ g + S
M
, (2.4)
where p is the pressure, g is the gravitational acceleration and S
M
 is the
momentum source term due to the presence of particles in the flow. 
2.2.2 Turbulence modelling
The RANS approach was adopted and the SST k-w model was chosen. In order to reproduce the
transient cavitation effects that occur due to the compressibility of a liquid-vapour mixture, the
turbulent viscosity μ
t
 needs to be modified according to Reboud et al. [52] as follows:
μ
t
=
f ( ρ ) k
ω
1
m ax
[
1
α
❑
,
S F
2
a
1
ω ]
,
(2.5)
where k is the turbulence kinetic energy, ω is the specific rate of dissipation of the
turbulence kinetic energy, α
❑
 is the damping coefficient, S is the strain rate magnitude,
F
2
 is the second blending function and a
1
 is the model constant equaling 0.31.
Correction of turbulent viscosity is achieved by modifying the mixture density function f ( ρ )
such that a lower turbulent viscosity is achieved in the presence of both the liquid and vapour
phases
f ( ρ ) = ρ
v
+
(
ρ − ρ
v )
n
(
ρ
l
− ρ
v )
n − 1
, (2.6)
where n is an arbitrary exponent with a recommended value of 10.
This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF (User
Defined Function) [48].
2.2.3 Cavitation modelling
(2.6)
where n is an arbitrary exponent with a recommended value of 10.
This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF (User De-
fined Function) [48].
2.2.3 Cavitation modelling
To model the mass transfer between the liquid and the vapour phase, an additional equation for 
conservation of vapour mass is solved
Prediction of cavitation and particle erosion in a radial divergent test section 7
----------
To model the mass transfer between the liquid and the vapour phase, an additional equation for
conservation of vapour mass is solved 
∂
(
α ρ
v )
∂ t
= ∇ ∙ ( α ρ
v
u ) = ρ
v
S
α v
,
(2.7)
where S
α v
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used,
which states that mass transfer source term S
α v
 is
S
α
v
=
{
ρ
v
ρ
l
ρ
α
v (1 − α
v )
3
R
b √
2
3
( p
v
− p )
ρ
l
     if p < p
v
ρ
v
ρ
l
ρ
α
v (
1 − α
v )
3
R
b
√
2
3
( p − p
v )
ρ
l
     if p > p
v
, (2.8)
where the bubble radius is
R
b
=
(
α
v
n
0
4
3
π ( 1 − α
v ) )
1 / 3
.
(2.9)
For the bubble number density n
0
 a default value of 
1 ∙10
11
 was used. The vapour
pressure was set to  p
v
=1854 Pa .
2.2.3 Cavitation erosion modelling
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in
terms of work done by the surrounding liquid on a bubble. In general , for a cavity with volume
V
V
, the potential energy of that cavity equals 
E
pot
=
(
p
d
− p
v )
V
V
, (2.10)
where p
d
 is the pressure in the surrounding liquid and p
v
 is the vapour pressure inside
the bubble. 
The cavitation potential power P
p ot
 is introduced as the total derivative of the potential
energy E
pot
 
P
pot
=
D
(
p
d
− p
v )
Dt
V
V
+
D V
V
Dt
( p
d
− p
v
) .
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of
magnitude lower than the second term [15], and taking into account that p
v
= c on s t .
(isothermal flow), the cavitation potential power is written as
(2.7)

JET  49
Prediction of cavitation and particle erosion in a radial divergent test section
where S
α
v
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used, which 
states that mass transfer source term S
α
v
 is
Prediction of cavitation and particle erosion in a radial divergent test section 7
----------
To model the mass transfer between the liquid and the vapour phase, an additional equation for
conservation of vapour mass is solved 
∂
(
α ρ
v )
∂ t
= ∇ ∙ ( α ρ
v
u ) = ρ
v
S
α v
,
(2.7)
where S
α v
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used,
which states that mass transfer source term S
α v
 is
S
α
v
=
{
ρ
v
ρ
l
ρ
α
v (1 − α
v )
3
R
b √
2
3
( p
v
− p )
ρ
l
     if p < p
v
ρ
v
ρ
l
ρ
α
v (
1 − α
v )
3
R
b
√
2
3
( p − p
v )
ρ
l
     if p > p
v
, (2.8)
where the bubble radius is
R
b
=
(
α
v
n
0
4
3
π ( 1 − α
v ) )
1 / 3
.
(2.9)
For the bubble number density n
0
 a default value of 
1 ∙10
11
 was used. The vapour
pressure was set to  p
v
=1854 Pa .
2.2.3 Cavitation erosion modelling
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in
terms of work done by the surrounding liquid on a bubble. In general , for a cavity with volume
V
V
, the potential energy of that cavity equals 
E
pot
=
(
p
d
− p
v )
V
V
, (2.10)
where p
d
 is the pressure in the surrounding liquid and p
v
 is the vapour pressure inside
the bubble. 
The cavitation potential power P
p ot
 is introduced as the total derivative of the potential
energy E
pot
 
P
pot
=
D
(
p
d
− p
v )
Dt
V
V
+
D V
V
Dt
( p
d
− p
v
) .
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of
magnitude lower than the second term [15], and taking into account that p
v
= c on s t .
(isothermal flow), the cavitation potential power is written as
(2.8)
where the bubble radius is
Prediction of cavitation and particle erosion in a radial divergent test section 7
----------
To model the mass transfer between the liquid and the vapour phase, an additional equation for
conservation of vapour mass is solved 
∂
(
α ρ
v )
∂ t
= ∇ ∙ ( α ρ
v
u ) = ρ
v
S
α v
,
(2.7)
where S
α v
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used,
which states that mass transfer source term S
α v
 is
S
α
v
=
{
ρ
v
ρ
l
ρ
α
v (1 − α
v )
3
R
b √
2
3
( p
v
− p )
ρ
l
     if p < p
v
ρ
v
ρ
l
ρ
α
v (
1 − α
v )
3
R
b
√
2
3
( p − p
v )
ρ
l
     if p > p
v
, (2.8)
where the bubble radius is
R
b
=
(
α
v
n
0
4
3
π ( 1 − α
v ) )
1 / 3
.
(2.9)
For the bubble number density n
0
 a default value of 
1 ∙10
11
 was used. The vapour
pressure was set to  p
v
=1854 Pa .
2.2.3 Cavitation erosion modelling
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in
terms of work done by the surrounding liquid on a bubble. In general , for a cavity with volume
V
V
, the potential energy of that cavity equals 
E
pot
=
(
p
d
− p
v )
V
V
, (2.10)
where p
d
 is the pressure in the surrounding liquid and p
v
 is the vapour pressure inside
the bubble. 
The cavitation potential power P
p ot
 is introduced as the total derivative of the potential
energy E
pot
 
P
pot
=
D
(
p
d
− p
v )
Dt
V
V
+
D V
V
Dt
( p
d
− p
v
) .
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of
magnitude lower than the second term [15], and taking into account that p
v
= c on s t .
(isothermal flow), the cavitation potential power is written as
(2.9)
For the bubble number density n
0
 a default value 1∙10
11
 of was used. The vapour pressure was 
set to p
v
 = 1854 Pa.
2.2.4 Cavitation erosion modelling
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in 
terms of work done by the surrounding liquid on a bubble. In general, for a cavity with volume 
V
V
, the potential energy of that cavity equals
Prediction of cavitation and particle erosion in a radial divergent test section 7
----------
To model the mass transfer between the liquid and the vapour phase, an additional equation for
conservation of vapour mass is solved 
∂
(
α ρ
v )
∂ t
= ∇ ∙ ( α ρ
v
u ) = ρ
v
S
α v
,
(2.7)
where S
α v
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used,
which states that mass transfer source term S
α v
 is
S
α
v
=
{
ρ
v
ρ
l
ρ
α
v (1 − α
v )
3
R
b √
2
3
( p
v
− p )
ρ
l
     if p < p
v
ρ
v
ρ
l
ρ
α
v (
1 − α
v )
3
R
b
√
2
3
( p − p
v )
ρ
l
     if p > p
v
, (2.8)
where the bubble radius is
R
b
=
(
α
v
n
0
4
3
π ( 1 − α
v ) )
1 / 3
.
(2.9)
For the bubble number density n
0
 a default value of 
1 ∙10
11
 was used. The vapour
pressure was set to  p
v
=1854 Pa .
2.2.3 Cavitation erosion modelling
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in
terms of work done by the surrounding liquid on a bubble. In general , for a cavity with volume
V
V
, the potential energy of that cavity equals 
E
pot
=
(
p
d
− p
v )
V
V
, (2.10)
where p
d
 is the pressure in the surrounding liquid and p
v
 is the vapour pressure inside
the bubble. 
The cavitation potential power P
p ot
 is introduced as the total derivative of the potential
energy E
pot
 
P
pot
=
D
(
p
d
− p
v )
Dt
V
V
+
D V
V
Dt
( p
d
− p
v
) .
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of
magnitude lower than the second term [15], and taking into account that p
v
= c on s t .
(isothermal flow), the cavitation potential power is written as
(2.10)
where p
d
 is the pressure in the surrounding liquid and p
v
 is the vapour pressure inside the bubble.
The cavitation potential power p
pot
 is introduced as the total derivative of the potential energy 
E
pot
Prediction of cavitation and particle erosion in a radial divergent test section 7
----------
To model the mass transfer between the liquid and the vapour phase, an additional equation for
conservation of vapour mass is solved 
∂
(
α ρ
v )
∂ t
= ∇ ∙ ( α ρ
v
u ) = ρ
v
S
α v
,
(2.7)
where S
α v
 is the source of the vapour phase. The Schnerr-Sauer cavitation model was used,
which states that mass transfer source term S
α v
 is
S
α
v
=
{
ρ
v
ρ
l
ρ
α
v (1 − α
v )
3
R
b √
2
3
( p
v
− p )
ρ
l
     if p < p
v
ρ
v
ρ
l
ρ
α
v (
1 − α
v )
3
R
b
√
2
3
( p − p
v )
ρ
l
     if p > p
v
, (2.8)
where the bubble radius is
R
b
=
(
α
v
n
0
4
3
π ( 1 − α
v ) )
1 / 3
.
(2.9)
For the bubble number density n
0
 a default value of 
1 ∙10
11
 was used. The vapour
pressure was set to  p
v
=1854 Pa .
2.2.3 Cavitation erosion modelling
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in
terms of work done by the surrounding liquid on a bubble. In general , for a cavity with volume
V
V
, the potential energy of that cavity equals 
E
pot
=
(
p
d
− p
v )
V
V
, (2.10)
where p
d
 is the pressure in the surrounding liquid and p
v
 is the vapour pressure inside
the bubble. 
The cavitation potential power P
p ot
 is introduced as the total derivative of the potential
energy E
pot
 
P
pot
=
D
(
p
d
− p
v )
Dt
V
V
+
D V
V
Dt
( p
d
− p
v
) .
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of
magnitude lower than the second term [15], and taking into account that p
v
= c on s t .
(isothermal flow), the cavitation potential power is written as
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of magni-
tude lower than the second term [15], and taking into account that p
v
 = const. (isothermal flow), 
the cavitation potential power is written as
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
P
pot
=
D V
V
Dt
(
p
d
− p
v )
. (2.12)
The driving pressure of the surrounding liquid p
d
 is taken as a time-averaged value in each
computational grid cell as proposed in [11]
p
d
=
1
t
❑ ∫
0
t
❑
p ( t ) dt , (2.13)
where p ( t ) is an instantaneous pressure that is then time-averaged over the current
simulated time t
❑
.
For the numerical computation of the cavitation potential power, it is convenient to calculate
the instantaneous change of the volume specific potential energy [11] or cavitation erosion
potential 
´ e
po t
=
P
pot
V
cel l
, (2.14)
where V
cell
 is the computational mesh cell volume. The total derivative of the vapour
volume in equation (2.12) can be calculated using three different but mathematically equivalent
formulations:
1
V
cel l
D V
V
Dt
=
{
∂ α
∂ t
+ u ∙ ∇ α
ρ
ρ
l
− ρ
v
∇ ∙ u
ρ
ρ
l
S
α
v
(2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it
was found that most accurate results were obtained using the third formulation from equation
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
D V
V
Dt
 with the
cavitation source term S
α v
. Finally, the cavitation erosion potential is given by
´ e
pot
=
(
p
d
− p
v )
ρ
ρ
l
S
α
v
. (2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This
accumulation and delayed release of energy is described by the additional conservation
equation for collapse induced kinetic energy ε
(2.12)
The driving pressure of the surrounding liquid p
d
 is taken as a time-averaged value in each com-
putational grid cell as proposed in [11] predicting the instantaneous surface impact power of 
collapsing cavities from the potential energy hypothesis (see Hammitt, 1963; Vogel and Lauter-
born, 1988
 Prediction of cavitation and particle erosion in a radial divergent test section 7 
   
---------- 
S
�
�
=
⎩
⎪
⎨
⎪
⎧𝜌𝜌 �
𝜌𝜌 �
𝜌𝜌 𝛼𝛼 �
(1 − 𝛼𝛼 �
)
3
𝑅𝑅 �
�
2
3
(𝑝𝑝 �
− 𝑝𝑝 )
𝜌𝜌 �
     if 𝑝𝑝 < 𝑝𝑝 �
𝜌𝜌 �
𝜌𝜌 �
𝜌𝜌 𝛼𝛼 �
(1 − 𝛼𝛼 �
)
3
𝑅𝑅 �
�
2
3
(𝑝𝑝 − 𝑝𝑝 �
)
𝜌𝜌 �
     if 𝑝𝑝 > 𝑝𝑝 �
, (2.8)
where the bubble radius is 
𝑅𝑅 �
= �
𝛼𝛼 �
𝑛𝑛 �
4
3
𝜋𝜋 (1 − 𝛼𝛼 �
)
�
� � ⁄
. 
(2.9)
For the bubble number density 𝑛𝑛 �
 a default value of 1 ∙ 10
��
 was used. The vapour pressure was 
set to  𝑝𝑝 �
= 1854 Pa. 
 
2.2.3 Cavitation erosion modelling 
Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in 
terms of work done by the surrounding liquid on a bubble. In general, for a cavity with volume 
𝑉𝑉 �
, the potential energy of that cavity equals  
𝐸𝐸 ���
= (𝑝𝑝 �
− 𝑝𝑝 �
)𝑉𝑉 �
, (2.10)
where 𝑝𝑝 �
 is the pressure in the surrounding liquid and 𝑝𝑝 �
 is the vapour pressure inside the bubble.  
The cavitation potential power 𝑃𝑃 ���
 is introduced as the total derivative of the potential energy 
𝐸𝐸 ���
  
𝑃𝑃 ���
=
𝐷𝐷 (𝑝𝑝 �
− 𝑝𝑝 �
)
𝐷𝐷𝐷𝐷 𝑉𝑉 �
+
𝐷𝐷 𝑉𝑉 �
𝐷𝐷𝐷𝐷
(𝑝𝑝 �
− 𝑝𝑝 �
). 
(2.11)
After neglecting the first term in equation (2.11), as it was found to be at least an order of 
magnitude lower than the second term [15], and taking into account that 𝑝𝑝 �
= const. (isothermal 
flow), the cavitation potential power is written as 
𝑃𝑃 ���
=
𝐷𝐷 𝑉𝑉 �
𝐷𝐷𝐷𝐷 (𝑝𝑝 �
− 𝑝𝑝 �
). (2.12)
The driving pressure of the surrounding liquid 𝑝𝑝 �
 is taken as a time-averaged value in each 
computational grid cell as proposed in [11] 
𝑝𝑝 �
=
1
𝐷𝐷 ∗
� 𝑝𝑝 (𝐷𝐷 )𝑑𝑑𝐷𝐷 �
∗
�
, (2.13)
where 𝑝𝑝 (𝐷𝐷 ) is an instantaneous pressure that is then time-averaged over the current simulated 
time 𝐷𝐷 ∗
. 
For the numerical computation of the cavitation potential power, it is convenient to calculate the 
instantaneous change of the volume specific potential energy [11] or cavitation erosion potential  
(2.13)
where p(t) is an instantaneous pressure that is then time-averaged over the current simulated 
time t
*
.

50 JET 50 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
For the numerical computation of the cavitation potential power, it is convenient to calculate the 
instantaneous change of the volume specific potential energy [11] or cavitation erosion potential
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
P
pot
=
D V
V
Dt
(
p
d
− p
v )
. (2.12)
The driving pressure of the surrounding liquid p
d
 is taken as a time-averaged value in each
computational grid cell as proposed in [11]
p
d
=
1
t
❑ ∫
0
t
❑
p ( t ) dt , (2.13)
where p ( t ) is an instantaneous pressure that is then time-averaged over the current
simulated time t
❑
.
For the numerical computation of the cavitation potential power, it is convenient to calculate
the instantaneous change of the volume specific potential energy [11] or cavitation erosion
potential 
´ e
po t
=
P
pot
V
cel l
, (2.14)
where V
cell
 is the computational mesh cell volume. The total derivative of the vapour
volume in equation (2.12) can be calculated using three different but mathematically equivalent
formulations:
1
V
cel l
D V
V
Dt
=
{
∂ α
∂ t
+ u ∙ ∇ α
ρ
ρ
l
− ρ
v
∇ ∙ u
ρ
ρ
l
S
α
v
(2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it
was found that most accurate results were obtained using the third formulation from equation
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
D V
V
Dt
 with the
cavitation source term S
α v
. Finally, the cavitation erosion potential is given by
´ e
pot
=
(
p
d
− p
v )
ρ
ρ
l
S
α
v
. (2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This
accumulation and delayed release of energy is described by the additional conservation
equation for collapse induced kinetic energy ε
(2.14)
where V
cell
 is the computational mesh cell volume. The total derivative of the vapour volume in 
equation (2.12) can be calculated using three different but mathematically equivalent formula-
tions:
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
P
pot
=
D V
V
Dt
(
p
d
− p
v )
. (2.12)
The driving pressure of the surrounding liquid p
d
 is taken as a time-averaged value in each
computational grid cell as proposed in [11]
p
d
=
1
t
❑ ∫
0
t
❑
p ( t ) dt , (2.13)
where p ( t ) is an instantaneous pressure that is then time-averaged over the current
simulated time t
❑
.
For the numerical computation of the cavitation potential power, it is convenient to calculate
the instantaneous change of the volume specific potential energy [11] or cavitation erosion
potential 
´ e
po t
=
P
pot
V
cel l
, (2.14)
where V
cell
 is the computational mesh cell volume. The total derivative of the vapour
volume in equation (2.12) can be calculated using three different but mathematically equivalent
formulations:
1
V
cel l
D V
V
Dt
=
{
∂ α
∂ t
+ u ∙ ∇ α
ρ
ρ
l
− ρ
v
∇ ∙ u
ρ
ρ
l
S
α
v
(2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it
was found that most accurate results were obtained using the third formulation from equation
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
D V
V
Dt
 with the
cavitation source term S
α v
. Finally, the cavitation erosion potential is given by
´ e
pot
=
(
p
d
− p
v )
ρ
ρ
l
S
α
v
. (2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This
accumulation and delayed release of energy is described by the additional conservation
equation for collapse induced kinetic energy ε
(2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it 
was found that most accurate results were obtained using the third formulation from equation 
(2.15) [15], [48], which expresses the total derivative of the vapour volume 
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022) 
  Issue 2 
---------- 
𝑒𝑒 ̇ ���
=
𝑃𝑃 ���
𝑉𝑉 ����
, (2.14)
where 𝑉𝑉 ����
 is the computational mesh cell volume. The total derivative of the vapour volume in 
equation (2.12) can be calculated using three different but mathematically equivalent 
formulations: 
1
𝑉𝑉 ����
𝐷𝐷𝑉𝑉
�
𝐷𝐷𝐷𝐷
=
⎩
⎪
⎨
⎪
⎧
𝜕𝜕𝜕𝜕
𝜕𝜕𝐷𝐷
+ 𝒖𝒖 ∙ ∇𝜕𝜕 𝜌𝜌 𝜌𝜌 �
− 𝜌𝜌 �
∇∙𝒖𝒖 𝜌𝜌 𝜌𝜌 �
S
�
�
 (2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it 
was found that most accurate results were obtained using the third formulation from equation 
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
��
�
��
 with the cavitation 
source term S
�
�
. Finally, the cavitation erosion potential is given by 
𝑒𝑒 ̇ ���
= (𝑝𝑝 �
− 𝑝𝑝 �
)
𝜌𝜌 𝜌𝜌 �
S
�
�
. 
(2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model 
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This 
accumulation and delayed release of energy is described by the additional conservation equation 
for collapse induced kinetic energy 𝜀𝜀 
𝜕𝜕𝜀𝜀
𝜕𝜕𝐷𝐷
+ ∇ ∙ (𝒖𝒖 𝒊𝒊 𝜀𝜀 ) = −𝑒𝑒 ̇ ���
(𝐷𝐷 ), 
(2.17)
where 𝒖𝒖 𝒊𝒊 is the collapse induced velocity, and 𝑒𝑒 ̇ ���
(𝐷𝐷 ) is the instantaneous radiated power of the 
pressure wave released in the cavity centre at the final stage of collapse. 
A model for transport equation (2.17) is then written in explicitly discretized form: 
𝜀𝜀 |
�� ∆�
= 𝜀𝜀 |
�
+
𝜕𝜕𝜀𝜀
𝜕𝜕𝐷𝐷
�
�
∆𝐷𝐷 = (1 − 𝛽𝛽 |
�
)� (𝐾𝐾 − 1)𝑒𝑒 ̇ ���
∆𝐷𝐷 − 𝜀𝜀 (𝑃𝑃 �
− 1)� �
�
. (2.18)
𝑃𝑃 �
 in equation (2.18) is the projection operator, which ensures that the specific kinetic energy of 
the collapse accumulates on the inside of the vapour-liquid interface, 𝐾𝐾 is the global energy 
conservation parameter and  𝛽𝛽 is the energy release criterion. These are further defined as: 
𝑃𝑃 �
= 𝑚𝑚𝑚𝑚𝑚𝑚 �
𝒖𝒖 ∙ ∇𝜀𝜀 |𝒖𝒖 ||∇𝜀𝜀 |
, 0� , (2.19)
𝐾𝐾 = 
∫
𝜀𝜀 ∆𝐷𝐷 �
𝑃𝑃 �
𝑑𝑑𝑉𝑉 ∫ 𝑒𝑒 ̇ ���
�
𝑑𝑑𝑉𝑉 , (2.20)
 with the cavita-
tion source term S
α
v
. Finally, the cavitation erosion potential is given by
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
P
pot
=
D V
V
Dt
(
p
d
− p
v )
. (2.12)
The driving pressure of the surrounding liquid p
d
 is taken as a time-averaged value in each
computational grid cell as proposed in [11]
p
d
=
1
t
❑ ∫
0
t
❑
p ( t ) dt , (2.13)
where p ( t ) is an instantaneous pressure that is then time-averaged over the current
simulated time t
❑
.
For the numerical computation of the cavitation potential power, it is convenient to calculate
the instantaneous change of the volume specific potential energy [11] or cavitation erosion
potential 
´ e
po t
=
P
pot
V
cel l
, (2.14)
where V
cell
 is the computational mesh cell volume. The total derivative of the vapour
volume in equation (2.12) can be calculated using three different but mathematically equivalent
formulations:
1
V
cel l
D V
V
Dt
=
{
∂ α
∂ t
+ u ∙ ∇ α
ρ
ρ
l
− ρ
v
∇ ∙ u
ρ
ρ
l
S
α
v
(2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it
was found that most accurate results were obtained using the third formulation from equation
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
D V
V
Dt
 with the
cavitation source term S
α v
. Finally, the cavitation erosion potential is given by
´ e
pot
=
(
p
d
− p
v )
ρ
ρ
l
S
α
v
. (2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This
accumulation and delayed release of energy is described by the additional conservation
equation for collapse induced kinetic energy ε
(2.16)
Based on the previous studies [12], [13] which suppresses the buoyant pressure gradient that 
otherwise deteriorates the sphericity of the bubbles. We measure the radius of the rebound 
bubble and estimate the shock energy as a function of the initial bubble radius (2-5.6mm Schen-
ke, Melissaris and van Terwisga proposed a model [14] to describe the focusing of initial cavity 
potential energy into the cavity collapse centre. This accumulation and delayed release of energy 
is described by the additional conservation equation for collapse induced kinetic energy ε
Prediction of cavitation and particle erosion in a radial divergent test section 9
----------
∂ ε
∂ t
+ ∇ ∙
(
u
i
ε
)
= − ´ e
rad
( t ) , (2.17)
where u
i
 is the collapse induced velocity, and ´ e
rad
( t ) is the instantaneous radiated
power of the pressure wave released in the cavity centre at the final stage of collapse.
A model for transport equation (2.17) is then written in explicitly discretized form:
ε ∨ ❑
t
+
∂ ε
∂ t
|
t
∆ t = ( 1 − β |
t
) [
( K − 1 ) ´ e
pot
∆ t − ε ( P
u
− 1 ) ] |
t
.
ε ∨ ❑
t + ∆ t
= ¿
(2.18)
P
u
 in equation (2.18) is the projection operator, which ensures that the specific kinetic
energy of the collapse accumulates on the inside of the vapour-liquid interface, K is the
global energy conservation parameter and  β is the energy release criterion. These are
further defined as:
P
u
= m ax
[
u ∙ ∇ ε
| u | | ∇ ε |
,0
]
, (2.19)
K =
∫
V
❑
ε
∆ t
P
u
dV
∫
V
❑
´ e
pot
dV
, (2.20)
β =
{
1, if p > p
d
∧ α
v
= 0
0, els e .
(2.21)
T h e e n e r g y r e l e a s e d a s a c o n s e q u e n c e o f t h e f i n i s h e d c o l l a p s e o f a c a v i t y a t t i m e l e v e l
t + ∆ t , where ∆ t is the time step, is expressed as the specific power of the pressure
wave 
´ e
ra d |
t + ∆ t
=
( βε ) ❑
t
∆ t
(2.22)
An infinite wave propagation speed is assumed in this model, therefore the energy received by a
surface element S is expressed as the instantaneous surface specific impact power
´ e
S
=
1
4 π
∫
V
❑
´ e
rad
[
(
x
cel l
− x
S )
∙ n
| x
cel l
− x
S |
3
]
dV , (2.23)
where x
cel l
 is the position vector of the computational cell centre, x
S
 is the position
v e c t o r o f t h e w a l l s u r f a c e e l e m e n t c e n t r e , a n d n is the surface element normal.
(2.17)
where u
i
 is the collapse induced velocity, é
rad
(t) and is the instantaneous radiated power of the 
pressure wave released in the cavity centre at the final stage of collapse.
A model for transport equation (2.17) is then written in explicitly discretized form:
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022) 
  Issue 2 
---------- 
𝑒𝑒 ̇ ���
=
𝑃𝑃 ���
𝑉𝑉 ����
, (2.14)
where 𝑉𝑉 ����
 is the computational mesh cell volume. The total derivative of the vapour volume in 
equation (2.12) can be calculated using three different but mathematically equivalent 
formulations: 
1
𝑉𝑉 ����
𝐷𝐷𝑉𝑉
�
𝐷𝐷𝐷𝐷
=
⎩
⎪
⎨
⎪
⎧
𝜕𝜕𝜕𝜕
𝜕𝜕𝐷𝐷
+ 𝒖𝒖 ∙ ∇𝜕𝜕 𝜌𝜌 𝜌𝜌 �
− 𝜌𝜌 �
∇∙𝒖𝒖 𝜌𝜌 𝜌𝜌 �
S
�
�
 (2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it 
was found that most accurate results were obtained using the third formulation from equation 
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
��
�
��
 with the cavitation 
source term S
�
�
. Finally, the cavitation erosion potential is given by 
𝑒𝑒 ̇ ���
= (𝑝𝑝 �
− 𝑝𝑝 �
)
𝜌𝜌 𝜌𝜌 �
S
�
�
. (2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model 
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This 
accumulation and delayed release of energy is described by the additional conservation equation 
for collapse induced kinetic energy 𝜀𝜀 
𝜕𝜕𝜀𝜀
𝜕𝜕𝐷𝐷
+ ∇ ∙ (𝒖𝒖 𝒊𝒊 𝜀𝜀 ) = −𝑒𝑒 ̇ ���
(𝐷𝐷 ), 
(2.17)
where 𝒖𝒖 𝒊𝒊 is the collapse induced velocity, and 𝑒𝑒 ̇ ���
(𝐷𝐷 ) is the instantaneous radiated power of the 
pressure wave released in the cavity centre at the final stage of collapse. 
A model for transport equation (2.17) is then written in explicitly discretized form: 
𝜀𝜀 |
�� ∆�
= 𝜀𝜀 |
�
+
𝜕𝜕𝜀𝜀
𝜕𝜕𝐷𝐷
�
�
∆𝐷𝐷 = (1 − 𝛽𝛽 |
�
)� (𝐾𝐾 − 1)𝑒𝑒 ̇ ���
∆𝐷𝐷 − 𝜀𝜀 (𝑃𝑃 �
− 1)� �
�
. (2.18)
𝑃𝑃 �
 in equation (2.18) is the projection operator, which ensures that the specific kinetic energy of 
the collapse accumulates on the inside of the vapour-liquid interface, 𝐾𝐾 is the global energy 
conservation parameter and  𝛽𝛽 is the energy release criterion. These are further defined as: 
𝑃𝑃 �
= 𝑚𝑚𝑚𝑚𝑚𝑚 �
𝒖𝒖 ∙ ∇𝜀𝜀 |𝒖𝒖 ||∇𝜀𝜀 |
, 0� , (2.19)
𝐾𝐾 = 
∫
𝜀𝜀 ∆𝐷𝐷 �
𝑃𝑃 �
𝑑𝑑𝑉𝑉 ∫ 𝑒𝑒 ̇ ���
�
𝑑𝑑𝑉𝑉 , (2.20)
(2.18)
P
u
 in equation (2.18) is the projection operator, which ensures that the specific kinetic energy of 
the collapse accumulates on the inside of the vapour-liquid interface, K is the global energy con-
servation parameter and β is the energy release criterion. These are further defined as:
Prediction of cavitation and particle erosion in a radial divergent test section 9
----------
∂ ε
∂ t
+ ∇ ∙
(
u
i
ε
)
= − ´ e
rad
( t ) , (2.17)
where u
i
 is the collapse induced velocity, and ´ e
rad
( t ) is the instantaneous radiated
power of the pressure wave released in the cavity centre at the final stage of collapse.
A model for transport equation (2.17) is then written in explicitly discretized form:
ε ∨ ❑
t
+
∂ ε
∂ t
|
t
∆ t = ( 1 − β |
t
) [
( K − 1 ) ´ e
pot
∆ t − ε ( P
u
− 1 ) ] |
t
.
ε ∨ ❑
t + ∆ t
= ¿
(2.18)
P
u
 in equation (2.18) is the projection operator, which ensures that the specific kinetic
energy of the collapse accumulates on the inside of the vapour-liquid interface, K is the
global energy conservation parameter and  β is the energy release criterion. These are
further defined as:
P
u
= m ax
[
u ∙ ∇ ε
| u | | ∇ ε |
,0
]
, (2.19)
K =
∫
V
❑
ε
∆ t
P
u
dV
∫
V
❑
´ e
pot
dV
, (2.20)
β =
{
1, if p > p
d
∧ α
v
= 0
0, els e .
(2.21)
T h e e n e r g y r e l e a s e d a s a c o n s e q u e n c e o f t h e f i n i s h e d c o l l a p s e o f a c a v i t y a t t i m e l e v e l
t + ∆ t , where ∆ t is the time step, is expressed as the specific power of the pressure
wave 
´ e
ra d |
t + ∆ t
=
( βε ) ❑
t
∆ t
(2.22)
An infinite wave propagation speed is assumed in this model, therefore the energy received by a
surface element S is expressed as the instantaneous surface specific impact power
´ e
S
=
1
4 π
∫
V
❑
´ e
rad
[
(
x
cel l
− x
S )
∙ n
| x
cel l
− x
S |
3
]
dV , (2.23)
where x
cel l
 is the position vector of the computational cell centre, x
S
 is the position
v e c t o r o f t h e w a l l s u r f a c e e l e m e n t c e n t r e , a n d n is the surface element normal.
(2.19)

JET  51
Prediction of cavitation and particle erosion in a radial divergent test section
8 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022) 
  Issue 2 
---------- 
𝑒𝑒 ̇ ���
=
𝑃𝑃 ���
𝑉𝑉 ����
, (2.14)
where 𝑉𝑉 ����
 is the computational mesh cell volume. The total derivative of the vapour volume in 
equation (2.12) can be calculated using three different but mathematically equivalent 
formulations: 
1
𝑉𝑉 ����
𝐷𝐷𝑉𝑉
�
𝐷𝐷𝐷𝐷
=
⎩
⎪
⎨
⎪
⎧
𝜕𝜕𝜕𝜕
𝜕𝜕𝐷𝐷
+ 𝒖𝒖 ∙ ∇𝜕𝜕 𝜌𝜌 𝜌𝜌 �
− 𝜌𝜌 �
∇∙𝒖𝒖 𝜌𝜌 𝜌𝜌 �
S
�
�
 (2.15)
However, the numerical error introduced varies between the three formulations [15], [48] and it 
was found that most accurate results were obtained using the third formulation from equation 
(2.15)[15], [48], which expresses the total derivative of the vapour volume 
��
�
��
 with the cavitation 
source term S
�
�
. Finally, the cavitation erosion potential is given by 
𝑒𝑒 ̇ ���
= (𝑝𝑝 �
− 𝑝𝑝 �
)
𝜌𝜌 𝜌𝜌 �
S
�
�
. (2.16)
Based on the previous studies [12], [13] Schenke, Melissaris and van Terwisga proposed a model 
[14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This 
accumulation and delayed release of energy is described by the additional conservation equation 
for collapse induced kinetic energy 𝜀𝜀 
𝜕𝜕𝜀𝜀
𝜕𝜕𝐷𝐷
+ ∇ ∙ (𝒖𝒖 𝒊𝒊 𝜀𝜀 ) = −𝑒𝑒 ̇ ���
(𝐷𝐷 ), 
(2.17)
where 𝒖𝒖 𝒊𝒊 is the collapse induced velocity, and 𝑒𝑒 ̇ ���
(𝐷𝐷 ) is the instantaneous radiated power of the 
pressure wave released in the cavity centre at the final stage of collapse. 
A model for transport equation (2.17) is then written in explicitly discretized form: 
𝜀𝜀 |
�� ∆�
= 𝜀𝜀 |
�
+
𝜕𝜕𝜀𝜀
𝜕𝜕𝐷𝐷
�
�
∆𝐷𝐷 = (1 − 𝛽𝛽 |
�
)� (𝐾𝐾 − 1)𝑒𝑒 ̇ ���
∆𝐷𝐷 − 𝜀𝜀 (𝑃𝑃 �
− 1)� �
�
. (2.18)
𝑃𝑃 �
 in equation (2.18) is the projection operator, which ensures that the specific kinetic energy of 
the collapse accumulates on the inside of the vapour-liquid interface, 𝐾𝐾 is the global energy 
conservation parameter and  𝛽𝛽 is the energy release criterion. These are further defined as: 
𝑃𝑃 �
= 𝑚𝑚𝑚𝑚𝑚𝑚 �
𝒖𝒖 ∙ ∇𝜀𝜀 |𝒖𝒖 ||∇𝜀𝜀 |
, 0� , (2.19)
𝐾𝐾 = 
∫
𝜀𝜀 ∆𝐷𝐷 �
𝑃𝑃 �
𝑑𝑑𝑉𝑉 ∫ 𝑒𝑒 ̇ ���
�
𝑑𝑑𝑉𝑉 , (2.20)
(2.20)
 Prediction of cavitation and particle erosion in a radial divergent test section 9 
   
---------- 
𝛽𝛽 = �
1,   if 𝑝𝑝 > 𝑝𝑝 �
 and 𝛼𝛼 �
= 0
0, else.
 (2.21)
The energy released as a consequence of the finished collapse of a cavity at time level 𝑡𝑡 +∆ 𝑡𝑡 , 
where ∆𝑡𝑡 is the time step, is expressed as the specific power of the pressure wave  
𝑒𝑒 ̇ ���
|
�� ∆�
=
(𝛽𝛽𝛽𝛽 )|
�
∆𝑡𝑡 
(2.22)
An infinite wave propagation speed is assumed in this model, therefore the energy received by a 
surface element 𝑆𝑆 is expressed as the instantaneous surface specific impact power 
𝑒𝑒 ̇ �
=
1
4𝜋𝜋 � 𝑒𝑒 ̇ ���
�
(𝒙𝒙 𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 − 𝒙𝒙 𝑺𝑺 ) ∙ 𝒏𝒏 |𝒙𝒙 𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 − 𝒙𝒙 𝑺𝑺 |
�
�𝑑𝑑𝑑𝑑 ,
�
 (2.23)
where 𝒙𝒙 𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 is the position vector of the computational cell centre, 𝒙𝒙 𝑺𝑺 is the position vector of the 
wall surface element centre, and 𝒏𝒏 is the surface element normal. Instantaneous surface specific 
impact energy with units of (
�
�
�
) is obtained by integrating equation (2.23) over the simulated 
time (𝑡𝑡 ): 
𝑒𝑒 �
= � 𝑒𝑒 ̇ �
�
�
𝑑𝑑𝑡𝑡 . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a 
UDF[48].  
 
2.2.4 Particle motion in a fluid 
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by Newton's 
second law. In a Lagrangian frame particle velocity is defined as 
𝒗𝒗 𝑷𝑷 =
𝑑𝑑 𝒙𝒙 𝑷𝑷 𝑑𝑑𝑡𝑡 , (2.25)
where 𝒙𝒙 𝑷𝑷 is the particle position in the global inertial frame of reference. We assume the spherical 
particle, therefore particles are assumed to undergo only translation without rotation. Newton's 
second law is then written as 
𝜌𝜌 �
𝑑𝑑 �
�
6
𝑑𝑑 𝒗𝒗 𝑷𝑷 𝑑𝑑𝑡𝑡 = 𝑭𝑭 , 
(2.26)
where the particle density is 𝜌𝜌 �
= 1700 kg/m
3
, the particle diameter is 𝑑𝑑 �
=5 μm and 𝑭𝑭 is the 
resultant force acting on the particle.  
The resultant force can be expressed as a sum of different contributing forces 
𝑭𝑭 = 𝑭𝑭 𝑫𝑫 + 𝑭𝑭 𝑩𝑩 + 𝑭𝑭 𝑮𝑮 + 𝑭𝑭 𝑷𝑷𝑮𝑮 + 𝑭𝑭 𝑽𝑽𝑽𝑽
, (2.27)
where 𝑭𝑭 𝑫𝑫 is the drag force, 𝑭𝑭 𝑩𝑩 is the buoyancy force, 𝑭𝑭 𝑮𝑮 is the force of gravity, 𝑭𝑭 𝑷𝑷𝑮𝑮 is the 
pressure gradient force and 𝑭𝑭 𝑽𝑽𝑽𝑽
 is the virtual mass force. 
(2.21)
The energy released as a consequence of the finished collapse of a cavity at time level t + ∆t, 
where ∆t is the time step, is expressed as the specific power of the pressure wave
 Prediction of cavitation and particle erosion in a radial divergent test section 9 
   
---------- 
𝛽𝛽 = �
1,   if 𝑝𝑝 > 𝑝𝑝 �
 and 𝛼𝛼 �
= 0
0, else.
 (2.21)
The energy released as a consequence of the finished collapse of a cavity at time level 𝑡𝑡 +∆ 𝑡𝑡 , 
where ∆𝑡𝑡 is the time step, is expressed as the specific power of the pressure wave  
𝑒𝑒 ̇ ���
|
�� ∆�
=
(𝛽𝛽𝛽𝛽 )|
�
∆𝑡𝑡 
(2.22)
An infinite wave propagation speed is assumed in this model, therefore the energy received by a 
surface element 𝑆𝑆 is expressed as the instantaneous surface specific impact power 
𝑒𝑒 ̇ �
=
1
4𝜋𝜋 � 𝑒𝑒 ̇ ���
�
(𝒙𝒙 𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 − 𝒙𝒙 𝑺𝑺 ) ∙ 𝒏𝒏 |𝒙𝒙 𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 − 𝒙𝒙 𝑺𝑺 |
�
�𝑑𝑑𝑑𝑑 ,
�
 (2.23)
where 𝒙𝒙 𝒄𝒄𝒄𝒄𝒄𝒄𝒄𝒄 is the position vector of the computational cell centre, 𝒙𝒙 𝑺𝑺 is the position vector of the 
wall surface element centre, and 𝒏𝒏 is the surface element normal. Instantaneous surface specific 
impact energy with units of (
�
�
�
) is obtained by integrating equation (2.23) over the simulated 
time (𝑡𝑡 ): 
𝑒𝑒 �
= � 𝑒𝑒 ̇ �
�
�
𝑑𝑑𝑡𝑡 . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a 
UDF[48].  
 
2.2.4 Particle motion in a fluid 
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by Newton's 
second law. In a Lagrangian frame particle velocity is defined as 
𝒗𝒗 𝑷𝑷 =
𝑑𝑑 𝒙𝒙 𝑷𝑷 𝑑𝑑𝑡𝑡 , 
(2.25)
where 𝒙𝒙 𝑷𝑷 is the particle position in the global inertial frame of reference. We assume the spherical 
particle, therefore particles are assumed to undergo only translation without rotation. Newton's 
second law is then written as 
𝜌𝜌 �
𝑑𝑑 �
�
6
𝑑𝑑 𝒗𝒗 𝑷𝑷 𝑑𝑑𝑡𝑡 = 𝑭𝑭 , 
(2.26)
where the particle density is 𝜌𝜌 �
= 1700 kg/m
3
, the particle diameter is 𝑑𝑑 �
=5 μm and 𝑭𝑭 is the 
resultant force acting on the particle.  
The resultant force can be expressed as a sum of different contributing forces 
𝑭𝑭 = 𝑭𝑭 𝑫𝑫 + 𝑭𝑭 𝑩𝑩 + 𝑭𝑭 𝑮𝑮 + 𝑭𝑭 𝑷𝑷𝑮𝑮 + 𝑭𝑭 𝑽𝑽𝑽𝑽
, (2.27)
where 𝑭𝑭 𝑫𝑫 is the drag force, 𝑭𝑭 𝑩𝑩 is the buoyancy force, 𝑭𝑭 𝑮𝑮 is the force of gravity, 𝑭𝑭 𝑷𝑷𝑮𝑮 is the 
pressure gradient force and 𝑭𝑭 𝑽𝑽𝑽𝑽
 is the virtual mass force. 
(2.22)
An infinite wave propagation speed is assumed in this model, therefore the energy received by a 
surface S element is expressed as the instantaneous surface specific impact power
Prediction of cavitation and particle erosion in a radial divergent test section 9
----------
∂ ε
∂ t
+ ∇ ∙
(
u
i
ε
)
= − ´ e
rad
( t ) , (2.17)
where u
i
 is the collapse induced velocity, and ´ e
rad
( t ) is the instantaneous radiated
power of the pressure wave released in the cavity centre at the final stage of collapse.
A model for transport equation (2.17) is then written in explicitly discretized form:
ε ∨ ❑
t
+
∂ ε
∂ t
|
t
∆ t = ( 1 − β |
t
) [
( K − 1 ) ´ e
pot
∆ t − ε ( P
u
− 1 ) ] |
t
.
ε ∨ ❑
t + ∆ t
= ¿
(2.18)
P
u
 in equation (2.18) is the projection operator, which ensures that the specific kinetic
energy of the collapse accumulates on the inside of the vapour-liquid interface, K is the
global energy conservation parameter and  β is the energy release criterion. These are
further defined as:
P
u
= m ax
[
u ∙ ∇ ε
| u | | ∇ ε |
,0
]
, (2.19)
K =
∫
V
❑
ε
∆ t
P
u
dV
∫
V
❑
´ e
pot
dV
, (2.20)
β =
{
1, if p > p
d
∧ α
v
= 0
0, els e .
(2.21)
T h e e n e r g y r e l e a s e d a s a c o n s e q u e n c e o f t h e f i n i s h e d c o l l a p s e o f a c a v i t y a t t i m e l e v e l
t + ∆ t , where ∆ t is the time step, is expressed as the specific power of the pressure
wave 
´ e
ra d |
t + ∆ t
=
( βε ) ❑
t
∆ t
(2.22)
An infinite wave propagation speed is assumed in this model, therefore the energy received by a
surface element S is expressed as the instantaneous surface specific impact power
´ e
S
=
1
4 π
∫
V
´ e
rad
[
(
x
cel l
− x
S )
∙ n
| x
cel l
− x
S |
3
]
dV , (2.23)
where x
cel l
 is the position vector of the computational cell centre, x
S
 is the position
v e c t o r o f t h e w a l l s u r f a c e e l e m e n t c e n t r e , a n d n is the surface element normal.
(2.23)
where x
cell
 is the position vector of the computational cell centre, x
S
 is the position vector of the 
wall surface element centre, and n is the surface element normal. Instantaneous surface specific 
impact energy with units of (
 J 
m²
) is obtained by integrating equation (2.23) over the simulated 
time (t):
10 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
Instantaneous surface specific impact energy with units of (
J
m
2
) is obtained by integrating
equation (2.23) over the simulated time ( t ):
e
S
=
∫
0
t
´ e
S
dt . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a
UDF[48]. 
2.2.4 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by
Newton's second law. In a Lagrangian frame particle velocity is defined as
v
P
=
d x
P
dt
, (2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the
spherical particle, therefore particles are assumed to undergo only translation without rotation.
Newton's second law is then written as
ρ
P
d
P
3
6
d v
P
dt
= F ,
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
= 5 μm and
F is the resultant force acting on the particle. 
The resultant force can be expressed as a sum of different contributing forces
F = F
D
+ F
B
+ F
G
+ F
PG
+ F
V M
, (2.27)
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity,
F
PG
 is the pressure gradient force and F
V M
 is the virtual mass force.
The drag force acting on the particle is written as
F
D
=
18 μd
P
6
C
D
ℜ
P
24
(
u − v
P )
, (2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
C
D
= c
1
+
c
2
ℜ
P
+
c
3
ℜ
p
2
, (2.29)
(2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via 
a UDF[48].
2.2.5 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by Newton’s 
second law. In a Lagrangian frame particle velocity is defined as
10 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
Instantaneous surface specific impact energy with units of (
J
m
2
) is obtained by integrating
equation (2.23) over the simulated time ( t ):
e
S
=
∫
0
t
´ e
S
dt . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a
UDF[48]. 
2.2.4 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by
Newton's second law. In a Lagrangian frame particle velocity is defined as
v
P
=
d x
P
dt
, (2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the
spherical particle, therefore particles are assumed to undergo only translation without rotation.
Newton's second law is then written as
ρ
P
d
P
3
6
d v
P
dt
= F ,
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
= 5 μm and
F is the resultant force acting on the particle. 
The resultant force can be expressed as a sum of different contributing forces
F = F
D
+ F
B
+ F
G
+ F
PG
+ F
V M
, (2.27)
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity,
F
PG
 is the pressure gradient force and F
V M
 is the virtual mass force.
The drag force acting on the particle is written as
F
D
=
18 μd
P
6
C
D
ℜ
P
24
(
u − v
P )
, (2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
C
D
= c
1
+
c
2
ℜ
P
+
c
3
ℜ
p
2
, (2.29)
(2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the spherical 
particle, therefore particles are assumed to undergo only translation without rotation. Newton’s 
second law is then written as
10 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
Instantaneous surface specific impact energy with units of (
J
m
2
) is obtained by integrating
equation (2.23) over the simulated time ( t ):
e
S
=
∫
0
t
´ e
S
dt . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a
UDF[48]. 
2.2.4 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by
Newton's second law. In a Lagrangian frame particle velocity is defined as
v
P
=
d x
P
dt
, (2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the
spherical particle, therefore particles are assumed to undergo only translation without rotation.
Newton's second law is then written as
ρ
P
d
P
3
6
d v
P
dt
= F ,
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
= 5 μm and
F is the resultant force acting on the particle. 
The resultant force can be expressed as a sum of different contributing forces
F = F
D
+ F
B
+ F
G
+ F
PG
+ F
V M
, (2.27)
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity,
F
PG
 is the pressure gradient force and F
V M
 is the virtual mass force.
The drag force acting on the particle is written as
F
D
=
18 μd
P
6
C
D
ℜ
P
24
(
u − v
P )
, (2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
C
D
= c
1
+
c
2
ℜ
P
+
c
3
ℜ
p
2
, (2.29)
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
 = 5 μm and is the result-
ant force F acting on the particle.
The resultant force can be expressed as a sum of different contributing forces
10 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
Instantaneous surface specific impact energy with units of (
J
m
2
) is obtained by integrating
equation (2.23) over the simulated time ( t ):
e
S
=
∫
0
t
´ e
S
dt . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a
UDF[48]. 
2.2.4 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by
Newton's second law. In a Lagrangian frame particle velocity is defined as
v
P
=
d x
P
dt
, (2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the
spherical particle, therefore particles are assumed to undergo only translation without rotation.
Newton's second law is then written as
ρ
P
d
P
3
6
d v
P
dt
= F ,
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
= 5 μm and
F is the resultant force acting on the particle. 
The resultant force can be expressed as a sum of different contributing forces
F = F
D
+ F
B
+ F
G
+ F
PG
+ F
V M
, (2.27)
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity,
F
PG
 is the pressure gradient force and F
V M
 is the virtual mass force.
The drag force acting on the particle is written as
F
D
=
18 μd
P
6
C
D
ℜ
P
24
(
u − v
P )
, (2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
C
D
= c
1
+
c
2
ℜ
P
+
c
3
ℜ
p
2
, (2.29)
(2.27)

52 JET 52 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity, F
PG
 is the pressure 
gradient force and F
VM
 is the virtual mass force.
The drag force acting on the particle is written as
10 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
Instantaneous surface specific impact energy with units of (
J
m
2
) is obtained by integrating
equation (2.23) over the simulated time ( t ):
e
S
=
∫
0
t
´ e
S
dt . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a
UDF[48]. 
2.2.4 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by
Newton's second law. In a Lagrangian frame particle velocity is defined as
v
P
=
d x
P
dt
, (2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the
spherical particle, therefore particles are assumed to undergo only translation without rotation.
Newton's second law is then written as
ρ
P
d
P
3
6
d v
P
dt
= F ,
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
= 5 μm and
F is the resultant force acting on the particle. 
The resultant force can be expressed as a sum of different contributing forces
F = F
D
+ F
B
+ F
G
+ F
PG
+ F
V M
, (2.27)
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity,
F
PG
 is the pressure gradient force and F
V M
 is the virtual mass force.
The drag force acting on the particle is written as
F
D
=
18 μd
P
6
C
D
ℜ
P
24
(
u − v
P )
, (2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
C
D
= c
1
+
c
2
ℜ
P
+
c
3
ℜ
p
2
, (2.29)
(2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
10 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
Instantaneous surface specific impact energy with units of (
J
m
2
) is obtained by integrating
equation (2.23) over the simulated time ( t ):
e
S
=
∫
0
t
´ e
S
dt . (2.24)
The described model to predict cavitation was previously implemented inside ANSYS Fluent via a
UDF[48]. 
2.2.4 Particle motion in a fluid
Particles are tracked in a Lagrangian frame, while the motion of particles is governed by
Newton's second law. In a Lagrangian frame particle velocity is defined as
v
P
=
d x
P
dt
, (2.25)
where x
P
 is the particle position in the global inertial frame of reference. We assume the
spherical particle, therefore particles are assumed to undergo only translation without rotation.
Newton's second law is then written as
ρ
P
d
P
3
6
d v
P
dt
= F ,
(2.26)
where the particle density is ρ
P
= 1700 kg/m
3
, the particle diameter is d
P
= 5 μm and
F is the resultant force acting on the particle. 
The resultant force can be expressed as a sum of different contributing forces
F = F
D
+ F
B
+ F
G
+ F
PG
+ F
V M
, (2.27)
where F
D
 is the drag force, F
B
 is the buoyancy force, F
G
 is the force of gravity,
F
PG
 is the pressure gradient force and F
V M
 is the virtual mass force.
The drag force acting on the particle is written as
F
D
=
18 μd
P
6
C
D
ℜ
P
24
(
u − v
P )
, (2.28)
where the drag coefficient is given by a correlation by Morsi and Alexander [53]
C
D
= c
1
+
c
2
ℜ
P
+
c
3
ℜ
p
2
, (2.29)
(2.29)
c
1
, c
2 
and c
3 
are constants of the c
D
 ̶ ℜ
P
 curve fit and particle Reynolds number is defined as
Prediction of cavitation and particle erosion in a radial divergent test section 11
----------
c
1
, c
2
 and c
3
 are constants of the C
D
  ̶  ℜ
P
 curve fit and particle Reynolds
number is defined as
ℜ
P
=
ρ d
P
| u − v
P
|
μ
.
(2.30)
The force resulting from buoyancy and gravity is
F
B
+ F
G
=
d
P
3
6
g
(
ρ
P
− ρ
)
,
(2.31)
the pressure gradient force is 
F
PG
=
d
P
3
6
ρ
D u
Dt
,
(2.32)
and the virtual mass force is 
F
V M
=
1
2
d
P
3
6
ρ
(
D u
Dt
−
d v
P
dt
)
. (2.32)
Two-way coupling between the continuous and discrete phase is achieved by including the 
momentum source term in equation 2.4. 
2.2.5 Particle-wall interaction
To resolve the interaction between the particle and the wall, a simplified hard sphere approach
is adopted, where the loss of particle momentum is accounted for via restitution coefficients in
the normal and tangential directions with respect to the observed wall. This scenario is
presented in Fig. 3.
Figure 3: The interaction between the particle and the wall.
In accordance with the hard sphere approach, the normal particle velocity after rebounding 
from the wall is
v
n ,2
= e
n
v
n ,1
, (2.34)
and the tangential particle velocity after rebounding from the wall is 
(2.30)
The force resulting from buoyancy and gravity is
Prediction of cavitation and particle erosion in a radial divergent test section 11
----------
c
1
, c
2
 and c
3
 are constants of the C
D
  ̶  ℜ
P
 curve fit and particle Reynolds
number is defined as
ℜ
P
=
ρ d
P
| u − v
P
|
μ
.
(2.30)
The force resulting from buoyancy and gravity is
F
B
+ F
G
=
d
P
3
6
g
(
ρ
P
− ρ
)
,
(2.31)
the pressure gradient force is 
F
PG
=
d
P
3
6
ρ
D u
Dt
,
(2.32)
and the virtual mass force is 
F
V M
=
1
2
d
P
3
6
ρ
(
D u
Dt
−
d v
P
dt
)
. (2.32)
Two-way coupling between the continuous and discrete phase is achieved by including the 
momentum source term in equation 2.4. 
2.2.5 Particle-wall interaction
To resolve the interaction between the particle and the wall, a simplified hard sphere approach
is adopted, where the loss of particle momentum is accounted for via restitution coefficients in
the normal and tangential directions with respect to the observed wall. This scenario is
presented in Fig. 3.
Figure 3: The interaction between the particle and the wall.
In accordance with the hard sphere approach, the normal particle velocity after rebounding 
from the wall is
v
n ,2
= e
n
v
n ,1
, (2.34)
and the tangential particle velocity after rebounding from the wall is 
(2.31)
the pressure gradient force is
Prediction of cavitation and particle erosion in a radial divergent test section 11
----------
c
1
, c
2
 and c
3
 are constants of the C
D
  ̶  ℜ
P
 curve fit and particle Reynolds
number is defined as
ℜ
P
=
ρ d
P
| u − v
P
|
μ
.
(2.30)
The force resulting from buoyancy and gravity is
F
B
+ F
G
=
d
P
3
6
g
(
ρ
P
− ρ
)
,
(2.31)
the pressure gradient force is 
F
PG
=
d
P
3
6
ρ
D u
Dt
,
(2.32)
and the virtual mass force is 
F
V M
=
1
2
d
P
3
6
ρ
(
D u
Dt
−
d v
P
dt
)
. (2.32)
Two-way coupling between the continuous and discrete phase is achieved by including the 
momentum source term in equation 2.4. 
2.2.5 Particle-wall interaction
To resolve the interaction between the particle and the wall, a simplified hard sphere approach
is adopted, where the loss of particle momentum is accounted for via restitution coefficients in
the normal and tangential directions with respect to the observed wall. This scenario is
presented in Fig. 3.
Figure 3: The interaction between the particle and the wall.
In accordance with the hard sphere approach, the normal particle velocity after rebounding 
from the wall is
v
n ,2
= e
n
v
n ,1
, (2.34)
and the tangential particle velocity after rebounding from the wall is 
(2.32)
and the virtual mass force is
Prediction of cavitation and particle erosion in a radial divergent test section 11
----------
c
1
, c
2
 and c
3
 are constants of the C
D
  ̶  ℜ
P
 curve fit and particle Reynolds
number is defined as
ℜ
P
=
ρ d
P
| u − v
P
|
μ
.
(2.30)
The force resulting from buoyancy and gravity is
F
B
+ F
G
=
d
P
3
6
g
(
ρ
P
− ρ
)
,
(2.31)
the pressure gradient force is 
F
PG
=
d
P
3
6
ρ
D u
Dt
,
(2.32)
and the virtual mass force is 
F
V M
=
1
2
d
P
3
6
ρ
(
D u
Dt
−
d v
P
dt
)
. (2.32)
Two-way coupling between the continuous and discrete phase is achieved by including the 
momentum source term in equation 2.4. 
2.2.5 Particle-wall interaction
To resolve the interaction between the particle and the wall, a simplified hard sphere approach
is adopted, where the loss of particle momentum is accounted for via restitution coefficients in
the normal and tangential directions with respect to the observed wall. This scenario is
presented in Fig. 3.
Figure 3: The interaction between the particle and the wall.
In accordance with the hard sphere approach, the normal particle velocity after rebounding 
from the wall is
v
n ,2
= e
n
v
n ,1
, (2.34)
and the tangential particle velocity after rebounding from the wall is 
(2.33)
Two-way coupling between the continuous and discrete phase is achieved by including the mo-
mentum source term in equation 2.4.
2.2.6 Particle-wall interaction
To resolve the interaction between the particle and the wall, a simplified hard sphere approach is 
adopted, where the loss of particle momentum is accounted for via restitution coefficients in the 
normal and tangential directions with respect to the observed wall. This scenario is presented 
in Fig. 3.
Figure 3: The interaction between the particle6 and the wall

JET  53
Prediction of cavitation and particle erosion in a radial divergent test section
In accordance with the hard sphere approach, the normal particle velocity after rebounding from 
the wall is
Prediction of cavitation and particle erosion in a radial divergent test section 11
----------
c
1
, c
2
 and c
3
 are constants of the C
D
  ̶  ℜ
P
 curve fit and particle Reynolds
number is defined as
ℜ
P
=
ρ d
P
| u − v
P
|
μ
.
(2.30)
The force resulting from buoyancy and gravity is
F
B
+ F
G
=
d
P
3
6
g
(
ρ
P
− ρ
)
,
(2.31)
the pressure gradient force is 
F
PG
=
d
P
3
6
ρ
D u
Dt
,
(2.32)
and the virtual mass force is 
F
V M
=
1
2
d
P
3
6
ρ
(
D u
Dt
−
d v
P
dt
)
. (2.32)
Two-way coupling between the continuous and discrete phase is achieved by including the 
momentum source term in equation 2.4. 
2.2.5 Particle-wall interaction
To resolve the interaction between the particle and the wall, a simplified hard sphere approach
is adopted, where the loss of particle momentum is accounted for via restitution coefficients in
the normal and tangential directions with respect to the observed wall. This scenario is
presented in Fig. 3.
Figure 3: The interaction between the particle and the wall.
In accordance with the hard sphere approach, the normal particle velocity after rebounding 
from the wall is
v
n ,2
= e
n
v
n ,1
, (2.34)
and the tangential particle velocity after rebounding from the wall is 
(2.34)
and the tangential particle velocity after rebounding from the wall is
12 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
v
t ,2
= e
t
v
t ,1
(2.35)
where e
n
 is the normal coefficient of restitution, e
t
 is the tangential coefficent of 
restitution, v
n ,1
 is the particle normal velocity before interaction with the wall, and v
t ,1
 
is the particle tangential velocity before interaction with the wall.
Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the 
impact angle. In this work model by Grant and Tabakoff [41], the formulae 
e
n
= 0.993 − 1.76 θ + 1.56 θ
2
− 0.49 θ
3
, (2.36)
e
t
= 0. 9 88 − 1 . 66 θ + 2. 11 6 θ
2
− 0 . 67 θ
3
, (2.37)
are used, where θ is the particle wall impact angle.
2.2.6 Particle erosion modelling
For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered
to be steel. Erosion rate of a surface element with units of [kg/m
2
] is obtained by totalling the
empirical expression for single particle erosion over all of the particles that have impacted the
observed surface element
ER =
∑
N
p
m
p
ρ
fa ce
S
fa ce
A v
1
b
f ( θ ) , (2.38)
where m
p
 is particle mass, ρ
fa ce
 i s t h e s u r f a c e d e n s i t y ( f o r s t e e l , a v a l u e o f
7800 kg / m
3
 was used), S
f a ce
 is the surface element area, A is the empirical
constant with the value of 2∙10
-9
, v
1
 is the particle impact velocity as shown in Fig. 3, b
is the material dependent velocity exponent with the value of 2.6, and f ( θ ) is the
dimensionless function of the particle impact angle. For steels, the dimensionless function is
given by
f ( θ ) = ∑
i = 1
8
( − 1 )
i + 1
B
i
α
i
, (2.39)
where the model constants B
i
 are given in Table 2.
Table 2: DNV dimensionless impact angle function model constants.
B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8
9.370 42.295
110.86
4
175.804 170.137 98.398 31.211 4.170
2.3 Boundary conditions and simulation setup
(2.35)
where e
n
 is the normal coefficient of restitution, e
t
 is the tangential coefficent of restitution, v
n,1
 
is the particle normal velocity before interaction with the wall, and v
t,1
 is the particle tangential 
velocity before interaction with the wall.
Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the 
impact angle. In this work model by Grant and Tabakoff [41], the formulae
12 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
v
t ,2
= e
t
v
t ,1
(2.35)
where e
n
 is the normal coefficient of restitution, e
t
 is the tangential coefficent of 
restitution, v
n ,1
 is the particle normal velocity before interaction with the wall, and v
t ,1
 
is the particle tangential velocity before interaction with the wall.
Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the 
impact angle. In this work model by Grant and Tabakoff [41], the formulae 
e
n
= 0.993 − 1.76 θ + 1.56 θ
2
− 0.49 θ
3
, (2.36)
e
t
= 0. 9 88 − 1 . 66 θ + 2. 11 6 θ
2
− 0 . 67 θ
3
, (2.37)
are used, where θ is the particle wall impact angle.
2.2.6 Particle erosion modelling
For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered
to be steel. Erosion rate of a surface element with units of [kg/m
2
] is obtained by totalling the
empirical expression for single particle erosion over all of the particles that have impacted the
observed surface element
ER =
∑
N
p
m
p
ρ
fa ce
S
fa ce
A v
1
b
f ( θ ) , (2.38)
where m
p
 is particle mass, ρ
fa ce
 i s t h e s u r f a c e d e n s i t y ( f o r s t e e l , a v a l u e o f
7800 kg / m
3
 was used), S
f a ce
 is the surface element area, A is the empirical
constant with the value of 2∙10
-9
, v
1
 is the particle impact velocity as shown in Fig. 3, b
is the material dependent velocity exponent with the value of 2.6, and f ( θ ) is the
dimensionless function of the particle impact angle. For steels, the dimensionless function is
given by
f ( θ ) = ∑
i = 1
8
( − 1 )
i + 1
B
i
α
i
, (2.39)
where the model constants B
i
 are given in Table 2.
Table 2: DNV dimensionless impact angle function model constants.
B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8
9.370 42.295
110.86
4
175.804 170.137 98.398 31.211 4.170
2.3 Boundary conditions and simulation setup
(2.36)
12 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
v
t ,2
= e
t
v
t ,1
(2.35)
where e
n
 is the normal coefficient of restitution, e
t
 is the tangential coefficent of 
restitution, v
n ,1
 is the particle normal velocity before interaction with the wall, and v
t ,1
 
is the particle tangential velocity before interaction with the wall.
Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the 
impact angle. In this work model by Grant and Tabakoff [41], the formulae 
e
n
= 0.993 − 1.76 θ + 1.56 θ
2
− 0.49 θ
3
, (2.36)
e
t
= 0. 9 88 − 1 . 66 θ + 2. 11 6 θ
2
− 0 . 67 θ
3
, (2.37)
are used, where θ is the particle wall impact angle.
2.2.6 Particle erosion modelling
For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered
to be steel. Erosion rate of a surface element with units of [kg/m
2
] is obtained by totalling the
empirical expression for single particle erosion over all of the particles that have impacted the
observed surface element
ER =
∑
N
p
m
p
ρ
fa ce
S
fa ce
A v
1
b
f ( θ ) , (2.38)
where m
p
 is particle mass, ρ
fa ce
 i s t h e s u r f a c e d e n s i t y ( f o r s t e e l , a v a l u e o f
7800 kg / m
3
 was used), S
f a ce
 is the surface element area, A is the empirical
constant with the value of 2∙10
-9
, v
1
 is the particle impact velocity as shown in Fig. 3, b
is the material dependent velocity exponent with the value of 2.6, and f ( θ ) is the
dimensionless function of the particle impact angle. For steels, the dimensionless function is
given by
f ( θ ) = ∑
i = 1
8
( − 1 )
i + 1
B
i
α
i
, (2.39)
where the model constants B
i
 are given in Table 2.
Table 2: DNV dimensionless impact angle function model constants.
B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8
9.370 42.295
110.86
4
175.804 170.137 98.398 31.211 4.170
2.3 Boundary conditions and simulation setup
(2.37)
are used, where is the particle wall impact angle.
2.2.7 Particle erosion modelling
For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered 
to be steel. Erosion rate of a surface element with units of [kg/m
2
] is obtained by totalling the 
empirical expression for single particle erosion over all of the particles that have impacted the 
observed surface element
12 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
v
t ,2
= e
t
v
t ,1
(2.35)
where e
n
 is the normal coefficient of restitution, e
t
 is the tangential coefficent of 
restitution, v
n ,1
 is the particle normal velocity before interaction with the wall, and v
t ,1
 
is the particle tangential velocity before interaction with the wall.
Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the 
impact angle. In this work model by Grant and Tabakoff [41], the formulae 
e
n
= 0.993 − 1.76 θ + 1.56 θ
2
− 0.49 θ
3
, (2.36)
e
t
= 0. 9 88 − 1 . 66 θ + 2. 11 6 θ
2
− 0 . 67 θ
3
, (2.37)
are used, where θ is the particle wall impact angle.
2.2.6 Particle erosion modelling
For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered
to be steel. Erosion rate of a surface element with units of [kg/m
2
] is obtained by totalling the
empirical expression for single particle erosion over all of the particles that have impacted the
observed surface element
ER =
∑
N
p
m
p
ρ
fa ce
S
fa ce
A v
1
b
f ( θ ) , (2.38)
where m
p
 is particle mass, ρ
fa ce
 i s t h e s u r f a c e d e n s i t y ( f o r s t e e l , a v a l u e o f
7800 kg / m
3
 was used), S
f a ce
 is the surface element area, A is the empirical
constant with the value of 2∙10
-9
, v
1
 is the particle impact velocity as shown in Fig. 3, b
is the material dependent velocity exponent with the value of 2.6, and f ( θ ) is the
dimensionless function of the particle impact angle. For steels, the dimensionless function is
given by
f ( θ ) = ∑
i = 1
8
( − 1 )
i + 1
B
i
α
i
, (2.39)
where the model constants B
i
 are given in Table 2.
Table 2: DNV dimensionless impact angle function model constants.
B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8
9.370 42.295
110.86
4
175.804 170.137 98.398 31.211 4.170
2.3 Boundary conditions and simulation setup
(2.38)
where m
P
 is particle mass, ρ
face
 is the surface density (for steel, a value 7800 kg/m
3
 of was used), 
S
face
 is the surface element area, A is the empirical constant with the value of 2∙10
-9
, v
1
 is the par-
ticle impact velocity as shown in Fig. 3, b is the material dependent velocity exponent with the 
value of 2.6, and f( θ) is the dimensionless function of the particle impact angle. For steels, the 
dimensionless function is given by
12 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
v
t ,2
= e
t
v
t ,1
(2.35)
where e
n
 is the normal coefficient of restitution, e
t
 is the tangential coefficent of 
restitution, v
n ,1
 is the particle normal velocity before interaction with the wall, and v
t ,1
 
is the particle tangential velocity before interaction with the wall.
Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the 
impact angle. In this work model by Grant and Tabakoff [41], the formulae 
e
n
= 0.993 − 1.76 θ + 1.56 θ
2
− 0.49 θ
3
, (2.36)
e
t
= 0. 9 88 − 1 . 66 θ + 2. 11 6 θ
2
− 0 . 67 θ
3
, (2.37)
are used, where θ is the particle wall impact angle.
2.2.6 Particle erosion modelling
For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered
to be steel. Erosion rate of a surface element with units of [kg/m
2
] is obtained by totalling the
empirical expression for single particle erosion over all of the particles that have impacted the
observed surface element
ER =
∑
N
p
m
p
ρ
fa ce
S
fa ce
A v
1
b
f ( θ ) , (2.38)
where m
p
 is particle mass, ρ
fa ce
 i s t h e s u r f a c e d e n s i t y ( f o r s t e e l , a v a l u e o f
7800 kg / m
3
 was used), S
f a ce
 is the surface element area, A is the empirical
constant with the value of 2∙10
-9
, v
1
 is the particle impact velocity as shown in Fig. 3, b
is the material dependent velocity exponent with the value of 2.6, and f ( θ ) is the
dimensionless function of the particle impact angle. For steels, the dimensionless function is
given by
f ( θ ) = ∑
i = 1
8
( − 1 )
i + 1
B
i
α
i
, (2.39)
where the model constants B
i
 are given in Table 2.
Table 2: DNV dimensionless impact angle function model constants.
B 1 B 2 B 3 B 4 B 5 B 6 B 7 B 8
9.370 42.295
110.86
4
175.804 170.137 98.398 31.211 4.170
2.3 Boundary conditions and simulation setup
(2.39)
where the model constants B
i
 are given in Table 2.
Table 2: DNV dimensionless impact angle function model constants
B
1
B
2
B
3
B
4
B
5
B
6
B
7
B
8
9.370 42.295 110.864 175.804 170.137 98.398 31.211 4.170

54 JET 54 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
2.3 Boundary conditions and simulation setup
Boundaries of the computational domain are shown and labelled in Fig. 4, while detailed bound-
ary conditions are presented in Table 3. The direction of gravitational acceleration was assigned 
in the z-direction.
For the Lagrangian tracking of particle, a Dispersed Phase Model (DPM) was used, where indi-
vidual particles can be grouped into parcels to reduce the computational effort. However, this 
behaviour was overwritten by assigning only one particle to each parcel. Along with the inlet 
velocity, presented in Table 3, it is also necessary to prescribe the particle mass flow rate, from 
which the number of particles entering the domain is calculated. Particle mass flow rate of 3.897 
kg/s was chosen to achieve the particle mass concentration of 0.05 kg/m
3
 which was reported by 
Gregorc [54] to be the insoluble sediment mass concentration of the river Drava.
Figure 4: Boundaries of the computational domain
Table 3: Boundary conditions
Boundary Boundary condition
Homogeneous mixture phase Dispersed phase
Inlet Normal velocity = 31 m/s Normal velocity = 31 m/s
Outlet Static pressure = 10.1 bar Opening (to leave the domain)
Wall No-slip wall Rebound; Grant & Tabakoff
Symmetry Symmetry Rebound; Grant & Tabakoff
The Navier-Stokes equations were solved using the SIMPLE algorithm. To evaluate the gradients, 
the Least Squares Cell-Based method was selected, while the pressure cell face values were in-
terpolated by the PRESTO! interpolation scheme. For spatial discretization, QUICK spatial discre-
tization scheme was selected for all equations. For temporal discretization of the unsteady terms 
in equations, the Bounded Second-Order Implicit time integration scheme was selected.

JET  55
Prediction of cavitation and particle erosion in a radial divergent test section
A time step of 1∙10
-6
 was chosen based on the previous cavitation study by the authors. The time 
step was checked and was found to be smaller than the particle relaxation time expressed as
14 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
A time step of 1∙10
-6
 was chosen based on the previous cavitation study by the authors. The time
step was checked and was found to be smaller than the particle relaxation time expressed as
τ
r
=
ρ
p
d
p
2
18 μ
24
C
D
ℜ
p
(2.40)
A scaled residuals convergence criterion of 1∙10
-3 
was achieved before the imposed limit of 100
iterations per time step. Overall, 0.016815 s of physical time were simulated.
3 RESULTS AND DISCUSSION
Erosion due to the cavitation and particles is evaluated on the bottom wall (when the geometry
is positioned as shown in Fig. 4) of the radial divergent test section. Figure 5. presents the
erosion prediction in terms of coloured contours, where in Fig. 5 a) coloured contours represent
the particle erosion results from equation 2.38 and in Fig. 5 b) coloured contours represent the
cavitation erosion results from equation 2.24. 
The main zone where cavitation erosion occurs is located between 19 mm and 32 mm from the
axis of symmetry, which agrees with the experimental results [28] and numerical simulation
results from different authors [31], [32], [51]. Since in these studies particles were not present in
the flow, we conclude that in this study particles do not change the flow behaviour to the extent
that it would change the cavitation erosion zone. It must be noted; however , that this is not
necessarily true for particles of a different diameter, density and volume fraction. 
Maximum particle erosion occurs downstream from the cavitation erosion zone. This can be
explained by observing the flow streamlines as shown in Fig. 6 a) and Fig. 6 b). Vapour cavities
disturb the flow in the downstream direction, multiple vortices can be seen in Fig 6. a) and Fig 6.
b). The relaxation time of the particles considered in this study is smaller than the flow
relax ation time, therefore particles respond quickly to any changes in the flow , such as the
change of flow direction due to a vortex. The vortices observed in Fig. 6 a) and Fig. 6. b) redirect
particles toward the bottom wall, causing them to impact the wall at a higher impact angle. In
DNV erosion model the dimensionless function of the particle impact angle (equation 2.39)
rapidly approaches zero as the impact angle approaches zero. From Fig. 6 a) and Fig. 6 b), it can
be seen that the streamlines are virtually parallel to the bottom wall in the region of cavitation,
therefore there is almost no particle erosion predicted (at least an order of magnitude smaller
than in the main zone between 40 mm and 52 mm from the symmetry axis).  
(2.40)
A scaled residuals convergence criterion of 1∙10
-3 
was achieved before the imposed limit of 100 
iterations per time step. Overall, 0.016815 s of physical time were simulated.
3 RESULTS AND DISCUSSION
Erosion due to the cavitation and particles is evaluated on the bottom wall (when the geome-
try is positioned as shown in Fig. 4) of the radial divergent test section. Figure 5. presents the 
erosion prediction in terms of coloured contours, where in Fig. 5 a) coloured contours represent 
the particle erosion results from equation 2.38 and in Fig. 5 b) coloured contours represent the 
cavitation erosion results from equation 2.24.
The main zone where cavitation erosion occurs is located between 19 mm and 32 mm from the 
axis of symmetry, which agrees with the experimental results [28] and numerical simulation re-
sults from different authors [31], [32], [51]. Since in these studies particles were not present in 
the flow, we conclude that in this study particles do not change the flow behaviour to the extent 
that it would change the cavitation erosion zone. It must be noted; however, that this is not nec-
essarily true for particles of a different diameter, density and volume fraction.
Maximum particle erosion occurs downstream from the cavitation erosion zone. This can be 
explained by observing the flow streamlines as shown in Fig. 6 a) and Fig. 6 b). Vapour cavities 
disturb the flow in the downstream direction, multiple vortices can be seen in Fig 6. a) and Fig 6. 
b). The relaxation time of the particles considered in this study is smaller than the flow relaxation 
time, therefore particles respond quickly to any changes in the flow, such as the change of flow 
direction due to a vortex. The vortices observed in Fig. 6 a) and Fig. 6. b) redirect particles toward 
the bottom wall, causing them to impact the wall at a higher impact angle. In DNV erosion model 
the dimensionless function of the particle impact angle (equation 2.39) rapidly approaches zero 
as the impact angle approaches zero. From Fig. 6 a) and Fig. 6 b), it can be seen that the stream-
lines are virtually parallel to the bottom wall in the region of cavitation, therefore there is almost 
no particle erosion predicted (at least an order of magnitude smaller than in the main zone be-
tween 40 mm and 52 mm from the symmetry axis).

56 JET 56 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
Figure 5: Erosion prediction on the bottom wall: a) particle erosion, b) cavitation erosion

JET  57
Prediction of cavitation and particle erosion in a radial divergent test section
Figure 6: Flow dynamics: a), b) cavitation represented by the iso-surface of vapour volume frac-
tion of 0.2 and flow dynamics represented by streamlines coloured by the flow velocity at two 
different times, c) time evolution of total vapour volume, with selected times highlighted with a 
red square
4 CONCLUSIONS
In this paper, a transient simulation of the cavitating and particle-laden flow is presented. Risk 
of erosion due to both the cavitation and the particles is assessed by employing the appropriate 
erosion modelling approaches.
By using the erosion models, it is found that for the investigated particles, the flow is not suf-
ficiently affected by the presence of the particles to change the location of the main cavitation 
erosion zone, indicating that a one-way momentum coupling between the continuous phase and 
the discrete phase would be sufficient.
A separate zone of pure particle erosion is predicted downstream of the cavitation erosion zone. 
This is explained by observing the formation of vortices that redirect particles towards the wall 
downstream of the vapour clouds.
The particles investigated in this study could be particularly useful for an experimental validation 
as there is no overlap between the cavitation erosion zone and the particle erosion zone, which 
would require the ability to distinguish one type of erosion from the other.

58 JET 58 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
5 ACKNOWLEDGEMENTS
The authors wish to thank the Slovenian Research Agency (ARRS) for their financial support with-
in the framework of the Research Programme P2-0196 Research in Power, Process and Environ-
mental Engineering.
References
[1] Lord Rayleigh: VIII. On the pressure developed in a liquid during the collapse of a spher-
ical cavity, Philosophical Magazine Series 6, Vol. 34, No. 200, pp. 94–98, 1917, doi: 
http://dx.doi.org/10.1080/14786440808635681
[2] F. Hammit: Observations on cavitation damage in a flowing system, ASME Journal of 
Basic Engineering, Vol. 85, No. 3, 1963
[3] A. Vogel, W. Lauterborn, and R. Timm: Optical and acoustic investigations of the 
dynamics of laser-produced cavitation bubbles near a solid boundary, Journal of 
Fluid Mechanics, Vol. 206, No. September 1989, pp. 299–338, 1989, doi: 10.1017/
S0022112089002314
[4] M. S. Plesset and R. B. Chapman: Collapse of an Initially Spherical Vapor Cavity in the 
Neighborhood of a Solid Boundary, Journal of Fluid Mechanics, Vol. 47, part 2, pp. 283–
290, 1971
[5] F. Pereira, F. Avellan, and P. H. Dupont: Prediction of cavitation erosion: An energy 
approach, Journal of Fluids Engineering, Transactions of the ASME, Vol. 120, No. 4, pp. 
719–727, 1998, doi: 10.1115/1.2820729
[6] R. Fortes Patella, J.-L. Reboud, and L. Briancon Marjollet: A Phenomenological and nu-
merical model for scaling the flow agressiveness in cavitation erosion, 2004. [Online]. 
Available: https://www.researchgate.net/publication/281921326
[7] T. Melissaris, N. Bulten, and T. J. C. Van Terwisga: On the applicability of cavitation 
erosion risk models with a URANS solver, Journal of Fluids Engineering, Transactions of 
the ASME, Vol. 141, No. 10, 2019, doi: 10.1115/1.4043169
[8] A. Philipp and W . Lauterborn: Cavitation erosion by single laser-produced bubbles, Jour-
nal of Fluid Mechanics, Vol. 361, pp. 75–116, 1998, doi: 10.1017/S0022112098008738
[9] C. Leclercq, A. Archer, R. Fortes-Patella, and F. Cerru: Numerical cavitation intensi-
ty on a hydrofoil for 3D homogeneous unsteady viscous flows, International Journal 
of Fluid Machinery and Systems, Vol. 10, No. 3, pp. 254–263, 2017, doi: 10.5293/
IJFMS.2017.10.3.254
[10] A. Van Oosterom and J. Strackee: The Solid Angle of a Plane Triangle, IEEE Transac-
tions on Biomedical Engineering, Vol. BME-30, No. 2, pp. 125–126, 1983, doi: 10.1109/
TBME.1983.325207
[11] S. Schenke and T. J. C. van Terwisga: An energy conservative method to predict the 
erosive aggressiveness of collapsing cavitating structures and cavitating flows from nu-

JET  59
Prediction of cavitation and particle erosion in a radial divergent test section
merical simulations, International Journal of Multiphase Flow, Vol. 111, pp. 200–218, 
2019, doi: 10.1016/j.ijmultiphaseflow.2018.11.016
[12] M. Tinguely, D. Obreschkow, P. Kobel, N. Dorsaz, A. De Bosset, and M. Farhat: Energy 
partition at the collapse of spherical cavitation bubbles, Physical Review E–Statistical, 
Nonlinear, and Soft Matter Physics, Vol. 86, No. 4, pp. 1–6, 2012, doi: 10.1103/Phys-
RevE.86.046315
[13] Y. C. Wang and C. E. Brennen: Numerical computation of shock waves in a spherical 
cloud of cavitation bubbles, Journal of Fluids Engineering, Transactions of the ASME, 
Vol. 121, No. 4, pp. 872–880, 1999, doi: 10.1115/1.2823549
[14] S. Schenke, T. Melissaris, and T. J. C. Van Terwisga: On the relevance of kinematics for 
cavitation implosion loads, Physics of Fluids, Vol. 31, No. 5, 2019, doi: 10.1063/1.5092711
[15] T. Melissaris, S. Schenke, N. Bulten, and T. J. C. van Terwisga: On the accuracy of pre-
dicting cavitation impact loads on marine propellers, Wear, Vol. 456–457, No. June, p. 
203393, 2020, doi: 10.1016/j.wear.2020.203393
[16] M. Pezdevsek, L. Kevorkijan, and I. Bilus: Cavitation Erosion Modelling – Comparison of 
Two Solid Angle Projection Approaches, International Journal of Simulation Modelling, 
Vol. 21, No. 2, pp. 249–260, 2022, doi: 10.2507/ijsimm21-2-600
[17] P . J. Zwart, A. G. Gerber, and T. Belamri: A two-phase flow model for predicting cavita-
tion dynamics, International Conference on Multiphase Flow, No. January 2004, p. 152, 
2004
[18] G. H. Schnerr and J. Sauer: Physical and Numerical Modeling of Unsteady Cavitation 
Dynamics, 4th International Conference on Multiphase Flow, No. June, pp. 1–12, 2001
[19] M. P. Kinzel, J. W. Lindau, and R. F. Kunz: A unified homogenous multiphase CFD mod-
el for cavitation, American Society of Mechanical Engineers, Fluids Engineering Divi-
sion (Publication) FEDSM, Vol. 1B-2017, No. January 2018, 2017, doi: 10.1115/FED-
SM2017-69363
[20] I. Bilus, M. Morgut, and E. Nobile: Simulation of sheet and cloud cavitation with ho-
mogenous transport models, International Journal of Simulation Modelling, Vol. 12, 
No. 2, pp. 94–106, 2013, doi: 10.2507/IJSIMM12(2)3.229
[21] E. Goncalves and R. F. Patella: Numerical simulation of cavitating flows with homo-
geneous models, Computers and Fluids, Vol. 38, No. 9, pp. 1682–1696, 2009, doi: 
10.1016/j.compfluid.2009.03.001
[22] M. Morgut and E. Nobile: Influence of the Mass Transfer Model on the Numerical Pre-
diction of the Cavitating Flow Around a Marine Propeller, Second International Sympo-
sium on Marine Propulsors smp’11, No. June, pp. 1–8, 2011
[23] A. Niedzwiedzka, G. H. Schnerr, and W. Sobieski: Review of numerical models of cavi-
tating flows with the use of the homogeneous approach, Archives of Thermodynamics, 
Vol. 37, No. 2, pp. 71–88, 2016, doi: 10.1515/aoter-2016-0013
[24] M. Pezdevšek, I. Biluš, and G. Hren: COMPARISON OF CAVITATION MODELS FOR THE 
PREDICTION OF CAVITATION AROUND A HYDROFOIL, Journal of Energy Technology, Vol. 
14, No. 1, pp. 41–55, 2021

60 JET 60 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
[25] T. Schenke, Sören; van Terwisga: Finite Mass Transfer Effects in Cavitation Modelling, 
Proceedings of the 19th Numerical Towing Tank Symposium, pp. 1–6, 2015
[26] W. Jian, M. Petkovšek, L. Houlin, B. Širok, and M. Dular: Combined numerical and 
experimental investigation of the cavitation erosion process, Journal of Fluids Engineer-
ing, Transactions of the ASME, Vol. 137, No. 5, 2015, doi: 10.1115/1.4029533
[27] J. B. Carrat, R. Fortes-Patella, and J. P . Franc: Experimental and numerical investigation 
of the erosive potential of a leading edge cavity, International Journal of Fluid Machin-
ery and Systems, Vol. 12, No. 2, pp. 136–146, 2019, doi: 10.5293/IJFMS.2019.12.2.136
[28] N. Berchiche, J. P. Franc, and J. M. Michel: A cavitation erosion model for ductile ma-
terials, Journal of Fluids Engineering, Transactions of the ASME, Vol. 124, No. 3, pp. 
601–606, 2002, doi: 10.1115/1.1486474
[29] M. Dular and M. Petkovšek: On the mechanisms of cavitation erosion–Coupling high 
speed videos to damage patterns, Experimental Thermal and Fluid Science, Vol. 68, No. 
June 2015, pp. 359–370, 2015, doi: 10.1016/j.expthermflusci.2015.06.001
[30] I. Biluš, M. Hočevar, M. Dular, and L. Lešnik: Numerical Prediction of Various Cavi-
tation Erosion Mechanisms, Journal of Fluids Engineering, Vol. 142, No. 4, 2020, doi: 
10.1115/1.4045365
[31] A. Peters: Numerical Modelling and Prediction of Cavitation Erosion Using Euler-Euler 
and Multi-Scale Euler-Lagrange Methods, Universität Duisburg-Essen, 2019
[32] M. H. Arabnejad, U. Svennberg, and R. E. Bensow: Numerical assessment of cavitation 
erosion risk using incompressible simulation of cavitating flows, Wear, Vol. 464–465, 
No. November 2020, p. 203529, 2021, doi: 10.1016/j.wear.2020.203529
[33] I. Finne: Erosion of surfaces, Wear, Vol. 3, pp. 87–103, 1960, doi: 10.1016/0043-
1648(60)90055-7
[34] J. G. A. Bitter: A Study of Erosion Phenomena Part I, Wear, Vol. 6, pp. 169–190, 1963
[35] J. G. A. Bitter: A study of erosion phenomena. Part II, Wear, Vol. 6, No. 3, pp. 169–190, 
1963, doi: 10.1016/0043-1648(63)90073-5
[36] Y. I. Oka, K. Okamura, and T. Yoshida: Practical estimation of erosion damage caused 
by solid particle impact: Part 1: Effects of impact parameters on a predictive equation, 
Wear, Vol. 259, No. 1–6, pp. 95–101, 2005, doi: 10.1016/j.wear.2005.01.039
[37] Y. I. Oka and T. Yoshida: Practical estimation of erosion damage caused by solid parti-
cle impact: Part 2: Mechanical properties of materials directly associated with erosion 
damage, Wear, Vol. 259, No. 1–6, pp. 102–109, 2005, doi: 10.1016/j.wear.2005.01.040
[38] K. R. Ahlert: Effects of Particle Impingement Angle and Surface Waiting on solid Partilce 
Erosion of AISI 1018 Steel, 1994
[39] Det Norske Veritas: Erosive wear in piping systems, 2007
[40] A. Forder, M. Thew, and D. Harrison: A numerical investigation of solid particle erosion 
experienced within oilfield control valves, Wear, Vol. 216, No. 2, pp. 184–193, 1998, doi: 
10.1016/S0043-1648(97)00217-2

JET  61
Prediction of cavitation and particle erosion in a radial divergent test section
[41] G. Grant and W. Tabakoff: An Experimental Investigation of the Erosive Characteristics 
of 2024 Aluminum Alloy, Springfield, 1973
[42] J. Sato, K. Usami, and T. Okamura: Basic Study of Coupled Damage Caused by 
Silt Abrasion and Cavitation Erosion, JSME International Journal, Vol. 34, No. 
3, pp. 292–297, 1991, [Online]. Available: https://www.jstage.jst.go.jp/article/
bpb1993/17/11/17_11_1460/_pdf/-char/ja
[43] K. Su, J. Wu, and D. Xia: Dual role of microparticles in synergistic cavitation–particle ero-
sion: Modeling and experiments, Wear, Vol. 470–471, p. 203633, 2021, doi: 10.1016/j.
wear.2021.203633
[44] F. Innings, E. Hultman, F. Forsberg, and B. Prakash: Understanding and analysis of wear 
in homogenizers for processing liquid food, Wear, Vol. 271, No. 9–10, pp. 2588–2598, 
2011, doi: 10.1016/j.wear.2011.01.084
[45] L. Kevorkijan, J. Ravnik, and I. Biluš: Numerično modeliranje kavitacijske erozije in 
abrazije delcev na profilu krila NACA 0015, in Kuhljevi dnevi 2021, 2021, No. September
[46] P. J. Dunstan and S. C. Li: Cavitation enhancement of silt erosion: Numerical studies, 
Wear, Vol. 268, No. 7–8, pp. 946–954, 2010, doi: 10.1016/j.wear.2009.12.036
[47] L. A. Teran, S. Laín, and S. A. Rodríguez: Synergy effect modelling of cavitation and hard 
particle erosion: Implementation and validation, Wear, Vol. 478–479, No. December 
2020, 2021
[48] L. Kevorkijan, L. Lešnik, and I. Biluš: Cavitation Erosion Modelling on a Radial Divergent 
Test Section Using RANS, Strojniski Vestnik/Journal of Mechanical Engineering, Vol. 68, 
No. 2, pp. 71–81, 2022, doi: 10.5545/sv-jme.2021.7364
[49] D. Greif: Kavitacijska erozija v vbrizgalnih komponentah z uporabo polidisperznega ka-
vitacijskega modela, 2012
[50] M. S. Mihatsch, S. J. Schmidt, and N. A. Adams: Cavitation erosion prediction based 
on analysis of flow dynamics and impact load spectra, Physics of Fluids, Vol. 27, No. 10, 
2015, doi: 10.1063/1.4932175
[51] S. Mouvanal, D. Chatterjee, S. Bakshi, A. Burkhardt, and V. Mohr: Numerical predic-
tion of potential cavitation erosion in fuel injectors, International Journal of Multiphase 
Flow, Vol. 104, pp. 113–124, 2018, doi: 10.1016/j.ijmultiphaseflow.2018.03.005
[52] J.-L. Reboud, B. Stutz, and O. Coutier-Delgosha: Two-phase flow structure of cavita-
tion: experiment and modelling of unsteady effects, 1998
[53] S. A. Morsi and A. J. Alexander: An investigation of particle trajectories in two-phase 
flow systems, Journal of Fluid Mechanics, Vol. 55, No. 2, pp. 193–208, 1972, doi: 
10.1017/S0022112072001806
[54] B. Gregorc: VPLIV TRDNO-KAPLJEVITIH ZMESI NA OBRATOVALNE KARAKTERISTIKE, Uni-
verza v Mariboru, 2011

62 JET 62 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
22 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
μ
t
turbulent viscosity
k turbulence kinetic energy
ω specific rate of dissipation of the turbulence kinetic energy
α
❑
damping coefficient
S strain rate magnitude
F
2
second blending function
a
1
SST turbulent model constant
f ( ρ ) mixture density function
n mixture density function modifying exponent
S
α
v
source of the vapour phase
R
b
bubble radius
n
0
bubble number density
p
v
vapour pressure
V
V
vapour (cavity) volume
E
pot
potential energy of a vapour cavity
p
d
pressure in the surrounding liquid
P
pot
cavitation potential power
p ( t ) instantaneous pressure
t
❑
current simulated time
´ e
p ot
cavitation erosion potential
V
cel l
computational mesh cell volume
ε collapse induced kinetic energy
u
i
collapse induced velocity
´ e
rad
( t ) instantaneous radiated power of the pressure wave
P
u
projection operator
K global energy conservation parameter
β energy release criterion
∆ t time step
Prediction of cavitation and particle erosion in a radial divergent test section 21
----------
[49] D. Greif: Kavitacijska erozija v vbrizgalnih komponentah z uporabo polidisperznega
kavitacijskega modela, 2012.
[50] M. S. Mihatsch, S. J. Schmidt, and N. A. Adams: Cavitation erosion prediction based on
analysis of flow dynamics and impact load spectra, Physics of Fluids, Vol. 27, No. 10,
2015, doi: 10.1063/1.4932175.
[51] S. Mouvanal, D. Chatterjee, S. Bakshi, A. Burkhardt, and V. Mohr: Numerical prediction
of potential cavitation erosion in fuel injectors, International Journal of Multiphase Flow,
Vol. 104, pp. 113–124, 2018, doi: 10.1016/j.ijmultiphaseflow.2018.03.005.
[52] J.-L. Reboud, B. Stutz, and O. Coutier-Delgosha: Two-phase flow structure of cavitation:
experiment and modelling of unsteady effects, 1998.
[53] S. A. Morsi and A. J. Alexander : An investigation of particle trajectories in two-phase
flow systems, Journal of Fluid Mechanics, Vol. 55, No. 2, pp. 193–208, 1972, doi:
10.1017/S0022112072001806.
[54] B. Gregorc: VPLIV TRDNO-KAPLJEVITIH ZMESI NA OBRATOVALNE KARAKTERISTIKE,
Univerza v Mariboru, 2011.
N o m en cl at u r e
(Symbols) (Symbol meaning)
ρ mixture density
μ mixture viscosity
α vapour volume fraction
ρ
l
liquid density
ρ
v
vapour density
μ
l
liquid viscosity
μ
v
vapour viscosity
t time
u mixture velocity field
p pressure
g gravitational acceleration
S
M
momentum source term

JET  63
Prediction of cavitation and particle erosion in a radial divergent test section
22 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
μ
t
turbulent viscosity
k turbulence kinetic energy
ω specific rate of dissipation of the turbulence kinetic energy
α
❑
damping coefficient
S strain rate magnitude
F
2
second blending function
a
1
SST turbulent model constant
f ( ρ ) mixture density function
n mixture density function modifying exponent
S
α
v
source of the vapour phase
R
b
bubble radius
n
0
bubble number density
p
v
vapour pressure
V
V
vapour (cavity) volume
E
pot
potential energy of a vapour cavity
p
d
pressure in the surrounding liquid
P
pot
cavitation potential power
p ( t ) instantaneous pressure
t
❑
current simulated time
´ e
p ot
cavitation erosion potential
V
cel l
computational mesh cell volume
ε collapse induced kinetic energy
u
i
collapse induced velocity
´ e
rad
( t ) instantaneous radiated power of the pressure wave
P
u
projection operator
K global energy conservation parameter
β energy release criterion
∆ t time step
Prediction of cavitation and particle erosion in a radial divergent test section 23
----------
´ e
S
instantaneous surface specific impact power
x
cell
position vector of the computational cell centre
x
S
position vector of the wall surface element centre
n surface element normal
e
S
instantaneous surface specific impact energy
v
P
particle velocity
x
P
particle position in the global inertial frame of reference
ρ
P
particle density
d
P
particle diameter
F resultant force acting on the particle
F
D
drag force
F
B
buoyancy force
F
G
force of gravity
F
PG
pressure gradient force
F
V M
virtual mass force
C
D
drag coefficient
c
1
first constant of the C
D
  ̶  ℜ
P
 curve fit
c
2
second constant of the C
D
  ̶  ℜ
P
 curve fit
c
3
third constant of the C
D
  ̶  ℜ
P
 curve fit
ℜ
P
particle Reynolds number
v
n ,1
particle normal velocity before interaction with the wall
v
n ,2
particle normal velocity after rebounding from the wall
e
n
normal coefficient of restitution
v
t ,1
particle tangential velocity before interaction with the wall
v
t ,2
particle tangential velocity after rebounding from the wall
e
t
tangential coefficient of restitution
θ particle wall impact angle
ER erosion rate

64 JET 64 JET
Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
Prediction of cavitation and particle erosion in a radial divergent test section 23
----------
´ e
S
instantaneous surface specific impact power
x
cell
position vector of the computational cell centre
x
S
position vector of the wall surface element centre
n surface element normal
e
S
instantaneous surface specific impact energy
v
P
particle velocity
x
P
particle position in the global inertial frame of reference
ρ
P
particle density
d
P
particle diameter
F resultant force acting on the particle
F
D
drag force
F
B
buoyancy force
F
G
force of gravity
F
PG
pressure gradient force
F
V M
virtual mass force
C
D
drag coefficient
c
1
first constant of the C
D
  ̶  ℜ
P
 curve fit
c
2
second constant of the C
D
  ̶  ℜ
P
 curve fit
c
3
third constant of the C
D
  ̶  ℜ
P
 curve fit
ℜ
P
particle Reynolds number
v
n ,1
particle normal velocity before interaction with the wall
v
n ,2
particle normal velocity after rebounding from the wall
e
n
normal coefficient of restitution
v
t ,1
particle tangential velocity before interaction with the wall
v
t ,2
particle tangential velocity after rebounding from the wall
e
t
tangential coefficient of restitution
θ particle wall impact angle
ER erosion rate
24 Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš JET Vol. 15 (2022)
Issue 2
----------
m
p
particle mass
ρ
fa ce
density of the wall surface
S
f ac e
surface element area
A particle erosion empirical constant
v
1
particle impact velocity
b particle erosion material dependent velocity exponent
f ( θ ) dimensionless function of the particle impact angle
B
i
DNV particle erosion model constants
τ
r
particle relaxation time