A robust and accurate finite element framework for cavitating flows with fluid-structure interaction
AA robust and accurate finite element framework for cavitating flows withfluid-structure interaction
Suraj R. Kashyap a , Rajeev K. Jaiman a, ∗ a Department of Mechanical Engineering, The University of British Columbia, Vancouver, BC V6T 1Z4
Abstract
In the current work, we present a unified variational mechanics framework for the cavitating turbulentflow and the structural motion via a stabilized finite element formulation. To model the finite mass trans-fer rate in cavitation phenomena, we employ the homogenous mixture-based approach via phenomeno-logical scalar transport differential equations given by the linear and nonlinear mass transfer functions.Stable linearizations of the finite mass transfer terms for the mass continuity equation and the reactionterm of the scalar transport equations are derived for the robust and accurate implementation. Thelinearized matrices for the cavitation equation are imparted a positivity-preserving property to addressnumerical oscillations arising from high-density gradients typical of two-phase cavitating flows. The pro-posed formulation is strongly coupled in a partitioned manner with an incompressible 3D Navier-Stokesfinite element solver, and the unsteady problem is advanced in time using a fully-implicit generalized- α time integration scheme. We first verify the implementation on the benchmark case of Rayleigh bub-ble collapse. We demonstrate the accuracy and convergence of the cavitation solver by comparing thenumerical solutions with the analytical solutions of the Rayleigh-Plesset equation for bubble dynamics.We find our solver to be robust for large time steps and the absence of spurious oscillations/spikes in thepressure field. The cavitating flow solver is coupled with a hybrid URANS-LES turbulence model witha turbulence viscosity corrected for the presence of vapor. We validate the coupled solver for a very highReynolds number turbulent cavitating flow over a NACA0012 hydrofoil section. Finally, the proposedmethod is solved in an Arbitrary Lagrangian-Eulerian framework to study turbulent cavitating flow overa pitching hydrofoil section and the coupled FSI results are explored for the characteristic features ofcavitating flows such as re-entrant jet and periodic cavity shedding. Keywords:
Cavitation, Fluid-Structure Interaction, Homogeneous-Mixture, Stabilized Finite Element,Partitioned iterative, Pitching Hydrofoil
1. Introduction
Cavitation is ubiquitous in natural and industrial industrial systems such as hydrofoils, nozzles,pumps, underwater vehicles and marine propellers. While cavitation can be a major source of noise,vibration and material erosion in these systems as unwanted effects [10, 34], useful applications of cav-itation have been developed in underwater cleaning [55, 12], ultrasonic cleaning [11, 60], biomedicalprocedures such as lithotripsy [5, 43] and enhanced drug delivery [56, 27], etc. The phenomenon ofcavitation involves the phase change of liquid into vapor and a highly complex interaction between thevapor and the liquid phases. Most liquids have inherent points of weaknesses in the form of entrainedmicroscopic gas bubbles, suspended particles, crevices along the shared boundaries with solid structures,and ephemeral voids created by the thermal motions of the liquid. When the flowing liquid encounters aregion of low pressure, it has a propensity to rupture at these locations of weaknesses due to the tensileforces [7]. This rupturing creates cavities in the fluid, which can be filled with the original entrained gasor by vapor generated by evaporation at the cavity interface. These cavities can then be convected bythe flow until they encounter a region of high pressure, where they can undergo violent collapse. Cavi-tation is often encountered in marine propeller operations where fluid acceleration over propeller blades ∗ Corresponding author
Email addresses: [email protected] (Suraj R. Kashyap), [email protected] (Rajeev K. Jaiman)
Preprint submitted to Computers & Mathematics with Applications February 22, 2021 a r X i v : . [ phy s i c s . f l u - dyn ] F e b enerates low-pressure regions near the blade surface. At these locations, nuclei present in seawater areprone to rupturing and forming vapor/gas-filled cavities. The cavities can then remain attached to or beshed from the blade surface in the form of cavitating structures with varying energy content [45]. Theonset of cavitation influences the fluid-structure dynamics of propeller operation and can have severalundesirable effects including performance degradation, vibration, material erosion and noise emission[7]. Traditional approaches to model the coupled fluid-structure dynamics in cavitating flows have oftenfocused on a one-way coupling between a flow solver for the cavitation hydrodynamics and a separatefinite element solver for structural deformation.The reduction of noise emission in marine vessels is of interest both from an industrial and a marine-environmental perspective. For example, [58] showed that in the cavitating regime, propeller noisedominates all other sources of self-noise from ships, including electrical noise, machinery noise andboundary layer noise. Recently, [10] provided an excellent review of noise from cavitating propellers andidentified two broad categories: (i) a broadband noise component resulting from the sudden collapse ofcavities and vortices, and (ii) tonal noise components from periodic fluctuations in the cavity volumes.In a classical work, [34] identified that tonal noise emission in propellers was a result of blade vibrationdue to irregular cavitation and vortex-shedding dynamics. Several experimental studies [10] have sinceconfirmed this. In addition, [3, 39, 4] showed that resonance in cavity-filled vortices shed from the bladetip can also emit intense tonal noise. Thus, a numerical study of propeller cavitation noise needs toconsider this complex hydrodynamic interplay between cavitation and vortex shedding, and their FSIeffects with the propeller blades. Figure 1 demonstrates some of these predominant noise-generatingmechanisms in hydrodynamic cavitation, with Ω f representing the fluid domain, Ω s the solid structuraldomain and Γ fs the fluid-structure interface.Ω f Ω s Γ fs collapsing cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge partial cavities & vortices U ∞ Figure 1: Representative schematic demonstrating some of the prominent noise sources from cavitating blades. Ω f and Ω s are the fluid and solid domains, Γ fs represents the interface between them. The emitted noise can be decomposed into twomain components: (i) tonal noise originating from blade vibration and resonance in tip vortex cavities, and (ii) broadbandnoise resulting from violently collapsing cavities and vortices. Cavitation in propellers manifests in the form of cavitating structures that exist across multiple ordersof spatial and temporal scales [4, 7], making the study of marine cavitation a challenging task. Thus,each of the different approaches developed for the numerical modeling of cavitating flows is generallycomputationally feasible only for a select set of flow configurations. A popular approach is to representthe fluid as a homogeneous mixture of liquid and vapor. The homogeneous mixture-based approachesdiffer in the estimation of the density field and can be broadly classified into two categories.The first category assumes equilibrium flow theory and the density is calculated using equations ofstate (EoS) [48, 52]. For isothermal flows, a barotropic equation of state is used to represent the densityfield as a function of the pressure [16, 13]. The equilibrium flow approach has the advantage of easierimplementation. It also does not require the use of empirical coefficients for modeling and relies on2ell-established equations of state. However, these models are generally used with the compressibleEuler equations and solved using density-based solvers. In [21], the authors compared a compressibledensity-based method and an incompressible pressure-correction method to study cavitating flows withthe barotropic EoS and reported better correlations with the compressible approach. Thus, this approachrequires very small time-steps to capture the pressure wave propagations in the compressible fluid [19].Another limitation of the equilibrium flow models based on barotropic EoS is that the gradients of thedensity and the pressure fields are parallel. Thus, the baroclinic torque which is proportional to the crossproduct ( ∇ ρ × ∇ p ) of these gradients is zero. This can result in inaccurate estimates of the vorticityproduction, which is an important feature of cavitating flows particularly in the closure region of attachedcavities [22]. Hence this approach may not be appropriate for modeling the physics of cavitating flowsover hydrofoils.The second category of homogeneous mixture-based approaches assumes the pure liquid and vaporphases to be incompressible. The two-phase mixture density is interpolated based on a phase indicator.This phase indicator is generally in the form of the local phase fraction of either the liquid or the vaporphase. The phase indicator in the computational domain is obtained as the solution of a scalar transportequation. We refer to these transport-equation based models as TEMs in the rest of the article. TheTEMs employ a source term, which is indicative of a finite mass transfer rate between the two phasesby the process of cavitation. TEMs generally vary in the formulation of the source term based onphenomenological arguments. Merkle et al. [41] used dimensional arguments for bubble clusters to relatethe source term to the local pressure and phase fraction of liquid. Kunz et al. [35] considered a similarsource term as [41] to model the evaporation process, but modified the condensation term employing asimplified form of the Ginzburg-Landau potential. Later models by Schnerr-Sauer [47], Zwart et al. [61]and Singhal et al. [54] assumed the cavities to be present in the form of clusters of spherical bubblesand used directly a simplified form of the Rayleigh-Plesset equation [7] for spherical bubble dynamics tomodel the mass transfer rate. These models differ in multiplier terms that were derived using differentphenomenological arguments for the underling bubble-bubble interaction. TEMs have been applied tothe study of several cavitating flow configurations, including hydrodynamic cavitation over hydrofoils[49, 20, 29].Despite wide applicability, one limitation of the TEM approach is the use of semi-empirical coefficientsin the source terms that need to be tuned for particular flow configurations. In addition, since thedensity is interpolated using the phase indicator, large spatial gradients in the density field exist acrossthe cavity interface. This can result in unphysical oscillations/spikes in the pressure in the vicinity of theinterface. For incompressible flow simulations, these pressure spikes can propagate rapidly throughoutthe computational domain, leading to numerical instability. These numerical artifacts have been reportedin several works such as [19], and special treatments are required for discontinuity capturing across theinterface. [50] suggested the ability to handle these spurious pressure spikes as one of the conditions fora robust solver for cavitating flows. The last decade has seen advances in the numerical study of FSI effects in non-cavitating marinepropeller operations. Simpler approaches have used methods based on inviscid potential flow theories.[37] and [40] used a coupled potential theory-based boundary element method (BEM) and finite elementmethod (FEM) approach for studying the hydro-elastic response of flexible marine propellers. [30]presented an optimization methodology for propellers considering FSI of the fluid and propeller bladesusing the panel method (PM), which is derived from BEM. A loosely coupled PM-FEM approach isused for the fluid-structure interaction between the hull wake and the propeller. While inviscid modelshave been demonstrated to be effective for making general design decisions for propellers because of lowcomputational cost, they are unable to capture the vortex-shedding process which is an important featureof cavitating flows. In addition, multiple modeling decisions have to be made on a case-by-case basisfor different propeller geometry and inflow conditions. Within the purview of viscous flow modeling,[36] used a tightly coupled CFD-FEM solver to study high Reynolds number flow over a flexible bladeundergoing vortex-induced vibration. The authors presented the requirement of tightly-coupled FSIsolvers for studying large-amplitude 3D vibrations at high Reynolds numbers.Although hydrodynamic studies of cavitating flows over propellers abound in literature, only lim-ited examples of studies that consider FSI have been found. [25] used a one-way coupled commercialcomputational fluid dynamics (CFD) solver to determine hydrodynamic loads on a cavitating hydrofoil3ndergoing prescribed pitching motion. [2] used a loose hybrid coupling to couple a commercial 2DURANS solver with a 2DOF hydrofoil model to study the effect of cavitation on the hydroelastic stabil-ity of hydrofoils. [59] used a similar approach to study the cavity shedding dynamics and flow-inducedvibration over a hydrofoil section. Reasonable agreements were reported with experimental measure-ments. However, there is a need for a robust and accurate unified framework for the strongly-coupledFSI studies of cavitating flows over propellers.The finite element method lends suitably to this purpose of modeling the caviating flows with fluid-structure interaction. It has long been staple for the study of structural deformations, and has alsobeen successfully applied to fluid flow studies [32, 33]. However, not much work has yet been doneon the modeling of cavitation using the finite element method. [6] presented a variational frameworkfor cavitating flows applied to prescribed motions of hydrokinetic turbines. A monolithic approachwas taken for the coupling of the cavitation TEM and flow equations. However, not much discussionhas been made on the cavity collapse pressures, shedding dynamics or the numerical stability at highliquid-vapor density ratios seen in marine flows. Recently, Joshi and Jaiman [32] presented a variationalfinite element framework to study the 3D hydroelastic response of marine riser in turbulent marineflows, undergoing high-amplitude vortex-induced vibrations. A strongly-coupled partitioned approachwas taken to solve the governing equations for the fluid flow, the structural deformation and the transportof the eddy viscosity. Further, Joshi and Jaiman [33] presented a positivity preserving variational (PPV)scheme for the numerical solution of two-phase flow of immiscible fluids. A discrete upwind operatorwas used locally in the interface region demonstrating oscillations. This acted in the form of addeddiffusion, ensuring positivity of the underlying element-level matrices. The scheme was demonstratedto be effective in reducing spurious pressure oscillations in fluids with high-density ratios. Numericalsolutions of cavitating flows using TEMs can benefit from similar treatment. The current work extendsthis framework by introducing the physics involved in the modeling of cavitating flows.
In the current work, we propose a novel finite element formulation for the numerical modeling ofcavitating flows. The objective is to integrate numerical modeling of cavitation into the framework ofthe stabilized finite element methods. Two cavitation TEMs based on homogeneous flow theory are usedto model finite mass transfer rates between the liquid and vapor phases by the process of cavitation. Thecavitation TEMs used are by Merkle et al. [41] and Schnerr-Sauer [47]. We present stable linearizationsof the two models for our proposed variational formulation. The linearized formulations are implementedin a variational finite element framework and the elemental matrices are imparted a positivity-preservingproperty for numerical stability in regions dominated by large density ratios. The fluid flow and cavitationsolvers are coupled in a staggered partitioned manner for versatility, and predictor-corrector iterationsare used for convergence stability and accuracy. A fully-implicit generalized- α time integration scheme[28] with user-controlled high-frequency damping is used to advance the solution in time, allowing fornumerical stability at relatively coarse spatial and temporal discretizations. Consistent with the sugges-tions in [50], the efficacy of the proposed method is demonstrated with respect to the requirements fora robust computational method for cavitating flows. The requirements can be summarized as (i) theaccurate prediction of the pressure field, (ii) the ability to handle large density ratios of the order of100-1000, and (iii) absence of spurious pressure spikes across the cavity interface.In the sections that follow, the governing equations and the proposed formulation are first presented.Next, a benchmark test of spherical vaporous bubble collapse is used to verify the numerical implementa-tion. The implementation is then validated on a turbulent cavitating flow configuration over a hydrofoilsection. Before concluding, the case of a pitching hydrofoil is considered to explore the compatibility ofthe formulation with FSI and its ability to predict flow features characteristic of cavitating flows. Finally,we close our paper with some concluding remarks.
2. Numerical formulation and implementation
In this section, we present a positivity preserving variational finite element implementation for thenumerical solution of TEMs applied to cavitation modeling. In particular, we consider two types ofTEM-based homogeneous mixture models classified as follows: • Model A: Cavitation model with nonlinear mass transfer rate by
Schnerr and Sauer [47]4
Model B: Cavitation model with linear mass transfer rate by
Merkle et al. [41]In the rest of manuscript, we shall refer to the models as A and B respectively. We design stablelinearizations of models A and B for the numerical implementation within the proposed variational finiteelement framework.
The strong forms of the governing equations are presented before introducing the variational formulation.We consider the fluid physical domain Ω f ( x f , t ) with an associated fluid boundary Γ f ( t ), where x f and t represent the spatial and temporal coordinates. The working fluid, consisting of the liquid and vaporphases, is assumed to be present in the form of a continuous homogeneous mixture. The phase indicator φ f ( x f , t ) is used to represent the phase fraction of the liquid phase at any coordinate ( x f , t ) in thehomogeneous two-phase liquid-vapor mixture. The fluid density ( ρ f ) and dynamic viscosity ( µ f ) aretaken as linear combinations of φ f ρ f = ρ l φ f + ρ v (cid:16) − φ f (cid:17) , (1) µ f = µ l φ f + µ v (cid:16) − φ f (cid:17) , (2)where ρ l and ρ v are the densities of the pure liquid and vapor phases, respectively. µ l and µ v are thedynamic viscosities of the liquid and the vapor phases. In the TEM approach, φ f is obtained as the solution of a scalar transport equation, which can be writtenin the conservative form in the ALE framework as: ∂ φ f ∂t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) χ + φ f ∇ · u f + (cid:0) u f − u m (cid:1) · ∇ φ f = ˙ mρ l , on ( x f , t ) ∈ Ω f (3)where χ is the referential coordinate system, u f = u f ( x f , t ) is the fluid velocity at each spatial point x f ∈ Ω f and u m is the relative velocity of the spatial coordinates x f with respect to the referentialcoordinate system χ . The source term ˙ m is representative of a finite mass transfer rate that governs therates of destruction and production of liquid by the process of cavitation. Cavitation TEMs vary in theway ˙ m is modeled. Cavitation Model A
For model A by [47], ˙ m is given as a non-linear function of φ f and p f ˙ m A = 3 ρ l ρ v ρ f R B (cid:115) ρ l | p f − p v | (cid:20) C c φ f (1 − φ f ) max (cid:0) p f − p v , (cid:1) + C v φ f (1 + φ nuc − φ f ) min (cid:0) p f − p v , (cid:1) (cid:21) (4)Model A attempts to relate the finite mass transfer rate to the rate of growth/collapse of an equiv-alent spherical bubble under an external pressure field, using a simplification of the Rayleigh-Plessetequation. Cavitation is assumed to initiate from nucleation sites present in the flow by a heterogeneousnucleation process [7]. The initial concentration of nuclei per unit volume ( n ) with an associated nucleidiameter( d nuc ) and is assumed to be a constant. It is also assumed that only vaporous cavitation occurs,and the effect of non-condensable gases is not considered. R B ( x f , t ) in Eq. (4) is representative of theequivalent radius of the vapor volume at the coordinates ( x f , t ), while φ nuc is the phase fraction of theinitial nucleation sites in an unit volume. These are calculated as R B = (cid:18) πn φ nuc − φ f φ f (cid:19) / and φ nuc = πn d nuc
61 + πn d nuc C c and the evaporationcoefficient C v , which require calibration for specific flow configurations. It has been applied to the studyof different cavitating flow configurations, including the collapse of vaporous bubbles [19] and cavitatingflow over hydrofoils [29].We see from Eq. (4) that ˙ m A depends on the local pressure p f as˙ m A ∝ p f − p v (cid:112) | p f − p v | p f ( x f , t ) ∈ IR − { } (6)It is observed that ˙ m A is not defined when p f = p v and can lead to numerical instability. In the currentwork, we take ˙ m A | p f = p v = lim p f → p v p f − p v (cid:112) | p f − p v | = 0 (7) Cavitation model B
For model B, ˙ m is a linear function of both φ f and the local pressure p f ˙ m B = C dest ρ l ρ v φ f12 ρ l U ∞ t ∞ min (cid:0) p f − p v , (cid:1) + C prod ρ l − φ f12 ρ l U ∞ t ∞ max (cid:0) p f − p v , (cid:1) (8)where p v is the saturation vapor pressure, U ∞ is the free-stream velocity and t ∞ is the mean flow time-scale. For hydrofoils, t ∞ is taken as t ∞ = C/U ∞ , where C is the chord length. C dest and C prod aresemi-empirical coefficients influencing the rates of destruction and production of liquid by the processof cavitation. In the current work, we assume C dest = 1 and C prod = 80 [51, 25] . Model B has beenderived using dimensional arguments for clusters of large bubbles and has been applied to macro-scalecavitation over hydrofoils [25, 6]. The unsteady Navier-Stokes equations for the fluid momentum and mass conservation can be written inan ALE framework as ρ f ∂ u f ∂t (cid:12)(cid:12)(cid:12)(cid:12) χ + ρ f (cid:0) u f − u m (cid:1) · ∇ u f − ∇ · σ = f f , on ( x f , t ) ∈ Ω f , (9) ∂ρ f ∂t (cid:12)(cid:12)(cid:12)(cid:12) χ + ρ f ∇ · u f + (cid:0) u f − u m (cid:1) · ∇ ρ f = 0 , on ( x f , t ) ∈ Ω f , (10)where f f is the body force applied on the fluid and σ = σ f + σ des (11)where σ f and σ des are the Cauchy stress tensor for a Newtonian fluid and the turbulent stress tensorrespectively, given by σ f = − p f I + µ f ( ∇ u f + ( ∇ u f ) T ) , (12) σ des = µ T ( ∇ u f + ( ∇ u f ) T ) , (13)where p f denotes the fluid pressure and µ T is the turbulent viscosity. σ des is modeled using the Boussinesqapproximation and in the current work a hybrid URANS-LES turbulence model is applied. The detailsof the turbulence model implementation can be found in Joshi and Jaiman [32].6 .1.3. Convective form of TEM and local fluid compressibility In the present work, the conservative form of the transport equation is re-arranged in the form aconvection-reaction equation. Taking the material derivative of Eq. (1) in the ALE framework, weobtain ∂ρ f ∂t (cid:12)(cid:12)(cid:12)(cid:12) χ + (cid:0) u f − u m (cid:1) · ∇ ρ f = ( ρ l − ρ v ) ∂ φ f ∂t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) χ + (cid:0) u f − u m (cid:1) · ∇ φ f (14)Combining equations (3), (10) and (14), the following forms of the mass continuity equation and thephase indicator transport equation are obtained, which are used in the current implementation. ∇ · u f = (cid:18) ρ l − ρ v (cid:19) ˙ m, on ( x f , t ) ∈ Ω f , (15) ∂ φ f ∂t (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) χ + (cid:0) u f − u m (cid:1) · ∇ φ f = ρ f ρ l ρ v ˙ m, on ( x f , t ) ∈ Ω f , (16)It is observed that the divergence of the velocity field is no longer zero, and local dilation effects areintroduced that are governed by the finite mass transfer rate. This local compressibility exists onlywithin the two-phase mixture and the pure phases are incompressible, since for the cavitation modelschosen no mass transfer occurs when φ f equals 0 or 1. The URANS implementations were developed for incompressible turbulent flows, and are not fully ca-pable of accounting for the local compressibility induced by the presence of vapor. For cavitating flowsover hydrofoils, this leads to an over-prediction of the turbulent stresses at the cavity closure and hasbeen reported by [15, 25]. This reduces the momentum of the characteristic re-entrant jet in cavitatingflows and does not allow it to penetrate the cavity. Thus, the periodic cavity breakup and shedding isundermined and the cavity remains attached to the surface in a quasi-steady state. This has also beenobserved by the authors in preliminary investigations, although not presented here. Hence, in the currentwork the turbulent viscosity is modified in the presence of vapor. An approach similar to that adoptedin [15] is taken. The turbulent viscosity µ T is modified as µ T mod = f ( ρ f ) ρ f µ T (17) f ( ρ f ) = ρ v + (cid:18) ρ v − ρ f ρ v − ρ l (cid:19) n ( ρ l − ρ v ) , n (cid:29) n = 10 is used in the current work, which has been shown to agree well with experimentalobservations of cavity shedding [15]. In Section 5, we study cavitating flow over a pitching hydrofoil. We briefly review the FSI boundaryconditions and the ALE mesh motion in the continuum setting. The modeling of FSI requires thesatisfaction of the velocity continuity and traction equilibrium at the fluid-structure boundary Γ fs . Letus consider a structural domain Ω s ⊂ R d with an associated structural boundary Γ s (0) at time t = 0.Let the mapping ϕ s ( x s , t ) map the deformation of the structure from its initial configuration Ω s to adeformed configuration Ω s ( t ) at time t , where x s denote the material coordinates. We denote the initialfluid-structure interface at t = 0 by Γ fs (0) = Γ f (0) ∩ Γ s (0). At time t the interface will then be deformedas Γ fs ( t ) = ϕ s (cid:0) Γ fs , t (cid:1) . The following kinematic and dynamic conditions are satisfied on Γ fs u f ( ϕ s ( x s , t ) , t ) = u s ( x s , t ) , ∀ x s ∈ Γ fs (19) (cid:90) ϕ s ( γ,t ) σ f · n f d Γ + (cid:90) γ σ s · n s d Γ = 0 , ∀ γ ⊂ Γ fs (20)where u s is the velocity of the structural domain, n f and n s are the unit normals to the deformed fluidelements ϕ s ( γ, t ) and their corresponding structural elements γ on the interface Γ fs respectively. Thestructural stress tensor σ s is modeled depending on the type of material.7he fluid spatial coordinates are updated to conform to the structural deformation. The motionof the coordinates which are not at Γ fs is modeled as an elastic material in equilibrium and the meshequation is solved as ∇ · σ m = , on Ω f , (21) η f = η f D , ∀ x f ∈ Γ m D (22)where σ m = (1 + k m ) (cid:104) ∇ η f + (cid:0) ∇ η f (cid:1) T + (cid:0) ∇ · η f (cid:1) I (cid:105) is the stress experienced at the fluid spatial coordi-nates due to the strain induced by the deformation of the interface, η f is the displacement of the fluidspatial coordinates. The amount of deformation of the spatial coordinates is controlled using the localstiffness parameter k m . Dirichlet conditions for the fluid mesh displacement η f D are satisfied on theboundary Γ m D . Both the fluid mass and momentum equations and the phase indicator transport equation are discretizedin time using a generalized- α predictor-corrector time integration method [14, 28]. For linear problems,the generalized- α method can be second-order accurate and unconditionally stable. This also enablesthe use of a single parameter called the spectral radius ρ ∞ , to achieve user-controlled high-frequencydamping. Let ∂ t φ f , n+ α m be the temporal derivative of φ f at time t n+ α m . Using the generalized- α method, φ f issolved for at the n + α time-level and integrated in time t ∈ (cid:2) t n , t n+1 (cid:3) according to the rules φ f , n+1 = φ f , n + ∆ t∂ t φ f , n + γ ∆ t ( ∂ t φ f , n+1 − ∂ t φ f , n ) ,∂ t φ f , n+ α m = ∂ t φ f , n + α m ( ∂ t φ f , n+1 − ∂ t φ f , n ) , (23) φ f , n+ α = φ f , n + α ( φ f , n+1 − φ f , n ) , where ∆ t is the time step size and ∂ t φ f , n+ α m is the temporal derivative of φ f at the n + α m time level. α m , α and γ are generalized- α parameters based on ρ ∞ [14].The transport equation is arranged in the following convection-reaction form for numerical solving G ( ∂ t φ f , n+ α m , φ f , n+ α ) = ∂ t φ f , n+ α m + (cid:0) u f − u m (cid:1) · ∇ φ f , n+ α + sφ f , n+ α − f = 0 , (24)where u f − u m is the convection velocity. The reaction coefficient s and the source term f for model Aare given by s A = − R B (cid:115) ρ l | p f − p v | (cid:20) C c (1 − φ f ) max (cid:0) p f − p v , (cid:1) + C v (1 + φ nuc − φ f ) min (cid:0) p f − p v , (cid:1) (cid:21) , (25) f A = 0 . (26)while for model B, we have s B = ρ f ρ v (cid:20) − ρ l ρ v C dest ρ l U ∞ t ∞ min ( p − p v ,
0) + C prod ρ l U ∞ t ∞ max ( p − p v , (cid:21) , (27) f B = ρ f ρ v (cid:20) C prod ρ l U ∞ t ∞ max ( p − p v , (cid:21) . (28) The Navier-Stokes equations are also discretized using the generalized- α time integration for consistencywith the phase indicator transport equation. The variational formulation employs the following relations8or time integration during t ∈ (cid:2) t n , t n+1 (cid:3) u f , n+1 = u f , n + ∆ t∂ t u f , n + γ ∆ t ( ∂ t u f , n+1 − ∂ t u f , n ) , (29) u f , n+ α = u f , n + α ( u f , n+1 − u f , n ) , (30) ∂ t u f , n+ α m = ∂ t u f , n + α m ( ∂ t u f , n+1 − ∂ t u f , n ) . (31)The semi-discrete forms of the Navier-Stokes equations are written as ρ f ∂ t u f , n+ α m (cid:12)(cid:12) χ + ρ f (cid:0) u f , n+ α − u m , n+ α (cid:1) · ∇ u f , n+ α − ∇ · σ n+ α − f n+ α = 0 , (32) ∇ · u f , n+ α − (cid:18) ρ l − ρ v (cid:19) ˙ m = 0 , (33) Next, we present the stabilized variational statements of the governing equations. We first present thepositivity preserving variational form of the cavitation model, followed by the Navier-Stokes equations.The stabilized variational form for the incompressible Navier-Stokes has been discussed in detail in otherwork [32, 31], which we shall not reproduce here. However, the presence of the non-zero divergence ofthe fluid velocity introduces additional terms to the formulation, which merits a brief discussion.
The fluid computational domain Ω f is spatially discretized into n el number of elements such thatΩ f = ∪ n el e=1 Ω f , e and ∅ = ∩ n el e=1 Ω f , e . The trial solutions are taken from from the space S h , which equal thegiven Dirichlet boundary condition at the boundary Γ D . The test functions are taken from the space V h ,which vanish on the Dirichlet boundary. We state the variational form of the phase indicator transportequation to find φ fh ( x f , t n+ α ) ∈ S h such that ∀ w h ∈ V h , (cid:90) Ω f (cid:18) w h ∂ t φ f , n+ α m h + w h (cid:0) u f , n+ α − u m , n+ α (cid:1) · ∇ φ f , n+ α h + w h sφ f , n+ α h − w h f (cid:19) dΩ f + n el (cid:88) e=1 (cid:90) Ω f , e (cid:18)(cid:0) (cid:0) u f , n+ α − u m (cid:1) · ∇ w h (cid:1) τ φ (cid:0) ∂ t φ f , n+ α m h + (cid:0) u f , n+ α − u m (cid:1) · ∇ φ f , n+ α h + sφ f , n+ α h − f (cid:1)(cid:19) dΩ f , e + n el (cid:88) e=1 (cid:90) Ω f , e χ | R φ ||∇ φ fh | k add s ∇ w h · (cid:18) (cid:0) u f , n+ α − u m (cid:1) ⊗ (cid:0) u f , n+ α − u m (cid:1) | ( u f , n+ α − u m ) | (cid:19) · ∇ φ f , n+ α h dΩ f , e (34)+ n el (cid:88) e=1 (cid:90) Ω f , e χ | R φ ||∇ φ f , n+ α h | k add c ∇ w h · (cid:18) I − (cid:0) u f , n+ α − u m (cid:1) ⊗ (cid:0) u f , n+ α − u m (cid:1) | ( u f , n+ α − u m ) | (cid:19) · ∇ φ f , n+ α h dΩ f , e = 0 , where the first line represents the standard Galerkin finite element terms. The second line consists oflinear stabilization terms with the stabilization parameter τ φ given by [53] τ φ = (cid:20)(cid:18) t (cid:19) + (cid:0) u f − u m (cid:1) · G (cid:0) u f − u m (cid:1) + s (cid:21) − / , (35)where G is the element contravariant metric tensor, which is defined as G = ∂ ξ T ∂ x f ∂ ξ ∂ x f , (36)where x f and ξ are the physical and parametric coordinates respectively. R φ is the residual of the phaseindicator transport equation given as R φ = ∂ t φ f , n+ α h + (cid:0) u f , n+ α − u m (cid:1) · ∇ φ f , n+ α h + sφ f , n+ α h − f. (37)The linear stabilization terms are used to address spurious oscillations in the solution when convectionand reaction effects are dominant. However, separate treatment is required to address oscillations in9he solution near the regions of high gradients [33]. Across the cavity interface, the transition from theliquid to the vapor phase takes place over the span of a few elements, and the spatial gradient of φ f isvery high. As seen in Eq. (1) and Eq. (2), physical properties of the fluid such as the density and theviscosity are obtained as weighted linear interpolations of φ f . Unbounded oscillations in φ f can result innegative values of ρ f and µ f , which are unphysical and can induce numerical instability. To address this,we introduce additional non-linear stabilization terms that impart a positivity property to the underlyingelement-level matrices. The third and fourth lines of Eq. (34) contain the positivity preserving nonlinearstabilization terms in the streamwise and crosswind directions respectively[33]. These essentially act asadded diffusion in the region of high gradients and ensure that the element-level matrix is an M -matrix.The PPV parameters χ , k add s and k add c for the phase indicator transport equation are obtained as: χ = 2 | s | h + 2 | ( u f − u m ) | , (38) k add s = max (cid:26) || (cid:0) u f − u m (cid:1) | − τ φ | (cid:0) u f − u m (cid:1) | s | h − τ φ | (cid:0) u f − u m (cid:1) | + sh , (cid:27) , (39) k add c = max (cid:26) | (cid:0) u f − u m (cid:1) | h sh , (cid:27) , (40)where | (cid:0) u f − u m (cid:1) | is the magnitude of the convection velocity and h is the characteristic element length[31]. We next present the Navier-Stokes equations and its variational form for the two-phase solver. The trial solution is taken from the function space S h , the values of which satisfy the Dirichlet boundarycondition at the Dirichlet boundary Γ D . Let V h be the space of test functions which vanish on Γ D .We formulate the variational statement for the fluid flow equations to find [ u n+ α h , p n+1h ] ∈ S h such that ∀ [ ψ h , q h ] ∈ V h , (cid:90) Ω f ψ f h · (cid:18) ρ f ∂ t u f , n+ α m h (cid:12)(cid:12)(cid:12) χ + ρ f (cid:16) u f , n+ αh − u m (cid:17) · ∇ u f , n+ αh (cid:19) d Ω+ (cid:90) Ω f ∇ ψ f h : σ n+ αh d Ω + (cid:90) Ω f ( t n+1 ) q h (cid:16) ∇ · u f , n+ αh (cid:17) d Ω+ n el (cid:88) e =1 (cid:90) Ω f τ m ρ f (cid:16) ρ f (cid:16) u f , n+ αh − u m (cid:17) · ∇ ψ f h + ∇ q h (cid:17) · R m d Ω e + n el (cid:88) e =1 (cid:90) Ω e ∇ · ψ f h τ c ρ f R c d Ω e (41) − n el (cid:88) e=1 (cid:90) Ω e τ m ψ h · ( R m · ∇ u n+ α h )dΩ e − n el (cid:88) e=1 (cid:90) Ω e ∇ ψ h ρ ( φ ) : ( τ m R m ⊗ τ m R m )dΩ e = (cid:90) Ω f ψ f h · f f , n+ α d Ω + (cid:90) Γ f N ψ f h · h f , n+ α d Γ + (cid:90) Ω f q h (cid:18) ρ l − ρ v (cid:19) ˙ m f , n+ α d Ω (cid:124) (cid:123)(cid:122) (cid:125) where the first and second lines contain the standard Galerkin finite element terms of the momentumand the continuity equation. The third line contains the Galerkin Least Squares stabilization terms forthe momentum and mass continuity equations. The fourth line contains stabilization terms based onthe multi-scale argument [26, 24]. The fifth line contains the Galerkin terms for the body force and theNeumann boundary in the momentum equation. The term in under-braces corresponds to the Galerkinprojection of the term in the mass conservation equation dependent on p f , and has been introduced inthis formulation. The element-wise residual of the momentum and the continuity equations denoted by R m and R c respectively are given by R m ( u f , p f ) = ρ f ∂ t u n+ α m h + ρ f (cid:0) u f , n+ α − u m (cid:1) · ∇ u f , n+ α h − ∇ · σ f , n+ α h − f f , n+ α h , (42) R c ( u f , p f , φ f ) = ∇ · u f , n+ α h − (cid:18) ρ l − ρ v (cid:19) ˙ m f , n+ α . (43)10 m and τ c are stabilization parameters [8, 57, 18] defined as τ m = (cid:20)(cid:18) t (cid:19) + u h · Gu h + C I (cid:18) µ ( φ ) ρ ( φ ) (cid:19) G : G (cid:21) − / , (44) τ c = 1tr( G ) τ m , (45)where C I is a constant derived from the element-wise inverse estimate [23] and tr( G ) denotes the traceof G . The Navier-Stokes equations are coupled with the phase fraction transport equation in a partitionedmanner as shown in Eq. (46) and Eq. (47). The equations are linearized with a Newton-Raphsontechnique to solve for the incremental velocity, pressure and phase indicator which are advanced in timeusing the Generalized- α time integration method. The linearized matrices for the fluid flow and thephase indicator transport equations are arranged in the form (cid:20) K Ω f G Ω f − G T Ω f C Ω f (cid:21) (cid:26) ∆ u f ,n + α ∆ p f ,n +1 (cid:27) = (cid:26) R m R c (cid:27) (46)[ K φ ] { ∆ φ f ,n + α } = { R φ } (47)where K Ω f is the stiffness matrix of the momentum equation, G Ω f is the discrete gradient operator and G T Ω f is the divergence operator. C Ω f consists of the pressure-pressure stabilization term and the termsin the mass continuity equation having dependency on the pressure. K φ is the stiffness matrix for thephase indicator transport equation, consisting of the transient, convection, reaction, linear stabilizationterms and the non-linear PPV stabilization terms. ∆ u f ,n + α , ∆ p f ,n +1 and ∆ φ f ,n + α are the incrementsin velocity, pressure and the phase indicator respectively.Linearization of the phase indicator transport equation requires the derivative of the terms in Eq. (2.2.1)with respect to φ f ,n + α . ∂G∂φ f ,n + α = ∂ (cid:0) ∂ t φ f , n+ α m (cid:1) ∂φ f ,n + α + (cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:24)(cid:58) ∂ (cid:0)(cid:0) u f − u m (cid:1) · ∇ φ f , n+ α (cid:1) ∂φ f ,n + α + ∂ (cid:0) sφ f , n+ α (cid:1) ∂φ f ,n + α − (cid:26)(cid:26)(cid:26)(cid:26)(cid:26)(cid:62) ∂f∂φ f ,n + α (48)= α m αγ ∆ t + ∂ (cid:0) sφ f , n+ α (cid:1) ∂φ f ,n + α (49)where the derivative of the transient term reduces to a constant obtained from the generalized- α relationsin Eq. (23). Owing to the symmetrical property of second order cross-derivatives, the convective termalso reduces to zero. The derivative of the reaction term ∂ (cid:0) sφ f , n+ α (cid:1) ∂φ f , n+ α needs to be calculated. Similarly,in the linearization of the mass continuity equation we encounter the term ∂ ˙ m f , n+ α ∂p f , n+1 . Care must be takenin forming the terms ∂ ˙ m∂p f , n+1 and ∂ (cid:0) sφ f , n+ α (cid:1) ∂φ f , n+ α for a stable linearization of the matrices C Ω f and K φ respectively. In model A, ˙ m varies non-linearly with both p and φ f , along with a discontinuity at p = p v .Model B presents a simpler linearization with ˙ m being a linear function of both p and φ f . We proposelinearlizations for the said terms in Tables. (1) and (2) respectively. The presented linearlizations havebeen tested to be stable over a range of temporal and spatial discretizations on two different cavitatingflow configurations. Algorithm 1 shows the algorithm for the staggered partitioned coupling of the implicit Navier-Stokesand cavitation solvers. This provides two unique advantages compared to monolithic couplings - (i)simplicity of linearization and implementation, and (ii) the flexibility to couple additional governing11 odel A Model B p n+1 h > p v C c φ f , n+ αh (1 − φ f , n+ αh ) ρ l ρ v ρ f R B (cid:115) ρ l (cid:0) p n+1 h − p v (cid:1) C prod ρ l − φ f , n+ αh ρ l U ∞ t ∞ p n+1 h < p v C v φ f , n+ αh (1 + φ nuc − φ f , n+ αh ) ρ l ρ v ρ f R B (cid:115) ρ l (cid:0) p v − p n+1 h (cid:1) C dest ρ l ρ v φ f , n+ αh ρ l U ∞ t ∞ p n+1 h = p v
0, using lim p n+1 h → p v p n+1 h − p v (cid:113)(cid:12)(cid:12) p n+1 h − p v (cid:12)(cid:12) = 0 0 Table 1: ∂ ˙ m f , n+ α ∂p f , n+1 in linearized matrix C Ω f Model A Model B p n+1 h > p v − C c (1 − φ f , n+ αh ) 3 R B (cid:115) (cid:0) p n+1 h − p v (cid:1) ρ l ρ f ρ v C prod ρ l U ∞ t ∞ (cid:0) p n+1 h − p v (cid:1) p n+1 h < p v C v (1 + φ nuc − φ f , n+ αh ) 3 R B (cid:115) (cid:0) p v − p n+1 h (cid:1) ρ l − ρ f ρ l ρ v C dest ρ l U ∞ t ∞ (cid:0) p n+1 h − p v (cid:1) p n+1 h = p v
0, using lim p n+1 h → p v p n+1 h − p v (cid:113)(cid:12)(cid:12) p n+1 h − p v (cid:12)(cid:12) = 0 0 Table 2: ∂ (cid:0) sφ f , n+ α (cid:1) ∂φ f , n+ α in linearized matrix K φ equations as demanded by the case being studied. For example, in Sections 4 and 5, the additionalequations for turbulence modeling and the ALE mesh update are similarly coupled using a staggeredpartitioned approach. For stable and accurate convergence, non-linear predictor-multicorrector iterationsare performed within each time step. Let us consider the time level t n and the associated velocity u f ( x , t n ), pressure p f ( x , t n ) and phase indicator φ f ( x , t n ) fields. Each non-linear predictor-correctoriteration consists of one pass through the steps [A] , [B] , [C] and [D] . In step [A] of the predictor-corrector iteration k, a predictor fluid velocity u f , n+1(k+1) and pressure p f , n+1(k+1) is obtained from the solutionof the Navier-Stokes equations(Eq. (46)). These are passed to the cavitation solver in step [B] . Anupdated phase indicator φ f , n+1(k+1) is obtained in step [C] by solving Eq. (47), which is then passed tothe Navier-Stokes solver in step [D] . This updated phase indicator value is used to interpolate the fluiddensity and dynamic viscosity, and prepare the Navier-Stokes matrices for the next iteration k + 1. Thiscyclic process is continued till the solvers have achieved the convergence criteria. Then the coupled solveris advanced to the next time level t n+1 . 12 lgorithm 1 Partitioned coupling of implicit Navier-Stokes and cavitation solvers
Input: u f , , p f , , φ f , for n ← n last do Predict solution (cid:20) u f , n+1(0) p f , n+1(0) φ f , n+1(0) (cid:21) ← (cid:20) u f , n p f , n φ f , n (cid:21) Interpolate density and viscosity fields ρ f ( φ f , n+1(0) ), µ f ( φ f , n+1(0) ) for k ← k max do[A] Navier-Stokes Solver [C] Cavitation SolverNon-linear[D] u n+1(k+1) , p n+1(k+1)
1. Interpolate solution2. Solve3. Correct solution4. Update solution 1. Interpolate solution2. Solve3. Correct solution4. Update solution u f , n+ α (k+1) ← u f , n + α ( u f , n+1(k) − u f , n ) p f , n+1(k+1) ← p f , n+1(k) ∆ u f , n+ α and ∆ p f , n+1 in Eq. (46) u f , n+ α (k+1) ← u f , n+ α (k+1) + ∆ u f , n+ α p f , n+1(k+1) ← p f , n+1(k+1) + ∆ p f , n+1 u f , n+1(k+1) ← u f , n + 1 α ( u f , n+ α (k+1) − u f , n ) p f , n+1(k+1) ← p f , n+1(k+1) + ∆ p f , n+1 φ f , n+ α (k+1) ← φ f , n + α ( φ f , n+1(k) − φ f , n )∆ φ f , n+ α in Eq. (47) φ f , n+ α (k+1) ← φ f , n+ α (k+1) + ∆ φ f , n+ α φ f , n+1(k+1) ← φ f , n + 1 α ( φ f , n+ α (k+1) − φ f , n ) iterations φ n+1(k+1) [B]end doend do A similar partitioned iterative coupling is used between the ALE fluid-cavitation and the structuraldisplacement. In the current study, the structural displacement is dictated by the prescribed motion.Alternatively, for fully-coupled FSI studies the structural displacement can be obtained by solving thegoverning equations for the structural mechanics as presented in Joshi and Jaiman [32]). At time t n ,let us consider a predictor structural displacement η s ( x s , t n ) obtained from the structural update inthe Lagrangian reference frame. These structural displacements are then passed to the fluid solvers,while satisfying kinematic conditions at the fluid-structure interface Γ fs as follows. At the time t n+1 , thedisplacement η f of the fluid spatial coordinates (i,e., Eulerian mesh coordinates) at the wetted boundaryΓ fs are equated to the structural displacement η s . η f , n+1 = η s , on Γ fs (50)In addition, the velocity continuity is satisfied at the interface is satisfied at the time t n+ α as u f , n+ α = u m , n+ α , on Γ fs (51)13ith the mesh velocity at the interface evaluated as u m , n+ α = η f , n+1 − η f , n ∆ t , on Γ fs (52)Away from the interface Γ fs and any Dirichlet conditions on η f on the Dirichlet boundary Γ m D , Eq. (21)is solved for updating the fluid spatial coordinates. The flow-cavitation equations are then solved withthe updated kinematic boundary conditions and the displaced fluid spatial coordinates. The updatedhydrodynamic forces are passed to the structural solver to correct the deformations. This sequence isrepeated iteratively till converged solutions are obtained.For the finite element discretization of the variables u f , p f and φ f , we consider equal-order interpo-lations. The Harwell-Boeing sparse matrix format is used to form and store the matrices for the linearsystem of equations. A Generalized Minimal RESidual (GMRES)[46] algorithm is used to solve the linearsystem. We observe that 2-3 non-linear iterations(as described in Algorithm (1)) are sufficient to obtaina converged solution at each time-step. The solver uses communication protocols based on standardmessage passing interface [1] for parallel computing on distributed memory clusters.
3. Numerical verification and convergence
A stabilized numerical method for the solution of cavitating flows has been presented. In this section weverify the implementation and discuss its convergence and efficacy.
The numerical implementation is verified by comparison with the analytical solution of the Rayleigh-Plesset equation for spherical bubble dynamics. It has influenced several works on the study of cavitationand bubble dynamics, and has seen many adaptations over the last century. We present here only salientfeatures relevant to the current study and interested readers are directed to comprehensive texts such as[7, 17].We consider the case of an iso-thermal collapse of a spherical vaporous bubble. A bubble of radius R is initialized in an infinite domain of liquid with a constant far-field pressure p ∞ . It is assumed thatthe bubble contents are homogeneous and consist only of saturated vapor and no non-condensable gases.It is also assumed that the pure liquid is incompressible. In the absence of any non-condensable gases,the uniform pressure p B ( t ) inside the bubble equals p v . Fig. (2) shows the schematic of the domain.Spherical symmetry is assumed and a one-dimensional analysis is performed along the radial direction. r R ( t ) p ( r, t ) u ( r, t ) p B ( t ) = p v Liquid Liquid − vapor interf aceV apor r ∞ , p ∞ Figure 2: Representative schematic for iso-thermal collapse of a vaporous spherical bubble
Using mass conservation in the incompressible liquid outside the bubble, the radial outward velocity u ( r, t ) at a radial distance r from the center of the bubble can be shown to follow an inverse square lawof the form u ( r, t ) ∝ /r (53)14ssuming that the liquid at the interface is in thermal equilibrium with the vapor at the saturationtemperature, no mass transfer in the form of evaporation or condensation occurs at the interface. Thus,at the interface the velocity u ( R, t ) = dR/dt . Thus Eq. (53) can be written as u ( r, t ) = R r dRdt (54)where R ( t ) is the bubble radius. Momentum conservation in the r -direction is given by − ρ l ∂p∂r = ∂u∂t + u ∂u∂r − ν l (cid:20) r ∂∂r (cid:18) r ∂u∂r (cid:19) − ur (cid:21) (55)Substituting the expression for u from Eq. (54) the following relation is obtained − ρ l ∂p∂r = 2 (cid:18) dRdt (cid:19) (cid:20) Rr − R r (cid:21) + R r d Rdt (56)Assuming surface tension and viscous effects to be negligible, and the absence of any mass transfer atthe interface, p r = R = p B ( t ) = p v . Using this assumption, spatial integration of Eq. (56) between r = R and r → ∞ , gives the well-known Rayleigh equation [44] p v − p ∞ ρ l = R d Rdt + 32 (cid:18) dRdt (cid:19) (57)Eq. (57) can be integrated to obtain the interface velocity during the collapse dRdt = − (cid:115) p v − p ∞ ρ l (cid:18) − R R (cid:19) (58)The interface acceleration is obtained as d Rdt = p v − p ∞ ρ l R R (59)Using the kinematic relations in Eq. (58) and Eq. (59), Eq. (56) can be integrated between r (in theliquid outside the bubble) and r → ∞ to obtain the pressure p ( r, t ) at the radial distance r from thecentre of the bubble p ( r, t ) − p ∞ p ∞ − p v = − (cid:18) − R R (cid:19) (cid:18) Rr − R r (cid:19) − R R r (60)In the absence of thermal effects and non-condensable gas content, we encounter the special case of avaporous bubble collapsing to zero volume. The total time t tc taken by such a bubble to collapse froman initial radius R is presented by Rayleigh’s relation [44]: t tc = 0 . R (cid:18) ρ l p ∞ − p v (cid:19) / (61) As we discussed in Section 1, the numerical modeling of cavitation requires careful consideration ofthe cavitation model based on the spatial-temporal scales of the particular flow configuration. We usecavitation model A for the numerical study of this case of micro-scale bubble collapse because of itsorigins in the Rayleigh-Plesset equation. In addition, model A has been previously used to study thisconfiguration in [19], where an implementation based on the finite volume method was used. For com-parison of the efficacy of the present variational finite element implementation, we set up our numericalcase using similar geometrical parameters and model coefficients as in [19].For the numerical study we consider a 3D spherical domain of R ∞ = 0 . R = 4 × − m is initialized in the centre of the domain.The spherical domain is discretized using 840264 hexahedral elements, with 49 nodes in the radialdirection resolving the initial bubble. The initial phase indicator and pressure inside the bubble are set15o 0 .
01 and p v = 2320Pa respectively, corresponding to the vapor phase. Outside the bubble, the phaseindicator is set to 1. The pressure in the liquid is initialized according to Eq. (60). Fig. (3) shows theinitial conditions for the phase indicator. A far-field pressure P ∞ = 1 × Pa is weakly enforced overthe outer surface boundary boundary of the domain using a traction boundary condition. The densitiesof the pure liquid and vapor phases are taken as ρ l = 1000 kg m − and ρ v = 0 . − , giving adensity ratio ρ l ρ v ≈ µ l = 0 .
001 kg m − s − and µ v = 10 − kg m − s − . The model parameters n and d nuc are assumed to be 10 and 10 − mrespectively. The condensation and evaporation coefficients are set to C c = C v = 100. Figure 3: Sliced octant of domain showing the initial phase indicator values. Inlay showing a close-up of the bubble at thecentre of the domain. Inside the bubble, the pressure is initialized to the vapor pressure p v , with φ f = 0 .
01 correspondingto the vapor phase. Outside the bubble, the initial pressure has a continuous distribution according to Eq. (60) and aconstant φ f = 1 corresponding to the liquid phase. The results of the numerical study are first compared to the exact solution of Eq. (58). Fig. (4a) showsthe evolution of the bubble radius R ( t ) with time. In the numerical study, we consider the bubble radiusto be the effective radius of the volume of vapor in the domain. R = (cid:18) π (cid:90) Ω f (cid:0) − φ f (cid:1) dΩ f (cid:19) / (62)Since no evaporation occurs outside the bubble, and no convection of vapor occurs though the boundaryat R ∞ (cid:18) ∂φ f ∂r = 0 (cid:19) , the only vapor content in the domain is in the bubble. This is also confirmed laterin Fig. (6a). Thus, Eq. (62) is seen to represent the bubble radius well. The bubble radius and thesolution time are non-dimensionalized using the initial bubble radius R and the Rayleigh collapse time t tc respectively. We consider predictions obtained using four different values of the numerical time-step∆ t . It can be observed from Fig. (4a) that solutions obtained using ∆ t ≤ × − s are able to capturethe evolution of the bubble radius well. To proceed, we sample the error between the numerical and exactbubble radius at four instances during the collapse process. These sampling instances correspond to thephysical times 0 . t tc , 0 . t tc , 0 . t tc and 0 . t tc . We compute the Root Mean Square Percentage16rror (RMSPE) using the sampled results asRMSPE = (cid:118)(cid:117)(cid:117)(cid:116) n · n (cid:88) i =1 ∆ R ,i · R rel ,i = R num,i R exact,i − R num,i is the numerically obtained bubble radius and R exact,i is the exact solution and i denotesthe index of the sampling time. Figure (4b) shows the calculated RMSPE for different values of ∆ t . Weobserve that the RMSPE is below 2 .
5% for ∆ t ≤ × − s, indicating good agreement with the exactsolution. Thus, for the rest of this study we present results obtained with ∆ t = 5 × − s. Remark . The solution is observed to be convergent and stable across a range of time-step sizes, evenat large values of the order of ∆ t = 1 × − s. No numerical oscillations are seen in the solution fieldsin the domain. [19] reported the presence of spurious pressure pulses in the domain for ∆ t > × − s.We note the efficacy of the present implementation in suppressing these spurious oscillations.We next look at the predicted pressure field, which is another quantity of interest in cavitating flows.It is well known that the cavity collapse process can result in high pressures in the domain, several timesthe magnitude of the collapse driving pressure p ∞ . Thus, for any effective study of the fluid-structureinteraction and noise effects in cavitating flows, it is important that the cavitation solver is able toaccurately capture the pressure field. Figure (5a) shows the non-dimensional pressure Π = p f − p ∞ p ∞ − p v in the domain at different times during the bubble collapse, compared to the analytical solution ofEq. (60). The pressure p f is taken in the radial direction along the z -coordinate axis, with r = 0 at thecenter of the bubble. The numerically obtained pressures are observed to be in good agreement withthe analytical solution. The peak pressures are close to the expected values, and an improvement inaccuracy is obtained compared to results presented in [19]. No spurious spikes or numerical oscillationsare observed in the pressure.However, it is observed that at the final stages of the collapse process, the pressure inside the bubble ishigher than p v . This deviates from the assumption that the bubble pressure p B ( t ) equals p v at all times.Increasing the coefficients C c and C v appears to resolve this. Fig. (5b) compares the pressure profile at t = 0 . τ R for three values of the evaporation and condensation coefficients. At C c = C v = 175, thepressure profile nearly matches the analytical solution. The pressure inside the bubble is also observedto approach closer to the theoretical value of the vapor pressure. On further increasing the C c = C v to 250, the vapor pressure inside the bubble is recovered. However, the peak pressure is observed to beunder-predicted. Thus, a systematic tuning of the semi-empirical coefficients is required, which is one ofthe limitations of such phenomenological models. However, further tuning of these parameters is beyondthe scope of the current work.Unlike in [19] where spherical symmetry was assumed in the numerical simulation, in the currentstudy we consider the full 3D spherical domain. Numerical discretization and the weakly enforced far-field pressure boundary condition can introduce asymmetry in the solution field. To determine theability of the solver to preserve the symmetrical nature of the solution, the non-dimensional pressure Πis plotted along the three cartesian directions. Fig. (5c) shows the comparison of Π (at t = 0 . t tc , C c = C v = 100) in the x , y and z directions as a function of the distance from the centre of the bubble.It is observed that the 3D solver is able to naturally preserve the symmetry of the collapsing bubble. Remark . We note the ability of the present implementation to accurately predict the pressure in thedomain, and the absence of spurious pressure oscillations across the bubble interface. The vapor pressureinside the cavity is also recovered, although it requires careful calibration of the model coefficients. It ispossible that an adaptive calibration (e.g., deep learning models [42, 9])) of these coefficients can aid ingeneralizing the model to multiple flow configurations. This can be explored in future studies.Next, we investigate the predicted values of the phase indicator φ f in the domain. As discussed previously,presence of large spatial gradients of φ f across the bubble interface can result in unbounded numericaloscillations in the density field, and in turn the pressure field. Fig. (6a) shows the distribution of φ f along17 a)(b) Figure 4: Numerical assessment of spherical bubble collapse problem: (a) comparison of bubble radius against analyticalsolution of the Rayleigh equation (b) root mean square percentage error in bubble radius calculation at different values ofthe numerical time-step size ∆ t a)(b) (c) Figure 5: Spherical bubble collapse problem: (a) comparison of non-dimensional pressure Π with the analytical solutionat different time steps, (b) sensitivity of non-dimensional pressure Π on the coefficients C c and C v , and (c) assessment ofspherical symmetry. the radial direction, while Fig. (6b) shows the gradient (cid:18) ∂φ f ∂r (cid:19) of φ f with respect to the radial distance r . We observe no oscillations in the solution in the vicinity of the bubble surface. The solution is seento be bounded within the range φ f ∈ [0 , ∂φ f ∂r is ≥ C c = C v = 175 is observed to best representthe expected distribution of φ f in the domain. Remark . We note the ability of the present implementation to recover monotone and bounded solutionsacross the bubble interface. The absence of spurious oscillations in φ f in the vicinity of the interface19 a)(b) Figure 6: Spherical bubble collapse problem: Distribution of φ f and (cid:18) ∂φ f ∂r (cid:19) in the vicinity of the bubble surface allows the solver to be stable at large values of the density ratio ρ l ρ v . It is worth emphasizing, in thecurrent study, the density ratio is taken ≈ . Turbulent cavitating flow over a hydrofoil: A validation study In this section, we validate the proposed numerical method on the case of turbulent cavitating flow overa hydrofoil. This is an often-encountered scenario in marine propellers where fluid acceleration over thehydrofoil surface can result in very low pressures and cavity inception near the blade leading edge. Forthis study, we use cavitation model B because of its origins in flows involving large bubble clusters. It iscomputationally less expensive, and has been previously applied to the study of macro-scale cavitationover hydrofoils. The turbulence model is validated first on non-cavitating flow before proceeding to thecase of cavitating flow. Fig. (7) shows the general schematic of the computational domain used in thesections to follow. C is the hydrofoil chord length, α is the angle of attack of the incoming flow, H is thechannel height and ν T is the kinematic turbulence viscosity. Specific details are given in the respectivecase descriptions. u f · n f = 0 , σ · n f = 0 u f = [ U ∞ , , u f · n f = 0 , σ · n f = 0 Outf low ∇ φ f · n f = 0 Inf low SlipSlip . C α . Cφ f = 1 σ · n f = 05 . C Hy x ∇ ν T · n f = 0 N o - slip u f = 0 , ν T = 0 Figure 7: Representative computational domain and associated boundary conditions for cavitating flow over hydrofoil
In the current study, we employ a hybrid URANS-LES model to model turbulence. We validate theturbulence model on non-cavitating flow over a hydrofoil section. A NACA66 hydrofoil with chordlength C = 0 . m and a span equal to 0 . C is used in the study. The height of the channel H = 1 . C .Fig. (8) shows the computational grid. The grid consists of 76364 hexahedral elements (eight-nodebricks) resolving the cross-section. 30 nodes resolve the spanwise direction. A target y + = yu τ /ν = 1was maintained in the discretization of the hydrofoil boundary layer, where y is the height of the firstnode from the wall, u τ is the friction velocity and ν is the kinematic viscosity of the single phase liquid.The flow Reynolds number is 8 × , with a free-stream velocity U ∞ = 5 .
333 m s − . Liquid waterat 25 ◦ C is taken as the working fluid, with a single-phase density of ρ l = 999 .
19 kg m − . A Dirichletvelocity condition equal to U ∞ is set at the inlet with a natural traction-free outflow condition at the flowexit. A symmetric boundary condition is used on the top and bottom surfaces with periodic conditionson the spanwise surfaces. The time-step ∆ t of the numerical simulation is set to 1 . × − ( T ref / T ref = C/U ∞ ).For our validation, the angle of attack α of the hydrofoil is varied between 0 ◦ − ◦ for the non-cavitatingflow condition, and the time-averaged lift ( C L ) and drag ( C D ) coefficients were monitored. Fig. (9) showsthe comparison of C L and C D obtained from the numerical simulation with the experimental results in[38]. The predicted numerical results are observed to agree well with the experimental values in thenon-cavitating regime, and are within the uncertainties for C L (∆ C L = 0 . C D (∆ C D = 0 . igure 8: Computational mesh for NACA66 hydrofoil. Inlay showing mesh in the vicinity of the hydrofoil (a) (b) Figure 9: NACA66 hydrofoil problem: Predicted C L and C D in the non-cavitating range compared with the experimentalvalues from [38] To proceed further, the case of turbulent cavitating flow over a hydrofoil is studied using our numericalsolver. A NACA0012 hydrofoil with α = 1 ◦ is consider with a Reynolds number of 2 × . The chordlength C = 1 with the height of the channel H = 6 C . Figure (10) shows the computational grid usedin the study. The domain is discretized with 90340 hexahedral elements and a 2D periodic boundarycondition is applied in the spanwise direction. A target y + equal to 1 is enforced at the hydrofoilsurface. The fluid domain is initialized with a liquid phase fraction of φ f = 1. A freestream velocity of U ∞ = 1m s − is applied at the inlet as a dirichlet boundary condition. A traction-free outflow boundarycondition is used, weakly setting p ∞ = 0. The phase fraction is set to 1 at the inlet, along with aNeumann boundary condition at the outflow. A density ratio ( ρ l /ρ v ) of 1000 is used in the study. Thecavitation number of the flow is defined as σ = p ∞ − p v . ρ l U ∞ , where p ∞ is the free-steam hydrostatic pressure.In the current study, the cavitation number of the flow is set to 0 .
42. The vapor pressure is set to meetthe cavitation number of the flow. A time-step of ∆ t = 1 × − s ( t ref / C p = p f − p ∞ . ρ l U ∞ ) on the22 igure 10: Computational mesh for NACA0012 hydrofoil. Inlay showing mesh in the vicinity of the hydrofoilFigure 11: Comparison of predicted pressure coefficient C p with results presented in [49] suction surface of the hydrofoil is compared against the numerical results of [49] based on Kunz et al.[35] model. A good agreement can be seen between these two numerical studies. The cavity pressurehas been captured consistently, which is important for the prediction of hydrodynamic loads on thehydrofoil surface. A sharp gradient in the pressure over the hydrofoil suction surface can be observed,23orresponding to the cavity closure location. Inlay shows the partial sheet cavity on the suction surfaceof the hydrofoil, marked by an iso-contour of φ f = 0 .
95. Notably, it is observed to stabilize over time toform a thin attached partial cavity on the hydrofoil surface. A slight shift in the cavity closure locationis observed which is attributed to the difference in the cavitation and turbulence models used in thetwo studies. In addition, the model coefficients used in the current study are taken from[25] and [51]where the geometries studied were different. No cavity separation and shedding are observed, and a fullyattached turbulent flow exists over the hydrofoil surface. The absence of a re-entrant jet is consistentwith the observations of [22] for thin cavities.
5. Application to fluid-structure interaction of caviating hydrofoil
In this closing section, we explore the ability of the proposed implementation to model freely movingcavitating hydrofoil and to predict some key features of cavitating flow over hydrofoils. We also test thecompatibility of the solver on configurations with moving solid boundaries. Before proceeding to ourfully-coupled cavitation and FSI demonstration, the primary objective of this section is to evaluate thefeasibility of the implementation for turbulent cavitating flows over a stationary hydrofoil section. Weconsider the same NACA0012 hydrofoil geometry as in Sec. (4.2) with the following modifications. Theheight of the channel H is reduced to 1 . C . The computational grid is also modified to a hybrid gridsystem comprising hexahedral (8-node brick) and prism (6-node wedge) elements. This is to preventthe formation of highly skewed elements (in fully hexahedral grids) resulting from mesh deformationsat high angles of attack. Figure (12) shows the representative computational grid in the vicinity of thehydrofoil, demonstrating deformation to α = 10 ◦ . Cavitation model B is used for the following studies.The values of the model coefficients and fluid properties defined in Section (4.2) are used. The cavitationnumber σ is increased to 1 . Re is set to 1 × , with U ∞ = 1m s − . Atraction-free outflow is used, setting the pressure weakly to 0 at the outflow boundary. Figure 12: Hybrid computational grid for NACA0012 hydrofoil deformed to α = 10 ◦ Leading-edge cavitation over hydrofoils can behave in different ways based on flow conditions definedby the flow Reynolds number Re , the cavitation number σ and the angle of attack α of the incomingflow. Under certain combinations of these flow parameters, periodic cavity growth and shedding canbe observed. Integral to this cavity shedding process is the formation of a periodic re-entrant jet alongthe hydrofoil surface. This jet periodically flows along the hydrofoil surface from the cavity closurelocation towards the leading edge, detaching the attached cavity. The capturing of this flow phenomenais of interest to the study of propeller vibration and noise because of two reasons. First, the periodicshedding of the cavity leads to periodic fluctuations in the hydrodynamic loading on the propeller blade.If the frequency of this loading is close to the natural frequency of the blade, it can result in structuralexcitation - leading to vibration and tonal noise emission. Second, the shed cavities exist in the form of24louds of vaporous bubbles. These bubbles are prone to collapsing near the trailing edge of the bladewith localized high-amplitude water hammer impacts. This can contribute to broadband underwaternoise emission. We investigate here the turbulent cavitating flow over a hydrofoil section at α = 10 ◦ .Figure (13) shows the results of the study at different instances during one cavity shedding cycle. Aniso-contour of φ f = 0 .
95 (in red) is used to represent the cavity surface. Also shown are the contoursof the vorticity in the direction z out of the plane of the figure. Fig. (13a) shows the inception of aleading-edge cavity. The cavity is observed to grow, primarily collocated with a leading-edge vortex(LEV), to the extent of the hydrofoil chord. In Fig. (13d) the clockwise rotating LEV interacts witha counter-clockwise trailing edge vortex (TEV). We observe that the interaction, along with a reversepressure gradient, leads to the formation of a re-entrant jet along the hydrofoil suction surface. Thisdetaches the cavity which is shed in the form of pockets (clouds) and is convected with the mean flow.Fig. (14) shows the streamlines in the domain at the beginning of the shedding process. The formationof the re-entrant jet can be observed to originate at the intersection of the LEV and TEV. Remark . We note the ability of the present implementation to capture some select physics of interestin turbulent cavitating flow over hydrofoils. Only preliminary results are presented for demonstration.The accuracy of the shedding frequency needs to be investigated, and can require careful calibration ofcoefficients in the cavitation model and the modified turbulent viscosity. Detailed investigation is beyondthe scope of the current work, and will be explored in future studies.
We next investigate turbulent cavitating flow over the NACA0012 section subject to a prescribedperiodic pitching motion. By using our ALE-based FSI framework, we linearly ramp the angle of attackof the hydrofoil between 0 ◦ and 15 ◦ . Fig. (15) shows the first four cycles of the prescribed motion.The frequency f pitch of the motion is 2 Hz. The rest of the study parameters are kept the same as inSection (5.1). Figure (16) shows the cavities in the domain marked by the iso-contour φ f = 0 .
95. Alsoplotted are the contours of z-vorticity. At the high pitching frequency f pitch = 2 Hz and the associatedhydrofoil acceleration, cavities are observed to originate on both surfaces because of the low pressuresduring the pitching motion. The primary cavity generation is at the hydrofoil leading edge. The cavityshedding frequency is low compared to the pitching frequency and the hydrofoil suction surface is seento be perennially covered by cavities that are continuously being shed and convected with the mean flow.The collocation between the vortices and the cavities can be discerned. Detailed investigations on theinfluence of the pitching frequency on the cavity and vortex shedding frequency are reserved for futurework.We note the compatibility of the cavitation and the flow solvers with structural deformation. Theversatility provided by the partitioned coupling was exploited to couple the additional solvers for turbu-lence modeling and ALE mesh update. The numerical solution obtained is stable and 2 −
6. Conclusion
A robust and accurate variational finite element formulation for the numerical study of cavitating flowshas been presented for stationary and moving hydrofoils. We introduced novel stabilized linearizationsof two cavitation transport equation models based on the two-phase homogeneous mixture theory. Thenumerical implementation has been employed to study two cavitating flow configurations with vastlydifferent temporal and spatial scales. An initial verification study on the micro-scale collapse of aspherical vaporous bubble has been shown to maintain stability across a range of time steps. Accuratesolutions of the pressure field were obtained, devoid of spurious numerical oscillations. The solutionof the phase indicator is demonstrated to be bounded, and the solver is numerically stable at largedensity ratios. The implementation has been validated on fully turbulent cavitating flow over a hydrofoilwith good agreement with previous numerical and experimental studies. We also explored the ability ofthe implementation to predict select characteristic features of macro-scale cavitating flows, including re-entrant jets, periodic cavity shedding and cavity-vortex interaction. Further, we examined the versatilityof the implementation to be coupled with solvers for studying fluid-structure interaction, demonstratingthe case of a pitching hydrofoil. In future work, the authors plan to extend the implementation to25 a) t (b) t + 0 . t ∞ (c) t + 0 . t ∞ (d) t + 0 . t ∞ (e) t + 0 . t ∞ (f) t + 0 . t ∞ (g) t + 0 . t ∞ (h) t + 0 . t ∞ Figure 13: Contours of Z-vorticity (positive is clockwise, negative is anti-clockwise) during one shedding cycle. Cavity (inred) marked by iso-contour of φ f = 0 .
95. First figure marked by time t for reference. Subsequent figures marked in termsof t and t ∞ = C/U ∞ igure 14: Vortex interaction and formation of re-entrant jet at trailing edge of hydrofoil. Cavity (in red) marked atiso-contour of φ f = 0 .
95. Figure 15: Representative prescribed pitching motion study fully-coupled FSI studies in cavitating flows. One potential application is cavitating flow-inducedvibrations taking into account the hydroelastic response of structures. Another potential application isthe study of material erosion resulting from cavitation bubble collapse using fully-Eulerian fluid-solidformulations.
Acknowledgements a) (+)5 ◦ (h) ( − )5 ◦ (b) (+)10 ◦ (g) ( − )10 ◦ (c) (+)12 . ◦ (f) ( − )12 . ◦ (d) (+)14 ◦ (e) ( − )14 ◦ Figure 16: Contours of Z-vorticity (positive is clockwise, negative is anti-clockwise) during one pitching cycle. Cavity (inred) marked by iso-contour of φ f = 0 .
95. (+) marks the upward pitching stroke and (-) marks the downward stroke. eferences [1] MPI: A message-passing interface version 3.1. ( . Technical report, 2015.[2] Deniz Tolga Akcabay and Yin Lu Young. Influence of cavitation on the hydroelastic stability ofhydrofoils. Journal of Fluids and Structures , 49:170–185, 2014.[3] VH Arakeri, H Higuchi, and REA Arndt. A model for predicting tip vortex cavitation characteristics.
Journal of Fluids Engineering , 110(2):190–193, 1988.[4] R Arndt, P Pennings, J Bosschers, and T Van Terwisga. The singing vortex.
Interface focus , 5(5):20150025, 2015.[5] Michael R Bailey, Robin O Cleveland, Tim Colonius, Lawrence A Crum, Andrew P Evan, James ELingeman, James A McAteer, Oleg A Sapozhnikov, and JC Williams. Cavitation in shock wavelithotripsy: the critical role of bubble activity in stone breakage and kidney trauma. In
IEEESymposium on Ultrasonics, 2003 , volume 1, pages 724–727. IEEE, 2003.[6] A Bayram and A Korobenko. Variational multiscale framework for cavitating flows.
ComputationalMechanics , 66(1):49–67, 2020.[7] Christopher Earls Brennen.
Cavitation and Bubble Dynamics . Cambridge University Press, 2013.doi: 10.1017/CBO9781107338760.[8] Alexander N Brooks and Thomas JR Hughes. Streamline upwind/petrov-galerkin formulations forconvection dominated flows with particular emphasis on the incompressible navier-stokes equations.
Computer methods in applied mechanics and engineering , 32(1-3):199–259, 1982.[9] Sandeep Reddy Bukka, Rachit Gupta, Allan Ross Magee, and Rajeev Kumar Jaiman. Assessmentof unsteady flow predictions using hybrid deep learning based reduced-order models.
Physics ofFluids , 33(1):013601, 2021.[10] John Carlton.
Marine propellers and propulsion . Butterworth-Heinemann, 2018.[11] G. L. Chahine, A Kapahi, J.-K. Choi, and C.-T. Hsiao. Modeling of surface cleaning by cavitationbubble dynamics and collapse.
Ultrasonics Sonochemistry , 29:528–549, 2016.[12] Georges L Chahine, Andrew F Conn, Virgil E Johnson Jr, et al. Cleaning and cutting with self-resonating pulsed water jets. In . Citeseer, 1983.[13] Yongliang Chen and SD Heister. Modeling hydrodynamic nonequilibrium in cavitating flows.
Journalof Fluids Engineering , 118(1):172–178, 1996.[14] Jintai Chung and GM Hulbert. A time integration algorithm for structural dynamics with improvednumerical dissipation: the generalized- α method. Journal of Applied Mechanics , 1993.[15] O Coutier-Delgosha, JL Reboud, and Y Delannoy. Numerical simulation of the unsteady behaviourof cavitating flows.
International journal for numerical methods in fluids , 42(5):527–548, 2003.[16] Manish Deshpande, Jinzhang Feng, and Charles L Merkle. Cavity flow predictions based on theeuler equations.
Journal of Fluids Engineering , 1994.[17] Jean-Pierre Franc and Jean-Marie Michel.
Fundamentals of cavitation , volume 76. Springer science& Business media, 2006.[18] Leopoldo P Franca and S´ergio L Frey. Stabilized finite element methods: Ii. the incompressiblenavier-stokes equations.
Computer Methods in Applied Mechanics and Engineering , 99(2-3):209–233, 1992.[19] Ebrahim Ghahramani, Mohammad Hossein Arabnejad, and Rickard E Bensow. A comparativestudy between numerical methods in simulation of cavitating bubbles.
International Journal ofMultiphase Flow , 111:339–359, 2019. 2920] A Gnanaskandan and Krishnan Mahesh. A numerical method to simulate turbulent cavitating flows.
International Journal of Multiphase Flow , 70:22–34, 2015.[21] Eric Goncalves, Maxime Champagnac, and Regiane Fortes Patella. Comparison of numerical solversfor cavitating flows.
International Journal of Computational Fluid Dynamics , 24(6):201–216, 2010.[22] Shridhar Gopalan and Joseph Katz. Flow structure and modeling issues in the closure region ofattached cavitation.
Physics of fluids , 12(4):895–911, 2000.[23] Isaac Harari and Thomas JR Hughes. What are c and h?: Inequalities for the analysis and design offinite element methods.
Computer Methods in Applied Mechanics and Engineering , 97(2):157–192,1992.[24] M-C Hsu, Yuri Bazilevs, Victor M Calo, Tayfun E Tezduyar, and Thomas JR Hughes. Improvingstability of stabilized and multiscale formulations in flow simulations at small time steps.
ComputerMethods in Applied Mechanics and Engineering , 199(13-16):828–840, 2010.[25] Biao Huang, Antoine Ducoin, and Yin Lu Young. Physical and numerical investigation of cavitatingflows around a pitching hydrofoil.
Physics of Fluids , 25(10):102109, 2013.[26] Thomas JR Hughes and Garth N Wells. Conservation properties for the galerkin and stabilisedforms of the advection–diffusion and incompressible navier–stokes equations.
Computer methods inapplied mechanics and engineering , 194(9-11):1141–1159, 2005.[27] Ghaleb A Husseini, Mario A Diaz De La Rosa, Eric S Richardson, Douglas A Christensen, andWilliam G Pitt. The role of cavitation in acoustically activated drug delivery.
Journal of ControlledRelease , 107(2):253–261, 2005.[28] Kenneth E Jansen, Christian H Whiting, and Gregory M Hulbert. A generalized- α method forintegrating the filtered navier–stokes equations with a stabilized finite element method. Computermethods in applied mechanics and engineering , 190(3-4):305–319, 2000.[29] B Ji, XW Luo, Roger EA Arndt, Xiaoxing Peng, and Yulin Wu. Large eddy simulation and theo-retical investigations of the transient cavitating vortical flow structure around a naca66 hydrofoil.
International Journal of Multiphase Flow , 68:121–134, 2015.[30] Jingwei Jiang, Haopeng Cai, Cheng Ma, Zhengfang Qian, Ke Chen, and Peng Wu. A ship pro-peller design methodology of multi-objective optimization considering fluid–structure interaction.
Engineering Applications of Computational Fluid Mechanics , 12(1):28–40, 2018.[31] Vaibhav Joshi and Rajeev K Jaiman. A positivity preserving variational method for multi-dimensional convection–diffusion–reaction equation.
Journal of Computational Physics , 339:247–284, 2017.[32] Vaibhav Joshi and Rajeev K Jaiman. A variationally bounded scheme for delayed detached eddysimulation: Application to vortex-induced vibration of offshore riser.
Computers & fluids , 157:84–111, 2017.[33] Vaibhav Joshi and Rajeev K Jaiman. A positivity preserving and conservative variational schemefor phase-field modeling of two-phase flows.
Journal of Computational Physics , 360:137–166, 2018.[34] W Kerr, JF Shannon, and RN Arnold. The problems of the singing propeller.
Proceedings of theInstitution of Mechanical Engineers , 144(1):54–90, 1940.[35] Robert F Kunz, David A Boger, Thomas S Chyczewski, D Stinebring, H Gibeling, and T Govindan.Multi-phase cfd analysis of natural and ventilated cavitation about submerged bodies. In
Proceedingsof the 3rd ASME-JSME Joint Fluids Engineering Conference , 1999.[36] Abe H Lee, Robert L Campbell, Brent A Craven, and Stephen A Hambric. Fluid–structure in-teraction simulation of vortex-induced vibration of a flexible hydrofoil.
Journal of Vibration andAcoustics , 139(4), 2017. 3037] Hyoungsuk Lee, Min-Churl Song, Jung-Chun Suh, and Bong-Jun Chang. Hydro-elastic analysisof marine propellers based on a bem-fem coupled fsi algorithm.
International journal of navalarchitecture and ocean engineering , 6(3):562–577, 2014.[38] Jean-Baptiste Leroux, Jacques Andr´e Astolfi, and Jean Yves Billard. An experimental study ofunsteady partial cavitation.
Journal of Fluids Engineering , 126(1):94–101, 2004.[39] B Maines and Roger EA Arndt. The case of the singing vortex.
Journal of Fluids Engineering ,1997.[40] Pieter Maljaars, Mirek Kaminski, and Henk Den Besten. Boundary element modelling aspects forthe hydro-elastic analysis of flexible marine propellers.
Journal of Marine Science and Engineering ,6(2):67, 2018.[41] Charles L Merkle. Computational modelling of the dynamics of sheet cavitation. In
Proc. of the3rd Int. Symp. on Cavitation, Grenoble, France, 1998 , 1998.[42] Tharindu P Miyanawala and Rajeev K Jaiman. An efficient deep learning technique for the navier-stokes equations: Application to unsteady wake flow dynamics. arXiv preprint arXiv:1710.09099 ,2017.[43] Yuriy A Pishchalnikov, Oleg A Sapozhnikov, Michael R Bailey, James C Williams Jr, Robin OCleveland, Tim Colonius, Lawrence A Crum, Andrew P Evan, and James A McAteer. Cavitationbubble cluster activity in the breakage of kidney stones by lithotripter shockwaves.
Journal ofendourology , 17(7):435–446, 2003.[44] Lord Rayleigh. Viii. on the pressure developed in a liquid during the collapse of a spherical cavity.
The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science , 34(200):94–98,1917.[45] Donald Ross and WA Kuperman. Mechanics of underwater noise, 1989. URL https://doi.org/10.1121/1.398685 .[46] Youcef Saad and Martin H Schultz. Gmres: A generalized minimal residual algorithm for solvingnonsymmetric linear systems.
SIAM Journal on scientific and statistical computing , 7(3):856–869,1986.[47] G¨unter H Schnerr and J¨urgen Sauer. Physical and numerical modeling of unsteady cavitationdynamics. In
Fourth international conference on multiphase flow , volume 1. ICMF New Orleans,2001.[48] G¨unter H Schnerr, Ismail H Sezal, and Steffen J Schmidt. Numerical investigation of three-dimensional cloud cavitation with special emphasis on collapse induced shock dynamics.
Physics ofFluids , 20(4):040703, 2008.[49] Inanc Senocak and Wei Shyy. Numerical simulation of turbulent flows with sheet cavitation. 2001.URL http://resolver.caltech.edu/cav2001:sessionA7.002 .[50] Inanc Senocak and Wei Shyy. A pressure-based method for turbulent cavitating flow computations.
Journal of Computational Physics , 176(2):363–383, 2002.[51] Inanc Senocak and Wei Shyy. Interfacial dynamics-based modelling of turbulent cavitating flows,part-1: Model development and steady-state computations.
International journal for numericalmethods in fluids , 44(9):975–995, 2004.[52] Ismail Hakki Sezal.
Compressible dynamics of cavitating 3-D multi-phase flows . PhD thesis, Tech-nische Universit¨at M¨unchen, 2009.[53] Farzin Shakib, Thomas JR Hughes, and Zdenˇek Johan. A new finite element formulation for compu-tational fluid dynamics: X. the compressible euler and navier-stokes equations.
Computer Methodsin Applied Mechanics and Engineering , 89(1-3):141–219, 1991.3154] Ashok K Singhal, Mahesh M Athavale, Huiying Li, and Yu Jiang. Mathematical basis and validationof the full cavitation model.
Journal of Fluids Engineering , 124(3):617–624, 2002.[55] WD Song, Hong MH, Lukyanchuk B, and Chong TC. Laser-induced cavitation bubbles for cleaningof solid surfaces.
Journal of Applied Physics , 95:2952—-2956, 2004.[56] Eleanor Stride and Constantin Coussios. Nucleation, mapping and control of cavitation for drugdelivery.
Nature Reviews Physics , 1(8):495–509, 2019.[57] Tayfun E Tezduyar, Sanjay Mittal, SE Ray, and R Shih. Incompressible flow computations with sta-bilized bilinear and linear equal-order-interpolation velocity-pressure elements.
Computer Methodsin Applied Mechanics and Engineering , 95(2):221–242, 1992.[58] Pieter Van Oossanen.
Calculation of performance and cavitation characteristics of propellers in-cluding effects on non-uniform flow and viscosity . PhD thesis, Delft University of Technology, 1974.URL http://resolver.tudelft.nl/uuid:daef4d65-e0cc-4796-88e0-6c6b2d2a54ac .[59] Qin Wu, Biao Huang, Guoyu Wang, and Shuliang Cao. The transient characteristics of cloudcavitating flow over a flexible hydrofoil.
International Journal of Multiphase Flow , 99:162–173,2018.[60] AG Zijlstra.
Acoustic surface cavitation . PhD thesis, University of Twente, 2011.[61] Philip J Zwart, Andrew G Gerber, Thabet Belamri, et al. A two-phase flow model for predictingcavitation dynamics. In