Numerical prediction of cavitation phenomena on marine vessel: Effect of the water environment profile on the propulsion performance

2.1 Shipping in NSR and Suez canal route (SCR)

One of the extreme areas currently possible to conquer is the Arctic Ocean in the North Pole. Arctic shipping routes are the maritime route connected entirely of both east–west land and the Atlantic–Pacific Ocean. Located in the north polar of the Earth, it is characterized as a zone with a cold climate and secured by ice. This environment has been untouchable and nearly undisturbed by human activities for a long time. Recently, the route through parts of the Arctic is conceivable since the worldwide warming impact has made a more critical portion of the Arctic melting within the summer. Mainly it happens in March, as shown in Figure 1 [27].

Figure 1

The total sea ice extent in March (left) and September 2012 (right). Purple lines reflect the median extent from 1979 to 2000 [29].

Maritime container cargo contains an extraordinary potential request for the NSR shipping route, which is currently shipped through the SCR between East Asia and Northwest Europe (Figure 2). Vehicles are one of the foremost critical potential cargoes. Furthermore, other prospective cargoes, such as Murmansk gas condensate, Hammerfest liquefied natural gas (LNG), and Kirkenes iron ore, had already been commercially exported since 2009 to East Asia [28]. As global warming has retreated from the Arctic seas, NSR trade in East Asia and Northwest Europe lately has started to gain momentum, leading to the direct result of a 40% reduction in sailing distances relative to the SCR. The NSR is expected to take 19.3 days to transport, which is 35.4% quicker than the SCR (30.4 days).

Figure 2

Comparison of shipping between NSR and SCR routes Hamburg to Yokohama [28].

The reduced effect of CO₂ pollution due to decreased distance is estimated between 13 and 35% for the NSR duty period of 105–225 days. This issue can entice the operators and shipowners based on the environmental benefits of shipping rather than the financial benefits [28]. Lindstad et al. [27] have investigated the combination of the emission obtained with the specified global warming potential (GWP) factors for the last 20 years. Twenty years of assessment based on relevancy to identify fuel and technology relations are used, resulting in the best climate impact reduction. It is compatible with the 2050 target of maintaining temperature stability at 2°C above pre-industrial levels.

It can be found that the positive values shown in Table 1 indicate a warming impact, and negative values have a cooling effect. Figure 3 exemplifies comparison data of trades through SCR and trades through NSR. This study uses Panamax dry bulkers for both shipping trades. The fuel options are light fuel oil (LFO), marine gas oil, and LNG. Each fuel option is divided into low power and high power. The bar shows that NSR yields the lowest cost compared to 18 USD/ton versus 29 USD/ton for NSR and SCR, respectively. However, compared to CO₂ equivalents (GWP₂₀), trades through NSR have a more significant climate effect per ton transported than trades through SCR. The LFO used in NSR trades, for example, is 65 kg CO₂-eq per ton transported, compared to 20 kg CO₂-eq in SCR trades.

Figure 3

Suez versus NSR trades cost and kg CO₂ equivalent (GWP 20) per ton transported per fuel type [27].

2.2 Water characteristic in ocean

The water characteristics in the ocean are dependent on the temperature. The natural water phase is liquid, but it could turn to a solid phase and gaseous/vapor phase determined from heat energy transferred by sunshine. Water has low thermal conductivity than air. It needs about four times the amount of energy to increase the temperature by 1°C, but it also means water is more difficult to release heat. Figure 4 depicts the heat maps of the ocean surface temperature of the Earth. The water near the equator has yellow until the red represents the warm zone and gets cold to the poles. The temperature of the world ocean surface imposes on Earth’s climate and weather. The temperature of the ocean also varies according to the depth. There is a surface layer categorized up to 200 m deep, with the same temperature. Below the surface layer is a thermocline that reaches 1,000 m, colder from top to bottom. However, in this case, the study of ship operation is only on the surface. So, it is clear that the scope of ocean water in the following discussion means surface water of the ocean.

Figure 4

Sea surface temperature [30].

Besides temperature, the salinity of seawater is an exciting topic in oceanography. Seawater has different salinity levels influenced by temperature, pressure, conductivity, etc. The spectrum of practical salinity S _p is in the range 2 ≤ S _p ≤ 42 [31]. The precise location of bulk material collection to prepare the standard seawater is not mentioned. The Standard Seawater Service issues guidance notes to ships collecting this bulk material, recommending that water be collected in deep water between longitudes 50 and 40°W during daylight hours. Millero et al. [32] created reference composition salinity S _R (or reference salinity for short) to be the best approximation in the mass fraction absolute salinity S _A of standard seawater. A seawater sample of reference composition whose practical salinity S _P is 35 has a reference salinity S _R of 35.16504 g kg⁻¹ with absolute uncertainty is ±0.007 g kg⁻¹. The difference between reference and practical salinities may be attributed to the early method of calculating salinity by evaporating water from seawater and weighing the remaining solid content. This method also evaporated certain volatile elements, accounting for the majority of the 0.16504 g kg⁻¹ salinity discrepancy.

