Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 UDK - UDC 621.224:532.528 Kratki znanstveni prispevek - Short scientific paper (1.03) Eksperimentalno numerična analiza kavitacijskega toka okoli lopatičnega profila Experimental and numerical analyses of the cavitational flows around a hydrofoil Ignacijo Biluš - Leopold Škerget - Andrej Predin - Matjaž Hriberšek Namen prispevka je predstaviti analizo kavitacijskih tokovnih razmer okoli lopatičnega profila NACA. Predstavljen je fizikalno-matematični model v obliki Navier-Stokesovih enačb, zapisanih za fizikalne lastnosti zmesi voda - para in na podlagi prenosne enačbe za ohranitev mase vodne pare, s katero je bil modeliran nastanek, razširjanje in izginjanje vodne pare v toku zmesi. Pri modeliranju fazne spremembe je bila uporabljena poenostavljena Rayleigh-Plessetova enačba, v kateri so upoštevani parametri, ki pomembneje vplivajo na dinamiko tokovnih pojavov v bližini osamljenega krogelnega parnega mehurčka v obdajajoči kapljevini. Matematično fizikalni model je bil vključen v programski paket CFX 5.6, rezultati pa primerjani z rezultati preizkusa, izvedenega v kavitacijskem tunelu. © 2005 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: lopatice turbinske, tok kavitacijski, analize toka, simuliranje numerično) In this paper we present an analysis of the cavitation flow conditions around a NACA hydrofoil. The mathematical model, in the form of Navier Stokes equations, based on the additional transport equation for vapour mass fraction inception, propagation and condensation is presented for the mixture’s (water -water vapour) properties. A simplified Rayleigh-Plesset equation is used when the phase change is modelled, where the parameters that influence the flow dynamics near the individual spherical bubble, surrounded by the liquid, are considered. Mathematical/physical model is included in the CFD code CFX 5.6. Simulation results are compared with experimental results from the cavitation tunnel, where the tested hydrofoil (blade) was placed. © 2005 Journal of Mechanical Engineering. All rights reserved. (Keywords: turbine blades, cavitation flow, water turbine flow analysis, numerical simulations) 0 UVOD 0 INTRODUCTION Kavitacija je pojav uparjanja vode in kondenzacije vodne pare v bližini lopatic v vodni turbini, ki je posledica neenakomernih razmer v tokovnem in temperaturnem polju tekočine. Opazovanje in razlage pojava kavitacije v vodnih turbinah so se začele že pred dvestopetdesetimi leti, ko je leta 1754 Euler prvi opisal omenjeni pojav. Od takrat se znanstveniki po svetu ukvarjajo s preučevanjem tega pojava, ki se pojavlja v hidravličnih strojih, ladijskih vijakih in hidravličnih napravah, vendar zaradi zapletenosti in številnih vplivnih parametrov kavitacija do danes fizikalno matematično še ni v popolnosti opisana in rešena. Cavitation is the water evaporation and condensation of water close to a turbine blade’s surface that is a consequence of unequal condi-tions in the flow and the temperature field. Cavita-tion observations of water turbines started about 250 years ago, when in 1754 Euler described the cavitation phenomenon. Since then, scientists have studied this phenomenon, which appears in hydrau-lic machinery, ships’ propellers and many hydrau-lic devices. However, the phenomenon is not clearly understood, either mathematically or physi-cally, and many parameters are connected to the phenomenon. 103 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 Pojav kavitacije je tesno povezan z lokalnimi časovno spremenljivimi tokovnimi lastnostmi, kot so lokalne hitrosti in termodinamični tlak, ki je nižji od pripadajočega uparjalnega tlaka kapljevine pri delovni temperaturi sistema. Eksperimentalne študije kažejo na več vrst nestabilnosti, med katerimi je najbolj izrazito pulziranje kavitacijskega oblaka [1]. Temeljni raziskovalni projekti s področja kavitacije so osredotočeni na opis mehanizmov dinamike dvofaznih tokov in medsebojne interakcije med kapljevito in plinasto fazo. Dodatno je moč obravnavo razširiti tudi na interakcijo tok - trdnina, saj imajo lokalna nihanja tlaka lahko velik vpliv na trdne površine sestavnih delov turbinskočrpalnih sistemov. S pospešenim razvojem merilnih metod in intenzivnim kopičenjem eksperimentalnih rezultatov je bilo v zadnjem času omogočeno tudi numerično modeliranje kavitacijskih tokov. Slednje izhaja s področja modeliranja večfaznih tokov, ki so bili razviti na temelju dvotekočinskega modela, torej na modelu Euler-Euler analize tokovnega polja v mirujoči prostorski točki, za vsako tekočinsko fazo posebej ([2] in [3]) . Zaradi velike zahtevnosti obravnavanega sistema, tako z vidika fizikalno-matematičnega modeliranja kakor tudi z vidika učinkovitega numeričnega izračuna, je bil glavni motiv preučevanja kavitacijskih tokov fizikalno-matematični model, ki bi z zadovoljivo natančnostjo opisal prenosne pojave v primeru kavitacije in omogočil razmeroma kratke računske čase na dostopni računalniški opremi. Tako so bili v raziskavi izmed številnih vplivov, ki se pojavljajo pri kavitaciji, obravnavani le pomembnejši, torej fazna sprememba (nukleacija), rast in velikost mehurčkov, turbulentni tok, snovske lastnosti obeh faz in geometrijska oblika obtekane trdne površine (lopatice). Deležnečistočin plinov, ki ne kondenzirajo, ter vplivi površinske napetosti in prenosnih pojavov na medfazni površini v prispevku niso obravnavani oziroma so obravnavani na makroskopskem integralnem nivoju. Prednost predstavljenega modela v primerjavi s podobnimi ([4] in [5]) je, da so vsi parametri eksperimentalno enostavno določljivi. Omenjeni motiv je bil dosežen z uporabo matematično-fizikalnega modela v obliki Navier-Stokesovih enačb, zapisanih in rešenih za fizikalne lastnosti zmesi voda - para in na podlagi prenosne enačbe za ohranitev mase vodne pare, s katero smo modelirali nastanek, razširjanje in izginjanje vodne pare v toku zmesi. The cavitation phenomenon is connected to locally time dependent flow properties, such as the local flow velocities and the thermodynamic pres-sure which is lower than vapour pressure at the given temperature. Experimental studies show a lot of in-stabilities sources, but the major one is in the strong cavitation cloud pulsation [1]. Basic research projects from the interest-ing area of cavitational flow deal with the dynamics of two-phase flow mechanisms that consider the interactions between the liquid and gas phases. It is possible to expand the analysed flow mechan-ics to the fluid–solid interaction studies, because local pressure oscillations have a large influence on the solid surfaces of parts of turbine/pump systems. The intensive development of measuring methods and the large number of experimental re-sults available, has enabled the numerical modelling of cavitational flows. The models are based on multi-phase-flow modelling principles, developed for two fluids. This means the Euler-Euler model for flow-field analyses in the fixed point of 3D flow, for each individual phase ([2] and [3]). Because of the analysed system’s complex-ity from the physical/mathematical modelling point of view, as well as from the powerful simulation development aspect, the main aim of the cavitational flow study was an efficient flow model for transient cavitating conditions with a short calculation time using the available computer resources. To satisfy these limitations only the most important or influen-tial parameters, i.e., the phase change (nucleation), bubble growth (bubble size), turbulent flow, fluid properties, and the geometry of the solid surface (blade surface) were considered. The impurities, the non-condensed gas and the surface tension’s influ-ence at the interfaces are not considered in this paper. The listed parameters are only considered from the macroscopic or integral point of view. Compar-ing to similar models ([4] and [5]), presented model has distinct advantage since it is easy to determine all included experimental parameters. The aim was achieved by using a math-ematical/physical model in the form of the Navier-Stokes equation, which was written and solved for the physical properties of water/vapour mixture based on the transient equation form for vapour continuity. Using this equation we modelled the ap-pearance, growth and collapse of the vapour bub-bles. 104 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 1 VODILNE ENAČBE ENOFAZNIH TOKOV 1.1 Zakon ohranitve mase Zakon ohranitve mase je izpeljan iz ugotovitve, da je masa sistema nespremenljiva veličina. Integralsko obliko zapišemo z naslednjo enačbo: 1 GOVERNING EQUATIONS FOR ONE-PHASE FLOW 1.1 Mass continuity equation The continuity equation results from the fundamental physical principle that mass is con-served. In integral form it can be written as: dt jrdV + jrvjnjdS = 0 (1). V 1.2 Zakon ohranitve gibalne količine Rezultirajoča sila okolice na prostornino je enaka časovnemu prirastku gibalne količine v prostornini in dotoku gibalne količine skozi njegovo površino. Integralska oblika zakona ohranitve gibalne količine je tako: drvi 1.2 Momentum equation The resulting force on the volume element is equal to the time increment of the momentum in the volume and the flux across the element surface. The momentum equation in integral form can be writ-ten as: dt kjer so: v hitrostno polje, fi prostorninska sila, p termodinamični tlak in r strižna napetost. 1.3 Zakon ohranitve turbulentne kinetične energije in raztrosa turbulentne kinetične energije Model turbulentne kinetične energije k in raztrosne hitrosti turbulentne kinetične energije e, oziroma model k-e, je najpomembnejši dvoenačbni turbulentni model, ki temelji na postopku turbulentne viskoznosti. Turbulentne napetosti (-r0v~v'j) izrazimo z Boussinesqueovim približkom: \rvidV + \virvjnjdS = \rfidV + \(-pöij+Tij)njdS (2), where v is the velocity flow field,/is the body force, p is the thermodynamic pressure, and ix. is the shear stress. 1.3 Conservation of turbulent kinetic energy and turbulent kinetic energy dissipation The two-equation model for the turbulent kinetic energy k and the dissipation of turbulent kinetic energy e, or the k - e model, is the most important two-equation turbulent model that is based on the turbulent viscosity principle. The turbulent stresses (-p^'~j ) are expressed with the Boussinesque approximation as follows: (-r0v~v'j) = r0PT kjer sta k povprečna turbulentna kinetična energija turbulentnih odstopanj in vT=(r,T/r) turbulentna viskoznost. Člen 2Öijk/3 je razširitev osnovne Boussinesqueove podmene in ga lahko prištejemo statičnemu tlaku. Značilne veličine so definirane z izrazi, npr. značilna hitrost je: raztrosna hitrost turbulentne kinetične energije e pa: dvi dvj ----+ , dx dx V j i % dij r0k (3), where k is the averaged kinetic energy of the turbulent fluctuations, and vT=(uT/r) the turbulent viscosity. The factor 2öijk/3 is an extension of the Boussinesque hypothesis that can be added to the static pressure. The characteristic properties are defined with the characteristic velocity 4k (4), and the dissipation velocity of the turbulent kinetic energy e : dvidvi dxdx (5) Eksperimentalno numerična analiza - Experimental and Numerical Analyses 105 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 podaja spremembo turbulentne energije toka v toplotno. Obe veličini k in e določimo iz dodatnih posameznih parcialnih diferencialnih enačb, ki vsebujejo nove stalnice in funkcije. Za k velja enačba: which defines the conversion of the turbulent energy into heat. Both the quantities k and e are deter-mined from additional, specific differential equations, which include new constants and functions. For k the following equation is valid dk _ dk d — + vi — = — dt dx dx P -e in podobno za e de _ de d — + vi — = — dt dx dx s uT | dk skJ dxi and similarly for e de dx ee C1e P -C2e kk (6) (7), kjer so stalnice modela C =0,09, sk=1,0, se=1,3, C1e =1,44 in C2e =1,92 [6]. m 2 MATEMATIČNI MODEL 2.1 Dinamika krogelnega mehurčka Osamljen krogelni mehurček je najpreprostejša pojavna oblika parne faze v obdajajoči kapljevini, saj ne upošteva deformacij mehurčka, ki se pojavijo zaradi nehomogenosti tokovnega polja, niti medsebojnega vpliva mehurčkov. Kljub omenjenim poenostavitvam je dobra osnova za modeliranje prenosnih pojavov v dvofaznih kavitacijskih tokovih, ki se pojavljajo v turbinskih strojih. Obravnavajmo iz tega razloga [7], krogelni mehurček polmera R(t) v obdajajoči kapljevini temperature T"«, in tlaka p(t). Vrednost p(t) naj bo znana, temperatura T., pa nespremenljiva. Na začetku predpostavimo nestisljivo kapljevino rL=konst. Vsebina mehurčka naj bo homogena, temperatura T in tlak p (t) v mehurčku pa enakomerna. Dinamična viskoznost kapljevine hL naj bo nespremenljiva. Radialno lego v tekočini podamo z razdaljo r od središča mehurčka (sl. 1a), tlak v poljubni točki T zunaj mehurčka označimo s p(r,t), radialna hitrost je u(r,t) in temperatura T(r,t). Iz zakona ohranitve mase izhaja ([7] in [8]), da je zaradi spremembe površine mehurčka s kvadratom polmera r, hitrost u(r, t) definirana kot: where the standard values of the constants are: C =0.09, sk=1.0, se=1.3, C1e =1.44 in C2e =1.92 [6]. 2 MATHEMATICAL MODEL 2.1 Bubble dynamic An isolated, spherical bubble represents the simplest form of vapour surrounded by a liquid. This situation does not consider the bubble deformation that appears as a consequence of an non-homogenous flow field or the interaction between the vapour bubbles. However, in spite of these simplifying conditions, the mathematical model can be considered as a good basis for numerical modelling of the transient flow conditions for the two-phase flow that appears in turbine machines. Consider a spherical bubble [7] of radius RB(t) (where t is the time) in an infinite domain of liquid whose temperature and pressure far from the bubble are T«, and p(t), respectively. The temperature is assumed to be constant and the pressure is assumed to be known. At the beginning it is also assumed that the fluid is incompressible (rL=const). The bubble contents are homogenous and the temperature T and pressure pB(t) within the bubble are uniform. The fluid dynamic viscosity hL is constant. The radial position within the liquid will be denoted by the distance r from the bubble center (Figure 1.a). Let the pressure at an arbitrary point T be represented by p(r,t), the radial outward velocity by u(r,t), and the temperature by T(r,t) From the mass continuity equation ([7] and [8]) it follows that the spherical bubble’s surface varies with the square of the radius r. Therefore, the flow velocity u(r,t) can be defined by: u(r,t) Zapišimo gibalno enačbo [7] za smer r: dt 2 rl\ (8). The momentum equation for the direction r [7] can be written as: 106 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 velikd uüdaljeriusl ud niehuiuka far from the bubble KAPLJEVINA LIQUID KAPLJEVINA LIQUID površina mehurčka bubble surface PLIN / PARA CAS/VAPOUR pnvrsina mehurčka bubble sur'ace 5 Sl. 1. Krogelni mehurček v tekočini (a) in izsek iz površine krogelnega mehurčka (b) [7] Fig. 1. Spherical bubble in a flow (a) and part of the spherical bubble surface (b) [7] 1 dp du du -----— = — + u — p dr dt dr 1 d ( 2du 2— r — rdr\ dr 2u (9), ker je u=F(t)/r2 dobimo z ločitvijo spremenljivk in z integracijo v mejah p^ RB) na površini mehurčka where u=F(t)/r2 can be defined with an integration between p^p(r=RB) on the bubble surface (r-+RB): rL p) 1dF 1F2 r dt 2r (10). Napetost na elementu površine krogelnega mehurčka podamo z izrazom: Za primer, ko ni masnega pretoka prek lupine, je t=0, od koder izhaja: The tension on the bubble surface element is 4hL dRB R dt 2s R (11). In the case without mass flow over the bub-ble sphere we can write t=0. From this we can derive the following: p(r=RB)=pB Z vstavitvijo izraza (11) v enačbo (9) dobimo: - 1(pb - 43LdRB - 2EL Pl{ RB dt RB 4hL dRB R dt 2s R (12). Combining Equations (12) and (10) we obtain: 1_dF 1_F2 r dt 2 r 4 ker je F = RB2(dr / dt) , velja pri r=RB: rL (pB-p) where F = RB2(dr / dt) for r=RB: 4udR 2s d2R RB dt rLRB + 3dRB dt 2 21 dt (13), (14). Enačba (14) predstavlja končno obliko Rayleigh-Plessetove enačbe, kjer so na levi strani gonilni, viskozni in člen površinske napetosti, na desni strani pa vztrajnostni člen. Z izpeljano Rayleigh-Plessetovo enačbo lahko zadovoljivo predstavimo dinamiko toka [8] v neposredni bližini osamljenega krogelnega parnega Equation (14) represents the final form of the Rayleigh–Plesset equation. The first term is the driving term, the second term is the viscous term and the last term on the left-hand side of the equation is the surface-tension term. The inertia term is on the right-hand side of the equation. With the derived Rayleigh-Plesset equation the flow dynamics [8] in the closest surroundings of Eksperimentalno numerična analiza - Experimental and Numerical Analyses 107 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 mehurčka v kapljeviti okolici, zaradi cesar je omenjena enačba podlaga za popis interakcije kapljevito -plinasto pri nastanku, razširjanju in izginjanju mehurčkov vodne pare. 2.2 Homogeni dvofazni tokovni model Definirajmo navidezno gostoto homogene mešanice voda - vodna para, ki je odvisna od masnega deleža parne faze (suhosti) f: an individual spherical vapour bubble in the liquid can be presented. In this case the equation describes the interaction between the liquid and gas phase by bubble growth and disappearance. 2.2 Homogenous two phase flow model The density of a homogenous mixture of water and water vapour can be defined with: 1 f 1-f =+ rrV rL (15), kjer je prostorninski delež parne faze: where f represents the vapour mass fraction, connected with the volume fraction by the equation a= f r rV (16). Za prenosno enačbo masnega deleža pare pišemo: With this notation the vapour-phase mass-fraction transport equation can be written as follows: dt kjer sta R in R izvirna člena oziroma stopnja uparjanja in stopnja kondenzacije, ki sta funkciji tokovnih veličin (tlak, hitrost) in snovskih lastnosti (gostote kapljevite in plinaste faze, uparjalnega tlaka, površinske napetosti), medtem ko je G kinematična difuzivnost spremenljivke f. Enačba (16) izhaja iz teorije homogenega dvofaznega toka, kar pomeni primerno poenostavitev iz naslednjih razlogov: - v inženirski praksi se kavitacija pojavlja v področjih nizkega tlaka, kjer so hitrosti razmeroma visoke in ni zdrsa med tekočo in plinasto fazo; - po navadi je faza pare v obliki drobnih mehurčkov. Ker pa je v okolici le-teh treba izbrati ustrezni fizikalni model za izračun velikosti (polmera mehurčka) in sile upora, kar pomeni velik problem, saj splošnega in zanesljivega modela še ni, je postopek po teoriji homogenega dvofaznega toka uporabna rešitev. Pričujoči model je omejen na izpeljavo izraza za fazni premeni R in R ([9] in [10]), z dodanim difuzivnim členom v prenosni enačbi (17). V kapljevinskem toku brez zdrsa na medfazni površini in brez upoštevanja vpliva plinov v mehurčku, dinamiko mehurčka podamo s spremenjeno Rayleigh-Plessetovo enačbo (14): R d2RB + 3(dRB B dt 2 2{ dt ( rf ) +V(rVfyv(rrvf)+Re-Rc (17), where Re and Rc represent the source terms or the evaporation/condensation rate. They are functions of the flow parameters (flow pressures and velocities) and the fluid properties (density, evaporation pressure, surface tension). G represents the kinematic diffusivity of the variable f. Equation (17) derives from the theory of homogenous two-phase flow, which represents a good simplifying approach for the following reasons: - In engineering practice cavitation appears in the low pressure areas, where the relative flow velocities are high, and therefore no slip between the liquid and gas phases exists. - The vapour phase usually has the form of small bubbles. In the bubble proximity we have to use a suitable models for bubble radius and drag determination, which causes the problem. Therefore, the homogenous two-phase model is the suitable solution. The presented homogenous two-phase flow model is limited to the derivation of the source terms (R and R ) ([9] and [10]), with the addition of the diffusion term (17). Bubble dynamics, without slip on the interfaces and without considering the gas effects in the bubble, can be defined with the modified Rayleigh-Plesset equation (14), written in the following form: rl 4uL R 2s rLRB (18). 108 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 Da bi izračunali neto stopnjo uparjanja R=R -R , To calculate the net evaporation rate R=R -R the zapišimo spremenjen zakon ohranitve mase za: e c mass continuity equation can be written as: e c - fazo pare: - vapour phase: R=DDt[arV] + arVV-v (19), - fazo kapljevine: - liquid phase: -R= DDt[( 1-a ) rL] + ( 1-a ) rLV.v (20), - homogeno mešanico: - homogenous mixture: 0 D r + rV-v (21). S kombinacijo enačb izpeljemo odvisnost: By combining the equations, the following depend- ency can be written: - (( 1-a ) rL+arV ) V-v= D----------- (22) ((1-a) rl +arV) oziroma: and: D ----(( 1-a ) rL +arV Rr = -r^___ ) (23). Dt (1-a)rL+arV Če ni spremembe gostote pare, niti If the vapour and liquid densities are con- spremembe gostote kapljevine in ob predpostavki, stant, and with the assumption that the volume phase da je sprememba prostorninskega deleža pare zaradi of vapour caused by convection can be neglected, konvekcije zanemarljiva, lahko zapišemo totalni the total derivation of the mixture density can be odvod gostote mešanice kot: written as: Drr- ( rL-rV ) da (24). Prostorninski delež plinaste faze je definiran The volume fraction of the gas phase with s številsko gostoto mehurčkov n in polmerom the bubble number density n and with the bubble mehurčka RB , in v kombinaciji z (21) izpeljemo: radius RB . Combining this with Equation (21), the following equation results: Dr =- ( rL-rV )( 4p n 13 3a )dRB (25). S kombinacijo Rayleigh-Plessetove enačbe With the combination of the Rayleigh-Plesset brez viskoznega člena in brez člena površinskih equation without the viscous term and the term of napetosti ter z uporabo zgoraj zapisanih zvez lahko, surface tension and (V ¦ v) a = 0, the following de- ob upoštevanju (V ¦ v)a = 0 in predpostavki pendency can be written (with the assumption Drv /Dt=0, izpeljemo odvisnosti: Drv/Dt=0), 12 (26). rVrL( 4p n 13 ( 3a2 3 če zanemarimo drugi odvod polmera f we neglect the second derivative of the mehurčka v zgornjem izrazu (pomemben samo v bubble radius in the equation above (important only začetni pospešeni fazi rasti mehurčka), velja in the first bubble-growing phase), the simplified Eksperimentalno numerična analiza - Experimental and Numerical Analyses 109 2 f pB-p I 2 d2RB rl J 3 dt Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 poenostavljena prenosna enačba (16) za paro v naslednji obliki: transport equation (16) for the vapour phase in the following form can be written as: 8t rf)+V-(rfv)=v(rTVf) + r 2\pbp 3 rl 1/2 (27). Drugi člen na desni strani zgornje enačbe pomeni uparjalni (izvirni) člen. Čeprav razmere ob uparjanju niso identične tistim ob kondenzaciji, poenostavljeno vzamemo, da je enačba enaka za oba postopka fazne premene. V zgornji enačbi so vsi členi, razen n, znane stalnice ali odvisne spremenljivke. Iz tega razloga bomo zapisali enačbo za izvirna člena v spremenjeni obliki [9]: Re=Cechrlrv s The second term on the right-hand side of the above equation is the vapour source term. It is known that the process during water evaporation is not equal to the water condensation process, but it can be simpli-fied, i.e. Equation (27) is the same for both processes. In the above Equation (27) all the terms ex-cept n are constants or dependent variables. This is the reason why the source terms are written in changed form [9]: 2pv V Rc =Cc ch s rLr kjer sta C in C stalnici, Vh pa določa lokalno relativno hitrost med kapljevino in paro in jo lahko ocenimo z enačbo Vh=4k . Zgornje enačbe temeljijo na naslednjih predpostavkah: - v mehurčastem toku je stopnja fazne spremembe sorazmerna Vc 2h, za večino praktičnih dvofaznih tokov pa lahko predpostavimo linearno odvisnost od hitrosti; - relativna hitrost med kapljevito in plinasto fazo je reda velikosti 1 do 10% povprečne hitrosti toka, kar v večini turbulentnih tokov ustreza redu velikosti lokalnih turbulentnih sprememb. Tako lahko za prvi približek zapišemo V = *Jk . 3 NUMERIČNA ANALIZA TOKA OKOLI LOPATE Numerično simulacijo obtekanje lopate smo s predstavljenim matematičnim modelom izvedli s programom CFX 5.6. Programski paket temelji na metodi končnih prostornin (MKP). Predstavljen matematično-fizikalni model smo vključili v običajni sistem Navier Stokesovih enačb v obliki dodatne konvektivno difuzivne prenosne enačbe za masni delež pare. Računsko območje je 3D prostornina velikosti 1200x1000x5 mm, z vseh strani zaprt s ploskvami. Območje računanja je razdeljeno na tetraedre. Take tetraedre imenujemo pretočni elementi, njihova oglišča pa so vozlišča. Opisane mreže imenujemo nestruktuirane. Lokacija vozlišča v 1/2 3rL 2 p-p 3rL (1"f) 1/2 f (28) (29), where C and C are constants, and Vh defines the local relative velocity between the liquid and the va-pour, and it can be assumed by Vh = 4k . Equations (28) and (29) are based on the following assumptions: - The rate of phase change is proportional to Vc 2h in bubbly flow, but the linear dependency on velocity can be predicted for conventional two-phase flows. - The relative velocity between the liquid and gas phases is in the range from 1 up to 10% of the average flow velocity, which is suitable for the local turbulent fluctuations in most turbulent flows. For the first approximation, V = *Jk can be used. 3 ANALYSES OF THE NUMERICAL FLOW AROUND THE BLADE The numerical simulation for the flow around the blade profile using the presented mathematical model was performed by CFX 5.6. The numerical code is based on the finite-volumes method (FVM). The presented mathematical/physical model is included in the conventional system of the Navier Stokes equation in the form of an additional convective-diffusion transport equation for the vapour mass fraction. The calculation area is a 3D volume with a size of 1200 x 1000 x 5 mm, closed from all sides. The area is divided into tetrahedrons. Such a tetrahedron is called a flow element, and its edges represent the calculating nodes. This type of mesh is an unstructured mesh. The location of the nodes is 110 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 prostoru je določena s kartezičnimi koordinatami x, y in z. V vsakem vozlu pretočnega elementa dobimo z numeričnim izračunom vrednosti odvisnih spremenljivk: tlaka, hitrosti, turbulentne kinetične energije, raztrosa turbulentne kinetične energije in masnega deleža pare. Testni izračun brez kavitacije smo izvedli pri različnih gostotah računske mreže okoli lopatičnega profila NACA 4418, prikazanega na sliki 2a Iz diagrama na sliki 2b je razviden potek tlaka pod in nad lopatico zaštiri različna števila vozlišč (preglednica 1). determined by the Cartesian coordinates x, y, and z. For each node of the flow element the dependent variables pressure, velocity, turbulent kinetic energy, dissipation of the turbulent kinetic energy and the vapour fraction are calculated. A test calculations without cavitation was performed for different calculation mesh densities and shown in Figure 2.a. The pressure distributions along the suction side and the pressure side of the blade for different mesh-refinement factors (Table 1) is evident in the diagram (Figure 2.b). a) b) 3,0E+05 2,0E+05 i 1,0E+05 0,0E+00 -1,0E+05 -2,0E+05 Algoritem izračuna: 1. Reševanje osnovnega sistema NS enačb za turbulentni tok. 2. Določitev R in R po enačbah (27) in (28). 3. Rešitev enačbe (16). Solver algorithm: 1. Solving the fundamental system of NS equations for turbulent flow 2. Determining the source terms R and R 3. Solving the additional transport equation. x [%] Sl. 2. Lopatični profil (a) in potek tlaka (b) vzdolž lopatice pri različnih gostotah računske mreže za režim a=0°, Re=1-106, c=3 Fig. 2. Blade profile (a) and pressure distribution (b) around the blade profile for different mesh-refinement factors for regime a=0°, Re=1-106, a=3 Preglednica 1. Podatki o računski mreži Table 1. Mesh data Zap. štev. Case no. Št. vozlišč Number of nodes 2846 NACA 4418 4394 3 13372 39824 4 29511 93061 Št. elementov Number of elements 8362 12608 Eksperimentalno numerična analiza - Experimental and Numerical Analyses 111 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 Sl. 3. Gostota računske mreže v neposredni bližini profila NACA 4418 Fig. 3. Mesh density near blade profile surface (profile NACA 4418) Iz diagrama je razvidno, da med primerom 3 in rezultati za gostoto računske mreže, označeno z zaporedno številko 4, ni bistvenih razlik. Z namenom hitrejšega računanja smo zato v nadaljevanju računali z mrežo s približno 40.000 elementi (zap. št. 3). Na sliki 3 sta prikazani računski mreži št.1 in št. 3 v neposredni bližini lopatičnega profila. Informacije o toku na mejnih ploskvah računskega območja smo definirali z naslednjimi robnimi pogoji: - stenski pogoj, - vstopni pogoj, - izstopni pogoj, - simetrični robni pogoj. Na steni ni bilo zdrsa, na vstopu smo predpisovali masni pretok, na izstopu pa termodinamični tlak. 3.1 Računski parametri in konvergenca Za ustaljeni izračun je bil uporabljen običajni turbulentni model k-e. Za diskretizacijo konvektivnega člena prenosne enačbe je bila uporabljena shema “upwind”, analiza pa je bila izvedena pri različnih vrednostih kinematične difuzivnosti G. Difuzivnost pomeni hitrost razširjanja skalarne veličine v primeru brez konvekcije in je v splošnem odvisna od lastnosti nosilne (kapljevite) faze in lastnosti faze, za katero rešujemo dodatno prenosno enačbo. Iz tega razloga smo numerično simulacijo ponovili pri različnih vrednostih G. Ciljni ostanek je bil r = 105, za odvisne spremenljivke (tlak, hitrost, turbulentna kinetična energija, raztros turbulentne kinetične energije in masni delež pare). From the diagram (Figure 2) it is evident that the calculation results do not deviate remarkably. From this fact it can be concluded that the density does not significantly affect the calculation results in this case. For this reason a mesh with approximately 40.000 elements (mesh type 3, table 1) was used. Figure 3 shows the mesh type 1 and 3 for the blade profile surface area. The information about the boundary layers of the calculating domain is de-fined by the following boundary conditions: - wall condition, - inlet condition, - outlet condition, - symmetry condition. No flow slip on the wall is considered. The mass flow rate is defined at the inlet. The thermo-dynamic pressure is defined at the outlet. 3.1 Calculation parameters and convergence The conventional k - e was used for the stationary calculation. For the convective term discretisation the “upwind” scheme was used. The calculation was performed for different values of the kinematic diffusion G. The diffusion represents the velocity of the scalar quantity propagation. In the case of no-convection, this in general depends on the liquid properties and/or on the properties of the phase that is currently calculated by the transport equation. For this reason the calculation is repeated for different values of G. The target residual was r = 10-5, for de-pendent values (pressure, velocity, turbulent kinetic energy, dissipation of the turbulent kinetic energy, and the vapour-phase mass fraction). 112 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 4 EKSPERIMENTALNA RAZISKAVA TOKA OKOLI LOPATE Eksperimentalna raziskava je bila izvedena v tunelu za preizkušanje Kaplanovih turbin v Turboinštitutu v Ljubljani [11]. Najpomembnejše notranje mere tunela so: B=150 mm in H=400 mm, kjer je ravni del tunela dolg L *=21L (sl. 4). Dolžina polirane lopate iz brona je L=150 mm in vpliva na vrednost Reynoldsovega števila: 4 EXPERIMENTAL INVESTIGATIONS OF FLOW AROUND THE BLADE PROFILE The experimental analysis was performed in the tunnel for the Kaplan turbine testing at the Turboinstitute in Ljubljana [11]. The general dimensions of the tunnel are as follows: width, B= 150 mm; and height, H=400 mm. The straight tunnel part is long, up to 21 times the blade profile lengths L*=21L, as shown in Figure 4. The length of the profile with a polished surface is made from brass L=150 mm. The Reynolds number is determined by: v Re = v-L u (30). Na vstopu v tunel je izveden konfuzor, na izstopu pa difuzor za upočasnitev toka in rekuperacijo kinetične energije. Tunel je opremljen z mehanizmom za spreminjanje nagibnega kota lopate ter oknom iz poliakrilnega stekla, za vizualno opazovanje. Potek in način meritev je podrobneje opisan v [11]. Spreminjanje termodinamičnega tlaka in s tem kavitacijskega koeficienta: For the tunnel intake an entrance tube was used, similar to the diffuser at the tunnel exit. At the exit part of the tunnel (in the diffuser) the flow ve-locities decrease; this allows recuperation. The tunnel was equipped with a mechanism to change the blade attack angle and a transparent Plexiglas window for the flow visualisation. The thermodynamic pressure variation was achieved with a vacuum pump. In this way the variation of the cavitation coefficient: s p- pV 12 rv2 (31) je omogočeno s sistemom z vakuumsko črpalko. V zgornji enačbi je p statični tlak, ki ga merimo v tunelu pred lopato na razdalji 350 mm pred vrtiščem lopate. is achieved. p in the Equation (31) presents the statical pressure measured in the cavitation tunnel 350 mm in front the blade rotating point. nffy.y"'"^: ^i^ Sl. 4. Kavitacijski tunel [11] Fig. 4. Cavitation tunnel [11] Eksperimentalno numerična analiza - Experimental and Numerical Analyses 113 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 Nagibni kot lopate a je mogoče spreminjati The blade inclination angle a can be changed v obe smeri, pri čemer pomeni kot a =0° vodoravno on both sides (positive and negative). The zero value a lego tetive profila, pozitivni koti pa dvigovanje =0° corresponds to the horizontal blade chord position vstopnega roba lopatice. and positive angles mean lifting of the leading edge. 5 PRIMERJAVA REZULTATOV 5 COMPARISON OF THE RESULTS Na slikah 5 do 8 je podana primerjava In Figures 5 to 8 is a results comparison for rezultatov za tri različne tokovne režime, definirane z three different flow regimes, defined by the non-di- brezrazsežnimašteviloma s in Re. mensional numbers s and Re. Na sliki 5 so razvidni zametki pare na sliki, ki In Figure 5 the starting cavitation is evident from prikazuje rezultate numerične simulacije in na the both results (experimental (b), and numerical (a)). The fotografiji preizkusa. Območje, kjer se pojavi para, je area where the cavitation occurs is larger in the case of the v primeru numerične simulacije nekoliko večje. numerical simulation than with the experimental result Na sliki 6 je razvidno, da pride do pojava In Figure 6 it is evident that the cavitation kavitacije tudi v primeru tokovnega režima phenomenon also occurs in the case of the flow de- definiranega z a=16°, Re=1106, s=2. Izračunan fined by, a=16°, Re=1106, s=2. The calculated cavita- Sl. 5. Primerjava eksperimentalnih (a) in numeričnih (b) rezultatov za spremenljivko f za režim a=10°, Re=8-105, s=3 Fig. 5. Comparison of the experimental (a) and numerical (b) results for value f in regime, a=10°, Re=8-105, s=3 Sl. 6. Primerjava eksperimentalnih (a) in numeričnih (b) rezultatov za spremenljivko f za režim a=16°, Re=1106, s=2 Fig. 6. Comparison of experimental (a) and numerical (b) results for value f in regime, a=16°, Re=1106, s=2 114 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 Sl. 7. Primerjava eksperimentalnih (a) in numeričnih (b) rezultatov za spremenljivko f za režim a=24°, Re=1106, s=1 Fig. 7. Comparison of experimental (a) and numerical (b) results for value f in regime, a=24°, Re=1106, s=1 Sl. 8. Primerjava numeričnih rezultatov za spremenljivko f pri različnih koeficientih difuzivnosti Fig. 8. Numerical result comparison for the value f with different diffusion coefficients kavitacijski oblak je nekoliko manjši od kavitacijskega oblaka, prikazanega na fotografiji, verjetno zaradi vpliva difuzivnosti G. Na sliki 7 je prikazana primerjava rezultatov za režim, pri katerem se je pri preizkusu kavitacijski oblak zelo razširil v tokovno brazdo. Iz primerjave je razvidno neujemanje oblike kavitacijskega oblaka. Razlika je verjetno posledica premajhne difuzivnosti, zaradi česar smo povečali vrednost na G=103. Iz rezultatov na sliki 8 je razviden vpliv difuzivnosti G na prenos masnega deleža pare f za primer a=24°, Re=1-106, s=1. Sklenemo lahko, da je ujemanje z rezultati preizkusa na sliki 7 v primeru večjega prispevka difuzivnega člena v enačbi (27) boljše. tion cloud is smaller than the cavitation cloud in the experiment. This is probably the effect of diffusion G. In Figure 7 is the comparison of the results in the operating regime where the cavitation cloud is propagated into the blade wake. A little disagreement of cavitation swirl form is evident from this comparison. The difference is probably caused by a too small diffusion coefficient. This is the reason why the diffusion coefficient is increased to G=10 3. From the results given in Figure 8 we can see the influence of diffusion on the mass flow part of the vapour phase f can be seen for the case, a=24°, Re=1106, s=1. It can be concluded that the agreement between the experimental and numerical results, given in Figure 8, is better when the diffusion part in Equation (27) is considered. Eksperimentalno numerična analiza - Experimental and Numerical Analyses 115 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 6 SKLEPI 6 CONCLUSIONS Na temelju slik 6 do 8 lahko zapišemo, da predstavljeni homogeni dvofazni model zadovoljivo popiše kavitacijske razmere pri obtekanju lopatičnega profila pri različnih natočnih kotih. Primeren je za napovedovanje pojava kavitacije in v razširjeni (konvektivno difuzivni) obliki zadovoljivo popiše tudi obliko kavitacijskega oblaka. V prenosni enačbi masnega deleža (suhosti) pare difuzivničlen vpliva na obliko kavitacijskega oblaka, saj z njim skušamo zajeti razširjanje kavitacijskega polja zaradi razlik v koncentracijah v mehurčku. Uporabljen običajni dvoenačbni turbulentni model k-e, je primeren za popis turbulentno kavitacijskih razmer pri obtekanju lopatičnih profilov. Predpostavljena nespremenljiva gostota pare in kapljevine ne pomeni prevelike napake. Konvekcija na medfazni površini ne vpliva na deleža pare. Predstavljeni homogeni kavitacijski model je dobro orodje za napovedovanje kavitacije v turbinskih strojih, saj so predstavljeni rezultati numerične simulacije časovno povprečene vrednosti, fotografije eksperimenta pa prikazujejo trenutno strukturo kavitacijskega oblaka, zaradi česar je primerjava omejena. Zahvala Based on Figures 6 to 8 we can conclude that the presented homogenous two-phase model gives good results for the flow pattern around the blade profile at different flow attack angles. The presented mathematical model is suit-able for the prediction of the cavitation cloud around a hydrofoil (shape, position and dimensions). The mass part of the vapour phase in the transportation equation influences the shape of the cavitation cloud. With this approach we try to con-sider the concentration change in the bubbles. The conventional two-equation turbulent model k - e is suitable for a determination of the turbulent cavitation flow conditions when the flow around blade profiles is calculated. The assumed constant density of the va-pour and the liquid do not cause important mis-takes. The convection at the bubble interface does not influence the vapour phase fraction. The presented homogenous cavitation model presents a good tool for cavitation prediction turbo machines, since averaged numerical simulation results agree to the instant cavitation cloud struc-ture at the photos. Acknowledgement Avtorji se zahvaljujejo osebju podjetja The authors would like to thank the per- Turboinštitut iz Ljubljane za vse rezultate meritev. sonnel of the Turboinštitut, Ljubljana, for all the measurement results. 7 SIMBOLI 7 SYMBOLS širina kavitacijskega tunela B cavitation tunnel width stalnica C constant stalnica C constant stalnica C constant stalnica C constant stalnica C constant f vapour mass fraction gostota prostorninske sile f body force funkcija F function višina kavitacijskega tunela H cavitational tunnel height turbulentna kinetična energija k turbulent kinetic energy dolžina lopatice L blade length številska gostota mehurčkov n bubble number density normala na površino r n, nj surface normal 116 Biluš I. - Škerget L. - Predin A. - Hriberšek M Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 nastanek turbulentne kinetične energije termodinamični tlak tlak v mehurčku uparjalni tlak razdalja od središča mehurčka neto fazna sprememba stopnja kondenzacije stopnja uparjanja polmer parnega mehurčka površina kontrolne prostornine čas temperatura obdajajoče kapljevine temperatura v mehurčku pare radialna hitrost značilna hitrost časovno povprečen vektor hitrosti nadzorna prostornina volumski delež parne faze Kroneckerjeva delta funkcija kinematična difuzivnost raztrosna hitrost turbulentne kinetične energije kinematična viskoznost turbulentna kinematična viskoznost dinamična viskoznost gostota mešanice gostota kapljevite faze gostota plinaste (parne) faze površinska napetost, kavitacijsko število napetostni tenzor napetost na elementu izseka iz krogelne površine tenzor strižnih napetosti P p pB p r R R R RB S t T« TB u V a d ij G e u uT h r rL rV s s ij t t ij turbulent kinetic energy production thermodynamic pressure pressure in the bubble vapour pressure distance from the bubble center net phase change condensation rate evaporisation rate bubble radius surface of control volume time surrounding fluid temperature temperature in bubble radial velocity characteristic velocity time averaged velocity vector control volume vapour volume fraction, Kronecker delta function kinematic difusivity disipation velocity of turbulent kinetic energy kinematic viscosity turbulent kinematic viscosity dynamic viscosity mixture density liquid density gas (vapour) density surface tension, cavitation number stress tensor tension on the bubble surface element shear stress tensor [ 1] [2] [3] [4] [5] [6] [7] [8] 8 LITERATURA 8 REFERENCES Hofmann, M. (2001) Ein Betrag zur Vermeidung des erosiven Potentials kavitierender Strömungen, Dissertation, Technische Universität Darmstadt, Darmstadt. Han, J., A. Alajbegocič (2002) Simulation of multiphase flows in complex geometry using a hybrid method combining the multifluid and the volume of fluid (VOF) approaches, Proceedings of ASME FEDSM‘02, Montreal, Canada, 1-6. Berg, E., M. Volmajer (2003) Validation of a CFD model for coupled simulation of nozzle flow, primary fuel jet break-up and spray formation, Proceedings of ICES03, Salzburg, Austria, 1-10. Sauer, J. (2000) Instationär kavitierende Strömungen - Ein neues Modell, basierend auf front Capturing (VoF) und Blasendynamik, Dissertation, Universität Karlsruhe, Karlsruhe. Lindau, J. W., R.F. Kunz (2002) High Reynolds number, unsteady, multiphase CFD modeling of cavitation flows, Journal of Fluids Engineering, Vol. 124, 607-616. Škerget, L. (1994) Mehanika tekočin, Univerza v Mariboru, Tehniška fakulteta. Brennen, C. E. (1995) Cavitation and bubble dynamics. Oxford University Press, London. Grist, E. (1999) Cavitation and the centrifugal pump. Taylor & Francis, London. Eksperimentalno numerična analiza - Experimental and Numerical Analyses 117 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)2, 103-118 [9] Singhal, A. K., M.M. Athavale (2002) Mathematical basis and validation of the full cavitation model. Journal of Fluids Engineering, Vol. 124, 617-624. [10] Dular, M., B. Širok, B. Stoeffel, B. Bachert (2003) Numerična simulacija kavitacije na osamljenem profilu v kavitacijskem tunelu, Slovensko društvo za mehaniko, Kuhljevi dnevi 2003, Zreče, 93-102. [11] Vujanič, V. (1992) Raziskava toka okoli lopate nameščene v kavitacijskem tunelu, magistrska naloga. Fakulteta za strojništvo, Ljubljana. Naslov avtorjev: dr. Ignacijo Biluš prof.dr. Leopold Škerget prof.dr. Andrej Predin prof.dr. Matjaž Hriberšek Univerza v Mariboru Fakulteta za strojništvo Smetanova ulica 17 2000 Maribor ignacijo.bilus@uni-mb.si leo@uni-mb.si andrej.predin@uni-mb.si matjaz.hribersek@uni-mb.si Authors’ Address: Dr. Ignacijo Biluš Prof. Dr. Leopold Škerget ProfDr. Andrej Predin ProfDr. Matjaž Hriberšek University of Maribor Faculty of Mechanical Eng. Smetanova 17 SI-2000 Maribor, Slovenia ignacijo.bilus@uni-mb.si leo@uni-mb.si andrej.predin@uni-mb.si matjaz.hribersek@uni-mb.si Prejeto: Received: 21.4.2004 Sprejeto: Accepted: 2.12.2004 Odprto za diskusijo: 1 leto Open for discussion: 1 year 118 Biluš I. - Škerget L. - Predin A. - Hriberšek M