Numerical Simulation of Cavitation Control around a Circular Cylinder Using Porous Surface by Volume Penalized Method

1. Introduction

The formation and subsequent collapse of vapor bubbles in a fluid can have detrimental effects on various engineering systems, including propellers, rudders, pumps, impeller and hydraulic turbines. Over the years, multiple numerical and experimental investigations have been conducted to unveil the mechanisms of cavitation around the cylinders, including the previous works by [1,2,3,4,5]. The phenomenon of cavitation can lead to significant damage and erosion on the solid boundaries, hydrodynamic efficiency reduction, and increasing of the noise and vibration on the immersible bodies [6,7]. To mitigate these undesirable effects, engineers have been exploring various methods to control and suppress the cavitation and its destructive effects. There have been several proposed methods for controlling the cavitation dynamics on the immersible bodies, which can be classified as either active or passive control method. Active control methods involve introducing additional force or momentum to the low-energy boundary layer, such as through the use of air injection or synthetic jets [8,9]. They employed synthetic jet actuators to effectively mitigate cavitation occurring on hydrofoils. Their results showed that the active control method can reduce the adverse effects of the cavitation during the operation of the hydrofoils.

By generation of controlled air or fluid jets through synthetic jet actuators, the cavitation dynamics and pressure pulsations around the hydrofoil can be manipulated at the cavitating regimes. Timoshevskiy et al. [10] found that tangential injection of liquid from the middle of a hydrofoil surface can effectively suppress cavitation. They presented that their technique is a promising way to suppress the partial unsteady cavity and cavitation cloud. In addition, the pressure fluctuations in the cloud cavitation regime was mitigated using the fluid injection. Despite the considerable advantages that active flow control methods, active methods need external energy sources and the implementation of air/water injection systems, adding complexity to the manufacturing process and limiting their applicability in long-distance navigation scenarios. It is worth noting that passive control methods for cavitation manipulation have attracted more attention due to their practicality and effectiveness compared to active approaches [11,12,13,14,15,16]. Kadivar et al. [17,18] studied the control of unsteady cloud and cavitation surge using different types of the vortex generators mounted on the surface of hydrofoils. They used wedge-type and hemispherical-type vortex generators to control the pressure pulsations due to the cavitation on the hydrofoils.

Modifying the surface properties, such as applying porous coatings to the body’s surface, has proven to be an effective and promising strategy for control of the cavitation. This approach involves altering the physical characteristics of the surface, introducing porous structures that can disrupt the formation and collapse of vapor bubbles in the cavitating flow [13]. The porous layer allows the fluid to penetrate into its structure, effectively reducing the pressure fluctuations and preventing the formation of the large-scale cavity structures. By introducing a porous layer, the cavitating flow is forced to pass through the interconnected pores, dissipating the energy of the flow and suppressing cavitation. Yu et al. [13] experimentally and numerically studied the dynamics of an unsteady cavitation around a hemisphere cylinder with and without porous layer. They concluded that the porous layer can suppress the cavitation inception and vortex structures around the cylinder. They indicated that different factors of the porous layer such as material properties, pore size, pore distribution, and porosity can influence the effectiveness of porous layers for control of the cavitation. Additionally, the thickness and surface roughness of the porous layer play a crucial role in determining the flow characteristics [14]. The ability to handle complexities in the geometry of porous layers, accurately and efficiently are the main issues in the computational fluid dynamics. To address these complexities, the Immersed Boundary Method (IBM) has emerged as an effective approach [19,20,21,22,23].

Within this methodological framework, the volume penalization technique, which is categorized as an immersed boundary method, plays a main role to simulate flow around complex geometries [21]. The volume penalization technique is particularly useful when dealing with flow problems that involve solid obstacles or porous media, where the grid generation or meshing can be challenging. In this technique, the flow domain is divided into two regions: the fluid region and the obstacle region. The flow equations, typically the Navier-Stokes equations, are solved over the entire domain, however a penalization term is introduced in the governing equations within the obstacle region [21]. This penalization term enforces the no-slip boundary condition and simulates the presence of the solid obstacle or porous medium. It essentially acts as a forcing term to account for the presence of the solid structure, ensuring that the flow behaves as if the porous medium is present [22].