Veronis measurements [33] revealed that after a temperature drops to −40°C, newly formed ice had a salinity of 10.16 g kg⁻¹. The ice salinity was as low as 5.64 g kg⁻¹ after the temperature dropped to −16°C. Because of the salt solution, the salinity of ice reduces with age. It stays between the crystals during freezing since it is denser than ice and therefore leaches. Temperature readings have a minimum precision of 0.01°C. Ocean temperatures range from –1.9 to 30°C. The lower boundary is defined by the accumulation of ice. However, supercooled water with a temperature below the freezing point may be encountered. Water temperature in the surface layer of a closed sea can be higher than 30°C in extreme conditions.

Cavitation phenomena cannot be separated from the evaporation and condensation process. This process depends on the temperature, pressure, water properties, and other thermodynamics properties in its specific environment. The water boiling point marks phase transformation from liquid to gaseous. Both freshwater and seawater have different boiling points because of their properties. The boiling point of water is determined by saturation pressure. In the case of this study, the saturation pressure of water is differentiated by freshwater and seawater. Figure 5 shows the difference between the boiling point and saturation point of both seawater and freshwater.

Figure 5

(a) Boiling point elevation of seawater based on salinity variation and (b) vapor pressure to temperature comparison relates to the salinity [34].

2.3 Previous cavitation research related to the temperature

Several researchers who studied the cavitation phenomena applied different techniques in approaching the cavitation phenomena. The effect of temperature on cavitation problems for various materials over a temperature range of 0–90°C has been studied by Plesset [35]. The findings showed that the optimum damage rate occurred at temperatures ranging from 40 to 50°C. Furthermore, Al-Arabi [36] used experimental observation to investigate the impact of water temperature on the efficiency and cavitation inception of a centrifugal pump. The section test was created to monitor the ratio of flow rate, suction pressure, rotational speed, and fluid temperature. The findings revealed that increasing the water temperature sped up cavitation and increased the net positive suction head.

Furthermore, an experimental test was conducted by Dular and Coutier-Delgosha [37] to research the effect of thermodynamics on a single cavitation bubble’s growth and collapse. Their research focused on studying temperature changes in liquid environments. Experiments were conducted in a water-filled cylinder at room temperature and atmospheric pressure. The findings revealed that the growth of the bubble was caused by the extension of an initial air bubble. The experimental findings were compared to the “thermal delay” model’s predictions. Temperature fluctuations are due to latent heat exchanges during the vaporization and condensation processes in this approach.

Tanaka [38] observed a centrifugal pump’s thermodynamic cavitation efficiency experimentally in 2011. The findings found that cavitation performance with liquid nitrogen was superior to cavitation performance with cold water, and the expected temperature drop due to the thermodynamic effect decreased as the flow coefficient decreased. Furthermore, at the same flow coefficient, the predicted temperature drop attributable to the thermodynamic impact on the low cavitation condition was greater than on the high cavitation efficiency impeller. Kim and Song [39] conducted experimental research on the influence of temperature on the critical cavitation number in a turbopump inducer. According to the data, nondimensional thermal parameters increased the onset to a lower cavitation number and decreased cavitation rate. However, at 340 K and 5,000 rpm, the critical cavitation number and revolving cavitation onset were independent of the nondimensional thermal parameter for more than 0.540. Yusvika et al. [40] used CFD-based numerical simulation to simulate propeller cavitation based on water physical properties. The simulation effects of propeller output characteristics, propeller cavitation patterns, and fluid state were compared at observed temperatures ranging from 0 to 50°C. The cavitation phenomenon was modeled using the RANS system with several reference frames (MRF). The obtained results revealed a cavitation pattern occurring more vigorously as water temperature increased.

3.1 Parameter of propeller performance

A parameter of propeller performances can be easily separated into open water and behind-hull properties, which have fundamental variations. This study observes the PPTC VP1304 propeller in open water by numerical simulation. The propeller thrust F (N) and torque Q ( Nm ) are expressed in a series of nondimensional characteristics, which are as follows [24]:

