SEPTEMBER 2022 www.fe.um.si/jet.html 2 JET JET 3 VOLUME 15 / Issue 2 Revija Journal of Energy Technology (JET) je indeksirana v bazah INSPEC© in Proquest’s Technology Research Database. The Journal of Energy Technology (JET) is indexed and abstracted in database INSPEC© and Proquest’s Technology Research Database. 4 JET JOURNAL OF ENERGY TECHNOLOGY Ustanovitelj / FOUNDER Fakulteta za energetiko, UNIVERZA V MARIBORU / FACULTY OF ENERGY TECHNOLOGY, UNIVERSITY OF MARIBOR Izdajatelj / PUBLISHER Fakulteta za energetiko, UNIVERZA V MARIBORU / FACULTY OF ENERGY TECHNOLOGY, UNIVERSITY OF MARIBOR Glavni in odgovorni urednik / EDITOR-IN-CHIEF Jurij AVSEC Souredniki / CO-EDITORS Bruno CVIKL Miralem HADŽISELIMOVIC Gorazd HREN Zdravko PRAUNSEIS Sebastijan SEME Bojan ŠTUMBERGER Janez USENIK Peter VIRTIC Ivan ŽAGAR Uredniško izdajateljski svet / PUBLISHING & EDITORIAL COUNCIL Dr. Anton BERGANT, Litostroj Power d.d., Slovenia Prof. dr. Marinko BARUKCIC, Josip Juraj Strossmayer University of Osijek, Croatia Prof. dr. Goga CVETKOVSKI, Ss. Cyril and Methodius University in Skopje, Macedonia Prof. dr. Nenad CVETKOVIC, University of Nis, Serbia Prof. ddr. Denis ÐONLAGIC, University of Maribor, Slovenia Doc. dr. Brigita FERCEC, University of Maribor, Slovenia JET 5 Prof. dr. Željko HEDERIC, Josip Juraj Strossmayer University of Osijek, Croatia Prof. dr. Marko JESENIK, University of Maribor, Slovenia Izr. prof. dr. Ivan Aleksander KODELI, Jožef Stefan Institute, Slovenia Prof. dr. Rebeka KOVACIC LUKMAN, University of Maribor, Slovenia Prof. dr. Milan MARCIC, University of Maribor, Slovenia Prof. dr. Igor MEDVED, Slovak University of Technology in Bratislava, Slovakia Izr. prof. dr. Matej MENCINGER, University of Maribor, Slovenia Prof. dr. Greg NATERER, Memorial University of Newfoundland, Canada Prof. dr. Enrico NOBILE, University of Trieste, Italia Prof. dr. Urška LAVRENCIC ŠTANGAR, University of Ljubljana, Slovenia Izr. prof. dr. Luka SNOJ, Jožef Stefan Institute, Slovenia Prof. Simon ŠPACAPAN, University of Maribor, Slovenia Prof. dr. Gorazd ŠTUMBERGER, University of Maribor, Slovenia Prof. dr. Anton TRNIK, Constantine the Philosopher University in Nitra, Slovakia Prof. dr. Zdravko VIRAG, University of Zagreb, Croatia Prof. dr. Mykhailo ZAGIRNYAK, Kremenchuk Mykhailo Ostrohradskyi National University, Ukraine Prof. dr. Marija ŽIVIC, Josip Juraj Strossmayer University of Osijek, Croatia 6 JET Tehnicni urednik / TECHNICAL EDITOR Sonja Novak Tehnicna podpora / TECHNICAL SUPPORT Tamara BRECKO BOGOVCIC Izhajanje revije / PUBLISHING Revija izhaja štirikrat letno v nakladi 100 izvodov. Clanki so dostopni na spletni strani revije – https://www.fe.um.si/jet.html / The journal is published four times a year. Articles are available at the journal’s home page – www.fe.um.si/en/jet.html. Cena posameznega izvoda revije(brez DDV) / Price per issue (VAT not included in price): 50,00 EUR Informacije o narocninah: https://www.fe.um.si/naro%C4%8Dnine.html / Subscription informa­tion: https://www.fe.um.si/en/subscriptions.html Lektoriranje / LANGUAGE EDITING TAIA INT d.o.o. Oblikovanje in tisk / DESIGN AND PRINT Tiskarna Koštomaj d.o.o. Naslovna fotografija / COVER PHOTOGRAPH Jurij AVSEC Oblikovanje znaka revije / JOURNAL AND LOGO DESIGN Andrej PREDIN Ustanovni urednik / FOUNDING EDITOR Andrej PREDIN Izdajanje revije JET financno podpira Javna agencija za raziskovalno dejavnost Republike Slovenije iz sredstev državnega proracuna iz naslova razpisa za sofinanciranje domacih znanstvenih periodicnih publikacij / The Journal of Energy Technology is co-financed by the Slovenian Research Agency. JET 7 Spoštovani bralci revije Journal of energy technology (JET) Podnebne spremembe, ki smo jim prica po vsem svetu, so obcutno zvišale temperaturo ozracja. V Sloveniji se je povprecna temperatura v zadnjih šestdesetih letih povišala za 2,4 stopinj Celzija. Pos­ledica višjih temperatur je tudi vrocinski stres, ki ga obcutijo tako ljudje kot tudi živali in rastline. Nekatere druge posledice globalnega segrevanja so tudi taljenje ledenikov, višanje gladine morja, pre­razporejanje padavin, zmanjšanje pridelka na kmetijskih površinah, ogrožen obstoj mnogih rastlinskih in živalskih vrst … Prav zaradi tega je na energetskem podrocju v Sloveniji in drugod po svetu treba preiti na obnovljive vire v mnogo vecji meri kot do sedaj. Poleg pogostejše uporabe obnovljivih virov bi bilo treba v en-ergetske sisteme veliko bolj izdatno vkljucevati alternativne energetske tehnologije, kot so na prim­er vodikove in metanolove tehnologije, toplotne crpalke in alternativne izvedbe vodnih ter vetrnih turbin. Pricujoca številka revije Journal of Energy Technology v nekaterih predstavljenih clankih obravnava prav to tematiko. Jurij AVSEC odgovorni urednik revije JET 8 JET Dear Readers of the Journal of Energy Technology (JET) The climate changes taking place in the world are having a significant effect on the increase in the temperature of the atmosphere. In Slovenia, the temperature has increased by 2.4 degrees celcius in the last sixty years alone. The result of the temperature increase is reflected in the heat stress felt by not just people, but also animals and plants. Apart from heat stress, the consequences of global warming include the melting of glaciers, rising sea levels, redistribution of precipitation, reduction of crops on agricultural land, and threats to the existence of many plant and animal species… It is for this reason that it is vital, both in Slovenia and across the world, to switch to the use of re­newable sources in the energy field on a much greater scale than ever before. Aside from the use of renewable resources, alternative energy technologies such as hydrogen technologies, methanol tech­nologies, heat pumps, and alternative versions of water and wind turbines must be included in energy systems in much greater quantities. This issue of the Journal of Energy Technology discusses this topic in some of the presented articles. Jurij AVSEC Editor-in-chief of JET JET 9 Table of Contents / Kazalo Hydrogen production using a thermochemical cycle Proizvodnja vodika s termokemicnim procesom Jurij Avsec, Urška Novosel, Dušan Strušnik. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Analytical and numerical analysis of trapped air pocket dynamic response due to pressure change in liquid-filled pipelines Analiticna in numericna analiza dinamicnega odziva zracnega mehurja na tlacno motnjo v cevovodih, napolnjenih z vodo Jošt Pekolj, Anton Bergant, Matjaž Perpar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Overview of GNSS systems and their operation Pregled GNSS sistemov in njihovo delovanje Boštjan Krošelj, Dalibor Igrec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Prediction of cavitation and particle erosion in a radial divergent test section Napoved kavitacijske erozije in erozije delcev v radialno divergentni testni sekciji Luka Kevorkijan, Luka Lešnik, Ignacijo Biluš . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Analysis of the influence of parameters when charging and discharging a capacitor using differential equations Analiza vpliva parametrov pri polnjenju in praznjenju kondenzatorja z uporabo diferencialnih enacb Matic Krašovic, Peter Virtic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Instructions for authors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 10 JET JET Volume 15 (2022) p.p. 11-20 Issue 2, 2022 Type of article: 1.01 www.fe.um.si/si/jet.html HYDROGEN PRODUCTION USING A THERMOCHEMICAL CYCLE PROIZVODNJA VODIKA S TERMOKEMICNIM PROCESOM Jurij Avsec1R, Urška Novosel1, Dušan Strušnik2 Keywords: Thermodynamics, energy, hydrogen production, solid phase, fluid phase Abstract Sustainable methods of clean fuel production are needed throughout the world due to depleting oil reserves and the need to reduce carbon dioxide emissions. The technology based on fuel cells for electricity production or the transport sector has already been developed. However, a key missing element is a large-scale method of hydrogen production. The copper-chlorine (CuCI) combined thermochemical cycle is a promising thermochemical cycle that can produce large amounts of cheap hydrogen. A particularly promising part of this process is its use in combi­nation with nuclear or thermal power plants. This paper focuses on a CuCl cycle and describes the models used to calculate thermodynamic and transport properties. This paper discusses the mathematical model for computing the thermodynamic properties for pure HCl and CuCl2. The mathematical model developed for the solid phase takes into account vibrations of atoms in molecules and intermolecular forces. This mathematical model can be used for the calculation of the thermodynamic properties of polyatomic crystals on the basis of the Einstein and Debye equations. The authors of this paper developed the model in the low temperature and high tem­perature region. All the analytical data have been compared with some experimental results and show a relatively good match. For the solid phase, the authors developed a model to calculate thermal conductivity based on electron and phonon contributions. R Corresponding author: Prof. Jurij Avsec, Ph.D., Tel: +386 (0)76 202 217, Fax: +386-2-620-2222, Mailing address: Hocevarjev trg 1, 8270 Krško, Slovenia. E-mail: jurij.avsec@um.si 1 Universitiy of Maribor, Faculty of Energy Technology, Hocevarjev trg 1, Krško 2 Energetika Ljubljana d.o.o., TE-TOL Unit, Toplarniška ulica 19, SI-1000 Ljubljana, Slovenia Povzetek Zaradi izcrpavanja zalog nafte in potrebe po zmanjšanju emisij ogljikovega dioksida so trajnostne metode proizvodnje cistega goriva potrebne v vseh državah sveta. Tehnologija na osnovi gorivnih celic je že razvita, vendar pa je kljucni manjkajoci element obsežna proizvodnja vodika. Kombini­rani termokemicni cikel baker-klor (CuCl) predstavlja obetaven termokemicni cikel za proizvodnjo poceni vodika v velikih kolicinah, kar je proces, ki je še posebej zanimiv v kombinaciji z jedrskimi elektrarnami ali termoelektrarnami. Ta clanek se osredotoca na cikel baker-klor (CuCl) in opisuje modele, kako izracunati termodinamicne in transportne lastnosti, pri cemer obravnava matema­ticni model za izracun termodinamicnih lastnosti za cisti HCl in CuCl2. Razvit matematicni model za trdno fazo upošteva vibracije atomov v molekulah in medmolekul­ske sile. Ta matematicni model, ki se lahko uporablja za izracun termodinamicnih lastnosti polia­tomskih kristalov na podlagi Einsteinovih in Debyejevih enacb, smo razvili v nizkotemperaturnem in visokotemperaturnem obmocju. Analiticne podatke smo primerjali z nekaterimi eksperimen­talnimi rezultati in kažejo dobro ujemanje. Za trdno fazo smo razvili model za izracun toplotne prevodnosti na podlagi prispevkov elektronov in fononov. 1 INTRODUCTION The ecological problems throughout the world are becoming significantly bigger by the year. One of the possible ways of reducing greenhouse gas emissions is the introduction of hydrogen technologies. However, the biggest problem is the production of large enough amounts of hydro­gen. Hydrogen can be obtained from hydrocarbons, but it can also be obtained from water. This article deals with the extraction of hydrogen from water, which, from an ecological perspective, is a much more acceptable process than the extraction of hydrogen from coal, oil and natural gas. Water decomposes naturally into hydrogen and oxygen at 2,500°C. Therefore, in recent decades, thermochemical processes in laboratories have improved and electrolysis processes are also be­ing improved. While the production of hydrogen in thermochemical processes mainly requires thermal energy, electricity is required for the electrolysis process. With the introduction of the new generation of nuclear reactors, where temperatures will be much higher than in existing reactors, the cogeneration of electricity and hydrogen [1-4] in connection with a nuclear power plant and thermochemical processes for hydrogen production can be a promising technology in the future. 2 THERMOCHEMICAL CYCLE FOR HYDROGEN PRODUCTION IN COMBINATION WITH A PWR NUCLEAR POWER PLANT Rather than deriving hydrogen from fossil fuels, a promising alternative is the thermochemical decomposition of water [1]. Electrolysis is a proven commercial technology that separates water into hydrogen and oxygen using electricity. Net electrolysis efficiencies (including both electricity and hydrogen generation) are ty pically about 32%. In contrast, thermochemical cycles to pro­duce hydrogen promise heat-to-hydrogen efficiencies of up to about 50%. Through thermochem­ical cycles, hydrogen is produced by chemical processes and heat supply at significantly lower temperatures than natural decomposition into hydrogen and oxygen. The temperatures usually necessary for thermochemical decomposition range between 750°C and 1,000°C (Table 1). Table 1: A short description of the most tested thermochemical cycles [1-9] Name of the cycle Maximum temperature (°C) Use efficiency (%) Advantages Reference Sulphur-iodine 823-900 42-51 Use efficiency beyond 60% is also planned [1] UT-3 Calcium-bromine 750 40-50 Lower maximum temperature [4] Vanadium-chlorine 925 40.5-42.5 [4] Hybrid Copper-chlorine cycle 550 46 Low maximum temperature [4] Hybrid copper-sulphur 827 68-73 High use efficiency [1,4] The aforementioned requirement for a relatively high temperature prevents the application of thermochemical cycles for conventional, water-cooled nuclear reactors. However, the use of high temperature gas-cooled reactors or fourth-generation reactors would be possible. The begin­nings of research into the thermochemical cycles of hydrogen production from water date back to around 1960. The key advantage of thermochemical cycles over others is the high use efficien­cy, which, in some cycles, goes beyond 50%. The next great advantage over other cycles is the circular process, where the environmental strain caused by CO2 is no longer a problem. In scientific literature, there are more than 200 thermochemical cycles [2], the majority of which have never been tested as pilot projects. Most of these cycles are being developed in the USA, Japan, Canada and France. The following table displays the most well known thermochemical cycles: of particular interest are the thermochemical hybrid cycles, which, apart from waste heat, require electrical energy from a nuclear power plant to operate. This article examines the thermophysical properties of a specific cycle called the copper-chlorine (CuCl) cycle, with particular relevance to nuclear-produced hydrogen. A conceptual schematic of the CuCl cycle is shown in literature [2,3]. In this cycle, water is decomposed into hydrogen and oxygen through intermediate CuCl compounds. Nuclear-based ‚water splitting‘ requires an in­termediate heat exchanger between a nuclear reactor and hydrogen plant, which transfers heat from the reactor coolant to the thermochemical cycle. An intermediate loop prevents exposure to radiation from the reactor coolant in the hydrogen plant, and also prevents corrosive fluids in the thermochemical cycle from entering the nuclear plant. Table 2 shows the most important components in the CuCl thermochemical process of hydrogen production. Table 2: Fundamental thermophysical properties for selected CuCl components HCl CuCl CuCl2 Cu2OCl2 CuO Hydrochloric acid Cuprous chloride Cupric chloride Copper oxychloride Cupric oxide Molecular weight [kg/kmol] 36.5 98.9 134.5 214 79.54 Figure 1: Total net efficiency for the production of hydrogen in comparison with electrolysis and a CuCl cycle The production price of hydrogen was calculated in relation to the PWR nuclear power plant. The results are shown in Table 3, namely the price of hydrogen produced through electrolysis or using a CuCl process. The CuCl process was based on the assumption that 1 MWh of electricity costs EUR 35, while in terms of hydrogen production through electrolysis, it was assumed that 79 kWh of power (51% efficiency) is required. The price of thermal energy from a nuclear power plant was calculated from the Rankine process based on the assumption of steam consumption at the low-pressure stage of the turbine. Table 3 and Figure 3 show the prices calculated for hydrogen production, excluding investment costs. Figure 2: CuCl process energy requirements for process [5] Table 3: Calculated prices for 1kg of hydrogen Electrolysis 2.76 EUR/kg H2 CuCl thermochemical process Electrical energy 0.3 EUR Price for heat at 80°C 0.233 EUR Price for heat at 400°C 0.217 EUR Price for heat at 600°C 0.315 EUR Total 1.065 EUR Figure 3: Comparison of hydrogen price THERMODYNAMIC PROPERTIES FOR SOLIDS In the theoretical formulation for solids it will be assumed that each form of motion of energy is independent of the others. Thus, the energy of the system of molecules can be written as a sum of the following individual contributions or decoupled forms of motion [2]: a) vibrational energy of molecules (Evib), due to the relative motion of atoms inside the mole­cules b) potential energy (Epot) of a system of molecules, which occurs due to the attractive or repul­sive intermolecular forces in a system of molecules c) energy of electrons (Eel), which is concentrated in the electrons or the electron shell of an atom or a molecule d) nuclear energy (Enuc), which is concentrated in the atom nucleus The partition function is defined as Z [6-10], which is applied to the system of particles with a certain volume V, temperature T and particle number N. Assuming that the energy spectrum is continuous, together with the other aforementioned assumptions, the canonical partition func­tion for the one-component system can be expressed in the following manner [9]: 1 . E . E . E . .. . .E . ... vib el nuc pot Z . ... exp .. . dp . dp ....dp ... exp .. . dr . dr ..dr (3.1) Nf ..... 12 N ....12N N!h kT kT . B .. B . The second term on the right-hand side of equation 3.1 is called the configurational integral, where f is the number of degrees of freedom of an individual molecule, p is momentum, r is the coordinate, and Evib, Eel, Enuc, and Epot represent the vibrational energy, electron energy and nuclear energy of individual molecules, and the potential energy between two molecules respectively. Similarly, the partition function Z for a multi-component system of indistinguishable molecules can be expressed as follows: 1 . Evib . Eel . Enuc . .. . . Epot .. .. . Z . ... exp .. . dp . dp ...dp . .. exp .. dr . dr ...dr (3.2) 12N 12N Ni!hNi fi .. ..kBT .... ..kBT .. . In equation 3.2, Ni is the number of molecules of the i-th component, while fi is the number of degrees of freedom of the i-th molecule. Using the canonical partition, the partition function Z of the one-component system as a product of partition functions becomes: Z = Zg Zvib Zel Znuc Zconf (3.3) The partition function Z is a product of terms of the ground state (g), translation (trans), vibration (vib), rotation (rot), internal rotation (ir), influence of electron excitation (el), influence of nuclei excitation (nuc) and the influence of the intermolecular potential energy (conf). Utilising the ca­nonical theory for computing the thermodynamic functions of the state leads to: .. lnZ . Pressure p . kT .. , .. V . T . ... k lnZ . . .. . . ln Z .. . V Entropy S . (3.4) T . ... , .. . T . . ... kT T . . . . ln Z T . ln Z V .. . .. ... . Enthalpy H V . ... , V T .. .. where T is temperature and V is the volume of the molecular system. The various derivatives and expressions of the fundamental equations shown in 3.4 have an important physical significance. This paper introduces expressions that are important in terms of energy exchange processes. The various derivatives below also have a practical significance [4-9]: U T H T ˆ...˜° - Heat capacity at constant volume per mole: C , (3.5) = v V TV . 2 .. ... . - Heat capacity at constant pressure per mole: C .. Cv . , (3.6) p . p where M is the molecular mass. The thermodynamic system consists of N particles associated by attractive forces. Atoms in a crystal lattice are not motionless, as they oscillate around their positions of equilibrium. At tem­peratures below the melting point, the motion of atoms is approximately harmonic [2,3]. This as­sembly of atoms has 3N-6 vibrational degrees of freedom. The six vibrational degrees of freedom have been omitted and the number of vibrational degrees of freedom marked with 3N. Through the knowledge of independent harmonic oscillators, the distribution function Z [7] can be derived as follows: ..h . .. 3N . exp .. .. .. 2k T .. B .. Z . (3.7) .. .h .. . 1 . exp .. .. .. . kT . .. B .. where . is the oscillation frequency of the crystal. The term h./k is the Einstein temperature. When comparing the experimental data for simple crystals, a relatively good match with the analytical calculations at higher temperatures is observed, whereas at lower temperatures the discrepancies are higher. This explains why Debye corrected Einstein’s model by taking account of the interactions between a numbers of quantised oscillators. The Debye approximation treats a solid as an isotropic elastic substance. Using the canonical distribution, the partition function may be written as: . /T 33 9 . D . .. D .. TD . ln Z .. N . 3N . ln . 1 . exp .. .... 3N d . (3.8) .3 8T .. T .. . D . .. 1 exp .. 0 In equation 3.8, .D is the Debye temperature. The following equation is ob­tained by expanding the third term in equation 3.8 into a series for a higher temperature range, . 3 21314 1 6 ..... .. .. ...... (3.9) exp ... 1 2 12 720 . Using equations 3.8 and 3.9 leads to the following expression: 34 5 .. 1 .. D . 1 .. D . 1 .. D . 3 .. ........ 9 . . .... . T .. 3 . T . 8 . T . 60 . T .. DD ln Z .. N . 3N . ln ..1 . exp .. .... 3N ..... 9 . (3.10) 8T .. T .. .. D .. 7 . 1 .. . 1 .. . DD .. ... ... .. . . 5040 . T . 272160 . T .. The Debye characteristic temperature was determined by means of the Grüneisen independent constant .: .. . D .CV , (3.11) where C is a constant, dependent on the material. This mathematical model can be used for the calculation of thermodynamic properties of polyatomic crystals. The derivations of the Einstein and Debye equations, outlined previously, apply specifically to monoatomic solids, namely those belonging to the cubic system. However, experiments have shown that the Debye equation also predicts the values of specific heat and other thermophysical properties for certain other mono-atomic solids, such as zinc, which crystallises in the hexagonal system. ANALYSIS AND RESULTS The results below show a comparison of the analytical model obtained using the mathematical model shown in the sections above. Figure 4 shows the results for CuCl2, a comparison between the calculated results and the results obtained using the Shomate equation. In the previous re­sults, the calculations of thermodynamic properties for solids are determined based on the fol­lowing Shomate equation. E Cp . A . BT . CT 2 . DT 3 . T2 , (4.1) H . C dT , S . C . (4.2) pp .. dT T Figure 5 shows the relative deviation for isobaric specific heat for Cu2OCl2 between the results obtained using the analytical model and the experimental results [12]. Figure 6 shows the results obtained using the mathematical model [6, 9] for isobaric specific heat of fluid and the Shomate equation. The results indicate a relatively good match between the mathematic model and the published results. Figure 4: Relative deviation for isochoric specific heat for CuCl2 between analytical data and the Shomate equation Figure 5: Relative deviation for isochoric specific heat between the analytical model and experi­mental results [14] Figure 6: Relative deviation for isobaric specific specific heat for HCl between analytical data and the Shomate equation References [1] D.A.J. Rand., R.M. Dell., Hydrogen energy, 2008, Royal Society of Chemistry, Cambridge. [2] K. Verfonden, Nuclear energy for hydrogen production, Energietechnik, Vol.58, 2007. [3] J. Avsec, U. Novosel, Application of alternative technologies in combination with nu­clear energy. Transactions of FAMENA, ISSN 1333-1124, 2016, Vol.40, spec. issue 1, pp.23-32. [4] J. Avsec, K. Watanabe, An approach to calculating thermodynamic properties of mix­tures including propane, n-butane and isobutene. Int. J. Thermophys., November 2005, Vol.26, No.6, pp.1,769-1,780. [5] J. Avsec., M. Oblak., Thermal vibrational analysis for simply supported beam and clamped beam. J. Sound Vib., Dec. 2007, Vol.308, Iss.3/5, pp.514-525. [6] J. Avsec, M. Marcic, Calculation of elastic modulus and other thermophysical properties for molecular crystals. J. Thermophys. Heat Transfer, July-September 2002, Vol.16, No.3, pp.463-468. [7] Z.P. Liu, Y.G. Li., J.F. Lu: Fluid Phase Equilibria, Vol.173, pp.189-209, 2000. [8] J. Avsec, Calculation of Transport Coefficients of R-32 and R-125 with the methods of Sta­tistical Thermodynamics and Kinetic Theories of Gas, Archives of Thermodynamics, Vol.24, No.3, pp.69-82. [9] J. Avsec, Marcic, M., Influence of Multipolar and Induction Interactions on the Speed of Sound. J. thermophys. heat transfer, October-December 2000, Vol.14, No.4, pp.496-503. [10] T.-H. Chung, M. Ajlan, L.L. Lee, K.E. Starling, Generalized Multiparameter Correlation for Nonpolar and Polar Fluid Transport properties, Ind. Eng. Chem. Fundam., Vol.23, No.1, 1984, pp.8-13. [11] C. Zamfirescu, I. Dincer, G.F. Naterer, Thermophysical properties of copper compounds in copper-chlorine thermochemical water splitting cycles, Int. J. Hyddrogen Energy, Vol.35, 2010, pp.4,389-4,852. [12] T. Parry, Thermodynamics and magnetism of Cu2OCl2, Master thesis, Birgham Young Uni­versity, Provo, Utah, USA, 2008. JET Volume 15 (2022) p.p. 21-30 Issue 2, 2022 Type of article: 1.01 www.fe.um.si/si/jet.html ANALYTICAL AND NUMERICAL ANALYSIS OF TRAPPED AIR POCKET DYNAMIC RESPONSE DUE TO PRESSURE CHANGE IN LIQUID-FILLED PIPELINES ANALITICNA IN NUMERICNAANALIZA DINAMICNEGA ODZIVA ZRACNEGA MEHURJA NA TLACNO MOTNJO V CEVOVODIH, NAPOLNJENIH Z VODO Jošt Pekolj1, Anton Bergant1,2R, Matjaž Perpar1 Keywords: pipelines, air pocket, transient pressures, analytical solution, method of characteris­tics Abstract This paper serves as an overview of the different approaches to modelling the filling and empty­ing of liquid-filled pipelines. Specifically, it compares an analytical and numerical analysis of the start-up of a liquid column in a one-dimensional pipeline with a trapped air pocket. The analytical analysis is based on the premise that there are no pressure waves travelling in the water column. For this reason, it is not suitable for pipeline systems with very small initial trapped air pockets R Corresponding author: Anton Bergant, PhD, Full time: Litostroj Power d.o.o., Litostrojska 50, 1000 Ljubljana, Slo­venia, anton.bergant@litostrojpower.eu; Part time: Faculty of Mechanical Engineering, University of Ljubljana, Ašk­erceva 6, 1000 Ljubljana, Slovenia 1 Faculty of Mechanical Engineering, University of Ljubljana, Aškerceva 6, 1000 Ljubljana, Slovenia 2 Litostroj Power d.o.o., Litostrojska 50, 1000 Ljubljana, Slovenia or very high reservoir pressures. In the pipelines with parameters, as shown in this paper, the results of the analytical analysis correlate sufficiently with those of the numerical analysis. Povzetek Prispevek predstavlja pregled dveh razlicnih pristopov k modeliranju polnjenja in praznjenja cevovodov s kapljevino. Primerja analiticno in numericno analizo polnjenja v enodimenzijskem cevovodu z ujetim zracnim mehurjem na enem koncu. Analiticna analiza temelji na predpostav­ki, da po vodnem stolpcu ne potujejo tlacni valovi. Zaradi tega ta pristop ni primeren za cevne sisteme z zelo majhnimi zacetnimi zracnimi mehurji ali zelo visokimi tlaki rezervoarja. Za predsta­vljene primere v tem delu se analiticni in numericni rezultati ujemajo. INTRODUCTION During filling or emptying of liquid-filled pipelines an air pocket, or pockets, of various sizes can get trapped at distinct positions (high point, pipe end). The presence of an air pocket influences the system’s response to a pressure disturbance. Pressure in an air pocket could increase by multiples of the original pressure change. The response of an air pocket can be simulated using different computational models. In this paper, the authors examine the one-dimensional analyt­ical and numerical models of air pocket dynamic response. Air pockets become trapped at the downstream end of the gravitational type of pipeline. The basis of the analytical analysis is to describe the behaviour of the system based on differential equations. The method of character­istics will be used for the numerical analysis, as it is convenient for tracking pressure waves in the water filled region of the pipe. The analytical and numerical results will be compared for different parameters of the hydraulic pipeline system during partial filling, as shown in Figure 1. At the leftmost part of the pipeline there is a reservoir with a pressure head Hrez, which is attached to the pipeline with a water column of length L(t) and an air pocket of length Lgas(t). The water and air regions are initially separated by a valve, which is located at x1. After the rapid valve opening, the water column moves towards the air pocket at velocity v(t). The response of an air pocket in similar pipelines has also been studied experimentally in [1], [2] and [3]. When comparing the two types of analysis, the numerical one shows better agreement with the experimental results, because it also includes the interaction of pressure waves and unsteady friction [2]. Figure 1: A hydraulic pipeline system with a trapped air pocket THEORETICAL MODELLING The following subchapters describe the most important equations used in analytical and numer­ical analysis. The equations are valid for the system shown in Figure 1. 2.1 Water region The water column in the water region of the pipe can be modelled as rigid or elastic. The one-di­mensional motion of the rigid water column is governed by the momentum conservation, as shown in equation (2.1): g· ( H -H ( t )) dv a ,rez a ,gas ( t )= + g· sin .- dt L ( t ) (2.1) --.. KvKe+1 ·v ( t )|v (t )|- ·v ( t )|v ( t )|- ·v2 ( t ) ·H ( v ( t )) . - 2 D 2 L ( t ) 2 L ( t ) Similar equations can be seen in [4] and [5]. Equation (2.1) is simplified by neglecting some of the terms; namely, friction factor ., coefficients of valve head losses Kv and reservoir exit losses Ke, while the slope term of the pipe will be set to zero (horizontal pipe). For the elastic water region, it is necessary to consider the motion of pressure waves in the water. This is based on the assumption that the waves are one-dimensional, which is justified, since the diameter of the pipe is substantially smaller than its length. The motion of the water is governed by the continuity and momentum conservation, as shown in equations (2.2) and (2.3) [6]: .H .H a 2 .v +v· -v· sin .+=0, (2.2) .t .x g.x .H .v .v .·v·|v| g· ++v· +=0 . (2.3) .x .t .x 2 D It is assumed that the velocity of pressure waves a is constant. Thereafter, the method of charac­teristics is used to transform equations (2.2) and (2.3) into four ordinary differential equations, which are solved using the finite difference method [6]. They are combined to derive a solution for the piezometric head Hi and volume flow rate Qi at a node in a characteristic grid, as shown in equations (2.4) and (2.5), where the friction head losses are neglected: C pBm + Cm B p Hi = , (2.4) B + B pm Cp-Cm Qi = . (2.5) B +B pm In order to calculate variables Hi and Qi at the boundary, a device-specific equation, or equations, replaces one of the water hammer compatibility equations (the reservoir and air pocket in this case) [2], [6]. 2.2 Air region The air pocket is modelled as a lumped body. Therefore, the presence of pressure waves is not included in the air region in the model. A similar description is illustrated in [6]. During the filling of the pipeline, the water column compresses the air pocket, which undergoes a polytropic pro­cess. The changes in pressure and volume of the air pocket are described using equation (2.6): H( t ) ·Ln ( t )=C A. (2.6) a , gas gas The gas constant CA is calculated by using the known initial conditions (pressure head and air pocket length) in the equation (2.6). From previous experimental investigations, such as [1-3], it was determined that the numerical model most closely matches the experimental data when the polytropic coefficient for the process in air equals 1.4. 2.3 Gas-water interface When using the rigid water column approach, equations (2.1) and (2.6) are coupled for the air and water regions. The result is a linear first-order ordinary differential equation. Equation (2.7) describes the motion of a gas-water interface [5]: 2 +( K +1 H dv . 2 a,rez C A +=2 g-+ sin .. (2.7) dL LD )v ( L ( x-L)n ) L It is assumed that the gas-water interface is always vertical (both in the rigid and elastic water approach). Equation (2.7) can be solved analytically (in this case) or numerically [5]. By using the elastic water column approach, the interface changes with the volume of the air pocket. ANALYTICAL AND NUMERICAL SOLUTIONS 3.1 Analytical solution The analytical method solves the differential equation for the gas-water interface (2.7) [5]. This equation was solved by using an integration factor. The solution is an expression for the speed of the water column in terms of its length (2.8): .....,.... ....*.·........*-...... .*. ·........*+ . . .. ..(....*). ....=......... ·.2....·. . (2.8) +sin............*... ·........* With the help of equation (2.6) for the air region, it was also possible to calculate how the pres­sure in the air pocket changes over time. 3.2 Numerical solution For the numerical analysis, the numerical model based on the method of characteristics [6] was used. The nodes in the numerical grid are fixed and are distributed over the length of the initial water column – from the leftmost node at the reservoir to the rightmost node at the valve. For the interior nodes of the water region characteristic grid, the piezometric head and volume flow rate are calculated using (2.4) and (2.5). At the leftmost node, at the boundary with the reservoir, the modified characteristic equations (2.4) and (2.5) are coupled with the boundary equation of the constant head reservoir [6]. For the node at the gas-water interface, the air region equa­tion coupled with the modified characteristic equations was used. When the water-air interface moves, it is assumed that the variables Hi and Qi at the interface and at the fixed node on the valve are the same. COMPARISONS OF ANALYTICAL AND NUMERICAL RESULTS The analytical and numerical solutions to the equations were calculated using different param­eters of the hydraulic pipeline system, as depicted in Figure 1. The parameters for the different pipelines are shown in Tables 1, 2 and 3. For all of the cases illustrated, the pressure loss coeffi­cient K and slope of the pipe .are equal to zero. Table 1: Pipeline system parameters for Case 1 [4], [5] Hrez [m] Ha, gas0 [m] L0 [m] Lgas0 [m] D [m] n [-] .v [kg/m3] a [m/s] 31 10.3 100 15 0.3 1.4 1000 1000 Table 2: Pipeline system parameters for Case 2 [1] Hrez [m] Ha, gas0 [m] L0 [m] Lgas0 [m] D [m] n [-] .v [kg/m3] a [m/s] 8.15 10.32 5.57 3.25 0.04 1.4 1000 400 Table 3: Pipeline system parameters for Case 3 [1] Hrez [m] Ha, gas0 [m] L0 [m] Lgas0 [m] D [m] n [-] .v [kg/m3] a [m/s] 12.23 10.28 5.57 3.25 0.04 1.4 1000 400 Figures 2 to 4 show the analytical solutions for the velocity of the water column in terms of its length for all the three cases, as illustrated in Tables 1, 2 and 3. The results shown in Figure 2 agree with those presented by Tijsseling et al. [5]. Figure 2: Water column velocity in terms of column length for Case 1 Figure 3: Water column velocity in terms of column length for Case 2 Figure 4: Water column velocity in terms of column length for Case 3 Figures 5 to 7 show the comparison of calculated pressure in the air pocket by using the analyt­ical and numerical methods (Chapter 3) for all three cases in Tables 1, 2 and 3. As can be seen from the figures, the air pocket pressure response for the discussed cases is similar to that of the analytical and numerical analysis. Figure 5: A comparison of analytical and numerical air pocket pressure for Case 1 Figure 6: A comparison of analytical and numerical air pocket pressure for Case 2 Figure 7: A comparison of analytical and numerical air pocket pressure for Case 3 CONCLUSION The analytical and numerical results match quite well, however there are still some differences. The reason for this is that the physical models are based on slightly different approaches, e.g. the rigid and elastic water column, and a different consideration of the water region between the valve and air pocket. The models match for those cases with large initial air pockets and low reservoir pressures. Since the analytical analysis does not include pressure waves in the water column and unsteady friction [2], it is not suitable for cases with small initial air pockets and high reservoir pressures. The numerical method could be improved by providing a better description of the gas-water interface. The part of the water column that extends over the valve will be in­cluded in the elastic water region in future research. References [1] L. Zhou, D. Liu, B. Karney, P. Wang: Phenomenon of White Mist in Pipelines Rapidly Fill­ing with Water with Entrapped Air Pockets, Journal of Hydraulic Engineering, Vol.139, Iss. 10, p.p.1041-1051, 2013 [2] A. Bergant, A. Tijsseling, Y. Kim, U. Karadžic, L. Zhou, M.F. Lambert, A.R. Simpson: Un­ steady Pressures Influenced by Trapped Air Pockets in Water-Filled Pipelines, Strojniški vestnik – Journal of Mechanical Engineering, Vol.64, Iss.6, pp.501-512, 2018 [3] L. Zhou, Y. Cao, B. Karney, A. Bergant, A. S. Tijsseling, D. Liu, P.Wang: Expulsion of Entrapped Air in a Rapidly Filling Horizontal Pipe, Journal of Hydraulic Engineering, Vol.146, Iss.7, p.p.1-16, 2020 [4] E. Cabrera, J. Abreu, R. Perez, A. Vela: Influence of Liquid Length Variation in Hydraulic Transients, Journal of Hydraulic Engineering, Vol.118, Iss.12, p.p.1639-1650, 1992 [5] A. S. Tijsseling, Q. Hou, Z. Bozkus: Analytical expressions for liquid-column velocities in pipelines with entrapped air, In: Proceedings of the ASME 2015 Pressure Vessels & Piping Division Conference, ASME, Paper PVP2015-45184, 2015 [6] F. B. Wylie, V. L. Streeter: Fluid Transients in Systems, Prentice Hall, 1993 Nomenclature A Area a Velocity of pressure waves Bm Proportionality constant Bp Proportionality constant CA Gas constant Cm Proportionality constant Cp Proportionality constant D Internal diameter of the pipe g Gravitational acceleration (g = 9.81 m/s2) H Heaviside step function Ha,gas Air pocket absolute pressure head Ha,gas0 Air pocket initial absolute pressure head Ha,rez Reservoir absolute pressure head Hi Pressure head at a node Hrez Reservoir pressure head K Combined pressure head loss coefficient Ke Pressure head loss coefficient for reservoir exit Kv Pressure head loss coefficient for the valve L Water column length L* Dummy length in integral L0 Initial water column length Lgas Air pocket length Lgas0 Initial air pocket length n Polytropic coefficient Qi Volumetric flow rate at a node t Time v Water column velocity x Space coordinate in direction of the pipeline x1 Location of the valve xL Full length of the pipeline . Friction factor . Slope of the pipe .v Water density Acknowledgments Acknowledgments The authors gratefully acknowledge the support of the Slovenian Research Agency (ARRS) con­ducted through research project L2-1825. JET Volume 15 (2022) p.p. 31-42 Issue 2, 2022 Type of article: 1.01 www.fe.um.si/si/jet.html OVERVIEW OF GNSS SYSTEMS AND THEIR OPERATION PREGLED GNSS SISTEMOV IN NJIHOVO DELOVANJE Boštjan Krošelj1R, Dalibor Igrec1 Keywords: GNSS, GPS, GLONASS, Galileo, BeiDou, navigation, satellite Abstract In this paper, we present the beginnings of navigation and what it was that drove its development forward. We will also present the development of global navigation satellite systems, known as GNSS, as well as all the most important GNSS systems and their development. Finally, the basic operation of such systems and the differences between them will be presented. Povzetek V nalogi bomo predstavili zacetke navigacije in kaj je poganjalo njen zgodovinski razvoj. Ogledali si bomo tudi razvoj globalnih navigacijskih satelitskih sistemov, ki jih poznamo pod kratico GNSS, se podrobneje poglobili v razvoj najpomembnejših med njimi in obravnavali njihovo osnovno delovanje ter razlike med njimi. INTRODUCTION Amongst the more important questions that humanity has asked itself throughout history is un­deniably ‘Where am I?’ We could also add to this question, ‘where is someone else?’ or ‘where is an object?’ These questions seem very simple at first, but humanity has been looking for a solution for a long time. The solution came with the introduction of global navigation satellite R Corresponding author: PhD student, Boštjan Krošelj, Faculty of Energy Technology, University of Maribor, Tel: +386 (0)41 436 901, Mailing address: Hocevarjev trg 1, Krško, Slovenia, E-mail address: bostjan.kroselj@gmail.com 1 University of Maribor, Faculty of Energy technology, Hocevarjev trg 1, Krško, Slovenia systems or GNSS for short. We use navigation systems every day. Often, we are not even aware of this, nor how much they make our daily tasks easier. Almost every smart device and every vehicle is already equipped with a GNSS signal receiver. Devices that show us the route and guide us to our destination have replaced traveling with a map in hand. We can see our location at any time, as well as how much longer it will take to reach our destination. HISTORY Man has always wanted to travel from point A to point B as fast and as safely as possible, ori­entating himself with the help of objects in nature. This; however, was much harder at sea. The Phoenicians navigated using the sun and the stars in 2000 BC. The distance travelled and the estimated speed of movement was measured using an hourglass. Even the Egyptians in 1500 BC relied on the sun and the stars. They also navigated using a simple depth gauge and wind as Vi­kings did between 900 and 1000 AD. It is clear that navigation began to develop at sea. Amongst the most important discoveries for navigation, especially by sea, was the compass. A lead cord at the end was also practical, used to sample the bottom and measure depth. The oldest known navigation system was used by the Greeks, who used a system of lighting bonfires on the tops of hills to sail across the Mediterranean Sea. These bonfires served as a point of reference for sailing at night. They started making and using maps with merely directions and ports. The problem was that they did not have a method that allowed the spherical surface of the Earth to be displayed. In around 1500, sailors began using devices to measure the angle of the sun and stars above the horizon, which they then used to determine latitude. Figure 1 shows a sextant precursor that was less accurate. Using the device, they measured the angle between the horizon and the sun or stars. [1] Figure 1: The predecessor of the sextant [1] Methods of measuring time were still very inaccurate. Even the best devices were losing about 10 minutes a day, which meant a huge loss over several weeks of sailing, and consequently an inaccurate location. Even speed measuring was improving very slowly. They were using the knot as a unit for measuring speed on water. It was not until the 18th century, when cartography de­veloped and it was thus possible to represent the earth in a spherical shape, that the compass came into play, allowing the determination of cardinal direction by accurately measuring time so that longitude could then be identified. In 1884, an agreement was signed confirming the Green­wich meridian as the Prime meridian. After the introduction of the universal radio signal from Greenwich, all that was needed was an ordinary clock and a radio on board to be able to calculate the exact location. In 1935, the British physicist Robert Watson developed the first useful radar system that was able to determine the location of objects outside the field of view using radio waves. Between 1940 and 1943, a navigation system called Loran was developed in the United States. It had limited signal coverage. [1] DEVELOPMENT OF GNSS There are currently four main types of GNSS systems in the world: the American GPS, the Russian GLONASS, the European Galileo and the Chinese BeiDou. In addition to the listed major GNSS, there are also regional navigation satellite systems known as RNSS. They complement GNSS by providing information on the accuracy, integrity and accessibility of the results provided by the GNSS system. 3.1 Development of GPS The development of GPS was encouraged by the launch of the Russian satellite Sputnik into orbit on October 4, 1957. Based on a signal emitted by Sputnik, a team of American experts be­gan their research. They theorised that the Doppler shift could be used to determine the exact location of a satellite. They then turned this idea upside down and put forward a new theory. If you know the location of the satellite, then you can determine the location of the navigator any­where on Earth. Based on this, the US Navy developed a system for locating submarines. In 1964, solar cells powered the Transit satellite. To determine location, the submarine had to expose the antenna to the surface within 10 to 16 minutes. The resulting location could then be provided with an accuracy of 25 metres in two dimensions. The deviation was much greater if the object was moving. Transit laid the foundations on which they later built the GPS system as we know it today. The US Air Force also conducted its satellite navigation programme in a three-dimensional space called the 621B. They focused their research on signal modulation, user data processing, receiver accuracy, error analysis, system costs, and Earth orbit. The greatest progress has been made in signal research. They also found the most effective signal, which they termed CDMA. The CDMA signal was scattered over a wide range of radio frequencies. With it, they achieved an accuracy of 5 metres in three-dimensional space. [2] With the merger of the Transit and 621B projects, a new project called NAVSTAR–GPS was cre­ated in 1973. The most important engineering challenges involved defining the specifics of the CDMA signal structure, developing space-based atomic clocks, achieving fast and accurate satel­lite orbit predictions, ensuring longer life of spacecraft, and developing custom GPS equipment. [3] In 1978, the first NAVSTAR BLOCK I development satellite was launched into orbit, followed by ten more across the years leading up to 1985. The system was designed for military purposes, which changed in 1983 when a Soviet fighter jet shot down a civilian passenger plane that ac­cidentally found itself in Soviet airspace. Following this event, US President Ronald W. Reagan announced that the NAVSTAR–GPS system would also be available for civilian use after it was fully operational. The initial operational capacity (IOC) was reached on 8 December 1993. BOCK I was replaced by nine BLOCK II satellites and 15 BLOCK IIA satellites. The system was declared fully operational (FOC) on 25 April 1995. An additional four BLOCK IIA satellites and 13 new BLOCK IIR satellites were later launched. [4], [5] 3.2 Development of GLONASS Development of GLONASS began in 1976, with the first satellite being launched into orbit in 1982. It formed part of the prototype generation Block I. The others were launched by 1985. The first generation was launched by 1990 and operated for four and a half years. Until the full establishment of the system, which required 24 operational satellites, the Russians had to launch 70 satellites. Many satellites were lost due to failed launches and many orbit failures. The system was declared FOC in 1996. The collapse of the Soviet Union greatly affected GLONASS. The sys­tem did not receive any investments and began to break down. This occurred to such an extent that by 2002, only seven satellites remained. It then began to slowly recover. The second-gener­ation satellite has remained in orbit since 2003, with an increased lifespan of up to seven years and an improved level of atomic clock stability. By 2011, the GLONASS system had reached FOC again. The second generation of satellites are equipped with an additional L2OF signal, which is intended for civilian use and allows the elimination of the ionosphere delay. [6] 3.3 Development of Galileo Galileo is a European GNSS that was not built for military purposes and is run by civil society. Europe wanted to offer an alternative to the American and Russian GNSS. They collaborated with the European Space Agency (ESA) on its development. Development began in 1999, when Germany, France, Italy and England merged, each with its own concept. A team of experts from all four countries then drew up a first phase plan, which was confirmed by the EU and ESA in 2003. [7] Galileo encountered financial problems because of issues with the US, which opposed the estab­lishment of the system and put pressure on the European Union to stop the project. Europe iden­tified its need for its GNSS and reached an agreement with the US concerning the security issues. The test satellite was launched in 2005. Due to financial difficulties, the project was taken over by the EU in 2007, while in the following year a second test satellite called GIOVE-B was launched into orbit, which together with the first GIOVE-A formed part of the first IOV phase. Based on the results, four IOV satellites were then launched in 2011 and 2012 for ground and space control segment testing purposes. In the second phase of the IOC, two satellites were launched into orbit in 2014, but they were used only for scientific purposes. By the end of 2016, the Galileo system was once again able to transmit a signal that could be used to determine location and time. There were 18 satellites in orbit, of which only 12 were operational. The Galileo system had 26 satellites in 2019, of which 22 were in operation. [6] 3.4 Development of BeiDou Eager for greater independence in order to increase its economic and social development, the Chinese government decided to establish its own GNSS. Development began in 1983 with the idea of two geostationary satellites, which were successfully tested in 1989, marking the begin­ning of a system called BeiDou. Experimental BDS-1 systems began to be developed in 1994. Two geostationary satellites were launched in 2000, while another was later launched in 2003. The system made it possible to locate and communicate via short messages in the area of China. The second phase of BDS-2 began in 2004. The system consisted of 14 satellites located at different altitudes in orbit that emitted signals at three different frequencies. The system allowed Asia and Pacific residents to determine location, navigation and exact time. The regional system was declared FOC on 27/12/2012 and was free of charge. The third phase of BDS-3 began in 2015 with five new satellites designed to perform tests of new signals. Exactly six years after the proc­lamation of the FOC, the system became global with 18 satellites in mid-Earth orbit and one in geostationary orbit. [6] GNSS BASIC OPERATION In order to determine the position of an object, GNSS uses the technique of accurately measur­ing the time required by the signal from the transmitter to the receiver. It is important that the transmitter is in a known location. The travel time of the signal is then multiplied by its speed. The result given is the distance from the transmitter to the receiver. This positioning technique is called geometric trilateration. The more accurate transmitter locations we have, the more accu­rate the receiver location will be. 4.1 Determining the position of the satellite The movement of satellites can be both predicted and accurately determined. The motion of sat­ellites is influenced by the gravity of the Earth, which is not completely round and does not have the same density. The forces of the Sun and Moon and other celestial bodies also influence the path of satellites. Satellites do not travel in a perfect vacuum and are subject to inhibitory force. The force of the sun’s radiation and the reflection of photons from the Earth also affects satel­lites. This occurs to such an extent that when moving towards the Sun, the satellites slow down and when moving away from the Sun, they accelerate. For this reason, ground base stations are used, with the help of which we are able to calculate the location of the satellite within 1 metre. Location data is divided into almanacs and ephemeris data. The almanac contains data on the path of satellites over a long period of time. Each satellite has an almanac for all other satellites in the constellation. The ephemeris data contain the exact location of the satellite for a short time in the future. Each satellite submits the location to itself, while the receiver must obtain informa­tion from each satellite in the field of view separately. With the location information, the receiver also obtains information about its validity. Receivers store data on ephemeris and almanacs in internal memory. Satellites have very accurate atomic clocks and simultaneously transmit a sig­nal, while the signal path defines the distance of the satellite from the receiver. With the help of a satellite signal, the receiver can determine the time. Determining the time is more accurate if there are multiple satellites. [7] 4.2 Two-dimensional location determination The measured distance of the satellite from the receiver is used for further calculation of the location. In the case of a single transmitter, such a result would mean that the receiver is located somewhere on the circle, distant for the result obtained. This is shown in Figure 2 a). If you want a more accurate result, you could use two transmitters with a known location. Hypothetically, we would determine the two points where the receiver is located. These two points define the intersections of the circles, as shown in Figure 2 b) as A and B. With the added third transmitter, our result is reduced to one point. This is the intersection of all three circles, as shown in Figure 2 c). It is important that all of the transmitters and the receiver have a precisely synchronised clock with each other. Figure 2: 2D distance determination with a) one, b) two or c) three transmitters 4.2 Three-dimensional location determination Let us start with one satellite with a known location and a receiver at an unknown location. Fig­ure 3 a) shows that the measurement result is a sphere around the satellite, while the location of the receiver is on the surface of this sphere. If we add a second satellite that emits a signal at the same time as the first, we see an additional sphere that is concentric to the second satellite. The receiver is also located on the surface of second sphere. Since the receiver must be located on the surface of both spheres at the same time, it may be on the perimeter of the shaded cir­cle shown by the arrow in Figure 3 b). It could also be at a point that touches both spheres, but this would only happen if the receiver was co-linear with the satellites. When the third satellite is added, the surface of its sphere intersects the shaded plane of Figure 3 b) at two points, as shown in Figure 3 c). Only one is right. Note that the points are mirrored with respect to the plane of the satellites. For a receiver on the ground, it makes sense that the right location is the one that is closer to the ground. [8] Figure 3: 3D distance determination with a) one, b) two or c) three transmitters [8] For the mathematical calculation of the positions of transmitters and receivers, it is necessary to choose a reference coordinate system. For the purpose of determining the position of satellites in orbit, it is convenient to use the Earth’s concentric inertial coordinate system ECI, where the starting point is at the centre of the Earth. The axes are fixed and directed in the direction of the stars, taking into account Newton’s laws of motion and gravity. It is necessary to determine the position of the earth at a given time. The ECI GPS coordinate system uses the orientation of the equatorial plane on 1.1.2000 at 12.00 and is known as the J2000 system. [8] 4.3 Measurement errors Even a small error in just one of the measurements can result in a difference of tens of metres. The most common sources of errors are atmospheric disturbances in the troposphere and ion­osphere. Disturbances are reflected in a decrease in the speed of light and have less impact at night. Precise atomic clocks in satellites, calibrated to 3 nanoseconds accurate, are very impor­tant. Other errors include electrical noise in the signal, signal path error as it can bounce off objects and prolong the path, incorrect satellite location, and satellite placement. The impact of errors can be reduced by increasing the number of satellites. 5 GPS GPS consists of a space, control and user segment. The space segment consists of 24 operational satellites arranged in six orbital planes, which are offset from each other by 60 ° and together form a complete 360 °. In addition to operational satellites, there are three backup satellites. All are located in the middle Earth orbit at an altitude of about 20,200 km above the earth’s surface. The satellites orbit the earth at an angle of 55 ° and an orbital time of 11 hours and 56 minutes. The radius of the satellites is 26,560 km. The satellites are positioned so that the receiver has a connection to at least four satellites anywhere on Earth at any time. The main task of the control segment is to monitor satellites and determine and predict their exact location, as well as to take care of the integrity of the system, the atomic clocks, atmosphere and satellite almanacs. The data collected by the control segments is then uploaded to the satellites. The control segment consists of 16 optimally positioned earth stations around the world that send data to satellites via 11 earth antennas. The main control station is located at the Schriever Air Force Base in Colora­do. The frequencies for sending the PRN signal are L1 (1575.42 MHz) and L2 (1227.60 MHz). The publicly available C/A code for civil use is implemented in the L1 signal, while the more precise P(Y) code, which is protected, is implemented in both frequencies. It uses WGS 84 for the GPS reference coordinate system. The user segment is usually called a GPS receiver, which processes signals sent from satellites to determine position, speed and time. A GPS receiver basically con­sists of an antenna, a receiver, a processor, a screen, and a power supply. It can be mounted on vehicles, ships or aircrafts, providing a precise position regardless of weather conditions. [6], [9] 5.1 GPS modernisation Later, L2C, L5 and L1C signals for civilian purposes were developed. The L2C has been available since 2005, when it was sent into orbit by Block IIR-M satellites. Due to the increased signal strength, it is more accessible in urban centres, cities with increased vegetation, and indoors. The simultaneous use of L2C and L1 C/A signals allows users to coordinate ionosphere shift. The L1C signal was developed with the Galileo system. The signal provides better reception in urban environments and has greater resistance to the reflection effect. For military purposes, the M code was developed. The L3 signal is reserved for the transmission of data from satellites to Earth stations and is used to detect nuclear explosions in the Earth’s atmosphere or near space. The L4 signal is used for additional research in the ionosphere. The L5 signal has been added to the system upgrade and is intended for a security system called Safety of Life. The signal is stronger than L2C and allows better reception, greater breakthrough and resistance. In combination with the L1 C/A and L2C signals, as well as the three-frequency technique, it enables more precise location determination up to 1 metre without the use of additional systems. [6] The GPS system satellite constellation consists of old and new. There is currently no longer a Block IIA generation satellite. There are 7 active Block IIR satellites, seven BLOCK IIR-M satellites, 12 Block IIF satellites and 4 GPS III / IIIF satellites. As of February 2, 2022, 30 satellites in the GPS constellation are operational. [10] GLONASS The space segment of the GLONASS system on 1 March 2022 consists of 25 satellites, of which 22 are operational. Two satellites are in the maintenance phase and one is in the test flight phase. The satellites are arranged in three MEOs with a spacing of 120 °. The satellites are spaced 45 ° apart in the orbital plane. They are located at an altitude of 19,100 km above the Earth’s surface. Each satellite needs 11 hours, 15 minutes and 44 seconds to make its way around the Earth. The layout allows constant coverage on the ground and altitudes up to 200 km. The constellation can be expanded to 64 satellites. Each satellite in the GLONASS constellation transmits the same PRN code, but at a different frequency. Unlike others, the GLONASS system uses the technique of differentiating FDMA signals and not CDMA. FDMA technology means multiple access with frequency distribution, while in CDMA technology satellites alternately send signals on one fre­quency without mixing the signals with each other. Signals are transmitted in two bands. L1 has a midrange of 1602 MHz and delays of 0.5625 MHz and L2 1246 MHz and delays of 0.4375 MHz. A constellation with 24 satellites needs 14 channels, which means that two of satellites located on opposite sides of the globe use the same channel. [6], [8] 6.1 GLONASS modernisation GLONASS satellites are known as Uragan. To this day, six satellite models have been launched, in­cluding a prototype and a zero-generation, belonging to Block I. Satellites of this generation were launched into orbit between 1982 and 1985 and had a lifetime of one year. The first generation included the Block IIa, Block IIb and Block IIv satellite models, which had their lifetime increased to four and a half years. They were launched between 1985 and 1990. The second generation of GLONASS-M satellites were into orbit from 2003 onwards. Their lifetime was extended to seven years. They also had improved atomic clock stability and were able to emit a second civil L2OF signal with the ability to eliminate ionosphere delay. The last satellite in the current constellation GLONASS was launched into orbit on October 25 2020 and is still in the test phase. Prior to this, two more satellites were launched in the same year, which are currently fully operational. They are all GLONASS-K generations, which will eventually be replaced by the GLONASS-K2 genera­tion. Although this should have started as early as 2020, to this day no satellite of this generation has been launched and they are still in development. At the same time, there is a plan to launch six GLONASS-B generation satellites into HEO in order to increase accuracy in urban centres, while GLONASS-KM satellites are also being developed for additional civilian signals. The ground control segment is also in the modernisation phase. The main control centre of the ground con­trol segment is located in Krasnoznamensk, near Moscow. The biggest problem is that almost all control stations were built on the territory of the former Soviet Republic, which means that satellites cannot be tracked and updated at all times. To solve this problem, two control stations were built, one in Brazil in 2014 and one in South Africa in 2017, which increased the accuracy of satellite locations in orbit, the accuracy of satellite atomic clocks, and improved signal coverage of the Earth. The construction of a station in India is also planned. Due to the use of FDMA signal, the manufacture of receivers was more complex and expensive. With the transition to CDMA signal, production is cheaper and more affordable. The ephemeris of the GLONASS system is given according to the ground reference coordinate system PZ-90. Its latest version PZ-90.11 was approved in 2014, is ITRF2014 compliant at the centimetre level, and allows better compatibility with other GNSS. [6] GALILEO The Galileo system is not yet at full operational capacity. When the system is complete, the sat­ellites will be divided into three orbital planes inclined at an angle of 56 ° to the equator. There will be one backup satellite in each plane in case of an operational failure. The satellite will take 14 hours to travel around the Earth. Galileo will be fully compatible with the American GPS sys­tem, which will further improve the accuracy of the system. Six to eight satellites are expected to be visible from most locations on Earth at any given time, which means an accuracy of a few centimetres. The satellites are located at an altitude of 23,222 km. There are currently 22 fully operational satellites and 4 test satellites in the constellation with a lifetime over 12 years. Full operational capacity is expected in 2023. [11] The ground segment consists of control and user. Control segment stations are set up around the world and allow you to obtain accurate information related to location, navigation and time. However, the majority of stations are within the European territory. The user segment has dif­ferent types of receivers. They are used for navigation in cars, mobile devices and other devices that require accurate location. [6] Galileo satellites emit 10 different navigation signals in the frequency bands E1, E5a, E5b, E6. At the same time, the satellites emit a SAR signal in the band intended for rescue services between 1544 MHz and 1545 MHz. The frequency ranges of Galileo signals are similar to the frequency ranges of the GPS system, so special modulation is required to avoid signal overlap. Galileo sig­nals provide a variety of services for both civilian and military use. Amongst the most famous is the Open Service, which is free for use primarily in vehicles and smart devices. It can be used by anyone who has a Galileo receiver anywhere in the world and allows for the accurate determi­nation of location. The system has up to 4-metre horizontal accuracy and up to 8-metre vertical accuracy when using two frequency receivers. This accuracy is slightly worse when using a single frequency receiver. High Accuracy Service is a paid closed service. The signal has improved ac­curacy and features simultaneous access to the signal on multiple frequencies. The Public Reg­ulated Service is intended for government use. Police, firefighters and similar services can use it. It is transmitted at two different frequencies to make the signal more resistant to intentional interference. The signal is intended for use in crisis situations. The SAR signal is designed to search for and rescue victims, allowing you to search around the world within 10 minutes and to be accurate to 5 km. A feedback signal is also introduced to inform the victim that their signal has been received and located. [6] 8 BEIDOU The BeiDou navigation system is also known by the acronym BDS. Satellites placed in three dif­ferent orbits represent the space segment. The satellites are located in GEO, IGSO and MEO orbits. Compared to other systems, BDS has the advantage that it has more satellites in higher orbits and better results in lower lying parts of the Earth. The earth segment consists of several earth stations run by the main control station. The clock is synchronised and data is uploaded onto the satellites by a special station, in addition to which there is an observation station, as well as management and control stations, that take care of the connection to the satellites. The BDS user segment offers a variety of products, systems and services that may be compatible with other GNSS. [12] The BDS system has a total of 35 satellites, 27 of which are intended for operational tasks and 3 of which will be spare. Like other systems, it is always fully operational 24 satellites. They will be divided into three orbital planes in the MEO with a delay of 120 °. The satellites will be at an altitude of 21,528 km and an inclination of 55 °. Satellites take 12 hours and 35 minutes to travel around the Earth. Five satellites will be located in GEO at an altitude of 35,786 km and the other three in IGSO with a delay of 120 °. The GEO and IGSO satellites are designed to cover China and the periphery of Asia, while the MEO satellites are designed to cover the rest of the Earth. [6] On 31 July 2020, the Chinese government declared the system fully operational, joining the US’s GPS, Russia’s GLONASS and Europe’s Galileo. The BeiDou-3 project was completed, with the sys­tem beginning to be used by more than 100 countries, mostly in Asia. With the new generation of BeiDou-3, China has started broadcasting a new civilian B1C signal at 1575.420 MHz, just like the L1 signal at GPS and the E1 signal at Galileo. This required modulation. Since open signals B2a and B2b are transmitted at the same frequencies as Galileo signals E5a and E5b, the pos­sibility of joint processing is open. The system enables both open and closed services, with the open system being for civilian use and free of charge. It provides accuracy of up to 5 metres in the Asia-Pacific region and up to 10 metres globally. A closed system for authorised users allows for more secure and accurate location, speed, time and communication. The control segment on Earth is limited, as all of its stations are located in China. The system uses its reference coordinate system called CGCS2000. [12] 9 CONCLUSION In this article, we have presented the global navigation satellite systems that are currently op­erational. The GPS system is still by far the most important and most developed, as well as the most widely used. Despite the fact that all but the Galileo system were initially built for military purposes, their use in civilian work has become so widespread that no one thinks of them as part of the military infrastructure anymore. In addition to GNSS, there are also smaller so-called regional systems that serve to obtain more accurate data in a certain area. Satellites of such sys­tems are placed only in geostationary orbit. Typical representatives of such systems include the Indian IRNSS and the Japanese QZSS. GNSS systems have become so common that we can hardly imagine life without them. Wher­ever we go, they accompany us and make our lives easier. We often forget their basic purpose, which is something we should fear. In wartime as well as in peacetime, such systems represent an excellent tool for espionage, monitoring the adversary’s military capabilities, as well as for co­ordinated attacks. We definitely want GNSS to be primarily used for civil and defence purposes. References [1] C. Riley: The History of Navigation, December 2021, Available: https://www.boatsafe. com/history-navigation/ [2] B. W. Parkinson, S. T. Powers: The origins of GPS and the pioneers who launched the system, North Coast Media, LLC, 2010 [3] G. W. Staff: The Origins of GPS, Fighting to Survive, 2022, Available: https://www. gpsworld.com/origins-gps-part-2-fighting-survive/ (15.1.2022) [4] E. Hall, History of the GNSS Industry and Milestones Ahead, North Coast Media, LLC, 2020 [5] M. Zrinjski, et al., Modernizacija GPS-a (GNSS-2), Geodetski list 59 (82), 2005 [6] M. Zrinjski, et al., Razvoj i modernizacija GNSS-a, Geodetski list 73 (96), 2019 [7] P. Mlakar, Globalni navigacijski satelitski sistemi, Fakulteta za racunalništvo in informa­tiko, Univerza v Ljubljani, 2010 [8] E. D. Kaplan, C. Hegarty, Understanding GPS/GNSS: Principles and Applications, Boston, Artech House, 2017 [9] P. Sirish Kumar, V. B. S. Srilatha Indira Dutt, The global positioning system: Popular accuracy measures, Materials today: proceedings 33: 4797-4801, 2020 [10] GPS.gov, GPS: The Global Positioning System, 2022, Available: https://www.gps.gov/ (1.3.2022) [11] EUSPA, Galileo is the European global satellite-based navigation system, 2022, Availa­ble: https://www.euspa.europa.eu/european-space/galileo/programme (27.02.2022) [12] B. N. S. System, System, 2022, Available: http://en.beidou.gov.cn/SYSTEMS/System/ (15.2.2022) 42 JET JET Volume 15 (2022) p.p. 43-64 Issue 2, 2022 Type of article: 1.01 www.fe.um.si/si/jet.html PREDICTION OF CAVITATION AND PARTICLE EROSION IN A RADIAL DIVERGENT TEST SECTION NAPOVED KAVITACIJSKE EROZIJE IN EROZIJE DELCEV V RADIALNO DIVERGENTNI TESTNI SEKCIJI Luka Kevorkijan1, Luka Lešnik1, Ignacijo Biluš1R Keywords: Cavitation, Particles ANSYS Fluent, Erosion, CFD, modelling, DPM Abstract The 3D unsteady, cavitating, particle-laden flow through a radial divergent test section was simulated with the homogeneous mixture model and Discrete Phase Model (DPM) within the commercial CFD code ANSYS Fluent. For turbulence, a RANS approach was adopted with the Reboud’s correction of turbulent viscosity in the k-. SST model. Cavitation erosion was predicted with the Schenke-Melissaris-Terwisga (SMT) model, while particle erosion was predicted with the Det Norske Veritas (DNV) model. Two distinct erosion zones were identified, one for pure cavitation erosion and one for pure particle erosion. The occurrence of the pure particle erosion zone downstream of the cavitation erosion zone was analysed. By observing the streamlines downstream of the cavitation structures, it was found that vortices form in the flow and redirect the particles towards the wall, causing a pure particle erosion zone on the wall. Particles under consideration in this study were not found to alter the flow to the extent that the cavitation ero­sion zone would be significantly altered compared with the results without solid particles which are reported in the literature. R Corresponding author: Associate Professor, Ignacijo Biluš, University of Maribor, Faculty of Mechanical Engineer­ing, Smetanova ulica 17, Maribor, Slovenia, Tel.: +386 (2) 220 7742, E-mail address: ignacijo.bilus@um.si 1 University of Maribor, Faculty of Mechanical Engineering, Smetanova ulica 17, Maribor, Slovenia Povzetek 3D nestacionaren, kavitirajoc tok z delci skozi radialno divergentno testno sekcijo je bil simuliran z modelom homogene zmesi in modelom diskretne faze (DPM) s komercialno RDT kodo ANSYS Fluent. Turbulenca je modelirana po pristopu RANS z Reboudovim popravkom turbulentne visko­znosti v modelu k-. SST. Izvedeni sta bili napoved kavitacijske erozije z modelom Schenke-Me-lissaris-Terwisga (SMT) in napoved erozije zaradi delcev z modelom Det Norske Veritas (DNV). Prepoznani sta bili dve razlicni erozijski coni, ena za zgolj kavitacijsko erozijo in ena za erozijo zgolj zaradi delcev. Analiziran je bil pojav cone ciste erozije zaradi delcev dolvodno od cone kavitacijske erozije. Z opazovanjem tokovnic dolvodno od kavitacijskih struktur je bilo ugotovljeno, da se v toku oblikujejo vrtinci in preusmerjajo delce proti steni, kar povzroca na steni cono erozije zgolj zaradi delcev. Ugotovljeno je bilo, da delci, obravnavani v tej študiji, ne spreminjajo toka do te mere, da bi se obmocje kavitacijske erozije znatno spremenilo v primerjavi z rezultati brez trdnih delcev, o katerih poroca literatura. INTRODUCTION Erosion of solid surfaces can occur in hydraulic systems as a result of multiple causes, including cavitation and solid particles. Cavitation is a phase change phenomenon that can arise when the pressure of a liquid drops to vaporisation pressure. As a result of this, vapour cavities (bubbles, or more generally vapour structures) form inside of the liquid. When vapour cavities are subjected to a pressure that ex­ceeds the vaporisation pressure, condensation takes place and they collapse. After the collapse of a vapour cavity, an instantaneous pressure peak occurs in the liquid. For spherical bubble collapse, Lord Rayleigh [1] derived the expression for the time evolution of the pressure in a liquid during the collapse. The pressure wave emitted after the collapse can cause the erosion of a solid surface as it impacts that surface. Hammit [2] introduced the hypothesis, based on experimental observations of cavitation damage in a Venturi channel, that erosion is related to the initial potential energy spectrum of vapour cavities. Vogel et al. [3] later derived the equation for potential energy of cavitation bubbles based on the observations of a single bubble collapse near a solid wall. However, when a bubble collapses in the vicinity of a solid wall, the collapse becomes asymmet­ric due to the presence of the wall on one side of the bubble. As a result of an asymmetric col­lapse, a liquid microjet develops and impacts the solid surface. Based on this observation, Plesset and Chapman [4] calculated the velocity of the impacting microjet and concluded that this can be the cause of solid surface erosion. Whether the pressure wave, microjet, or some combination of both causes cavitation erosion in larger flow scales (complex macroscopic cavitation structures as opposed to a single bubble) remains an open area of research. In engineering applications, erosion is undesired as it can reduce efficiency and ultimately cause failure of a given hydraulic system, for example hydraulic turbines, pumps, fuel injectors in in­ternal combustion engines, pipe elbows and valves. It is therefore important to incorporate the study of erosion risk in the design process of a hydraulic system. In the past, most of the insight into the cavitation erosion process was obtained by experiments. Recently, however, approaches based on computational fluid dynamics (CFD) have emerged. From the observations made by Hammit [2] and Vogel et al. [3], an energy cascade mechanism has been proposed by both Pereira et al. [5] and Fortes-Patella et al. [6] to explain the transfer of potential energy from large cavitation structures to the solid surface. Fortes-Patella et al. [6] introduced a cavitation erosion model based on an assumption that the erosive aggressiveness of large cavitation structures is proportional to their distance from the solid surface. Melissaris et al. [7] cavitation erosion prediction becomes more and more critical, as the requirements for more efficient propellers increase. Model testing is yet the most typical way a propeller designer can, nowadays, get an estimation of the erosion risk on the propeller blades. However, cavita­tion erosion prediction using computational fluid dynamics (CFD analysed the applicability of this cavitation erosion model, with a simplification that potential energy is calculated only in the first cell layer next to a solid surface. This finding is supported by previous observations by Philipp and Lauterborn [8], who found that cavitation structures in contact with the solid surface are the most aggressive. Leclercq et al. [9] which can be described as the fluid mechanical loa­ding leading to cavitation damage. The estimation of this quantity is a challenging problem both in terms of modeling the cavitating flow and predicting the erosion due to cavitation. For this purpose, a numerical methodology was proposed to estimate cavitation intensity from 3D un­steady cavitating flow simulations. CFD calculations were carried out using Code_Saturne, which enables U-RANS equations resolution for a homogeneous fluid mixture using the Merkle’s mo­del, coupled to a k-e turbulence model with the Reboud’s correction. A post-process cavitation intensity prediction model was developed based on pressure and void fraction derivatives. This model is applied on a flow around a hydrofoil using different physical (inlet velocities included all cavitation structures in cavitation erosion prediction by projecting the release of potential energy from the cells that are not in contact with a solid surface via the solid angle. While Leclercq et al. calculated solid angles for individual cells using a discrete form written by Van Osteroom and Strackee [10], Schenke and van Terwisga [11] instead introduced a continuous form of the solid angle. Based on the description of a global energy balance of a bubble in an infinite liquid [12], and the observed focusing of potential energy into the centre of a bubble cloud collapse [13], Schenke et al. [14] have introduced a model to describe the conversion of cavitation potential energy into kinetic energy, which is then released as a pressure wave. Melissaris et al. [15] have analysed three mathematical forms of a vapour volume fraction total derivative used in the cal­culation of cavitation potential energy, concluding that the least numerical error is introduced when using the form with mass source term. Melissaris et al. [15] and Pezdevšek et al. [16] also compared the two solid angle projection approaches. The above-described cavitation erosion models have been developed with an important limita­tion that both the vapour and the liquid phase are considered as a homogeneous mixture. This approach is widely used, and several cavitation models [17]–[19] exist and have been analysed [20]–[25]. However, an alternative approach is to track some if not all of the vapour phase as individual bubbles in a Lagrangian frame [26]–[32]. Hybrid (multiscale) approaches, where only sufficiently small cavitation structures are converted into bubbles and tracked in a Lagrangian frame, are particularly promising for cavitation erosion prediction. When the solid particles present in the flow impact a solid surface, they too can cause erosion of that surface. Finnie [33] identified key parameters that determine the mass loss of a solid sur­face that is impacted by solid particles. An expression that relates the solid particle impact angle to the erosion of the solid surface has been written. Similarly, Bitter [34], [35] later proposed a different impact angle dependency. There are many empirical expressions in the literature, some of the more often used are by Oka et al. [36], [37], Ahlert [38] and Det Norske Veritas (DNV) [39]. Empirical particle erosion models can be used in CFD simulations in combination with Lagrangian particle tracking, where particles are assumed to be point-masses. When particles hit a solid surface, they rebound from that surface. The rebound behaviour is described using empirical expressions that relate the particle velocity after the rebound with the particle impact angle [40], [41]. In the past, researches usually focused their attention on one phenomenon only, either the cav­itation erosion or the particle erosion. Few studies have been conducted, where both phenom­ena are present at the same time. These are mostly experimental studies [42]–[44], but some studies apply CFD simulations as well [45]–[47]. In this work a CFD simulation of cavitation in a particle-laden flow has been conducted on a geom­etry commonly used in cavitation erosion studies using a commercial CFD software ANSYS Fluent 2021R2. Both the cavitation and the particle erosion are predicted using the Schenke-Melissar-is-van Terwisga cavitation erosion model and the DNV particle erosion model respectively. The Schenke-Melissaris-van Terwisga cavitation erosion model was previously implemented within the ANSYS Fluent using User Defined Function (UDF), written in a C programming language [48]. The mesh was generated with a Fluent meshing module version 2020 R2. SIMULATION OF CAVITATION IN A PARTICLE-LADEN FLOW 2.1 Geometry and mesh Radial divergent test section [28] geometry (Fig. 1) was chosen, as it is often the choice of numer­ical cavitation erosion studies [31], [32], [48]–[51] and because the present authors previously analysed the implemented cavitation erosion model on this geometry [48]. Figure 1: Radial divergent test section geometry and dimensions The mesh was generated using ANSYS Fluent meshing module. To reduce the computational effort, only a 45° sector was meshed. The resulting poly-hexcore mesh is shown in Fig. 2 and the mesh statistics are presented in Table 1. The mesh nodalization study was conducted by the authors in the previous work [48]. Figure 2: Mesh of a 45° sector: a) whole domain, b) zoomed in on the refined part of the domain, c) zoomed in on the radius Table 1: Mesh statistics Mesh statistic Value Number of cells 1,396,703 Minimum orthogonal quality 0.28 Maximum aspect ratio 65 Minimum y+ dimensionless wall distance 0.174 Maximum y+ dimensionless wall distance 51.9 2.2 Governing equations 2.2.1 Homogeneous mixture Vapour and liquid phases are treated as a homogeneous mixture, the mixture density (.) and the mixture viscosity (µ) are determined via the mixing rule ....=.........+(1-....)..... (2.1) ....=.........+(1-....)..... (2.2) where ais the vapour volume fraction, the liquid and vapour densities of water are .l =998.85 kg/m³ and .v = 0.01389 kg/m³. The liquid and vapour viscosities of water are µl = 0.0011 Pa·s and µv = 9.63 · 10-6 Pa·s. The flow of mixture is governed by the mixture mass conservation equation .. +. ·( .u)=0 , (2.3) .t where u is the mixture velocity field, and by the mixture momentum conservation equation . ( .u ) +. · ( . uu)=- . p+. · [ µ (. u+. uT )]+ .g+ SM , (2.4) .t where p is the pressure, g is the gravitational acceleration and SM is the momentum source term due to the presence of particles in the flow. 2.2.2 Turbulence modelling The RANS approach was adopted and the SST k-. model was chosen. In order to reproduce the transient cavitation effects that occur due to the compressibility of a liquid-vapour mixture, the turbulent viscosity µt needs to be modified according to Reboud et al. [52] as follows: ....(....).... 1 .....= , .... ......... (2.5) .............1 ....*,.......... where k is the turbulence kinetic energy, . is the specific rate of dissipation of the turbulence kinetic energy, a* is the damping coefficient, S is the strain rate magnitude, F2 is the second blending function and a1 is the model constant equaling 0.31. Correction of turbulent viscosity is achieved by modifying the mixture density function f(.) such that a lower turbulent viscosity is achieved in the presence of both the liquid and vapour phases ( .-. )n f ( .)=.v + v )n- 1 , (2.6) -. ( .l v where n is an arbitrary exponent with a recommended value of 10. This modification of turbulent viscosity was implemented in ANSYS Fluent via a UDF (User De­fined Function) [48]. 2.2.3 Cavitation modelling To model the mass transfer between the liquid and the vapour phase, an additional equation for conservation of vapour mass is solved . ( a . v ) =. · ( a. u)=.vSv, (2.7) .t va where Sav is the source of the vapour phase. The Schnerr-Sauer cavitation model was used, which states that mass transfer source term Sav is ..( p-p ) vl 2 v a(1 -a) 3 if p< p vv v . Rb 3 .l Sa =, (2.8) v .v .l 2 ( p- p) v a (1 -a ) 3 if p> p vv v . Rb 3 .l {vv where the bubble radius is ( a ) 1/ 3 Rb =v. 4 (2.9) n0 p ( 1 -a) 3 v For the bubble number density n0 a default value 1·1011 of was used. The vapour pressure was set to pv = 1854 Pa. 2.2.4 Cavitation erosion modelling Hammit [7] and Lauterborn et al. [8] expressed the potential energy of a cavitation bubble in terms of work done by the surrounding liquid on a bubble. In general, for a cavity with volume VV, the potential energy of that cavity equals Epot =( pd- pv ) VV, (2.10) where pd is the pressure in the surrounding liquid and pv is the vapour pressure inside the bubble. The cavitation potential power ppot is introduced as the total derivative of the potential energy Epot D ( pd- pv ) DV V P= V +( p-p) . (2.11) pot Vdv Dt Dt After neglecting the first term in equation (2.11), as it was found to be at least an order of magni­tude lower than the second term [15], and taking into account that pv = const. (isothermal flow), the cavitation potential power is written as DV V P=( p-p) . (2.12) pot dv Dt The driving pressure of the surrounding liquid pd is taken as a time-averaged value in each com­putational grid cell as proposed in [11] predicting the instantaneous surface impact power of collapsing cavities from the potential energy hypothesis (see Hammitt, 1963; Vogel and Lauter-born, 1988 ..... = 1..*....(....)........ , (2.13) ....* where p(t) is an instantaneous pressure that is then time-averaged over the current simulated time t*. For the numerical computation of the cavitation potential power, it is convenient to calculate the instantaneous change of the volume specific potential energy [11] or cavitation erosion potential P e ´pot = pot , (2.14) V cell where Vcell is the computational mesh cell volume. The total derivative of the vapour volume in equation (2.12) can be calculated using three different but mathematically equivalent formula­tions: { .a +u·. a .t 1 DV V . =. ·u (2.15) V Dt .l-. cell v . Sa v .l However, the numerical error introduced varies between the three formulations [15], [48] and it was found that most accurate results were obtained using the third formulation from equation (2.15) [15], [48], which expresses the total derivative of the vapour volume ... .. with the cavita­ tion source term Sav. Finally, the cavitation erosion potential is given by e ´pot =( pd- pv) .Sa . (2.16) .l v Based on the previous studies [12], [13] which suppresses the buoyant pressure gradient that otherwise deteriorates the sphericity of the bubbles. We measure the radius of the rebound bubble and estimate the shock energy as a function of the initial bubble radius (2-5.6mm Schen­ke, Melissaris and van Terwisga proposed a model [14] to describe the focusing of initial cavity potential energy into the cavity collapse centre. This accumulation and delayed release of energy is described by the additional conservation equation for collapse induced kinetic energy e .e +. ·( ui e )=-e´rad ( t ) , (2.17) .t where ui is the collapse induced velocity, érad(t) and is the instantaneous radiated power of the pressure wave released in the cavity centre at the final stage of collapse. A model for transport equation (2.17) is then written in explicitly discretized form: ........ ....|....=....|.+.....=(1-....|.).(....-1)............-....(.....-1).... (2.18) ......... Pu in equation (2.18) is the projection operator, which ensures that the specific kinetic energy of the collapse accumulates on the inside of the vapour-liquid interface, K is the global energy con­servation parameter and ß is the energy release criterion. These are further defined as: u·. e Pu =max,0, (2.19) [|u||. e|] .... . ............. . ..... ....= , (2.20) .. ....... ........ .1, if....>.....and.....=0 ....=(2.21) 0, else. The energy released as a consequence of the finished collapse of a cavity at time level t + .t, where .t is the time step, is expressed as the specific power of the pressure wave (........)|. .......|....=(2.22) ..... An infinite wave propagation speed is assumed in this model, therefore the energy received by a surface S element is expressed as the instantaneous surface specific impact power ( x-xS) ·n 1 cell e´ S = .´ erad dV , (2.23) 4 p [] V |xcell -x S|3 where xcell is the position vector of the computational cell centre, xS is the position vector of the wall surface element centre, and n is the surface element normal. Instantaneous surface specific J impact energy with units of ( m²) is obtained by integrating equation (2.23) over the simulated time (t): t e=. ´ edt. (2.24) SS 0 The described model to predict cavitation was previously implemented inside ANSYS Fluent via a UDF[48]. 2.2.5 Particle motion in a fluid Particles are tracked in a Lagrangian frame, while the motion of particles is governed by Newton’s second law. In a Lagrangian frame particle velocity is defined as d xP v = , (2.25) P dt where xP is the particle position in the global inertial frame of reference. We assume the spherical particle, therefore particles are assumed to undergo only translation without rotation. Newton’s second law is then written as d 3 P dv P .=F, (2.26) P 6 dt where the particle density is .P = 1700 kg/m3, the particle diameter is dP = 5 µm and is the result­ant force F acting on the particle. The resultant force can be expressed as a sum of different contributing forces F=F+ F+F+F + F D BG PGVM , (2.27) where FD is the drag force, FB is the buoyancy force, FG is the force of gravity, FPG is the pressure gradient force and FVM is the virtual mass force. The drag force acting on the particle is written as 18 µdCR PD P F=( u-v P) , (2.28) D 6 24 where the drag coefficient is given by a correlation by Morsi and Alexander [53] cc C=c+ 2 + 3 , D 12 (2.29) R P R p c1, c2 and c3 are constants of the cD . RP curve fit and particle Reynolds number is defined as (2.30) The force resulting from buoyancy and gravity is (2.31) the pressure gradient force is (2.32) and the virtual mass force is (2.33) Two-way coupling between the continuous and discrete phase is achieved by including the mo­mentum source term in equation 2.4. 2.2.6 Particle-wall interaction To resolve the interaction between the particle and the wall, a simplified hard sphere approach is adopted, where the loss of particle momentum is accounted for via restitution coefficients in the normal and tangential directions with respect to the observed wall. This scenario is presented in Fig. 3. Figure 3: The interaction between the particle6 and the wall In accordance with the hard sphere approach, the normal particle velocity after rebounding from the wall is (2.34) and the tangential particle velocity after rebounding from the wall is v=ev t ,2 tt ,1 (2.35) where en is the normal coefficient of restitution, et is the tangential coefficent of restitution, vn,1 is the particle normal velocity before interaction with the wall, and vt,1 is the particle tangential velocity before interaction with the wall. Coefficients of restitution are calculated using empirical relations, expressed as polynoms of the impact angle. In this work model by Grant and Tabakoff [41], the formulae en =0.993 -1.76 .+1.56 .2 -0.49 .3 , (2.36) et =0.988 -1.66 .+2.116 .2 - 0.67 .3 , (2.37) are used, where is the particle wall impact angle. 2.2.7 Particle erosion modelling For particle erosion modelling, the DNV model [39] was used. The walls’ material is considered to be steel. Erosion rate of a surface element with units of [kg/m2] is obtained by totalling the empirical expression for single particle erosion over all of the particles that have impacted the observed surface element m ER=. p Avbf ( . ) , 1 (2.38) .S N p face face where mP is particle mass, .face is the surface density (for steel, a value 7800 kg/m3 of was used), Sface is the surface element area, A is the empirical constant with the value of 2·10-9, v1 is the par­ticle impact velocity as shown in Fig. 3, b is the material dependent velocity exponent with the value of 2.6, and f(.) is the dimensionless function of the particle impact angle. For steels, the dimensionless function is given by 8 i f ( . )=. ( -1 )i+1 Bia, (2.39) i=1 where the model constants Bi are given in Table 2. Table 2: DNV dimensionless impact angle function model constants B1 B2 B3 B4 B5 B6 B7 B8 9.370 42.295 110.864 175.804 170.137 98.398 31.211 4.170 2.3 Boundary conditions and simulation setup Boundaries of the computational domain are shown and labelled in Fig. 4, while detailed bound­ary conditions are presented in Table 3. The direction of gravitational acceleration was assigned in the z-direction. For the Lagrangian tracking of particle, a Dispersed Phase Model (DPM) was used, where indi­vidual particles can be grouped into parcels to reduce the computational effort. However, this behaviour was overwritten by assigning only one particle to each parcel. Along with the inlet velocity, presented in Table 3, it is also necessary to prescribe the particle mass flow rate, from which the number of particles entering the domain is calculated. Particle mass flow rate of 3.897 kg/s was chosen to achieve the particle mass concentration of 0.05 kg/m3 which was reported by Gregorc [54] to be the insoluble sediment mass concentration of the river Drava. Figure 4: Boundaries of the computational domain Table 3: Boundary conditions Boundary Boundary condition Homogeneous mixture phase Dispersed phase Inlet Normal velocity = 31 m/s Normal velocity = 31 m/s Outlet Static pressure = 10.1 bar Opening (to leave the domain) Wall No-slip wall Rebound; Grant & Tabakoff Symmetry Symmetry Rebound; Grant & Tabakoff The Navier-Stokes equations were solved using the SIMPLE algorithm. To evaluate the gradients, the Least Squares Cell-Based method was selected, while the pressure cell face values were in­terpolated by the PRESTO! interpolation scheme. For spatial discretization, QUICK spatial discre­tization scheme was selected for all equations. For temporal discretization of the unsteady terms in equations, the Bounded Second-Order Implicit time integration scheme was selected. A time step of 1·10-6 was chosen based on the previous cavitation study by the authors. The time step was checked and was found to be smaller than the particle relaxation time expressed as .d 2 pp 24 tr = (2.40) 18 µ CD R p A scaled residuals convergence criterion of 1·10-3 was achieved before the imposed limit of 100 iterations per time step. Overall, 0.016815 s of physical time were simulated. RESULTS AND DISCUSSION Erosion due to the cavitation and particles is evaluated on the bottom wall (when the geome­try is positioned as shown in Fig. 4) of the radial divergent test section. Figure 5. presents the erosion prediction in terms of coloured contours, where in Fig. 5 a) coloured contours represent the particle erosion results from equation 2.38 and in Fig. 5 b) coloured contours represent the cavitation erosion results from equation 2.24. The main zone where cavitation erosion occurs is located between 19 mm and 32 mm from the axis of symmetry, which agrees with the experimental results [28] and numerical simulation re­sults from different authors [31], [32], [51]. Since in these studies particles were not present in the flow, we conclude that in this study particles do not change the flow behaviour to the extent that it would change the cavitation erosion zone. It must be noted; however, that this is not nec­essarily true for particles of a different diameter, density and volume fraction. Maximum particle erosion occurs downstream from the cavitation erosion zone. This can be explained by observing the flow streamlines as shown in Fig. 6 a) and Fig. 6 b). Vapour cavities disturb the flow in the downstream direction, multiple vortices can be seen in Fig 6. a) and Fig 6. b). The relaxation time of the particles considered in this study is smaller than the flow relaxation time, therefore particles respond quickly to any changes in the flow, such as the change of flow direction due to a vortex. The vortices observed in Fig. 6 a) and Fig. 6. b) redirect particles toward the bottom wall, causing them to impact the wall at a higher impact angle. In DNV erosion model the dimensionless function of the particle impact angle (equation 2.39) rapidly approaches zero as the impact angle approaches zero. From Fig. 6 a) and Fig. 6 b), it can be seen that the stream­lines are virtually parallel to the bottom wall in the region of cavitation, therefore there is almost no particle erosion predicted (at least an order of magnitude smaller than in the main zone be­tween 40 mm and 52 mm from the symmetry axis). Figure 5: Erosion prediction on the bottom wall: a) particle erosion, b) cavitation erosion Figure 6: Flow dynamics: a), b) cavitation represented by the iso-surface of vapour volume frac­tion of 0.2 and flow dynamics represented by streamlines coloured by the flow velocity at two different times, c) time evolution of total vapour volume, with selected times highlighted with a red square CONCLUSIONS In this paper, a transient simulation of the cavitating and particle-laden flow is presented. Risk of erosion due to both the cavitation and the particles is assessed by employing the appropriate erosion modelling approaches. By using the erosion models, it is found that for the investigated particles, the flow is not suf­ficiently affected by the presence of the particles to change the location of the main cavitation erosion zone, indicating that a one-way momentum coupling between the continuous phase and the discrete phase would be sufficient. A separate zone of pure particle erosion is predicted downstream of the cavitation erosion zone. This is explained by observing the formation of vortices that redirect particles towards the wall downstream of the vapour clouds. The particles investigated in this study could be particularly useful for an experimental validation as there is no overlap between the cavitation erosion zone and the particle erosion zone, which would require the ability to distinguish one type of erosion from the other. ACKNOWLEDGEMENTS The authors wish to thank the Slovenian Research Agency (ARRS) for their financial support with­in the framework of the Research Programme P2-0196 Research in Power, Process and Environ­mental Engineering. References [1] Lord Rayleigh: VIII. On the pressure developed in a liquid during the collapse of a spher­ical cavity, Philosophical Magazine Series 6, Vol. 34, No. 200, pp. 94–98, 1917, doi: http://dx.doi.org/10.1080/14786440808635681 [2] F. Hammit: Observations on cavitation damage in a flowing system, ASME Journal of Basic Engineering, Vol. 85, No. 3, 1963 [3] A. Vogel, W. Lauterborn, and R. Timm: Optical and acoustic investigations of the dynamics of laser-produced cavitation bubbles near a solid boundary, Journal of Fluid Mechanics, Vol. 206, No. September 1989, pp. 299–338, 1989, doi: 10.1017/ S0022112089002314 [4] M. S. Plesset and R. B. Chapman: Collapse of an Initially Spherical Vapor Cavity in the Neighborhood of a Solid Boundary, Journal of Fluid Mechanics, Vol. 47, part 2, pp. 283– 290, 1971 [5] F. Pereira, F. Avellan, and P. H. Dupont: Prediction of cavitation erosion: An energy approach, Journal of Fluids Engineering, Transactions of the ASME, Vol. 120, No. 4, pp. 719–727, 1998, doi: 10.1115/1.2820729 [6] R. Fortes Patella, J.-L. Reboud, and L. Briancon Marjollet: A Phenomenological and nu­merical model for scaling the flow agressiveness in cavitation erosion, 2004. [Online]. Available: https://www.researchgate.net/publication/281921326 [7] T. Melissaris, N. Bulten, and T. J. C. Van Terwisga: On the applicability of cavitation erosion risk models with a URANS solver, Journal of Fluids Engineering, Transactions of the ASME, Vol. 141, No. 10, 2019, doi: 10.1115/1.4043169 [8] A. Philipp and W. Lauterborn: Cavitation erosion by singlelaser-produced bubbles, Jour­nal of Fluid Mechanics, Vol. 361, pp. 75–116, 1998, doi: 10.1017/S0022112098008738 [9] C. Leclercq, A. Archer, R. Fortes-Patella, and F. Cerru: Numerical cavitation intensi­ty on a hydrofoil for 3D homogeneous unsteady viscous flows, International Journal of Fluid Machinery and Systems, Vol. 10, No. 3, pp. 254–263, 2017, doi: 10.5293/ IJFMS.2017.10.3.254 [10] A. Van Oosterom and J. Strackee: The Solid Angle of a Plane Triangle, IEEE Transac­tions on Biomedical Engineering, Vol. BME-30, No. 2, pp. 125–126, 1983, doi: 10.1109/ TBME.1983.325207 [11] S. Schenke and T. J. C. van Terwisga: An energy conservative method to predict the erosive aggressiveness of collapsing cavitating structures and cavitating flows from nu­ merical simulations, International Journal of Multiphase Flow, Vol. 111, pp. 200–218, 2019, doi: 10.1016/j.ijmultiphaseflow.2018.11.016 [12] M. Tinguely, D. Obreschkow, P. Kobel, N. Dorsaz, A. De Bosset, and M. Farhat: Energy partition at the collapse of spherical cavitation bubbles, Physical Review E–Statistical, Nonlinear, and Soft Matter Physics, Vol. 86, No. 4, pp. 1–6, 2012, doi: 10.1103/Phys­RevE.86.046315 [13] Y. C. Wang and C. E. Brennen: Numerical computation of shock waves in a spherical cloud of cavitation bubbles, Journal of Fluids Engineering, Transactions of the ASME, Vol. 121, No. 4, pp. 872–880, 1999, doi: 10.1115/1.2823549 [14] S. Schenke, T. Melissaris, and T. J. C. Van Terwisga: On the relevance of kinematics for cavitation implosion loads, Physicsof Fluids, Vol. 31, No. 5, 2019, doi: 10.1063/1.5092711 [15] T. Melissaris, S. Schenke, N. Bulten, and T. J. C. van Terwisga: On the accuracy of pre­dicting cavitation impact loads on marine propellers, Wear, Vol. 456–457, No. June, p. 203393, 2020, doi: 10.1016/j.wear.2020.203393 [16] M. Pezdevsek, L. Kevorkijan, and I. Bilus: Cavitation Erosion Modelling – Comparison of Two Solid Angle Projection Approaches, International Journal of Simulation Modelling, Vol. 21, No. 2, pp. 249–260, 2022, doi: 10.2507/ijsimm21-2-600 [17] P. J. Zwart, A. G. Gerber, and T. Belamri: A two-phase flow model for predicting cavita­tion dynamics, International Conference on Multiphase Flow, No. January 2004, p. 152, 2004 [18] G. H. Schnerr and J. Sauer: Physical and Numerical Modeling of Unsteady Cavitation Dynamics, 4th International Conference on Multiphase Flow, No. June, pp. 1–12, 2001 [19] M. P. Kinzel, J. W. Lindau, and R. F. Kunz: A unified homogenous multiphase CFD mod­el for cavitation, American Society of Mechanical Engineers, Fluids Engineering Divi­sion (Publication) FEDSM, Vol. 1B-2017, No. January 2018, 2017, doi: 10.1115/FED­SM2017-69363 [20] I. Bilus, M. Morgut, and E. Nobile: Simulation of sheet and cloud cavitation with ho-mogenous transport models, International Journal of Simulation Modelling, Vol. 12, No. 2, pp. 94–106, 2013, doi: 10.2507/IJSIMM12(2)3.229 [21] E. Goncalves and R. F. Patella: Numerical simulation of cavitating flows with homo­geneous models, Computers and Fluids, Vol. 38, No. 9, pp. 1682–1696, 2009, doi: 10.1016/j.compfluid.2009.03.001 [22] M. Morgut and E. Nobile: Influence of the Mass Transfer Model on the Numerical Pre­diction of the Cavitating Flow Around a Marine Propeller, Second International Sympo­sium on Marine Propulsors smp’11, No. June, pp. 1–8, 2011 [23] A. Niedzwiedzka, G. H. Schnerr, and W. Sobieski: Review of numerical models of cavi­tating flows with the use of the homogeneous approach, Archives of Thermodynamics, Vol. 37, No. 2, pp. 71–88, 2016, doi: 10.1515/aoter-2016-0013 [24] M. Pezdevšek, I. Biluš, and G. Hren: COMPARISON OF CAVITATION MODELS FOR THE PREDICTION OF CAVITATION AROUND A HYDROFOIL, Journal of Energy Technology, Vol. 14, No. 1, pp. 41–55, 2021 [25] T. Schenke, Sen; van Terwisga: Finite Mass Transfer Effects in Cavitation Modelling, Proceedings of the 19th Numerical Towing Tank Symposium, pp. 1–6, 2015 [26] W. Jian, M. Petkovšek, L. Houlin, B. Širok, and M. Dular: Combined numerical and experimental investigation of the cavitation erosion process, Journal of Fluids Engineer­ing, Transactions of the ASME, Vol. 137, No. 5, 2015, doi: 10.1115/1.4029533 [27] J. B. Carrat, R. Fortes-Patella, and J. P. Franc: Experimental and numerical investigation of the erosive potential of a leading edge cavity, International Journal of Fluid Machin­ery and Systems, Vol. 12, No. 2, pp. 136–146, 2019, doi: 10.5293/IJFMS.2019.12.2.136 [28] N. Berchiche, J. P. Franc, and J. M. Michel: A cavitation erosion model for ductile ma­terials, Journal of Fluids Engineering, Transactions of the ASME, Vol. 124, No. 3, pp. 601–606, 2002, doi: 10.1115/1.1486474 [29] M. Dular and M. Petkovšek: On the mechanisms of cavitation erosion–Coupling high speed videos to damage patterns, Experimental Thermal and Fluid Science, Vol. 68, No. June 2015, pp. 359–370, 2015, doi: 10.1016/j.expthermflusci.2015.06.001 [30] I. Biluš, M. Hocevar, M. Dular, and L. Lešnik: Numerical Prediction of Various Cavi­tation Erosion Mechanisms, Journal of Fluids Engineering, Vol. 142, No. 4, 2020, doi: 10.1115/1.4045365 [31] A. Peters: Numerical Modelling and Prediction of Cavitation Erosion Using Euler-Euler and Multi-Scale Euler-Lagrange Methods, Universität Duisburg-Essen, 2019 [32] M. H. Arabnejad, U. Svennberg, and R. E. Bensow: Numerical assessment of cavitation erosion risk using incompressible simulation of cavitating flows, Wear, Vol. 464–465, No. November 2020, p. 203529, 2021, doi: 10.1016/j.wear.2020.203529 [33] I. Finne: Erosion of surfaces, Wear, Vol. 3, pp. 87–103, 1960, doi: 10.1016/0043­1648(60)90055-7 [34] J. G. A. Bitter: A Study of Erosion Phenomena Part I, Wear, Vol. 6, pp. 169–190, 1963 [35] J. G. A. Bitter: A study of erosion phenomena. Part II, Wear, Vol. 6, No. 3, pp. 169–190, 1963, doi: 10.1016/0043-1648(63)90073-5 [36] Y. I. Oka, K. Okamura, and T. Yoshida: Practical estimation of erosion damage caused by solid particle impact: Part 1: Effects of impact parameters on a predictive equation, Wear, Vol. 259, No. 1–6, pp. 95–101, 2005, doi: 10.1016/j.wear.2005.01.039 [37] Y. I. Oka and T. Yoshida: Practical estimation of erosion damage caused by solid parti­cle impact: Part 2: Mechanical properties of materials directly associated with erosion damage, Wear, Vol. 259, No. 1–6, pp. 102–109, 2005, doi: 10.1016/j.wear.2005.01.040 [38] K. R. Ahlert: Effects of Particle Impingement Angle and Surface Waiting on solid Partilce Erosion of AISI 1018 Steel, 1994 [39] Det Norske Veritas: Erosive wear in piping systems, 2007 [40] A. Forder, M. Thew, and D. Harrison: A numerical investigation of solid particle erosion experienced within oilfield control valves, Wear, Vol. 216, No. 2, pp. 184–193, 1998, doi: 10.1016/S0043-1648(97)00217-2 [41] G. Grant and W. Tabakoff: An Experimental Investigation of the Erosive Characteristics of 2024 Aluminum Alloy, Springfield, 1973 [42] J. Sato, K. Usami, and T. Okamura: Basic Study of Coupled Damage Caused by Silt Abrasion and Cavitation Erosion, JSME International Journal, Vol. 34, No. 3, pp. 292–297, 1991, [Online]. Available: https://www.jstage.jst.go.jp/article/ bpb1993/17/11/17_11_1460/_pdf/-char/ja [43] K. Su, J. Wu, and D. Xia: Dual role of microparticles in synergistic cavitation–particle ero­sion: Modeling and experiments, Wear, Vol. 470–471, p. 203633, 2021, doi: 10.1016/j. wear.2021.203633 [44] F. Innings, E. Hultman, F. Forsberg, and B. Prakash: Understanding and analysis of wear in homogenizers for processing liquid food, Wear, Vol. 271, No. 9–10, pp. 2588–2598, 2011, doi: 10.1016/j.wear.2011.01.084 [45] L. Kevorkijan, J. Ravnik, and I. Biluš: Numericno modeliranje kavitacijske erozije in abrazije delcev na profilu krila NACA 0015, in Kuhljevi dnevi 2021, 2021, No. September [46] P. J. Dunstan and S. C. Li: Cavitation enhancement of silt erosion: Numerical studies, Wear, Vol. 268, No. 7–8, pp. 946–954, 2010, doi: 10.1016/j.wear.2009.12.036 [47] L. A. Teran, S. Laín, and S. A. Rodríguez: Synergy effect modelling of cavitation and hard particle erosion: Implementation and validation, Wear, Vol. 478–479, No. December 2020, 2021 [48] L. Kevorkijan, L. Lešnik, and I. Biluš: Cavitation Erosion Modelling on a Radial Divergent Test Section Using RANS, Strojniski Vestnik/Journal of Mechanical Engineering, Vol. 68, No. 2, pp. 71–81, 2022, doi: 10.5545/sv-jme.2021.7364 [49] D. Greif: Kavitacijska erozija v vbrizgalnih komponentah z uporabo polidisperznega ka­vitacijskega modela, 2012 [50] M. S. Mihatsch, S. J. Schmidt, and N. A. Adams: Cavitation erosion prediction based on analysis of flow dynamics and impact load spectra, Physics of Fluids, Vol. 27, No. 10, 2015, doi: 10.1063/1.4932175 [51] S. Mouvanal, D. Chatterjee, S. Bakshi, A. Burkhardt, and V. Mohr: Numerical predic­tion of potential cavitation erosion in fuel injectors, International Journal of Multiphase Flow, Vol. 104, pp. 113–124, 2018, doi: 10.1016/j.ijmultiphaseflow.2018.03.005 [52] J.-L. Reboud, B. Stutz, and O. Coutier-Delgosha: Two-phase flow structure of cavita­tion: experiment and modelling of unsteady effects, 1998 [53] S. A. Morsi and A. J. Alexander: An investigation of particle trajectories in two-phase flow systems, Journal of Fluid Mechanics, Vol. 55, No. 2, pp. 193–208, 1972, doi: 10.1017/S0022112072001806 [54] B. Gregorc: VPLIV TRDNO-KAPLJEVITIH ZMESI NA OBRATOVALNE KARAKTERISTIKE, Uni-verza v Mariboru, 2011 Nomenclature (Symbols) . µ a .l . v µl µ v t u p g SM µt k . . aS F2 a 1 f ( .) n S av Rb n 0 pv VV (Symbol meaning) mixture density mixture viscosity vapour volume fraction liquid density vapour density liquid viscosity vapour viscosity time mixture velocity field pressure gravitational acceleration momentum source term turbulent viscosity turbulence kinetic energy specific rate of dissipation of the turbulence kinetic energy damping coefficient strain rate magnitude second blending function SST turbulent model constant mixture density function mixture density function modifying exponent source of the vapour phase bubble radius bubble number density vapour pressure vapour (cavity) volume Epot potential energy of a vapour cavity pressure in the surrounding liquid pd Ppot cavitation potential power p(t ) instantaneous pressure . t current simulated time e ´ pot cavitation erosion potential V computational mesh cell volume cell e collapse induced kinetic energy ui collapse induced velocity e´(t ) instantaneous radiated power of the pressure wave rad P projection operator u K global energy conservation parameter ß energy release criterion .t time step e´ S instantaneous surface specific impact power xcell position vector of the computational cell centre xS position vector of the wall surface element centre n surface element normal eS instantaneous surface specific impact energy vP particle velocity xP particle position in the global inertial frame of reference particle density .P dparticle diameter P F resultant force acting on the particle drag force FD Fbuoyancy force B Fforce of gravity G pressure gradient force FPG FVM virtual mass force CD drag coefficient R first constant of the C . curve fit c 1 DP R second constant of the C D . . curve fit P R third constant of the Ccurve fit c 3 DP R particle Reynolds number P vn ,1 particle normal velocity before interaction with the wall vn ,2 particle normal velocity after rebounding from the wall e normal coefficient of restitution n vt ,1 particle tangential velocity before interaction with the wall vt ,2 particle tangential velocity after rebounding from the wall etangential coefficient of restitution t . particle wall impact angle ER erosion rate mp particle mass .density of the wall surface face Sface surface element area A particle erosion empirical constant v1 particle impact velocity b particle erosion material dependent velocity exponent f ( . ) dimensionless function of the particle impact angle Bi DNV particle erosion model constants t particle relaxation time r JET Volume 15 (2022) p.p. 65-80 Issue 2, 2022 Type of article: 1.01 www.fe.um.si/si/jet.html ANALYSIS OF THE INFLUENCE OF PARAMETERS WHEN CHARGING AND DISCHARGING A CAPACITOR USING DIFFERENTIAL EQUATIONS ANALIZA VPLIVA PARAMETROV PRI POLNJENJU IN PRAZNJENJU KONDENZATORJA Z UPORABO DIFERENCIALNIH ENACB Matic KrašovicR, Peter Virtic1 Keywords: RC circuit, RLC circuit, state space, transfer function, differential equations, parame­ters R, L, C, damping Abstract Two different electrical circuits were analysed in this paper. The first was an RC circuit consisting of a resistor R, a capacitor C, and a DC voltage u, while the second was an RLC circuit consisting of a resistor R, a capacitor C, a DC voltage u and an inductor L. Both circuits were described using the state space mathematical model and, on this basis, graphs for charging and discharging a capacitor were plotted. Both circuits were described using differential equations for electrical current through a capacitor and voltage over a capacitor. Finally, different values of the R (re­sistance), L (inductance) and C (capacitance) parameters were taken, and graphs were plotted for voltage over capacitor. The purpose of the study was to explore how different values of parame­ters influence capacitor charging and discharging. R Corresponding author: Matic Krašovic, University of Maribor, Faculty of Energy Technology, Hocevarjev trg 1, 8270 Krško, Slovenia, Email: matic.krasovic@student.um.si 1 University of Maribor, Faculty of Energy Technology Povzetek V tem clanku sta bili analizirani dve razlicni elektricni vezji. Prvo je bilo RC vezje, sestavljeno iz upora R, kondenzatorja C in enosmerne napetosti u. Drugo je bilo RLC vezje, ki je bilo podobno prejšnjemu vezju, le da je bila dodana tuljava L. Obe vezji sta bili pretvorjeni v matematicni model v prostoru stanj in na podlagi tega so bili izrisani razlicni grafi za polnjenje in praznjenje konden­zatorja. Obe vezji sta opisani z diferencialnimi enacbami za elektricni tok skozi kondenzator in napetost na kondenzatorju. Na koncu so bile izbrane razlicne vrednosti parametrov R (uporno­sti), L (induktivnosti) in C (kapacitivnosti) in na podlagi tega so bili izrisani grafi za napetost na kondenzatorju. Namen je bil videti, kako razlicne vrednosti parametrov vplivajo na polnjenje in praznjenje kondenzatorja. INTRODUCTION Nowadays, different electrical circuits can be found almost everywhere. Analysing these circuits and describing them mathematically, as well as seeing what is actually going on inside them, can really help engineers to better understand such circuits. This was among the purposes of this article, i.e. to describe two different circuits in detail using differential equations and state space. A state space is a mathematical notation of a circuit in the form of matrices, where differential equations are pre-existing, and is a straightforward way to describe the response of the system. Below there are two simple equations that will be used later in this paper, which demonstrate that they have a derivative in it: .....=....... (1.1) (1.2) When these equations are used as a part of a system of equations necessary to obtain a solution, a differential equation is obtained. In [1] the authors discussed RLC circuits with resonant frequency and bandwidth, while in [2] the authors described RLC circuits using system theoretic properties with differential equations and matrices. In [3] the authors compared parameters and their influence on the streamer discharge process of a gas park switch in a series RC circuit, and in [4] the authors used a mathematical approach of optimised waveform relaxation for a small RC circuit. Numerous models in the form of matrices were also used. In [5] the authors introduced capacitor charging and discharging through the Ohm and Kirchhoff Laws using charge and capacitance. In comparison with the authors in [1] and [2], the authors of this paper used a similar approach, albeit an analytical and numerical one, including the use of differential equations and matrices for capacitor charging and discharging in an RLC circuit. In order to compare this paper with that of the authors in [3], one of the main themes is the influence of parameters; the difference being that the authors [3] used a slightly different approach and plotted different graphs. In compari­son to the authors in [4], the authors of this paper also used a mathematical approach, with more for voltage and electrical current through capacitor and less for waveform relaxation. Finally, the greatest similarity of this paper is to that of author [5], who used charge q and capacitance C to describe current through the capacitor in an RC circuit when being charged and discharged. This paper, however, has mostly electrical current i, voltage u and resistance R. In contrast to existing articles, the authors of this paper introduced RC and RLC circuit state space models in order to observe how the voltage over capacitor changes when different values of parameters are used. For instance, in the RLC circuit, when a bigger inductance L was used, the capacitor was charged with a higher overshoot and settling time. A more in-depth description of the impact of the parameters is found in the following sections. RC CIRCUIT ANALYSIS 2.1 Analytical method Figure 1 shows an RC circuit consisting of a resistor R, a capacitor C and a DC supply u. Figure 1: RC circuit with a resistor R, capacitor C and a DC supply u At the beginning of charging, the capacitor should be empty so that uc=0V. When the supply is turned on, after a short time the capacitor becomes partially charged. By writing a loop for the voltages: ...........+..... (2.1.1) and using a derivation, while knowing that: ................. (2.1.2) and: .............. (2.1.3) it is possible to solve a differential equation and get the electrical current through the capacitor: .·. .......·...... (2.1.4) Where i is the electrical current through the capacitor, u is the supply voltage, R is the resistance, e is the Euler constant, C is the capacitance and t is time. The similar applies to the voltage over capacitor and, after solving it, the following is obtained: .....=.....1-........ (2.1.5) The above was carried out analytically; the numerical method is described below. 2.2 Numerical method 2.2.1 MATLAB and State space First, a circuit needed to be modelled using the state space model, which was programmed in MATLAB. ..... =.................. (2.2.1.1) ....=................... (2.2.1.2) 1 [.....1]=.- .·[....1]...1.·.... (2.2.1.3) ................ y=[1].·[.....]..0·.... (2.2.1.4) The matrices represent a mathematical notation, called state space. A is called the state or system matrix, b is the input matrix, c is the output matrix and d is a feedthrough or feedforward matrix. Based on this notation, graphs were plotted at the given values: R=1000O and C=1000e-6F. Figure 2: Charging and discharging of a capacitor in an RC circuit Graphs represent the transfer function of the system. As can be seen, when it is charging the voltage is higher over time and lower when discharging. 2.2.2 MATLAB and Simulink A similar thing was done, again in MATLAB, but now with the addition of Simulink. Again, some calculations needed to be done. The loop, from section 1, had to be converted into a derivative of voltage over capacitor. After solving it, the following equation is obtained: ...... ...... (2.2.2.1) ..=.... Based on the previously derived differential equation, a block diagram was built in Simulink for this circuit. Figure 3: Simulink block diagram The code was written in MATLAB, connected to Simulink through the MATLAB ‘sim’ function, so that the code could work perfectly, based on the Simulink block diagram. The value of R was 1000O and that of C was 1000e-6F. The graphs that were plotted are: supply voltage u, voltage across capacitor and resistor and the electrical current through the capacitor. Figure 4: Voltages and current in the RC circuit As can be seen above, at some point the voltage across the capacitor is the same as the supply voltage, which means it reached the highest voltage value. Thereafter, however, the capacitor starts discharging, and the voltage slowly falls to zero, as does the current. Another point to men­tion is the time constant t (‘tao’). This constant represents the time when the capacitor is 63% charged, which can easily be calculated as t = R·C. In this case, R was 1000O and C was 1000e-6F. Therefore, when entered into the equation, t is equal to 1s. However, since it started from t=1s, the graph shows that t=2s, so the 1s needs to be subtracted from 2s, thus resulting in t=1s. RLC CIRCUIT ANALYSIS 3.1 Analytical method Figure 5 shows an RLC circuit consisting of a resistor R, a capacitor C, a DC-supply u and an indu­ctor L. Figure 5: RLC circuit with a resistor R, capacitor C, inductor L and a DC supply u Some equations needed to be written for the analytical approach. First, the loop for voltages: ...........+.....+..... (3.1.1) A similar approach was used, as in the case of the RC circuit, but uL was also used, which by de­finition is: .............. (3.1.2) In order to consider everything and put it into an equation, a second order differential equation for the electrical current through the capacitor was solved: .... ...+ . .......0 (3.1.3) ...+.... .. Where ....... is a second derivative of electrical current through capacitor, ..... is a first derivate of the electrical current through the capacitor and L is inductance. For the voltage over capacitor, the equation that was solved is a linear differential equation of the first order: ... ... .......... ........ (3.1.4) .. ...... .. Similar to the RC circuit, a numerical method was also used. 3.2 Numerical method 3.2.1 MATLAB and state space To begin, the RLC circuit had to be converted into a state space. ..... =.................. (3.2.1.1) ....=................... (3.2.1.2) -. -. ...... .. . ..... ..· ........... 1 = .... .·.... (3.2.1.3) .... 0 0 ......·.....y=[0 1]. ...0·.... (3.2.1.4) Matrices a, b, c and d again represent a mathematical notation called state space. Based on the matrices, graphs were plotted at values: R=10O, C=1000e-6F, L=1H. Figure 6: Charging and discharging of a capacitor in an RLC circuit As is clearly demonstrated above, the flow is not the same as in the RC circuit. This is because the damping is present due to the two energy sinks that are in the circuit, which are the capacitor and inductor. The higher the resistance R, the bigger the damping, therefore it was necessary to take the right value of resistance R, at R=10O. The response at charging is extremely fast, then damping does its job and the amplitude slowly lowers until it hits zero, then at some point it starts discharging, as is seen in the second graph. 3.2.2 MATLAB and Simulink A numerical method was also used in MATLAB and Simulink. First, an equation for voltage over capacitor needed to be calculated, even more so than the one from chapter 3.1, equation (3.1.4). Since in the circuit all three elements are in series, i.e. connected one after another, the electrical current is the same through all three elements, hence it can be said that: .................. ....... (3.2.2.1) Because all three currents are the same, instead of iL, which is the same as iC, it can be written as: .............. (3.2.2.2) In equation (3.1.4), instead of iL, the equation (3.2.2.2) is written, thus: ... ... .......... ........ (3.2.2.3) .. ...... .. After some work and solving the equation, again a derivative of uC is needed, and since it already ....... , is in the equation, when iL is written as .. , this is what is solved: .... ... ............. ........ (3.2.2.4) ..... .. .... ... Where ... is a second derivative of voltage over capacitor over time, .. is a first derivate of voltage over capacitor over time. The equation represents a second order differential equation for voltage over capacitor. Based on this equation, again a Simulink block diagram was constructed to represent the voltage over capacitor, so that the graphs could be plotted. Figure 7: Simulink block diagram Again, the code had to be written in MATLAB and connected to Simulink through the MATLAB ‘sim’ function, so that the code could work perfectly, based on the Simulink block diagram. The value of R was 10O, of C was 1000e-6F and L was 1H. The graphs that were plotted are: supply voltage, voltage over capacitor, voltage over resistor, voltage over inductor and electrical current through a capacitor. Figure 8: Voltages and current in the RLC circuit As can be clearly seen above, damping of the system occurred. When the voltage over capacitor reaches a constant value, it means that the capacitor is fully charged, then at some point it starts discharging and falls to zero. It is interesting to note that the voltage over capacitor and inductor are reversed. CHANGING THE PARAMETERS The main purpose of this article was to see how parameters, when being varied, influence the charging and discharging of a capacitor. Therefore, different variations of parameters in an RC and RLC circuit were made. After every subchapter, a few main observations were written. More detailed observations were made in the conclusion, so that the variations of the parameters could be compared. 4.1 RC circuit 4.1.1 Constant R, variation of C Initially, R was constant, and C was varied. Figure 9: Capacitor charging when C was varied Values of parameters: R C1;C2;C3 1000O 1000e-6F 1000O 1000e-5F 1000O 1000e-4F It is seen above that the lower the C, the more the capacitor was charged (step response). When the capacitance C was 1000e-4H, the capacitor was barely charged. 4.1.2 Constant C, variation of R Here, capacitance C was constant, while the resistance R was being varied. -uC1 - uC2 -uC3 -uC1 - uC2 -uC3 Values of parameters: R1;R2;R3 C 100O 1000e-6F 1000O 1000e-6F 10000O 1000e-6F It is noticeable that, as in the previous example when R was constant, when R=100O, i.e. at the lowest value, the capacitor is most highly charged. 4.2 RLC circuit 4.2.1 Constant C and R, variation of L Here, capacitance C and resistance R were constant, while induction L was varied. Values of parameters: R C L1;L2;L3 10O 1000e-6F 1H 10O 1000e-6F 5H 10O 1000e-6F 10H It seems that the higher the inductance, the more the capacitor charges, and at the lowest indu­ctance it charges with the highest – albeit the least – overshoot. 4.2.2 Constant L and R, variation of C Here, inductance L and resistance R were constant, and capacitance C was varied. -uC1 - uC2 -uC3 Values: R C1;C2;C3 L 10O 1000e-6F 1H 10O 1000e-5F 1H 10O 1000e-4F 1H The graph again shows that at the smallest capacitance, the capacitor filled up the most, whereas at the largest C, the charge behaved like a first order system, due to the damping being equal to 1, which is called critical damping. 4.2.3 Constant C and L, variation of R Finally, capacitance C and inductance L were constant, while resistance R was varied. -uC1 - uC2 -uC3 Values: R1;R2;R3 C L 100O 1000e-6F 1H 1000O 1000e-6F 1H 10000O 1000e-6F 1H These graphs illustrate the worst response of the system. The fastest response was at R=100O, while at R=10000O the response was almost linear. CONCLUSION With the aim of finding at which parameters the capacitor charges the most (which was also the main point of the article), different values of parameter and variations were taken and made in RC and RLC circuits, which were clearly illustrated in the graphs. For the RC circuit, it can be seen that when varying resistance R and capacitance C, the response barely changed, and both were almost the same. Therefore, at constant R, it is better to take the lower C, so that capacitor is most highly charged. At constant R, it is also better to take the lower, so that capacitor can be charged more rapidly. Different things happened in the RLC circuit. First, when the variation of inductance L was made, three concise waves were plotted, which indicated that at the highest inductance, the capacitor charged with the highest overshoot with the highest settling time, while at the lowest inductance, the capacitor charged with the lowest overshoot, which was the shortest settling time. Next, the variation of capacitance C was opposite to that of inductance. At the lowest capaci­tance, the capacitor charged with the highest overshoot and also had the fastest response. At the highest capacitance C, the response did not even exceed 1V and it behaved like a first order system (if only one capacitor/inductor was in the system, but the circuit had two), due to the damping being equal to 1, which is called critical damping. Finally, when the variation of resistance R was made, something different happened. A similar response occurred at R=100O and R=1000O. When R was 100O it was the fastest response, while at uC=1V it was fully charged. At R=1000O the capacitor was also fully charged at uC=1V, however the response was slightly slower and it took more time for the capacitor to charge to maximum. At the highest resistance R, i.e. 10000O, the response was the slowest and did not reach the vol­tage, as the previous two did, and was almost linear. Therefore, it can be concluded that taking the lowest resistance is the best option, so the capacitor is charged the fastest. The results of the findings and research at the RLC-circuit were seen when varying the inductance L and keeping R and C constant. The differences between RC and RLC-circuit are described and can be easily seen. References [1] Atsushi Sakurai, Bo Zhao, Zhuomin M. Zhang: Resonant frequency and bandwidth of metamaterial emitters and absorbers predicted by an RLC model, Journal of Quantita­tive Spectroscopy & Radiative Transfer 149 (2014) [2] Timo Reis, Fedor Glazov: System theoretic properties of linear RLC circuits, IFAC Paper-sOnLine (2021) [3] Shichao Zheng, Zhongjian Kang, Lei Li, Anqi Zhang, Kai Zhao, Yaxun Zhou: Influence of series RC circuit parameters on the streamer discharge process of gas park switch, Elsevier, Vacuum 193 (2021), 110518 [4] Mohammad Al-Khaleel, Martin J. Gander, Albert E. Ruehli: A mathematical analysis of optimized waveform relaxation for a small RC circuit, Elsevier, Applied Numerical Mathematics 75 (2014) [5] Roman Ya. Kezerashvili: Teaching RC and RL Circuits Using Computer-supported Exper­iments, Elsevier, IERI Procedia 2 (2012) Nomenclature (Symbols) (Symbol meaning) t time i current iL current through inductor iR current through resistor iC current through capacitor u supply voltage uC voltage over capacitor uR voltage over resistor uL voltage over inductor C capacitance L inductance R resistance e Euler constant duC/dt First derivate of voltage over capacitor over time diL/dt First derivate of current through inductor over time d2uC/dt2 Second derivate of voltage over capacitor over time d2iC/dt2 First derivate of current through capacitor over time 80 JET Author instructions ww .fe.um.si/en/je w ww.fe.um.si/jet.h tml MAIN TITLE OF THE PAPER SLOVENIAN TITLE Author1, Author 2, Corresponding authorR Keywords: (Up to 10 keywords) Abstract Abstract should be up to 500 words long, with no pictures, photos, equations, tables, only text. Povzetek (Abstract in Slovenian language) Submission of Manuscripts: All manuscripts must be submitted in English by e-mail to the editorial office at jet@um.si to ensure fast processing. Instructions for authors are also available online at http://www.fe.um.si/en/jet/author-instructions.html. Preparation of manuscripts: Manuscripts must be typed in English in prescribed journal form (MS Word editor). A MS Word template is available at the Journal Home page. A title page consists of the main title in the English and Slovenian language; the author(s) name(s) as well as the address, affiliation, E-mail address, telephone and fax numbers of author(s). Corre­sponding author must be indicated. Main title: should be centred and written with capital letters (ARIAL bold 18 pt), in first paragraph in English language, in second paragraph in Slovenian language. Key words: A list of 3 up to 6 key words is essential for indexing purposes. (CALIBRI 10pt) Abstract: Abstract should be up to 500 words long, with no pictures, photos, equations, tables, ­ text only. Povzetek: - Abstract in Slovenian language. Main text should be structured logically in chapters, sections and sub-sections. Type of letters is Calibri, 10pt, full justified. R Corresponding author: Title, Name and Surname, Organisation, Department, Address, Tel.: +XXX x xxx xxx, E-mail address: x.x@xxx.xx 1 Organisation, Department, Address 2 Organisation, Department, Address Authors names and surnames Units and abbreviations: Required are SI units. Abbreviations must be given in text when first men­tioned. Proofreading: The proof will be send by e-mail to the corresponding author in MS Word’s Track changes function. Corresponding author is required to make their proof corrections with accepting or rejecting the tracked changes in document and answer all open comments of proof reader. The corresponding author is responsible to introduce corrections of data in the paper. The Editors are not responsible for damage or loss of submitted text. Contributors are advised to keep copies of their texts, illustrations and all other materials. The statements, opinions and data contained in this publication are solely those of the individual authors and not of the publisher and the Editors. Neither the publisher nor the Editors can accept any legal responsibility for errors that could appear during the process. Copyright: Submissions of a publication article implies transfer of the copyright from the author(s) to the publisher upon acceptance of the paper. Accepted papers become the permanent property of “Journal of Energy Technology”. All articles published in this journal are protected by copyright, which covers the exclusive rights to reproduce and distribute the article as well as all translation rights. No material can be published without written permission of the publisher. Chapter examples: MAIN CHAPTER (Arial bold, 12pt, after paragraph 6pt space) 1.1 Section (Arial bold, 11pt, after paragraph 6pt space) 1.1.1 Sub-section (Arial bold, 10pt, after paragraph 6pt space) Example of Equation (lined 2 cm from left margin, equation number in normal brackets (section. equation number), lined right margin, paragraph space 6pt before in after line): Equation (1.1) Tables should have a legend that includes the title of the table at the top of the table. Each table should be cited in the text. Table legend example: Table 1: Name of the table (centred, on top of the table) Paper title Figures and images should be labelled sequentially numbered (Arabic numbers) and cited in the text – Fig.1 or Figure 1. The legend should be below the image, picture, photo or drawing. Figure legend example: Figure 1: Name of the figure (centred, on bottom of figure, photo, or drawing) References [1] N. Surname: Title, Journal Title, Vol., Iss., p.p., Year of Publication [2] N. Surname: Title, Publisher, Year of Publication [3] N. Surname: Title [online], Publisher or Journal Title, Vol., Iss., p.p., Year of Publication. Available: website (date accessed) Examples: [1] J. Usenik: Mathematical model of the power supply system control, Journal of Energy Technology, Vol. 2, Iss. 3, p.p. 29 – 46, 2009 [2] J. J. DiStefano, A.R. Stubberud, I. J. Williams: Theory and Problems of Feedback and Con­trol Systems, McGraw-Hill Book Company, 1987 [3] T. Žagar, L. Kegel: Preparation of National programme for SF and RW management taking into account the possible future evolution of ERDO [online], Journal of Energy Technol­ogy, Vol. 9, Iss. 1, p.p. 39 – 50, 2016. Available: http://www.fe.um.si/images/jet /Volume 9_Issue1/03-JET_marec_2016-PREPARATION_OF_NATIONAL.pdf (7. 10. 2016) Example of reference-1 citation: In text [1], text continue. Nomenclature (Symbols) (Symbol meaning) t time JET I Journal of Energy Technology I Vol. 15, Issue 2, September 2022 I University of Maribor, Faculty of Energy Technology ISSN 1855-5748 84 JET 9 771855 574008