Implementation and validation of an extended Schnerr-Sauer cavitation model for non-isothermal flows in OpenFOAM

Available online at www.sciencedirect.com Available online at www.sciencedirect.com Available online at www.sciencedirect.com

Availableonline onlineatatwww.sciencedirect.com www.sciencedirect.com Available

ScienceDirect Energy Procedia 00 (2017) 000–000 ScienceDirect Energy Procedia 00 (2017) 000–000 Energy Procedia (2017) 000–000 Energy Procedia 126 (201709) 58–65 Energy Procedia 00 00 (2017) 000–000

nd

www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia www.elsevier.com/locate/procedia

72 Conference of the Italian Thermal Machines Engineering Association, ATI2017, 6–8 nd September 2017, Lecce, Italy Association, ATI2017, 6–8 72 72nd Conference Conference of of the the Italian Italian Thermal Thermal Machines Machines Engineering Engineering Association, ATI2017, 6–8 September 2017, Lecce, Italy September 2017, Lecce, Italy Implementation and validation of an extended Schnerr-Sauer

Implementation and of extended Implementation andforvalidation validation ofonan an extended Schnerr-Sauer cavitation model non-isothermal flows in Schnerr-Sauer OpenFOAM The 15th International Symposium District Heating and Cooling cavitation model for non-isothermal flows in OpenFOAM a,∗ a cavitation model for non-isothermal flows in OpenFOAM Maria Grazia De Giorgi , Antonio Ficarella , Donato Fontanarosaa a,∗ a a Assessing the feasibility offorusing the demand-outdoor Maria Grazia De Ficarella ,, Donato Fontanarosa a,∗, Antonio a University of Salento, Dep. of Engineering Innovation, Via peraheat Monteroni, 73100-Lecce, Italy Maria Grazia De Giorgi Giorgi , Antonio Ficarella Donato Fontanarosa University of Salento, Dep. of Engineering for Innovation, district Via per Monteroni, 73100-Lecce, Italy temperature function for a long-term heat demand forecast University of Salento, Dep. of Engineering for Innovation, Via per Monteroni, 73100-Lecce, Italy a a a

Abstract

I. Andrića,b,c*, A. Pinaa, P. Ferrãoa, J. Fournierb., B. Lacarrièrec, O. Le Correc

Abstract In the present work cavitation in liquid hydrogen and nitrogen was investigated by using the open source toolbox OpenFOAM. Abstract a IN+ Center forperformed Innovation, by Technology anda Policy Research model, - Instituto Superior Técnico, Av. Roviscomixture Pais 1, 1049-001 Lisbon, Portugal Simulations of mass and transfer on the homogeneous approach in OpenFOAM. combination b In the presentwere work cavitation inmeans liquid hydrogen nitrogen wasbased investigated by using the open source toolbox Veolia Recherche & Innovation, 291 Avenue Dreyfous Daniel, 78520 Limay, France In the present work cavitation in liquid hydrogen and nitrogen was investigated by using the open source toolbox OpenFOAM. with the Volume ofperformed Fluid (VOF) fora the reconstruction the liquid-vapor Two additional transport in equations were c Simulations were by method means of mass transfer model, based on theinterface. homogeneous mixture combination Département Systèmes Énergétiques et Environnement - IMT Atlantique, rue Alfred Kastler, 44300approach Nantes, France Simulations were performed by means of advection a mass transfer model, based on the4homogeneous mixture approach in combination considered, i.e. the liquid(VOF) volume fraction and the temperature equation. The implementation of an extended Schnerrwith the Volume of Fluid method for the reconstruction the liquid-vapor interface. Two additional transport equations were with the Volume of Fluid method forofthethe reconstruction the in liquid-vapor interface. Two additional transport equations were Sauer model allowed for (VOF) the introduction thermal effects termsequation. of latent heatimplementation release/absorption convective heat considered, i.e. the liquid volume fraction advection and the temperature The of anand extended Schnerrconsidered, i.e.the theliquid-vapor liquid volume fractionAadvection and the temperature equation. The implementation ofthe an local extended Schnerrtransfer inside interface. set of Antoine-like equations relate the saturation conditions to conditions. Sauer model allowed for the introduction of the thermal effects in terms of latent heat release/absorption and convective heat Sauer model allowed for the introduction of the thermal effects in terms of latent heat release/absorption and convective heat Abstract transfer inside the liquid-vapor interface. A set of Antoine-like equations relate the saturation conditions to the local conditions. transfer the liquid-vapor A set of Antoine-like equations relate the saturation conditions to the local conditions. c 2017 inside  The Authors. Publishedinterface. by Elsevier Ltd.