where K T and K Q are the dimensionless coefficients of thrust and torque, D  (m) is the diameter of the propeller, J is the advance coefficient, ρ ( kg/m 3 ) is the fluid density, η is the open water efficiency, V _A (m/s) is the freestream fluid velocity, and n   ( rps ) is the propeller rotational speed.

3.2 Calculation of cavitation event

Cavitation flow conditions were organized based on the cavitation number σ n at the specified advance coefficient of testing. The pressure coefficient c _p and cavitation number σ n can be defined as follows [24]:

where P _v is the absolute vapor pressure, P is the local pressure, and P Ref is the pressure used for reference. Based on the theory of the cavitation inception criterion, cavitation will occur if σ n ≤ c _p.

3.3 Multiphase Reynolds-averaged Navier–Stokes (RANS) method

Cavitation flow analysis is gas-liquid two-phase mixture that combine water & water vapor into a single fluid mixed-phase state. The RANS equation can solve the multiphase flow if the combination of the two liquids is homogeneous. To approach turbulent flow, the RANS equation assumes that the fluid is incompressible and employs the eddy viscous model [23]. The following is the governing equation:

This equation represents the multiphase RANS equations. The volume continuity and momentum equation for the liquid–vapor mixture is shown on the left, and on the right, the volume fraction for the liquid phase. ρ l and ρ v are liquid and vapor density, respectively. P avg represents the average pressure, U indicates the average velocity, S is the additional momentum, and γ   indicates the water volume fraction. Then, α is the vapor volume fraction:

The multiphase flow is constructed from the mixing of liquid–vapor density as well as its dynamics viscosity. Each phase’s density is expressed by a scalar volume fraction, which is as follows:

equal to:

where ρ m represents the mixture density, m ̇ is the interphase transfer rate due to cavitation, and μ m is the dynamic viscosity of the mixture. Furthermore, the eddy viscous turbulent model was used to close the governing equation model and treat the two liquid–vapor phases as one phase. The standard fully k−ε turbulence model and shear stress transport (SST) turbulence models were used for each simulation in the test case. Meanwhile, for the observation part, only the SST turbulence model was used because it considers the validation in the test case, which has superior accuracy with the SST turbulence model to k−ε.

3.4 Bubble cavity inception process and mass transfer model

The first stage of cavitation begins with the formation of bubbles. It assumes that the cavitation bubble is a homogeneous phase and analyzes only one bubble during the flows, despite it being impossible in practice. The bubbles grow larger following its pressure at the local area. The Rayleigh–Plesset equation is used to model bubble growth. It represents the increasing diameter of bubbles in a fluid as a time function. The equation measures growth in a single bubble when assuming no thermodynamic flow or temperature decrease from bubble initiation to collapse. For test and observation scenarios, the equation has been used as a mass transfer model in CFX, and it is as follows [41]:

where R B represents the bubble radius, P v indicates the pressure of the bubble at the fluid reference temperature, and   P is the pressure around the bubble. ρ l and 2 σ reflect the fluid density and surface tension coefficient of water vapor and fluid, respectively. When the second derivative is removed, the equation is reduced to:

The rate of bubble volume growth is described as follows:

while the rate of change in mass of the bubble can be described as follows:

If the number of bubbles per unit volume is defined by N B , then the equation for the bubble vapor volume fraction α is

The total mass transfer rate per volume is then set as follows:

As the vaporization process occurs, the bubble volume fraction increases with increasing bubble diameter, and the nucleation density decreases. As a result, the volume fraction of water vapor produced during vaporization is a feature of r n ( 1 − α ) . The rate of vaporization and condensation bubble formation can be described as follows:

where r n represents the nucleation site of the vapor volume fraction, F v and F c reflect the empirical calibration coefficient of evaporation and condensation, respectively, and α represents the vapor volume fraction. The following model parameters have been shown to be effective for a wide range of fluids: R B = 10 – 6 m , r n = 5 × 10 – 4 , F v = 50 , F c = 0.01 . This model of cavitation flow is known as the Zwart mass transfer model or the Zwart cavitation model.

3.5 Kunz cavitation model

This mass transfer model parameter is set up in two ways to build m ̇ + and destroy m ̇ − of liquid. The liquid-to-vapor transition is measured as a function of how much pressure is below the vapor pressure. Otherwise, the conversion of vapor to liquid is dependent on a third-order polynomial function of volume fraction, C. The specific mass transfer rate is denoted by the formula m = (m ⁺ ) + (m ⁻).

where t ∞ = L/U denotes the mean flow time scale, U ∞ (m/s) represents the freestream velocity, L is the characteristic length scale, and C _dest and C _prod are two empirical coefficients, respectively, where C _dest = 100 and C _prod = 100 are assumed.