Although numerous studies were performed to control cavitation on hydrofoils using passive flow control methods, numerical simulation of the porous layer potential as a cavitation control mechanism in a typical test case such as a circular cylinder has not received extensive attention in previous researches. Consequently, there is a significant gap in the existing literature regarding a comprehensive numerical investigation of this particular cavitation control strategy. For this research topic, fluid dynamics, materials science, and computational modeling are integrated to simulate and quantify the effect of porous surfaces on the cavitation control around blunt bodies. In the scope of our present research, one of the primary objective was to conduct a comprehensive investigation of the influence of the permeability of a porous layer with a certain thickness on control of the cavitation around a circular cylinder. To achieve this aim, we employed numerical simulation techniques, allowing us to explore and understand the dynamics of cavitation phenomenon around the circular cylinder with and without porous layers. The rest of the present work is organized as follows: the solution methodology based on the penalized Navier–Stokes equations with the homogeneous mixture model and numerical procedure are presented in the Section 2. Validation and results are discussed in the Section 3. Finally, the conclusions are given in the Section 4.

2. Methodology

In this research, the penalized Navier–Stokes equations for porous media was utilized, as proposed by [23], to simulate two-phase flow in the presence of the porosity. Furthermore, our study incorporates the homogeneous mixture model [24], which assumes that the phase mixture behaves as a single compressible fluid. Within this model, a linear combination of liquid and vapor densities was considered while assuming isothermal conditions and local kinematic equilibrium between the phases:

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_j\right)=-\left(\frac{1}{\varnothing }-1\right)\chi \frac{\partial }{\partial x_j}\left[\rho \left(u_j\right)\right]$$

(1)

$$\frac{\partial \mathit{\rho }{\mathit{u}}_{\mathit{i}}}{\partial \mathit{t}}+\frac{\partial }{\partial {\mathit{x}}_{\mathit{j}}}\left(\mathit{\rho }{{\mathit{u}}_{\mathit{i}}\mathit{u}}_{\mathit{j}}\right)=\frac{\partial \mathit{p}}{\partial {\mathit{x}}_{\mathit{i}}}+\frac{\partial {\mathit{\tau }}_{\mathit{i}\mathit{j}}}{\partial {\mathit{x}}_{\mathit{j}}}-\mathit{\chi }\frac{\mathit{\mu }}{\mathit{K}}\left({\mathit{u}}_{\mathit{i}}\right)-{\mathit{u}}_{\mathit{i}}\left(\frac{1}{\mathit{\varnothing }}-1\right)\mathit{\chi }\frac{\partial }{\partial {\mathit{x}}_{\mathit{j}}}\left[\mathit{\rho }\left({\mathit{u}}_{\mathit{j}}\right)\right]$$

(2)

$$\frac{\partial \mathit{\rho }\mathit{Y}}{\partial \mathit{t}}+\frac{\partial }{\partial {\mathit{x}}_{\mathit{j}}}\left(\mathit{Y}{\mathit{u}}_{\mathit{j}}\right)=0$$

(3)

the last term in the right side of the continuity equation (equation 1) and last two terms in the right side of the momentum equations (equation 2) are a volume penalization of the flow inside the body. Here, $\varnothing$ is the porosity and $\mu$ is the mixture viscosity. The permeability denoted as $K$, is determined by the structural of porous materials. The permeability signifies the material’s capacity to allow the flow of fluids or gases through its porous structure. Following the homogenous equilibrium model for compressible two-phase gas-liquid systems [24], the mixture density $\rho$ is expressed as a linear combination of gas density ${\rho }_g$ and liquid density ${\rho }_l$in the following manner:

$$\rho =\left(1-\alpha \right){\rho }_l+\alpha {\rho }_g$$

(4)

The relationship between the void fraction $\alpha$ and the gas mass fraction $Y$ can be described as follows:

$$1-Y=\left(1-\alpha \right){\rho }_l/\rho$$

(5)

$$Y=\alpha {\rho }_g/\rho$$

(6)

In the case of a pure liquid phase, the equation of state is described as follows [25]:

$$p+p_c={\rho }_{l}k\left(T+T_c\right)$$

(7)

For the vapor phase, the ideal gas equation of state is assumed:

$$p={\rho }_{g}RT$$

(8)

where $p_c=1944.61MPa$, $T_c=3837K$, $k=472.27J/kgK$ and $R=462J/kgK$ are the pressure, temperature, liquid and gas constants, respectively. The mixture viscosity is expressed by [25]:

$$\mu =\left(1-\alpha \right)\left(1+2.5\alpha \right){\mu }_l+\alpha {\mu }_g$$

(9)

In the penalization technique, the mask function $\chi$ is defined as follows:

$$\chi =\left\{\begin{array}{c}1 inside the rigid body\\ 0 otherwises\end{array}\right.$$

(10)

When $\chi =0$, the ordinary Navier–Stokes equations are regenerated. For the spatial derivatives, the fourth-order compact finite-difference method [26] is applied and to avoid the oscillatory solutions, the high-order low-pass spatial filter [27] has been applied to the solution.

Two-phase flows, marked by significant sound speed variations due to the varying compressibility of the liquid and the two-phase mixture, pose a challenge in simulation. Demonstrating a solution, a preconditioned technique is established. The governing equations (1)-(3) with preconditioning are rewritten in the vector form as:

$$\Gamma \frac{\partial W}{\partial \tau }+\frac{\partial Q}{\partial t}+\frac{\partial F_j}{\partial x_j}=S$$

(11)

where, t is the physical time, τ is a pseudo-time, $W$ and $\Gamma$ are the primitive variable vector and Jacobian matrix, respectively. $Q$ is the conservative vector, $F$ is the flux vector and $S$ is the penalized term. For the two-dimensional flows, the vectors and matrix are defined as:

$$W=\left\{\begin{array}{c}P\\ \mathrm{u}\\ \mathrm{v}\\ \mathrm{Y}\end{array}\right\}Q=\left\{\begin{array}{c}\mathsf{\rho }\\ \mathsf{\rho }\mathrm{u}\\ \mathsf{\rho }\mathrm{v}\\ \mathsf{\rho }\mathrm{Y}\end{array}\right\}$$

$$F=\left\{\begin{array}{c}\mathsf{\rho }{u}_j\\ \mathsf{\rho }\mathrm{u}{u}_j+p{\delta }_{ij}+{\tau }_{1j}\\ \mathsf{\rho }\mathrm{v}{u}_j+p{\delta }_{2j}+{\tau }_{2j}\\ \mathsf{\rho }{u}_j\mathrm{Y}\end{array}\right\},S=\left\{\begin{array}{c}-\left(\frac{1}{\varnothing }-1\right)\chi \frac{\partial }{\partial x_j}\left[\rho \left(u_j\right)\right]\\ -\chi \frac{\mu }{K}\left(u\right)-u\left(\frac{1}{\varnothing }-1\right)\chi \frac{\partial }{\partial x_j}\left[\rho \left(u_j\right)\right]\\ -\chi \frac{\mu }{K}\left(v\right)-v\left(\frac{1}{\varnothing }-1\right)\chi \frac{\partial }{\partial x_j}\left[\rho \left(u_j\right)\right]\\ 0\end{array}\right\}$$

(12)

$$\Gamma =\left[\begin{array}{c}{\rho }_p\rho 0{\rho }_Y\\ u{\rho }_p\rho 0{u\rho }_Y\\ v{\rho }_p\rho 0{v\rho }_Y\\ Y{\rho }_p\rho 0{Y\rho }_Y+\rho \end{array}\right]$$

Here, ${\rho }_p=\frac{\partial \rho }{\partial p}{|}_T$, which represents inverse of the square of the isothermal speed of sound. Also ${\rho }_Y=\frac{\partial \rho }{\partial Y}{|}_T={\rho }^2\left(\frac{1}{{\rho }_l}-\frac{1}{{\rho }_g}\right)$. To ensure reliable accuracy, ${\rho }_p$ in $\Gamma$ is replaced by ${\acute{\rho }}_p=\frac{1}{u_{\infty }^2}$ and preconditioning matrix ${\Gamma }_P$ is introduced as follows [28]:

$${\Gamma }_P=\left[\begin{array}{c}{\acute{\rho }}_p\rho 0{\rho }_Y\\ u{\acute{\rho }}_p\rho 0{u\rho }_Y\\ v{\acute{\rho }}_p\rho 0{v\rho }_Y\\ Y{\acute{\rho }}_p\rho 0{Y\rho }_Y+\rho \end{array}\right]$$

(13)

The second-order backward differencing for the physical-time derivative and explicit Euler for the pseudo-time derivative [29] is applied.

$${\Gamma }_P\frac{W^{m+1}-W^m}{∆\tau }+\frac{3Q^{n+1}-4Q^n+Q^{n-1}}{2∆t}={RHS}^m$$

(14)

$${RHS}^m=S\left(Q^m\right)-\frac{\partial F_j\left(Q^m\right)}{\partial x_j}$$

(15)

In this paper, n denotes the index for physical time, while m represents the index for pseudo-time (see [29] for more details).

3. Results and Discussion

In this section, the flow field over an isolated circular cylinder at Re=200 and cavitation numbers $\sigma =1$ was governed with porous layer thickness $h=0.3D$ and porosity φ=0.98. The cavitation number $\sigma$ and Reynolds number Re was given by $\sigma =\left(p_0-p_v\right)/\left(0.5{\rho }_0{u}_0^2\right)$ and $Re={\rho }_0{u}_{0}D/{\mu }_0$. A small void fraction of ${\alpha }_o$ = 0.01 was used to start the cavitation. Considering that the cell diameter of the porous material was in the micrometer orders, the permeability $K$ typically spans from ${10}^{-12}$ to ${10}^{-9}{m}^2$ [30]. Indeed, as $K$ decreases, the fluid permeability to the porous layer diminishes, and it assumes a more solid-like behavior. A rectangular grid system with non-uniform spacing was applied for the mesh type. The non-reflective boundary conditions presented by Tam and Dong [31] are imposed for all the free boundaries. To remove the possible reflections from the boundaries, the sponge layer with the grid stretching was utilized. The computational domain is presented schematically in Figure 1.

3.1. Validation

In order to examine the accuracy and efficiency of the methodology, we firstly considered single phase flow around the circular cylinder with Re=150. The Mach number of uniform flow was set to 0.2. The cylinder was wrapped by a porous layer with a thickness$h=0.4D$. In order to compare the results with the available data, we assumed $1/\eta =\frac{\mu }{K}$ in Eq. (2) is constant. Three values of $\eta ={10}^{-3},$10 and $\eta ={10}^3,$ with porosity φ=0.98 are studied for the corresponding simulation. For better understanding the effect of permeability on the aerodynamic results, the instantaneous spanwise vorticity fields of cylinders covered by the porous layer of different $\eta$ are shown in Figure 2. The outer solid line represents the porous layer. It can be observed, when the permeability is very low ($\eta ={10}^{-3}$), the porous materials act like solid, and the vortex shedding happens from the surface of the porous coating. In contrast, when the permeability is quite high $\eta ={10}^3$, the porous materials act like pure fluid and the vortex shedding occurs around the inner solid cylinder. At $\eta =10$, shear flow around the body is slightly deviated towards the outer layer, resulting in an attenuation of the vortex strength.

In Figure 3, the time-averaged pressure contours for the smooth cylinder ($\eta ={10}^3$), the porous cylinder ($\eta =10$), and the larger cylinder ( $\eta ={10}^{-3}$) are displayed. It is evident that the presence of the porous layer results in a weakening of the low-pressure zones. In the cylinder with extremely low permeability, the low-pressure zone shifts towards the outer layer of the porous material. Figure 4 displays the mean pressure coefficient $C_p$ around the cylinder at the porous layer. Despite the increased windward pressure at η=10 when compared to the smooth cylinder, the presence of porosity leads to a recovery in the pressure coefficient on the leeward side of the cylinder. Table 1 shows the drag coefficient $C_d$, which is match well with the previous works [32] at the single-phase flow around the circular cylinder.

3.2. Grid-Independence Study

To ensure the grid independence, three different mesh types were created with varying grid resolutions: coarse, medium, and fine meshes. These meshes were designed for the comparison purposes. The grid numbers for each mesh type in the context of a smooth cylinder case are outlined in Table 2. The grid size dependency is represented by the drag coefficient ($C_d$) and the Strouhal number (St), reflecting the sensitivity of the hydrodynamic characteristics to variations of the grid resolution. The comparison between the results obtained from the medium mesh and the fine mesh indicated only a minor discrepancy. Therefore, the medium mesh was selected for subsequent numerical simulations proved to be adequate. This decision was made to optimize computational resources and reduce costs while still maintaining an acceptable level of accuracy in capturing the essential flow dynamics.

3.3. Results

In this section, the hydrodynamics and the results of the cavitating flow over an isolated circular cylinder are presented. Unlike the previous section, $\frac{\mu }{K}$ is not constant here because the viscosity of the mixture varies within the field. Therefore, the permeability changes are considered for values of $K={10}^{-12},{10}^{-11},{10}^{-10}$. Figure 5 displays the time-averaged streamwise velocity along the centerline. There is a negative velocity region downstream of the cylinder, which is actually a distinct recirculation zone. For the smooth cylinder, the velocity reaches its minimum value around $X/D=0.9$. Beyond this point, the flow velocity gradually begins to recover and increase along the downstream range at $X/D=0.9$. As observed in streamwise velocity profile, the porous covering with low permeability exhibits a shift in the minimum flow velocity from the cylinder surface compared to the smooth cylinder. At $K={10}^{-11}$, the minimum velocity shifts to $X/D=1.03$, while at $K={10}^{-12}$, it shifts to $X/D=1.5$. Figure 6 shows the streamwise velocity profiles along the vertical centerline at different $X/D$ locations. Due to the blockage of the porous matrix, the streamwise velocity decreases. Once the porous media is applied around the cylinder, the wake region surrounding the cylinder widens. Specifically, the velocity within the wake, particularly near the central area (at $Y/D=0$), initially experiences a reduction followed by an enhancement compared to the velocity observed around a smooth cylinder. Based on the results presented in the Figure 5 and Figure 6, the porous layer with a permeability of $K={10}^{-10}$ does not significantly affect the mean velocity profiles. Therefore, two other cases are further investigated in the following.

The time histories of the lift and drag coefficients are shown in Figure 7 for the cylinder coated with different permeability of $K={10}^{-12},{10}^{-11}$ and with the thickness of the porous layer. Here, a reduction around 15% for $K={10}^{-11}$ and 39% for $K={10}^{-12}$ in the amplitudes of the lift fluctuation are evident across the cylinder coated with porous materials compared to the uncoated smooth cylinder. The drag coefficient is also depicted in Figure 7. In the figure, it can be deduced that the drag coefficient decreased around 2% with applying the porous layer at $K={10}^{-11}$. However, as permeability diminishes further, the drag coefficient begins to increase more than 13%. The reason for this behavior is the formation of blockage against the flow due to the reduction in the permeability. As the permeability of the porous layer decreases excessively, it behaves like a solid body, resulting in significant less flow passing through the porous layer. In addition, Table 3 shows the comparison of the lift coefficient and Strouhal number with various permeability. The results reveals that the porous layer leads to a reduction in the frequency of the dominant vortex behind the cylinders.

With the addition of a porous layer on the cylinder surface, the strength of the separated shear layers tends to decrease. In the smooth cylinder, vorticity primarily spreads from the surface into the surrounding external flow.

However, with the porous layer, vorticity diffusion occurs not only into the external flow field but also penetrates into the internal flow field within the porous layer itself. The right column of Figure 8 shows the mean vorticity contours around the smooth cylinder and porous-coated cylinders. The solid line encircling the cylinder represents the boundary between the porous material and the fluid. By applying the porous layer, the results reveal that the vortex range within the wake region of the cylinder becomes weaker and more elongated. A secondary shedding vortex develops along the solid wall of the cylinder at $K={10}^{-12}$, but it was swiftly dissipated by the effects of the porous layer on the leeward side.

Figure 9 provides evidence that the existence of the porous layer leads to a noticeable reduction in the strength of cavitation regions. The presence of the porous layer weakens and the attached cavitation zones on the cylinder’s wall, particularly with the low permeability $K={10}^{-12}$. The porous materials altered the local pressure distribution around the surface of the cylinder. They act as pressure manipulators, smoothing out the pressure gradients that might otherwise lead to the cavitation inception. This pressure redistribution is often linked to the ability of porous media to diffuse and absorb gas or vapor, reducing the low-pressure zones. Particular for the case $K={10}^{-12}$, no strong shock waves due to the collapse of cavitation structures behind the cylinder was observed. It can be deduced that the application of the porous layer plays an important role in the weakening or mitigating of the destructive effects of the unsteady cavitation such as strong pressure waves behind the circular cylinder. Therefore, the porous layer can be suggested as a passive control method to control the unsteady cavitation in different applications.