Conference of most the Italian Thermal Machines Engineering Peer-review under responsibility of Elsevier the scientific the 72nd as heating networks are commonly addressed in theofliterature one of the effective solutions for decreasing the © 2017 Authors. Published by Ltd. committee cDistrict  2017 The The Authors. Published by Elsevier c  2017 The Authors. Published by Elsevier Ltd. nd Association. nd Conference Peer-review under responsibility of the scientific committee of the 72 of the Italian Thermal Machines Engineering greenhouse gas emissions from the building sector. These systems require high investments which are returned through the heat Peer-review under responsibility of the scientific committee of the 72nd Conference of the Italian Thermal Machines Engineering Peer-review under responsibility of the scientific committee of the 72 Conference of the Italian Thermal Machines Engineering Association sales. Due to the changed climate conditions and building renovation policies, heat demand in the future could decrease, Association. Keywords: Association.cavitation, thermo-sensitive fluids, thermal effects, Schnerr-Sauer, Volume Of Fluid, OpenFOAM prolonging the investment return period. Keywords: cavitation, thermo-sensitive fluids, thermal effects, Schnerr-Sauer, Volume Of Fluid, OpenFOAM cavitation, thermo-sensitive fluids, the thermal effects, of Schnerr-Sauer, Volume Of Fluid, OpenFOAM Keywords: The main scope of this paper is to assess feasibility using the heat demand – outdoor temperature function for heat demand forecast. The district of Alvalade, located in Lisbon (Portugal), was used as a case study. The district is consisted of 665 buildings that vary in both construction period and typology. Three weather scenarios (low, medium, high) and three district 1.renovation Introduction scenarios were developed (shallow, intermediate, deep). To estimate the error, obtained heat demand values were compared with results from a dynamic heat demand model, previously developed and validated by the authors. 1. Introduction 1. In Introduction turbomachinery cavitation has tochange be avoided since it causes erosion, noise and which areapplications coupled The results showed thatthe when only weather is considered, the margin of error could bevibrations, acceptable for some with a strong loss of performances and a severe deterioration of the operational lifetime of the machine. Thermo(the error in annual demand was lower than 20% for all weather scenarios considered). However, after introducing renovation In the cavitation has to be avoided since it causes erosion, noise vibrations, which are coupled In turbomachinery turbomachinery thecryogenic cavitation has to behot avoided since causes erosion, noise and and vibrations, whichproperties are coupled scenarios, the error value increased up to 59.5% (depending on itthe and renovation combination considered). sensitive fluids, such as fluids and fluids, present aweather strong dependence ofscenarios thermo-physical on with a strong loss of performances and a severe deterioration of the operational lifetime of the machine. Thermowith a strong loss of performances and a severe deterioration of the operational lifetime of the machine. ThermoThe value of slope coefficient increased on average within the range of 3.8% up to 8% per decade, that corresponds to the temperature and pressure. They experience a significant temperature drop inside the cavity due to phase change and sensitive fluids, such fluids and hot present aa strong dependence of properties on sensitive fluids, such as asofcryogenic cryogenic fluids hot fluids, fluids, present strong dependence ofonthermo-physical thermo-physical properties on decrease in the number heating hours ofand 22-139h during the heating season (depending the combination of weather and heat transfer phenomena, which play an essential role in the cavitating process. Consequently the prediction of the temperature and They experience a significant temperature drop the due to change and temperature and pressure. pressure. They experience significant temperature drop inside inside the cavity cavityper due to phase phase changeon and renovation considered). On theflows otheraishand, function intercept increased for 7.8-12.7% decade (depending the behavior ofscenarios non-isothermal cavitating of primary importance. heat transfer phenomena, which play an essential role in the cavitating process. Consequently the prediction of the coupled scenarios). The values suggested could be used to modify the function parameters for the scenarios considered, and heat transfer phenomena, which play an essential role in the cavitating process. Consequently the prediction of the Cavitation in thermo-sensitive fluids has been largely investigated. Stahl et al. [1] and Stepanoff [2] first studied behavior of cavitating flows improve the accuracy of heat demand estimations. behavior of non-isothermal non-isothermal cavitating flows is is of of primary primary importance. importance.

thermal effectsinduring cavitation in highhas temperature water flows by introducing the parameter B-factor, defined as Cavitation fluids been investigated. Stahl [1] and [2] Cavitation in thermo-sensitive thermo-sensitive fluids has volumes. been largely largely investigated. Stahl et etofal. al.the [1]latent and Stepanoff Stepanoff [2] first first studied studied the ratio between the vapor and the liquid Due to the equivalence heat of vaporization and thermal effects during cavitation in high temperature water flows by introducing the parameter B-factor, defined as © 2017 The Authors. Published by Elsevier Ltd. thermal effects during cavitation in high temperature water flows by introducing the parameter B-factor, defined as the sensible heat absorption, B-factor represents also the ratio between the actual temperature drop and the nominal the ratio between the vapor and the liquid volumes. Due to the equivalence of the latent heat of vaporization and Peer-review under responsibility of the Scientific Committee of The 15th International Symposium on District Heating and the ratio between the vapor and the liquid volumes. Due to the equivalence of the latent heat of vaporization and Cooling. the the sensible sensible heat heat absorption, absorption, B-factor B-factor represents represents also also the the ratio ratio between between the the actual actual temperature temperature drop drop and and the the nominal nominal