3.6 Thermodynamic parameter

The thermodynamic parameter is known as the change in water temperature and is used to calculate the behavior of bubble dynamics. Thermodynamic behavior ∑ can be denoted as [42]:

where L denotes the latent heat of vaporization, ρ represents the density, the specific heat at constant pressure is indicated by C p , and T is the temperature of the liquid.

4.1 The cavitation tests conducted in the Schiffbau-Versuchsanstalt (SVA) potsdam

In this work, the cavitation experiments were conducted using the cavitation tunnel K15A (Kempf and Remmers) of the SVA Potsdam [43]. A dimension with a length of 2593.5 mm and a cross-section of 850 × 850 mm was used for cavitation observation tests. Figure 6 depicts the cavitation tunnel in the potsdam cavitation test.

Figure 6

Cavitation tunnel in the potsdam cavitation test [43].

Throughout measurement, the propeller was mounted according to the configuration shown in Figure 6, with a 12° inclination of the propeller toward the inflow direction or the X-axis. In this case, the propeller in an inhomogeneous environment is powered by the cavitation tunnel’s homogeneous inflow. Tangential velocities work on one side of the propeller against and on the other side with the direction of blade rotation. Figure 7 shows the cavitation test section.

Figure 7

Cavitation test section [43].

For the simulation, this study has slight differences compared to the testing model. The simplification technique is applied, which only involves a propeller without a hub and shaft to achieve optimal element mesh due to resource limitation of computer capabilities. This method was valid because the hub behind the propeller would not disrupt the inflow. This technique was also used in simulation techniques by Morgut et al. [44], as shown in Figure 8.

Figure 8

Morgut et al. simulation model [44].

The inflow velocity in the numerical calculations can be determined with the advance coefficient. The pressure level in the numerical calculations can be resolved with the cavitation number. All evaluations shall be conducted with the propeller being in the zero-degree position. The origin of the geometry corresponds to the propeller plane position (x = 0 m). The calculations must be carried out according to the thrust identity (KT values are provided, non-cavitating), with the inflow velocity free to alter. The revolution rate is kept constant as it is used to define the cavitation number.

4.2 Propeller model and test case scenario

With the aim of testing and validation, SVA Potsdam developed a PPTC model VP1304 propeller model. The propeller consists of five blades with a controllable pitch propeller and a right-handed direction of rotation, as shown in Figure 9. The experimental data and geometries are available for viewing on the company website [43]. The VP1304 model was used in previous research for numerical simulation and validation results. Instead of experimental analysis, numerical simulation of the propeller output under cavitating flows using CFD-based software has been a good choice.

Figure 9

3D images of the Potsdam Test Case Propeller (PPTC1) VP1304 model.

The simulation began with the simulation of the cavitating flow for validation purposes under existing conditions and then continued with the observational scenario. Fluid properties of water and vapor are set according to the water condition in the experimental testing tank. The physical properties used in numerical simulations are assumed that water is freshwater and seawater. The saturation vapor pressure is a representation of the reference temperature. As shown in Table 2, this section consists of three operating cases.

Table 2 displays the operational parameter, i.e., the observation scenario as a case to be studied. The simulation uses J , σ n , and n as fixed parameters, while other parameters are temperature dependent and are the most important to observe. Water temperature varies at 1, 10, and 20°C. The current cavitating flow for the test case is subjected to fluid properties at 22, 21.9, and 23.2°C. New water material at specified water temperature is set to define water and vapor’s thermodynamic and transport properties. Furthermore, Table 3 offers a comprehensive description of the fluid’s properties. Using the CalcSteam online calculator, all fluids are measured based on the water temperature. In Table 3, all parameters with different temperatures are shown, with the new outlet boundary conditions as the static pressure choice determined using the equation of cavitation number.

4.3 Model validation

For validation purposes, the design of the computational domain used for the simulation test under consideration is a cylinder, as shown in Figure 10. The diameter D of the VP1304 propeller is 250 mm. The width and length of the rotating domain are 270 and 300 mm, respectively. The test section has an area of 850 mm × 850 mm perpendicular to the flow direction. The fixed domain is extended to 2593.5 mm overall. The simulations were run in 3D and assumed a steady-state condition. Because of the open water conditions, the propeller moved forward in a homogeneous uniform flow. The multiple reference frame (MRF) method was used to simulate propeller rotation, with the rotational speed remaining constant at n = 20 rps . The domain is divided into two parts: a rotating part and a fixed part. The MRF method analyzes the rotating domain in relation to other domains, while the general grid interface (GGI) method approximates a continuous surface [45].