∗ Corresponding author. Address: University of Salento, Dep. of Engineering for Innovation, Research Center for Energy and Environment Keywords: Heat demand; Forecast; Climate change (UNISALENTO-DII-CREA), Via per Monteroni, LECCE I-73100, Italy. Tel: +39 0832297759. ∗ Corresponding author. Address: of Salento, Dep. of Engineering for Innovation, Research Center for Energy and Environment ∗ Corresponding author. Address: University University of Salento,(Maria Dep. ofGrazia Engineering for Innovation, Research Center for Energy Environment E-mail addresses: [email protected] De Giorgi)., [email protected] (Antonioand Ficarella)., do(UNISALENTO-DII-CREA), Via per Monteroni, LECCE I-73100, Italy. Tel: +39 0832297759. (UNISALENTO-DII-CREA), per Monteroni, LECCE I-73100, Italy. Tel: +39 0832297759. [email protected] (Donato Fontanarosa). E-mail addresses: [email protected] (Maria Grazia De Giorgi)., [email protected] (Antonio Ficarella)., doE-mail  caddresses: 1876-6102 2017 [email protected] Authors.(Donato Published by Elsevier Ltd.(Maria Grazia De Giorgi)., [email protected] (Antonio Ficarella)., [email protected] Fontanarosa). nd [email protected] (Donato Fontanarosa). Peer-review under responsibility ofPublished the scientific committee 1876-6102 © 2017 The Authors. byElsevier Elsevier Ltd.of the 72 Conference of the Italian Thermal Machines Engineering Association. c 2017 The Authors. Published by 1876-6102  Ltd. c under 1876-6102  2017 The Authors. Published by Elsevier Ltd. of Thend15th International Symposium on District Heating and Cooling. Peer-review responsibility thescientific Scientific Committee Peer-review under responsibility ofofthe committee of the 72 Conference of the Italian Thermal Machines Engineering Association. Peer-review under responsibility of the scientific committee of the 72nd Conference of the Italian Thermal Machines Engineering Association. 1876-6102 © 2017 The Authors. Published by Elsevier Ltd. Peer-review under responsibility of the scientific committee of the 72nd Conference of the Italian Thermal Machines Engineering Association 10.1016/j.egypro.2017.08.057

2

Maria Grazia De Giorgi et al. / Energy Procedia 126 (201709) 58–65 De Giorgi et. al / Energy Procedia 00 (2017) 000–000

59