Figure 10

Design of computational domain.

The single rotational speed of the water was applied to obtain the value of advance coefficient J and cavitation number σ n . The advance coefficient J was used to measure the inlet velocity. The turbulent pressure boundary on the inlet was set at 1%. The pressure outlet for each run was determined using the cavitation number σ n equation. The pressure outlet was P out = P ref with a static pressure option for all simulations. The free slip boundary condition was added to the domain’s outer boundary, while the propeller or solid surface was set with a no-slip boundary condition. The domain motion on the rotating domain was set to rotate with angular velocity. Furthermore, the stationary motion was applied. First, SST with automated wall treatment was used to influence the turbulent transition flow. In this work, the standard k−ε is applied in the next step for the full turbulence model with a scalable wall function. The Rayleigh–Plesset equation was used to describe cavitation bubble growth. Table 4 shows additional configuration parameters.

4.4 Meshing strategy

The chosen meshing method is an unstructured grid and hybrid mesh, which provides high precision. The rotating and fixed domain parts of the domain were segregated in this analysis. The mesh was created with CFX’s solver mesh (Figure 11). The rotating domain was initialized automatically, and the fixed domain was defined using the tetrahedron method. The CFX solver will join the rotating and fixed interfaces meshing method using a simplified grid interface (GGI). The mesh sizing with the body of influence selection to approach more precision for the critical segment was used to reduce the number of mesh elements and maximize the mesh density. Then, 11 layers with inflation on the propeller surface to obtain the viscous effect as the turbulent boundary layer were placed. For all propeller surfaces, all meshes have the same y + value of approximately 10. The following equation can be used to calculate the y + value [24]:

where ν denotes the kinematic viscosity of the fluid, τ w represents the shear stress at the wall, ρ l denotes the density of the water, and y represents the normal distance from the wall.

Figure 11

Mesh generated in Ansys program: (a) mesh interface for the rotating domain; and (b) mesh interface for a fixed domain.

4.5 Solution and solver setting

The CFD platform ANSYS CFX 19 was used for simulation, which offers solutions for three-dimensional viscous scenarios. The propeller cavitation simulation was performed under the assumption of steady-state conditions. The “upwind” advection scheme was added to the turbulent numerical high-resolution transport equation. The dynamic control model used the velocity–pressure coupling and multiphase control models.

4.6 Verification and validation tests

Verification and validation tests are conducted to investigate the effects of fluid properties and temperature effects on the cavitation process. Following that, the effects of propeller output at three operating conditions: ( J = 1.019 , σ n = 2.024 ) , ( J = 1.269 , σ n = 1.424 ) , and ( J = 1.408 , σ n = 2.000 ) are compared. The output of thrust and torque coefficient is calculated in this section using CFX results by comparing two turbulence models of k−ε and SST. Following that, in each case, all nondimensional coefficients of output ( K T  and  K Q ) are determined. The Zwart cavitation and Kunz Cavitation models’ validation results under various operating conditions are shown in Tables 5 and 6, respectively.

Furthermore, the validation stage was extended by performing a mesh dependence analysis to achieve optimum meshing. Table 7 displays the meshing information for the rotating and fixed domains depending on the number of elements. The optimum meshing value was calculated by taking simulation time and data error precision into account. It can be found that the longer the running time, the more the total of mesh elements. The result of K T value was determined using Case 3 at the existing temperature, and the error value was determined by comparing the simulation results to the experimental data collected in SVA potsdam’s cavitation tunnel [24].

There are inconsistencies in the results of simulation data in experiments of both cavitation models with variable errors in Tables 5 and 6, assuming steady state and heat transfer during the process are set with total energy. The comparison of numerical results shows that the predicted cavitation at the three operating conditions by both turbulence models is in good agreement with the corresponding experimental results. For the K T  and  K Q coefficients, the agreement between predicted results based on both numerical and experiments is good. The relative error value differs between k−ε and SST turbulence models, with the k−ε turbulence model providing an almost higher relative error, while SST produces more variable results. Based on the simulation results in Tables 5 and 6, the SST turbulence model is superior to the k−ε turbulence model for modeling cavitation. As a result, the simulation in the second stage was entirely based on SST turbulence models. Figures 12–15 provide a comparative graphic of K T  and  K Q values based on observations and simulations of various turbulence models.

Figure 12

Numerical results of thrust coefficient (K _T ) of Zwart Cavitation Model.

Figure 13

Numerical results of torque coefficient ( 10 K Q ) of the Zwart cavitation model.

Figure 14

Numerical results of thrust coefficient ( K T ) of the Kunz cavitation model.