temperature drop due to the heat absorption. Hord et al. [3–5] extended the investigation to cryogenic cavitating flows of liquid hydrogen and liquid nitrogen in different cavitating conditions for different geometries. Concerning the numerics of cavitating flows, various computational approaches have been developed with different levels of complexity. Basically they can be classified into two categories: single phase methods ([6, 7]) and multiphase methods. Multiphase models treat both liquid and vapor phases by means of different approaches. Among them, mixture models have aroused increasing interest since they have shown promising results with a reduced computational effort. A special class of mixture models is represented by mass transfer models. They introduce a transport equation for the vapor volume fraction, and a cavitation model which computes vapor destruction/production source terms related to the phase change process. Different cavitation models have been used in a large number of numerical works dealing with a wide range of applications. Rodio et al. [8] investigated the influence of convective heat transfer modeling on the estimation of thermal effects in cryogenic cavitating flows. Chen et al. [9] studied the behavior of cavitating flows of water and fluoroketone in a wide range of free-stream temperatures and velocities. Zhu et al. [10, 11] studied the interactions of vortices, thermal effects and cavitation in liquid hydrogen cavitating flows. De Giorgi et al. [12] simulated with success cryogenic cavitation by introducing both inertial and heat transfer control bubble growth. Zhang et al. [13] focused on the numerical modeling of thermodynamic phase-change theory during cavitation of cryogenic flows. Goncalv`es [14] developed a 4-equations model for non isothermal cavitation; recently, this model was used by Goncalv`es and Zeidan [15] for numerical investigations concerning the unsteady cavitation in liquid hydrogen flows. Despite of the enhancements achieved, cavitation in thermo-sensitive fluids is still far from being fully modeled owing to the strong coupling between heat transfer mechanisms and phase change phenomena. In this context, the present work aims to provide an extended cavitation model for non-isothermal cryogenic flows with OpenFOAM, which is considered a good alternative to well-known commercial codes nowadays. In particular, a mass transfer model was implemented in combination with a Volume of Fluid (VOF) phase-fraction based interface capturing approach, and an extended Schnerr-Sauer cavitation model [16]. The Reynolds-Averaged Navier-Stokes (RANS) equations have been used by assuming of a homogeneous mixture of two incompressible and immiscible fluids. The advection of the vapor phase into the liquid bulk is modeled by means of the transport equation of the liquid volume fraction. Starting from the isothermal model, the conservation of the energy has been introduced in the form of a temperature equation. Therefore, two source terms have been added into the temperature equation so as to consider the latent heat of vaporization/condensation and the convective heat transfer. The saturation properties of the fluid have been expressed as a function of the local conditions by means of Antoine-like equations. 2. Numerical approach 2.1. Governing Equations The set of governing equations is composed by the conservation equations of the mixture mass, momentum and energy, shown in Eqs. (1)-(3), and the transport equation of the liquid volume fraction given in Eq. (4). The quantities denoted by the subscript m refer to the mixture, while the subscripts l, v and turb denote to the liquid, the vapor and the turbulent properties respectively. Um is the velocity vector field of the mixture, ρ is the density, µ and ν are the dynamic and the kinematic viscosities, g is the standard gravity, Fσ is the surface tension force, T is the temperature, c p is the specific heat capacity at constant pressure, λ is the thermal diffusivity, and hL is the latent heat of vaporization/condensation. t and x are the temporal and the spatial variables, while the subscripts (i, j,k) denote the directions of the Cartesian coordinates. The liquid volume fraction is defined as α = (VlV+Vl v ) . ∂ρm + ∇ · (ρm Um ) = 0 ∂t    ∂U ∂Um, j 2 ∂Uk  ∂(ρm Um,i ) ∂(ρm Um,i Um, j ) ∂p ∂ Q˙ conv m,i =− + + − δi j + ρm gi + Fσ,i + + (µm + µturb ) ∂t ∂x j ∂xi ∂x j ∂x j ∂xi 3 ∂xk (ρc p )m ∂T R h T L + ∇ · (Um T ) = [λm + λturb ]∇2 T + ∂t (ρc p )m  ρ  ∂α m + Um . ∇(α) = R ∂t ρl ρv

(1) (2) (3) (4)

Maria Grazia De Giorgi et al. / Energy Procedia 126 (201709) 58–65 De Giorgi et. al / Energy Procedia 00 (2017) 000–000

60

3

The mixture properties are defined by means of phase averaging  to obtain   so as   ρm = ρl α + ρv (1 − α), µm = k k µl α + µv (1 − α), (ρc p )m = ρl c p,l α + ρv c p,v (1 − α) and λm = ρc p = ρc p α + ρck p (1 − α). The turbulent thermal m l v   t , where νt is the turbulent kinematic viscosity and Prturb is the turbulent Prandtl diffusivity is defined as λturb = Prνturb

number, which has been has been set to 0.85. The mass transfer rate R is computed by means of the cavitation model presented in the subsection 2.2. The latent heat source term RT and the convective heat transfer term Q˙ conv are modeled by means of a proper thermal model, which is described in the subsection 2.3. 2.2. Cavitation model

The total mass transfer rate R is composed by the destruction term Re due to the evaporation, and production term Rc related to the condensation. Based on the Rayleight-Plesset equation of the dynamics of a spherical bubble, the Schnerr-Sauer model leads to the following relations: 3 ρv ρl Rc = +Cc α(1 − α) ρm Rb



2 |p − pvap (T )| 3 ρl

3 ρv ρl Rv = −Cv α(1 + αNUC − α) ρm Rb

and



2 |p − pvap (T )| 3 ρl

(5)

where the subscripts c and v refer to the condensation and the evaporation processes. In the equations above, Cc and Cv are two tuning parameters used to control condensation and evaporation independently. The vapor pressure pvap (T ) is derived from the saturation pressure p sat (T ) corrected by an estimation of the local values of the turbulent pressure fluctuations, as proposed by Singhal et al. [17] (see Eq. (6)). A numerical nucleation term αNUC is introduced in order to initialize the cavitation process. It is defined as a function of the diameter of nuclei dNUC by the Eq. (7). The radius of the vapor bubble Rb is related to the density of nuclei nb the by means of relations in Eq. (7). pvap (T ) = p sat (T ) + αvapor = (1 + αNUC − α) =

1

pturb , 2

4 πR3b nb 3 , + 43 πR3b nb

pturb = 0.39ρm k with

αNUC =

(6) 3 πdNUC nb 6

1+

3 πdNUC nb 6

(7)

2.3. Thermal model The thermal model deals with two terms, i.e. the latent heat source term RT and the convective heat transfer term Q˙ conv . RT is similar to the mass transfer rate R, differing just in the tuning coefficients Cc,T and Cv,T , which replace Cc and Cv . This modification has been introduced because of the strong coupling inside the liquid-vapor interface between the condensation and the convective heat transfer. The convective heat transfer term is subtracted from the latent heat source term. As a result, the global thermal source becomes Q˙ tot = Q˙ c + Q˙ v with   ρv ρl 2 |p − pvap (T )| − hb |T sat (p) − T | hL ρm 3 ρl    3 ρv ρl 2 |p − pvap (T )| − α) − hb |T sat (p) − T | hL Rb ρm 3 ρl

3 Q˙ c = (Rc,T hL − Q˙ conv,c ) = +Cc,T α(1 − α) Rb Q˙ v = (Rv,T hL − Q˙ conv,v ) = −Cv,T α(1 + αNUC



(8) (9)

where T sat (p) is the saturation temperature, expressed as a function of the local pressure. The convective heat transfer coefficient hb is computed locally by means of the relation by Christopher et al. [18], shown in Eq. (10). hb =

  ¯ 1/2 pvap (T ) 2σ ˆ ρv h2L  M 1− , with 2−σ ˆ T vap 2πNA T v 2hL ρv

T v = T sat,0 +

2σT sat,0 . hL ρm Rb

(10)

¯ is the molecular weight, T sat,0 is the saturation temperature In the expressions above, NA is the Avogadro number, M ˆ is the accommodation factor of the fluid, and σ is the surface tension of the liquid at at the upstream pressure p0 , σ the saturation conditions.

4

Maria Grazia De Giorgi et al. / Energy Procedia 126 (201709) 58–65 De Giorgi et. al / Energy Procedia 00 (2017) 000–000

61

Fig. 1. Computational domain and boundaries.

3. Numerical setup 3.1. Reference Test Cases, Computational domain and Boundary Conditions (BCs) c Version 3.0.1, based Numerical simulations were performed by using the open source CFD toolbox OpenFOAM on a Finite Volume formulation. Two different test cases have been used for the validation of the implemented cavitation model, i.e. the Hord’s experimental tests numbered 248C and 296B [4]. They refer to two different fluids, the liquid hydrogen (LH2) and the liquid nitrogen (LN2) respectively. All the settings are reported in Tab. 1, where Σ is the thermodynamic parameter defined by Brennen [19]. A set of extended Antoine equations relate saturation conditions to the local conditions, so as to introduce the dependence of the saturation pressure to the local temperature p sat (T ), as well as the dependence of the saturation temperature to the local pressure T sat (p). Table 1. Numerical setup.

Test Case

Fluid

T∞ [K]

[m/s]

U∞

kσ [-]

Cavity Length [m]

Re

248C

LH2

20.46

51.2

1.60

0.0139

1.76 × 107

296B

LN2

88.54

23.7

1.61

0.0127

Σ

[ms−3/2 ]

[-]

7

1.08 × 10

963652 202498

The computational domain and boundaries are shown in Fig. 1. The geometry reproduces the experimental domain used by Hord [4], whose results have been used for comparisons. The mesh is 2D and composed of 99956 structured cells. 530 elements on the hydrofoil. A maximum y+ equal to 12 was reached. Concerning the setup of the BCs, a no-slip condition was imposed on the hydrofoil and the upper wall. Each simulation was constrained by setting the inlet velocity and the pressure outlet, while the initial temperature field was supposed uniform. The database of the National Institute of Standard and Technology (NIST) [20] was used for the determination of the transport properties of fluids. 3.2. Numerics The governing equations were solved by using the PIMPLE algorithm [21], which is provided in OpenFOAM for simulating unsteady incompressible flows. It applies a prediction-correction routine by merging the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm, and the Pressure Implicit with Splitting of Operator (PISO) algorithm. The Multidimensional Universal Limiter for Explicit Solution (MULES) algorithm [22] was used in order to compute the field of the liquid volume fraction. A fixed time step equal to 1 × 10−5 permitted to enhance the numerical stability and achieve a quasi-steady solution. The non cavitating steady-state solution initialized unsteady computations. The k − ω Shear Stress Transport model (SST) by Menter [23] is used to take account of the turbulence. The coefficients of the model have been set to their standard values. Time derivatives were computed with the Euler scheme. The following numerical schemes were used in reference to the convective fluxes: the Self Filtered Central Differencing (SFCD) scheme for the velocity field, the van Leer scheme for the volume fraction of liquid, the Gauss upwind scheme for the temperature field, and the Total Variation Diminishing (TVD) scheme for all remaining fields.

Maria Grazia De Giorgi et al. / Energy Procedia 126 (201709) 58–65 De Giorgi et. al / Energy Procedia 00 (2017) 000–000

62

5

4. Results and discussion Several simulations were performed, each of that refers to the specific numerical setup reported in Tab. 2. Computations were stopped at 0.01 s, when a quasi-steady solution was reached. Table 2. Test matrix. Notes: Cc = Cv = 1.0.

Simulation #

Test Case

TEqn

˙ conv Q

Cc,T

Cv,T [-]

[-]

[m−3 ]

nb

dNUC

SIM1

248C

OFF

OFF

-

-

-

SIM2

248C

ON

ON

0.5

1.0

0.007

1.6 × 108

6.0 × 10−5

SIM3

248C

ON

OFF

0.6

1.0

-

5.4 × 108

6.0 × 10−5

[-]

[-]

[-]

σ ˆ

SIM4

248C

ON

ON

0.6

1.0

0.007

SIM5

296B

OFF

OFF

-

-

-

SIM6

296B

ON

OFF

1.0

1.0

-

SIM7

296B

ON

ON

0.6

1.0

0.007

SIM8

296B

ON

ON

0.7

1.0

0.007

SIM9

296B

OFF

OFF

-

-

-

[m]

5.4 × 108 5.4 × 10

8

4.8 × 10

8

4.8 × 10

8

6.0 × 10−5 6.0 × 10−5

2.4 × 108

6.0 × 10−6

4.8 × 108

6.0 × 10−6

6.0 × 10−6 6.0 × 10−6

4.8 × 108

6.0 × 10−6

The static pressure and temperature on the hydrofoil surface were used for validation. Figg. 2 and 3 show the comparison of numerical results with the experimental data by Hord, for test cases 248C and 296B respectively. 3.5

×10 5

20.8 SIM2 SIM3 SIM4 SIM1 Hord

3 2.5

20.6 20.4

T [K]

p [Pa]

20.2 2 1.5

20 19.8 19.6

1

SIM2 SIM3 SIM4 SIM1 Hord

19.4 0.5

19.2

0

(a)

19 0.03

0.04

0.05 x [m]

0.06

0.07

(b)

0.03

0.04

0.05 x [m]

0.06

0.07

Fig. 2. Validation of numerical results. Comparison with the experimental data: (a) static pressure; (b) temperature. Hord’s Test Case: 248C (LH2).

The entire investigation started from the analysis of the isothermal model, i.e. simulations SIM1 for the liquid hydrogen and SIM5 for the liquid nitrogen. The evaporation and condensation mass transfer rates were numerically balanced (Cv = Cc = 1), so that a sensitivity analysis was conducted just by tuning the nucleation parameters nb and dNUC . The best agreement with the experimental data was found by setting (nb = 1.6 × 108 m−3 , dNUC = 6.0 × 10−5 m) for the test case involving the liquid hydrogen, and (nb = 2.4 × 108 m−3 , dNUC = 6.0 × 10−6 m) for the liquid nitrogen test case. The pressure profiles in Figg. 2(a) and 3(a) show that the cavity length is in accordance with the experimental data, yet the pressure drop inside the vapor bubble is limited by the vapor pressure which is constant and referred to the inlet temperature. Thus, the activation of the latent heat sources without taking account of the convective heat transfer was introduced first. A new calibration of the nucleation parameters was required since thermal effects produced a contraction of the vapor cavity. In particular the nuclei diameter was preserved while the number of nuclei nb was gradually increased. The best agreement with the Hord’s data was obtained with the simulation SIM3 of the liquid hydrogen and the simulation SIM6 involving the liquid nitrogen. Concerning the extent of the increment of the number of nuclei, nb was doubled in case of the liquid hydrogen (SIM3), and it was more than tripled for the liquid nitrogen

6

De Giorgi et. al Procedia 00 (2017) 000–000 Maria Grazia De/ Energy Giorgi et al. / Energy Procedia 126 (201709) 58–65 9

#10 5

89 SIM7 SIM6 SIM8 Hord SIM9 SIM5

8 7

88.5

88 T [K]

6 p [Pa]

63

5 4

87.5

87

SIM7 SIM6 SIM8 Hord SIM9 SIM5

3 86.5

2 1

(a)

86 0.03

0.04

0.05 x [m]

0.06

0.07

(b)

0.03

0.04

0.05 x [m]

0.06

0.07

Fig. 3. Validation of numerical results. Comparison with the experimental data: (a) static pressure; (b) temperature. Hord’s Test Case: 296B (LN2).

simulation (SIM6). Our investigations showed that the pressure and temperature profiles have a similar shape for both the liquid hydrogen and the liquid nitrogen test cases when the thermal model was set at the same conditions. Simulation SIM3 and simulation SIM6 have a different thermal condensation coefficient, Cc,T . In fact, when the condensation is numerically balanced with respect to the evaporation as done in the simulation SIM6, the increased nucleating behaviour related to the higher nb produces a higher temperature drop, which numerically is coupled with a temperature increment at the closure of the cavity close to the hydrofoil surface. Thus, in order to control this numerical effect the thermal condensation was reduced by means of the coefficient Cc,T , as done in the simulation SIM3. Nevertheless, a good agreement with the temperature measurements by Hord just by varying Cc,T was far from being reached unless the convective heat transfer was introduced. Simulations SIM2 and SIM4 for the liquid hydrogen, as well as simulations SIM7 and SIM9 involving the liquid nitrogen, refer to the fully non-isothermal model. They differ in the value of the thermal condensation coefficient Cc,T which influences both the latent and the convective heat transfer in the condensation regions. The introduction of the convective heat transfer combined with the reduction of Cc,T produced results in agreement with the experimental data. as shown by simulations SIM4 and SIM2 for the liquid hydrogen, and simulations SIM8 ans SIM7 in case of the liquid nitrogen. The best matching was obtained with simulations SIM4 and SIM8, for the liquid hydrogen and the liquid nitrogen respectively. The evaluation of the effect of the coefficient Cc,T arises out of the comparison of simulations SIM4 and SIM2 for the liquid hydrogen, and simulations SIM8 and SIM7 in case of the liquid nitrogen. In general Cc,T weakly affects the pressure field because the lower temperature downstream in the wake produces a reduction of the vapor pressure which has no effect on the pressure field. On the contrary, the temperature field is very sensitive to Cc,T . In fact, it determines the temperature peak so as to influence the temperature field downstream in the wake. In particular the temperature peak rises as Cc,T increases due to the fact that it falls inside the evaporation region. The isothermal simulation SIM9 highlights that the cavity length extends over the entire hydrofoil when the best fitting non-isothermal nb is introduced into the isothermal models. This numerical behavior is represented in Fig. 4(b), where the contour of the vapor fraction α of the simulation SIM9 (plot in the center) shows the vapor overproduction in comparison with the simulation SIM5 (plot on the bottom). The same effect is observed in Fig. 3(a) by considering the green curve (SIM9) and the yellow curve (SIM5). Further outcomes come out from the analysis of the contour plots of the vapor cavity, the mass transfer rates, the temperature field and the heat sources. The isothermal model computes cavities having a similar shape for both liquid hydrogen and liquid nitrogen, as shown in Fig. 4 by comparing simulations SIM1 and SIM5. With respect to the liquid nitrogen, the liquid hydrogen is characterized by a greater numerical length of the cavity, in accordance with LH2 LN2 the experimental lengths by Hord (Lexp = 0.0139 m > Lexp = 0.0127 m). In case of fully non-isothermal modeling, the liquid nitrogen (simulation SIM8) shows a more squashed vapor cavity with respect to the liquid hydrogen (simulation SIM4) which has a larger cavitation region. Finally, the shape of the vapor cavities is strictly related to both the mass transfer rates displayed in Fig. 5 and the heat sources shown in Fig. 6. Furthermore, the simulation involving the liquid nitrogen exhibits mass transfer rates higher than the one with hydrogen, as displayed in Fig. 6(a),

Maria Grazia De Giorgi et al. / Energy Procedia 126 (201709) 58–65 De Giorgi et. al / Energy Procedia 00 (2017) 000–000

64

7

while the heat sources Q˙ have the same magnitude but a different spatial field which is in line with the vapor volume fraction field. alpha

alpha

(b)

(a)

Fig. 4. Field of the vapor volume fraction: (a) Test case 248C: SIM1, SIM4 and SIM3, with an experimental cavity length equal to 0.0139m; (b) Test case 296B: SIM5, SIM9 and SIM8, with an experimental cavity length equal to 0.0127m.

mdot

(a)

T

(b) Fig. 5. Contour plots for the test case 248C: (a) mass transfer rate; (b) temperature.

qdot

mdot

qdot

mdot

(a)

(b)

Fig. 6. Influence of the liquid properties on cavitation fields: (a) mass transfer rate; (b) heat transfer rate. Test case 248C: SIM1 and SIM4. Test case 296B: SIM7.

8

Maria Grazia De Giorgi et al. / Energy Procedia 126 (201709) 58–65 De Giorgi et. al / Energy Procedia 00 (2017) 000–000

65

5. Conclusions and future works The present work dealt with the numerical modeling of thermodynamic cavitation for a cryogenic flow on hydrofoils, by using the open source CFD toolbox OpenFOAM. Results showed a good fitting with experimental data for two test cases, one for the liquid hydrogen and the other one for the liquid nitrogen. The outcome is that condensation plays a significant role since it is strongly coupled with the convective heat transfer mechanism at the liquid-vapor interface. In particular results showed that the introduced numerical coefficient of condensation Cc,T , which is related to the thermal model, weakly affects on the pressure field; on the contrary the temperature field is very sensitive to it. The drawback of the proposed model is the introduction of several numerical parameters, which are related to each test case. Next steps will be to focus on the overcoming of this drawback. After a sensitivity analysis aimed at the quantification of the influence of the most relevant numerical parameters, an optimization strategy will be searched in order to get an optimum numerical setup for some specific applications. References [1] H. A. Stahl, A. J. Stepanoff, N. Phillipsburg, Thermodynamic aspects of cavitation in centrifugal pumps, ASME J. Basic Eng 78 (1956) 1691–1693. [2] A. Stepanoff, Cavitation properties of liquids, Journal of Engineering for Power 86 (2) (1964) 195–199. [3] J. Hord, L. Anderson, W. Hall, Cavitation in liquid cryogens, volume i: Venturi, nasa rept, Technical Report NASA CR-2054, National Aeronautics and Space Administration (May 1972). [4] J. Hord, Cavitation in liquid cryogens. 2: Hydrofoil, Technical Report NASA-CR-2156, National Aeronautics and Space Administration (January 1973). [5] J. Hord, Cavitation in liquid cryogens iiiogives, Technical Report NASA CR-2242, National Aeronautics and Space Administration (May 1973). [6] J. Feng, C. L. Merkle, Numerical modeling of the thermodynamic effects of cavitation, Journal of fluids engineering 119 (1997) 421. [7] B. Stutz, J.-L. Reboud, Two-phase flow structure of sheet cavitation, Physics of Fluids 9 (12) (1997) 3678–3686. [8] M. G. Rodio, M. De Giorgi, A. Ficarella, Influence of convective heat transfer modeling on the estimation of thermal effects in cryogenic cavitating flows, International Journal of Heat and Mass Transfer 55 (23) (2012) 6538–6554. [9] T. Chen, B. Huang, G. Wang, H. Zhang, Y. Wang, Numerical investigation of thermo-sensitive cavitating flows in a wide range of free-stream temperatures and velocities in fluoroketone, International Journal of Heat and Mass Transfer 112 (2017) 125–136. [10] J. Zhu, Y. Chen, D. Zhao, X. Zhang, Extension of the schnerr–sauer model for cryogenic cavitation, European Journal of Mechanics-B/Fluids 52 (2015) 1–10. [11] J. Zhu, D. Zhao, L. Xu, X. Zhang, Interactions of vortices, thermal effects and cavitation in liquid hydrogen cavitating flows, International Journal of Hydrogen Energy 41 (1) (2016) 614–631. [12] M. G. De Giorgi, A. Ficarella, Simulation of cryogenic cavitation by using both inertial and heat transfer control bubble growth, in: Proc. 39th AIAA Fluid Dynamics Conference, San Antonio, Texas, 2009, pp. 22–25. [13] X. Zhang, Z. Wu, S. Xiang, L. Qiu, Modeling cavitation flow of cryogenic fluids with thermodynamic phase-change theory, Chinese Science Bulletin 58 (4-5) (2013) 567–574. [14] E. Goncalv`es, Modeling for non isothermal cavitation using 4-equation models, International Journal of Heat and Mass Transfer 76 (2014) 247–262. [15] E. Goncalv`es, D. Zeidan, Numerical simulation of unsteady cavitation in liquid hydrogen flows, International Journal of Engineering Systems Modelling and Simulation 9 (1) (2017) 41–52. [16] G. H. Schnerr, J. Sauer, Physical and numerical modeling of unsteady cavitation dynamics, in: Fourth International Conference on Multiphase Flow, ICMF New Orleans, 2001. [17] A. K. Singhal, M. M. Athavale, H. Li, Y. Jiang, Mathematical basis and validation of the full cavitation model, Transactions-American Society of Mechanical Engineers Journal of Fluids Engineering 124 (3) (2002) 617–624. [18] D. M. Christopher, H. Wang, X. Peng, Numerical analysis of the dynamics of moving vapor bubbles, International Journal of Heat and Mass Transfer 49 (19) (2006) 3626–3633. [19] C. E. Brennen, Cavitation and bubble dynamics, Cambridge University Press, 2013. [20] E. Lemmon, M. McLinden, D. Friend, P. Linstrom, W. Mallard, Nist chemistry webbook, nist standard reference database number 69, National Institute of Standards and Technology, Gaithersburg. R [21] T. Holzmann, Mathematics, numerics, derivations and openfoam. [22] S. M. Dami´an, An extended mixture model for the simultaneous treatment of short and long scale interfaces, Dissertationsschrift, Universidad Nacional del Litoral, Argentinien. Krepper, E., Lucas, D., Frank, T., Prasser, H.-M., & Zwart, PJ (2008). The inhomogeneous MUSIG model for the simulation of polydispersed flows. Nuclear Engineering and Design 238 (7) (2013) 1690–1702. [23] F. R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA journal 32 (8) (1994) 1598–1605.