19th SVSFEM ANSYS Users' Group Meeting and Conference 2011

Proceedings of

19th SVSFEM ANSYS Users' Group Meeting and Conference 2011 th

th

19 – 20 of October 2011, Congress Center in Prague, Czech Republic

Sponzoři:

19th SVSFEM ANSYS Users' Group Meeting and Conference 2011

TechSoft Engineering & SVSFEM

© TechSoft s.r.o, SVS FEM s.r.o.

ISBN: 978-80-905040-0-4

http://aum.svsfem.cz

2



Content

VÝVOJ METODIKY SIMULACÍ PÁDOVÝCH ZKOUŠEK V LS-DYNA A EXPERIMENTÁLNÍ OVĚŘENÍ .................... 5 J. Holub1, M. Michálek1, W. Ganzert2, M. Angerbauer2 OPTIMALIZACE TALÍŘOVÉ PRUŽINY TURBODMYCHADLA ........................................................................... 13 Radek Jandora FSI ANALÝZA JEZOVÉ KLAPKOVÉ HRADÍCÍ KONSTRUKCE ............................................................................ 21 Jirsák V., Kantor M., Nowak P. METODIKA FSI APLIKACE PŘI SIMULACI VYNUCENÉHO KMITÁNÍ LOPATEK PARNÍ TURBÍNY ..................... 30 J. Jurenka1, D. Gallovič1, T. Míšek2 SOFTWARE ANSYS PRO NÁVRH A OPTIMALIZACI ELEKTRICKÝCH STROJŮ A ZAŘÍZENÍ, MOŽNOSTI MULTIFYZIKÁLNÍCH ANALÝZ ....................................................................................................................... 42 Jakub Hromádka, Jindřich Kubák A steady state approach to calculation of valve pressure rise rate characteristics .................................... 49 Ján Oravec ŘEŠENÍ VNĚJŠÍ AERODYNAMIKY KOLEJOVÝCH VOZIDEL .............................................................................. 58 Milan Schuster SIMULACE PODCHLAZENÉHO VARU VE FLUENTU 13 ................................................................................. 65 Ing. Ladislav Vyskočil Ph.D. VLIV ČÁSTEČNÉ SMÁČIVOSTI NA PROUDĚNÍ A KAVITAČNÍ OBLAST V LAVALOVĚ DÝZE .............................. 87 B. FRODLOVÁ1, P. RUDOLF2, L. ZAVADIL3, M. KOZUBKOVÁ4, J. RAUTOVÁ5 PROBLEMATIKA NUMERICKÉ SIMULACE FENOMÉNŮ SPOJENÝCH S ŽÁROVÝMI NÁSTŘIKY ........................ 95 Ing. Tomáš Zábranský, Ing. Jan Hamza CFD CODE AS A HELPER FOR THE DESIGN OF THE WATER JET CUTTING DEVICE ...................................... 102 ZDENĚK ŘÍHA NUMERICKÁ SIMULACE PROUDĚNÍ STUPNĚM S VYROVNÁVACÍMI ŠTĚRBINAMI..................................... 110 Bartoloměj Rudas, Zdeněk Šimka, Petr Milčák, Ladislav Tajč, Michal Hoznedl A STEADY STATE APPROACH TO CALCULATION OF VALVE PRESSURE RISE RATE CHARACTERISTICS ...... 118 Ján Oravec MODELOVANIE SPAĽOVANIA HNEDÉHO UHLIA V KOTLY A ĎALŠIE PRÍKLADY ......................................... 126 Ing. Rastislav Beniak, Ing. Csaba Latika


3



PARNÍ STROJ - VÝPOČTY S POHYBLIVÝMI SÍTĚMI ..................................................................................... 135 P. ZÁcha CFD SIMULACE TEPLOTNĚ-HYDRAULICKÉ CHARAKTERISTIKY NA MODELU PALIVOVÉ TYČI V OBLASTI DISTANČNÍ MŘÍŽKY .................................................................................................................................... 144 D. Lávička PREDICTION OF FLOW THROUGH THE BURNER IN THE INDUSTRIAL-SIZE COMBUSTION CHAMBER ....... 150 Jiří Vondál VYUŽITÍ 3D NUMERICKÉHO MODELOVÁNÍ PRO STANOVENÍ EFEKTU PŘEKÁŽEK V KANALIZAČNÍM POTRUBÍ .................................................................................................................................................... 160 Karel Kříž*, Jaroslav Pollert ml.*, Martin Kantor** MODELOVÁNÍ PŘEPADU VODY PŘES POHYBLIVOU KLAPKOVOU KONSTRUKCI ........................................ 169 V. Jirsák, M. Kantor, P. Sklenář FSI ANALÝZA BRZDOVÉHO KOTOUČE TRAMVAJE ...................................................................................... 180 Michal Moštěk THERMODYNAMIC DESIGN OF THE ROOF ................................................................................................. 198 Ing.arch. Ing. Milan Palko, PhD. WINDOW FOR PASSIVE HOUSES AND ITS LOCATION IN THE WALL ........................................................... 206 Ing.arch. Ing. Milan Palko, PhD. HEAT AND MOISTURE ANALYSIS OF FUNCTIONAL GAPS FOR WOODEN WINDOWS ................................ 213 Ing.arch. Ing. Milan Palko, PhD. ALUMINUM WINDOW FOR THE PASSIVE HOUSE ..................................................................................... 225 Ing.arch. Ing. Milan Palko,PhD.


4



VÝVOJ METODIKY SIMULACÍ PÁDOVÝCH ZKOUŠEK V LS-DYNA A EXPERIMENTÁLNÍ OVĚŘENÍ J. HOLUB1, M. MICHÁLEK1, W. GANZERT2, M. ANGERBAUER2 1

MBtech Bohemia, Praha; 2ANDREAS STIHL AG & Co. KG, Waiblingen

Abstract: The design of hand-held gasoline powered or battery driven forestry and gardening equipment is dominated by the demand of lightweight products. However an ergonomic design with a user friendly configuration, the unlimited operator security and compliance with the governmental regulations are a major goal. The sum of all leads to devices with diligently optimized components yielding sophisticated structures from the load and stress point of view. This paper describes the systematic approach that was undertaken towards developing a numerical drop test model. The lessons learned from the initial steps were incorporated in the final analysis model and resulted in an overall high agreement between the analytical and experimental results. Keywords: Drop test, FEM, experimental verification, various drop orientation, result composition. 1

Úvod

Ve snaze o zkrácení vývojového času nových produktů, ve snaze zvýšení jejich kvality a v neposlední řadě pro snížení nákladů vývoje nabývá na významu digital prototyping a proniká i do oblastí, kde doposud převládalo experimentální testovaní. Naším cílem bylo vyvinout a ověřit spolehlivou, robustní, ale i rychlou a levnou metodiku pro simulace pádových resp. nárazových zkoušek, která bude respektovat všechna specifika dané skupiny produktů. Cílem bylo také maximálně využít synergie s ostatními již vysoce rozvinutými výpočetními metodami (například pro řešení dynamických jevů) na stejných produktech, které jsou v současné době řešeny v FEM v programech ANSYS a PERMAS a v MBS v programu ADAMS (Ganzert, 2010). Pro řešení pádových a nárazových problémů jsme zvolili explicitní FE program LS-Dyna, produkt firmy Livermore Software Technology Corp. (HALLQUIST, 2007). Vývojová práce byla realizována na akumulátorových zahradních nůžkách STIHL HSA 85 (Obrázek 1) (www.stihl.cz), ale měla by být přenositelná na všechny podobné stoje. Pro AKU- produkty platí přísnější kriteria u pádových zkoušek než pro produkty se spalovacím motorem a také samotný vývoj akumulátorových produktů je obvykle rychlejší a nezbývá zde příliš prostoru pro experimentální iterace při vývoji pádové odolnosti. Akumulátorové nůžky sestávají z relativně malého počtu dílů, kde je pro pádové zkoušky dominantní skříň nůžek na kterou jsou dále navěšeny všechny ostatní části. Skříň je geometricky složitý díl se spoustou žeber a z výrobních důvodů s výrazně proměnnou tloušťkou stěny.


5



Obrázek 1 - Akumulátorové zahradní nůžky STIHL HSA 85

Protože jedním z velmi důležitých parametrů je hmotnost celého systému, jsou zde použity lehké materiály jako slitina magnézia pro skříň převodovky nebo celá paleta plastů: polyamid, polypropylen, ABS, POM. Konkrétně skříň nůžek je vyrobena z PA6-GF30I, jehož vlastnosti jsou velice citlivé na orientaci skelných vláken, rychlost zatěžování, teplotu a další technologické vlivy. Výzkum spojený s definicí materiálových modelů a vhodných optimalizačních kriterií není předmětem této publikace. 2

Volba typu elementů a jemnosti sítě

První otázkou bylo, které typy elementů jsou použitelné pro naší aplikaci, přičemž jsme se nechtěli omezit jen na v LS-Dyna spolehlivé SHELL a HEXAHEDRAL elementy, ale z praktických důvodu pro nás byly zajímavé i TETRAHEDRAL elementy. Vytvořili jsme sérii jednoduchých modelů, pro které bylo k dispozici i analytické řešení a porovnávaly jsme s ním různé typy elementu a vliv stupně diskretizace. Příklad jednoduchého modelu s velkou citlivostí na kvalitu elementu je vizualizován na Obrázku 2. Na těchto modelech jsme vyhodnocovali různé veličiny. Na Obrázku 3 je zobrazena například deformace a tedy tuhost pro různé elementy. Důležitá byla pochopitelně i výpočtová cena jednotlivých elementů. Při jednoduchém např. tahovém zatížení vykazovaly všechny elementy relativně dobrou spolehlivost. Při ohybovém namáhaní byla citlivost větší. Samozřejmě se ověřilo, že při dostatečně jemné diskretizaci dávají všechny typy elementů relativně přesné výsledky. Velice zajímavé pro svou přesnost byly TETRA elementy typu 4 a 16. U TET4 se ale projevila vysoka citlivost na kvalitu sítě u složitějších modelů.

Obrázek 2 - Základní model pro testovaní elementů


6



Obrázek 3 - Deformační vlastnosti jednotlivých typů elementů

3

Tvorba sítě

Pro reálnou geometrii skříně nůžek byla poté vytvořena FEM síť s TETRA a SHELL elementy (Obrázek 4). HEXAHEDRAL elementy nebyly pro svoji náročnost při generování sítě dále zkoumány. TETRA elementy mají velkou výhodu v geometricky věrném, velice rychlém a automatizovatelném síťovaní, které dostatečně dobře vystihuje vrubové efekty. Shellová sít bude naopak počítat velice rychle, bude správně vystihovat ohybovou tuhost, ale nezobrazí vhodně vrubové efekty, nevystihne věrně tloušťku stěny a její tvorba je podstatně náročnější. Vhodnost jednotlivých variant byla ověřena experimentem.

Obrázek 4 - Síť na skříni s TETRA a SHELL elementy


7


4


Pádová zkouška samostatného dílu

Obrázek 5 - Pádová zkouška samostatného dílu

Samotná skříň akumulátorových nůžek byla na 3 místech osazena tenzometry a jedním senzorem zrychlení a byla hozena z 1 metru na tuhou podložku. Pádová zkouška byla snímána vysokorychlostní kamerou (Obrázek 5). Naměřené hodnoty byly porovnány s FEM výpočty s různými typy elementů. Na porovnání s tenzometrem v dolní části skříně (Obrázek 6) je vidět, že všechny testované elementy vykazují dobrou schodu s experimentem. Když se podrobněji podíváme na deformace a napjatost v místě měření, je zřejmé, že se jedna o globální ohyb kde je přes gradient napětí velký počet elementů a dobrý výsledek je pochopitelný.

Obrázek 6 - Porovnání výsledků z měření s výpočty

V místech, kde byla napjatost podobná spíše ohybu desky, vykazují nejlepší výsledky modely z shell elementů a z kvadratických TET16 elementů. Protože ale SHELL elementy mají




spoustu víše zmíněných praktických nevýhod rozhodli jsme se dále pokračovat pouze s TETRA elementy. Tetrahedron element typu 13 je vlastně vylepšený TET10 téměř bez volumetric lockingu (shear locking ale zůstavá). Tento element jsme jako kompromis vybrali jako základní element pro rozsáhlé simulace. Jako záložní a ve speciálních situacích podstaně přesnější element zůstává 10-ti uzlový kvadratický tetrahedron typu 16, který je ale minimálně 6x výpočtově náročnější a ztrácí svou přesnost při větších deformacích. Ověření použitelnosti tetrahedronové sítě otevírá cestu k automatizaci generování sítí a k opravdu masivnímu nasazení těchto simulací v komerčním použití u komplikovaných geometrií. 5

FEM sestava HSA

Na základě předchozích zjištění byl připraven model celé sestavy akumulátorových nůžek. Všechny vnější a podstatné vnitřní díly byly diskretizovány solidovými tetrahedrony a celý model má 23 částí a asi 800000 elementů (Obrázek 7). Zbylé vnitřní části a spojovací části byly nahrazeny hmotnými body a navázány pomoci RBE elementů. Dopadový fundament je ideálně modelován jako Rigidwall Planar Forces. Pro ostatní kontakty byl použit Automatic_Single_Surface, pro některé pevné vazby TIDE kontakt. Tento model byl opět experimentálně ověřen.

Obrázek 7 - FEM model sestavy

6

Pádová zkouška HAS

Při experimentu bylo ponecháno stejné senzorové osazení a opět byl použit záznam vysokorychlostní kamerou (Obrázek 8). Naměřená zrychlení vykazují opět alespoň při první vlně relativně dobrou shodu s výpočtem (Obrázek 8). Na Obrázku 9 je zobrazeno porovnání s tenzometrem, který je více vzdálen od místa dopadu. Shoda je menší než u pádu samostatného dílu, ale jelikož jde o relativně jednoduchý, universální a přenositelný model, který simuluje poměrně komplikovaný jev, byla shoda považována za dostatečnou.




Obrázek 8 - Pádová zkouška celého stroje a porovnání naměřených zrychlení

Obrázek 9 - Porovnání naměřených přetvoření s výpočtem

7

Postprocesing - globální identifikace kritických mist

S takto verifikovaným modelem resp. metodikou je možné počítat libovolné pádové orientace. Na Obrázku 10 jsou zobrazeny základní pádové orientace, které byli realizovány v tomto projektu.




Obrázek 10 - Pádové orientace

Při postprocessingu v průmyslové praxi je kladen velký důraz na efektivitu, kvalitu a rychlost a proto byla základní metodika doplněna o spoustu automatizovaných kroků. Na Obrázku 11 je např. zobrazeno automaticky generované globální maximum optimalizačního kriteria ze všech časových kroků, které přehledně identifikuje kritická místa a je doplněno časovým průběhem pro vybranou oblast.

Obrázek 11 - Globální maximum z pádové orientace

V dalším kroku je pak automaticky generováno globální maximum pro všechny pádové orientace (Obrázek 12). Takováto vizualizace velice zjednoduší a zpřehlední prezentaci výsledků a projektovému týmu umožní efektivní cestou využít informace. Pokud je někde překročen daný limit optimalizačního kritéria, je možné pomocí zpětné identifikace učit při které pádové orientaci k tomu došlo (Obrázek 12 vpravo nahoře). Tento proces umožní například i definování nejhorších pádových orientací, pro které je pak proveden kontrolní experiment.




Obrázek 12 - Globální maxima ze všech pádových orientací

8

Závěr

Popsaná metodika podporuje v současné době vývoj celé řady produktů, je neustále porovnávána s prováděnými experimenty, s kterými dosahuje vysoké shody, a je dále kalibrována, vylepšována a automatizována. Takovéto simulace, které věrně popisují reálné chování, pomáhají při pochopení a objasnění příčin nežádoucích jevů a umožňují efektivní cestou optimalizaci technických systémů. Dlouhodobě se ukazuje, že při řešení komplexního problému jako jsou pádové a nárazové zkoušky nabývají explicitní výpočtové metody na velikém významu. Reference GANZERT, W., HOLUB, J.: Dynamic analysis of crankshafts in highly loaded hand-held power tools, Setkání Uživatelů, INFO DNY MSC.SOFTWARE 2010 HALLQUIST J. O.: LS-DYNA Theory Manual - Version 971, May 2007 www.stihl.cz




OPTIMALIZACE TALÍŘOVÉ PRUŽINY TURBODMYCHADLA RADEK JANDORA Honeywell, spol. s r.o. – HTS CZ o.z. Abstract: Durability tests of turbochargers with variable geometry identified wear which was caused by insufficient preload of a spring shroud. This preload should fix the nozzle assembly and help to avoid vibrations. Using simulations in Ansys Mechanical, conditions leading to this wear were identified. Then the shroud design was optimized to ensure sufficient preload in all regimes of the turbocharger. Keywords: turbocharger, optimization, spring shroud, thermal load. 1

Úvod

Při testování životnosti turbodmychadel s proměnnou geometrií bylo na lopatkovém věnci nalezeno opotřebení, které mělo takový tvar, že po obvodu plochy, kde se opírá přítlačná talířová pružina, bylo pravidelně rozmístěno dvanáct míst s vysokým opotřebením. Obdobné opotřebení bylo nalezeno i na pružině samotné (viz obrázek 1).

Obrázek 1 - Opotřebení na lopatkovém věnci a přítlačné talířové pružině

Při kontrolním měření a výpočtu vlastních tvarů odlehčené pružiny byl nalezen odpovídající tvar na výrazné vlastní frekvenci 10550 Hz (10723 Hz při výpočtu), viz obrázek 2. Toto je frekvence, která odpovídá buzení vyvolanému rotujícím turbínovým kolem.




a) b) c) Obrázek 2 - Modální analýza přítlačné pružiny: a) výpočet, b) měření vlastního tvaru, c) měření vlastní frekvence

Jako pravděpodobná příčina vzniku opotřebení byla takto odhalena ztráta předepnutí přítlačné pružiny a bylo zahájeno zkoumání, proč ke ztrátě dochází a jak by mu bylo možné zabránit. 2

Tepelně-mechanický výpočet předepnutí přítlačné pružiny

Za účelem rychlého zjištění předepnutí pružiny při provozní teplotě byl vytvořen zjednodušený rovinný model v řezu turbodmychadla. Sestava měla všechny části, které jsou v řetězci předepnutí pružiny, viz obrázek 3.

Obrázek 3 - Model pro tepelně-mechanický výpočet

Při provozní teplotě se pružina dotýká dvou částí, které jsou v naprosto protichůdných podmínkách. Střední skříň je chlazena vodou, aby byla udržována nízká teplota v ložiscích rotoru. Na povrchu skříně tak je cca. 230 °C. Naproti tomu turbínová skříň a lopatkový věnec jsou vystaveny přímo výfukovým plynům a teplota v nich dosahuje až 650 °C. V přítlačné pružině je tak teplotní gradient až 400 °C, viz obrázek 4a. Rozdílná teplota v pružině a rozdílné roztažení na různých průměrech vede k tomu, že pružina se tepelně deformuje a ohýbá směrem od opěrné plochy, viz obrázek 4b. U nezatížené pružiny se vněší průměr vzdálí o 0,37 mm od kontaktní plochy, u zatížené pružiny se to projeví snížením předepnutí o 55 %, a to při nominálních rozměrech. U kombinace nejnižší




pružiny a největší vzdálenosti opěrných ploch může dojít k úplné ztrátě předpětí, viz srovnání s optimalizovaným návrhem na obrázku 12. Tímto tedy byla příčina odlehčení pružiny potvrzena.

a) b) Obrázek 4 - Výsledky výpočtu: a) teploty v pružině, b) tepelné deformace

3

Optimalizace tvaru přítlačné pružiny

Za účelem odstranění ztráty předpětí přítlačné pružiny byla provedena optimalizace jejího tvaru. Možné změny tvaru pružiny byly omezeny návazností na okolní současti (viz obrázek 5): 

Plocha vnitřního disku pružiny (v kontaktu se střední skříní) musí být dostatečná, aby bylo možné namontovat pojistný kroužek zajišťující pružinu.



Horní plochá část pružiny se musí přiblížit až do doporučené vzdálenosti od turbínového kola, aby se minimalizovalo proudění za kolo.



Vnější průměr pružiny je omezen vnitřním průměrem lopatkového věnce.

Měnit se tedy mohly tyto parametry: 

Sklon a velikost šikmé přechodové části z vnitřního disku (v kontaktu se střední skříní) po konturu turbínového kola.



Výška od povrchu vnitřního disku (kontaktní plochy střední skříně) po vnější opěrnou plochu na lopatkovém věnci.



Tloušťka použitého plechu.

Byla provedena studie účinku každého z těchto parametrů a poté optimalizace tvaru z přihlédnutím k působení každého z nich.




Obrázek 5 - Optimalizovaná přítlačná pružina a okolní části

Vliv sklonu šikmé části Sklon šikmé části mění zejména poměr mezi velikostí přechodové části pružiny a plochou částí. Rozdíl je patrný ze dvou profilů v obrázku 6, kde varianta 80° rychle přejde z jedné ploché části do druhé, kdežto varianta 30° vytvoří spíše komolý kužel a plocha za kolem je minimální.

Obrázek 6 - Varianty tvaru podle sklonu: vlevo – 80°, vpravo – 30°

Obrázek 7 - Vliv sklonu šikmé části na předepnutí pružiny

Jak je patrné z grafu na obrázku 7, s rostoucím podílem kuželového profilu významně roste tuhost pružiny za pokojové teploty. Tuhost při zahřátí se ovšem tak významně nemění.




Naopak, zvětšuje se rozdíl mezi předepnutím za pokojové teploty a při vysoké teplotě, což je nežádoucí jev. Vliv výšky vnější opěrné plochy Další zkoumaný parametr, výška vnější opěrné plochy, modifikoval pouze vnější část pružiny. Spolu se změnou tvaru pružiny byl odpovídajícím způsobem změněn i lopatkový věnec, aby bylo dosaženo stejného předepnutí za pokojové teploty (viz obrázek 8).

Obrázek 8 - Varianty tvaru podle výšky vnější opěrné plochy: vlevo – 5.83 mm, vpravo – 3.5 mm

Obrázek 9 - Účinky výšky vnější opěrné plochy: nahoře – předepnutí, vlevo dole – tepelné deformace, vpravo dole – tuhost pružiny

Graf na obrázku 9 nahoře ukazuje, že při snižování výšky vnější přítlačné plochy se zmenšuje rozdíl mezi předepnutím za pokojové teploty a při zahřátí, což je cílem optimalizace.




Podrobnější zkoumání ostatních výsledků odhaluje příčinu tohoto chování. Se zmenšením výšky opěrné plochy se zvětšuje vnější šikmá přehodová část z kontury u kola po opěrnou plochu u lopatkového věnce. Při zahřátí pak roztažnost materiálu působí takovým způsobem, že tepelné deformace ve vnitřní šikmé části způsobují stlačení pružiny. Naopak tepelné deformace vracejí vnější plochu pružiny zpět k nedeformované poloze. Takto se při snížení výšky vnější opěrné plochy sníží tepelné deformace pružiny až na polovinu (viz obrázek 9 vlevo dole). Takováto pružina má ovšem vyšší tuhost (viz obrázek 9 vpravo dole), což je nutné kompenzovat. Vliv tloušťky plechu Tuhost pružiny lze ovlivňovat změnou tloušťky plechu použitého pro její výrobu, přičemž tuhost roste s tloušťkou plechu, viz obrázek 10.

Obrázek 10 - Účinek tloušťky plechu: vlevo – předepnutí, vpravo – tuhost pružiny

4

Výsledný tvar

Při zohlednění všech předchozích závěrů byla provedena optimalizace pružiny. Výsledný tvar má minimální možnou velikost vnitřní šikmé části s takovými rádiusy, aby bylo minimalizováno napětí v těchto místech. Vnější šikmá část byla naopak zvolena co nejdelší, aby byly omezeny tepelné deformace. Tloušťka plechu byla zvolena podle požadovaného rozsahu předepnutí. Optimalizovaný tvar pak má sníženy tepelné deformace o 33 % (viz obrázek 11). Nejvyšší napětí při konfiguraci nejvyšší pružiny a nejmenšího prostoru je srovnatelné vůči původnímu tvaru (charakter napjatosti viz obrázek 12, zatížení obou pružin je pro jiný případ – původní tvar je v nominálních rozměrech, optimalizovaný v největším zatíženích). Rozdíl mezi předepnutím pružiny při pokojové teplotě a při zahřátí byl snížen o 50 %. Celková charakteristika předepnutí pružiny je plošší, protože díky tenčímu plechu byla snížena tuhost pružiny, a proto u konfigurace s nejmenší výškou pružiny v největším prostoru nedochází ke ztrátě kontaktu, ale zůstává stále předepnutí, jaké má původní tvar při zahřátí v nominálních rozměrech.




Obrázek 11 - Tepelné deformace: vlevo – původní tvar, vpravo – optimalizovaný tvar

Obrázek 12 - Redukované napětí: vlevo – původní tvar, vpravo – optimalizovaný tvar. Pozn. Původní tvar je v nominálních rozměrech, optimalizovaný konfiguraci s maximálním stlačením.

Obrázek 13 - Srovnání závislosti obou tvarů na odchylku míry stlačení od nominální hodnoty

5

Závěr

Optimalizace tvaru přítlačné pružiny vedla k mírnému zvýšení předepnutí při nominálních rozměrech při pokojové teplotě o přibližně 10 %. Přitom dochází k mnohem menší ztrátě předepnutí při zahřátí díky tvaru, který lépe kompenzuje rozdílnou roztažnost materiálu. Díky nižší tuhosti pružiny pak je jak menší zatížení u nejvyšší pružiny v nejmenším prostoru, tak




je stále zachováno předepnutí v opačném případě, tedy u nejmenší pružiny v největším prostoru. Oba tvary byly dále srovnány v simulaci tepelného zatížení celého turbodmychadla, kde minimální předepnutí v průběhu celého teplotních cyklu bylo zvýšeno o 60 % (viz obrázek 14).

Obrázek 14 - Výsledky Předepnutí v průběhu tepelného zatěžovacího cyklu Kontaktní adresa: Tuřanka 100/1387, 627 00 Brno, [email protected]




FSI ANALÝZA JEZOVÉ KLAPKOVÉ HRADÍCÍ KONSTRUKCE JIRSÁK V., KANTOR M., NOWAK P. České vysoké učení technické v Praze, Fakulta stavební Abstract: This article deals with FSI analysis of an usual type of a weir flap gate. As a first step numerical model of a free surface fluid flow was successfully verified on a small experimental physical model and then applied on the real structure. Calculated hydrodynamic pressure field from CFD analysis was later used as a boundary condition for structural analysis of the weir flap gate. Keywords: FSI analysis, weir flap gate 1

Úvod

Klapkové hradící konstrukce spadají do kategorie pohyblivých jezů, jejichž hlavní funkcí je vzdutí vody nad jezem na požadovanou úroveň. Pohyblivé jezové konstrukce se většinou umísťují u vodních elektráren, kde nadržují horní vodu na nejvyšší možné úrovni a zajišťují tak maximální možné využití vodní energie. Klapky jsou poklopové jezy s osou otáčení na spodní stavbě jezu a jsou v současné hydrotechnické výstavbě nejpoužívanějším typem. Většinou je celá technologická část jezu (klapka s příslušnými ovládacími a upevňovacími prvky) vyrobena ve výrobním závodě a jako kompletní prefabrikát dopravena na staveniště, kde je přikotvena k vybudované spodní stavbě jezu. Vývoj klapek byl doprovázen poměrně rozsáhlým hydrotechnickým výzkumem, jehož výsledkem je návrh optimálního geometrického tvaru, a to jak z hlediska hydraulického, tak i z hlediska statického. Tímto způsobem vznikla tzv. dutá klapka, jejíž konstrukce je tvořena dvěma válcovými stěnami a tím je dána její vysoká prostorová tuhost a velká únosnost. Cílem této studie bylo provedení detailní FSI (Fluid-Structure Interaction) analýzy daného typu klapkové konstrukce, poukázání na slabá místa konstrukce a návrh vhodných doporučení. 2

Zadání

Úkolem bylo provedení detailní FSI (Fluid-Structure Interaction) analýzy dutého klapkového uzávěru. Byl vybrán běžný typ konstrukce používaný v Čechách. Výpočet byl proveden pro pět vybraných zatěžovacích stavů.




Nejprve byl sestaven numerický model přepadu vody přes jezovou konstrukci. Z tohoto fluidálního modelu bylo převzato vypočtené hydrodynamické tlakové zatížení a to bylo aplikováno na samotnou ocelovou konstrukci klapkového uzávěru pro jednotlivé zatěžovací stavy. Účinek dolní vody nebyl uvažován. Konstrukce klapky je uložena na dvou spodních bočních kulových ložiskách, jednom středním čepovém ložisku a dvou vrchních bočních kulových čepech pro uchycení ovládacích hydromotorů. Hrazená výška je 2,9 m, hrazená šířka jedné klapky je 12 m a celková přelivná šířka je 48 m.

Obrázek 1 - Příčný řez klapkou a spodní stavbou jezu

3

Výpočet hydrodynamického tlakového zatížení

V této práci se věnujeme komplexnímu výpočtu úlohy pomocí FSI numerické analýzy. Pro získání přesného průběhu a celkového rozložení hydrodynamického tlakového zatížení byla provedena numerická analýza proudění v programu ANSYS CFX. Jedná se o stacionární výpočet přepadu vodního paprsku přes těleso klapkového uzávěru, tedy o proudění o volné hladině. Výpočet byl proveden pomocí homogenního vícefázového modelu proudění se standardním typem turbulence k-e. Tento numerický model byl nejprve úspěšně ověřen na zmenšeném fyzikálním modelu a poté aplikován na uvažovanou skutečnou realizaci. Výpočet byl proveden již ověřeným postupem a to pro jednu sklopenou polohu klapky na výšku 0,6H od horní hladiny vody. Tato poloha zhruba odpovídá maximálnímu vyvozenému momentu na klapce. Vypočtené hydrodynamické tlakové zatížení bylo následně použito pro vstupní okrajovou podmínku strukturální analýzy samotné konstrukce klapkového uzávěru.




Obrázek 2 - Porovnání experimentálně naměřených tlaků a vypočítaných tlaků pomocí CFD na zmenšeném experimentálním modelu

Obrázek 3 - Zobrazení rozhraní voda-vzduch včetně zvýraznění oblastí se zvýšeným zavířením na numerickém modelu reálné konstrukce

4

Materiál a okrajové podmínky

Charakteristika materiálu: těleso klapky

ocelový svařenec

boční rámy

ocelový svařenec

Uvažované charakteristiky materiálu

model lineární, elastický, homogenní, isotropní

Modul pružnosti

210 000 MPa

Poissonův součinitel

0,3

Objemová hmotnost

7850 kg/m3



Vlastní hmotnost konstrukce


cca 13 000 kg

Zatěžovací stavy: Vzhledem k omezenému rozsahu analýzy byly vybrány pouze charakteristické zatěžovací stavy. ZS I – Kontrola reakcí - Tento zatěžovací stav slouží pouze pro kontrolu správnosti výpočtu vnitřních sil a reakcí. Jsou uvažovány shodné okrajové podmínky jako v provozním stavu - ZS III, pouze pro zjednodušení výpočtu není uvažována vlastní tíha konstrukce. ZS II – Montážní stav - Klapka je umístěna v montážní, tedy v plně vztyčené poloze. Je uchycena bočními výsuvnými čepy v kulových pouzdrech ještě před přivařením stolice středního ložiska a poloha je zajištěna aretačními táhly. Zatížena je pouze vlastní tíhou. Tento ZS slouží pro zjištění deformace v bodě středního ložiska. ZS III – Provozní stav - Klapka je sklopena do polohy zhruba odpovídající maximálnímu sklopnému momentu, tedy do polohy 0,6H od hladiny horní vody. Je zatížena jak vlastní tíhou tak i hydrodynamickým tlakovým zatížením převzatým z předchozího odpovídajícího numerického fluidálního modelu přepadu vody přes danou jezovou klapku. Uvažována jsou všechna uložení – dvě spodní boční kulová ložiska, jedno střední čepové ložisko a dva vrchní boční kulové čepy pro uchycení ovládacích hydromotorů. Vliv dolní vody není uvažován. ZS IV – Provozní stav + pokles střední ložiskové stolice - Tento ZS vychází z předchozího ZSII. Zde je navíc uvažován posun střední ložiskové stolice o 1 mm ve směrech proti osám Y a Z. Tyto posuny zhruba odpovídají deformaci v montážní poloze bez nadlehčení středního ložiska. ZS V – Havarijní stav (1 hydromotor) - V tomto ZS se uvažuje havárie jednoho z hydromotorů. Celé zatížení (stejné jako v ZSII) je tedy převzato jediným hydromotorem a spodními ložisky. MESH a výpočet Vygenerování sítě, nastavení okrajových podmínek a samotný výpočet byl proveden v programu ANSYS Workbench Mechanical 12.1. Pro vytvoření sítě konečných prvků byly použity prostorové prvky ve tvaru šestistěnu a čtyřstěnu typu Solid186 a Solid187. Kontakty stykových ploch mezi jednotlivými dílci byly modelovány kontaktními prvky typu Conta174 a Targe170. Kontakty v místech všech čepů a uchycení hydromotorů byly modelovány pomocí kloubových kontaktních prvků MPC184.




Obrázek 4 - Celkový pohled na výpočetní síť z návodní strany

Obrázek 5 - Detail výpočetní sítě

Obrázek 6 - Vektorové zobrazení hydrodynamického tlakového zatížení převzatého z CFD analýzy

5

Výsledky




Pro první zatěžovací stav (ZSI – Kontrola reakcí) byla vybrána částečně sklopená poloha klapky, pro kterou bylo zjištěno zatížení hydrodynamickým tlakem. Na základě vypočtených reakcí a hydrodynamického tlakového zatížení byla úspěšně provedena kontrola vypočtených hodnot na základě součtové i momentové podmínky. V druhém zatěžovacím stavu (ZSII – Montážní stav), který odpovídá fázi montáže po vysunutí bočních ložiskových čepů a zaaretování v plně vztyčené poloze před montáží středního ložiska byla zjištěna deformace celé konstrukce v místě středního ložiska pro pozdější posouzení vlivu poklesu ložiskové stolice na napětí v konstrukci. V tomto místě byla vypočtena výsledná totální deformace (vektorový součet všech tří složek deformace) o hodnotě 0,85 mm. Třetí zatěžovací stav (ZSIII – Provozní stav) odpovídá částečnému sklopení ocelové klapky, která je zatížena hydrodynamickým tlakem přepadajícího paprsku. V této poloze vyvozuje klapka maximální sklopný moment vzhledem k ose otáčení. Vypočtená maximální napětí jsou v tomto ZS na střední příčné výztuze. V ostatních místech konstrukce jsou napětí (von-Mises) poměrně nízká a nepřevyšují cca 130 MPa. Maximálního napětí je dosaženo na střední příčné výztuze v místě vybrání u spoje hlavního hradícího plechu a dolního krycího plechu ložisek. Zde je dosaženo maximálního tahového i tlakového napětí cca 230 MPa a ekvivalentní napětí von-Mises zde dosahuje hodnoty cca 225 MPa (špička cca 270 MPa). K maximální deformaci dochází na horním okraji hlavního výztužného plechu v rovině symetrie a to o hodnotě totální deformace cca 9,2 mm. Vypočtená deformace ve svislém směru je zde cca 8,5 mm.

Obrázek 7 - ZSIII - celkový pohled (po vodě) na rozložení ekvivalentního napětí von-Mises – max cca 90 MPa, v místě uchycení hydromotorů max cca 130 MPa

Obrázek 8 - ZSIII - celkový pohled (proti vodě) na rozložení ekvivalentního napětí von-Mises – na bočních příčných výztuhách max cca 110 MPa




Obrázek 9 - ZSIII - Detail maxima napětí (von-Mises) – max špička cca 270 MPa, v těsném okolí cca 220 MPa.

Obrázek 10 - Zobrazení rozložení totální deformace – max. 9,2 mm.

Ve čtvrtém zatěžovacím stavu (ZSIV - Provozní stav + pokles střední ložiskové stolice) se vypočtené hodnoty napětí i deformací se v tomto ZS oproti předchozímu ZSIII lehce snížila. Je to dáno poklesem podpory ve směru deformace konstrukce a tím je způsobeno odlehčení napětí v konstrukci. Naopak se zde nepatrně zvýšila hodnota reakce v uchycení střední ložiskové podpory. V posledním pátém zatěžovacím stavu (ZS V – Havarijní stav - 1 hydromotor) se uvažuje havárie jednoho z hydromotorů. Uložení konstrukce vychází ze ZS III a je pouze odebrán jeden boční hydromotor. Celé zatížení je tedy převzato jediným hydromotorem a spodními ložisky. Napětí (von-Mises) na hradícím plechu a střední příčné výztuze se oproti ZSIII nijak zásadně nezvýšila. Zvýšila se napětí na hlavním podélném výztužném plechu v prvním poli na hodnotu cca 100 MPa a k výraznému navýšení napětí došlo v hlavním příčném výztužném plechu u uchycení hydromotoru. Zde je dosaženo maximálního napětí cca 330 MPa. Špička napětí je pak na výztuhách prstence uchycení bočního čepu o hodnotě cca 600 MPa. Vypočítaná deformace hlavního hradícího plechu na straně hydromotoru 2 mm v příčném směru je poměrně zanedbatelná a odpadá tak obava zpříčení (opření) konstrukce klapky o boční pilíř.




Obrázek 11 - ZSV - celkový pohled (po vodě) na rozložení ekvivalentního napětí von-Mises – max cca 90 MPa, v místě uchycení hydromotorů max cca 210 MPa

Obrázek 12 - ZSV - celkový pohled (proti vodě) na rozložení ekvivalentního napětí von-Mises – na střední příčné výztuze max cca 250 MPa, v prvním poli na výztužném plechu max cca 100 MPa.

Obrázek 13 - ZSV - Deformace v příčném směru (směr osy spodních ložisek) – maximum na hradícím plechu cca 2 mm.

6

Závěr

Uvedená analýza ilustruje využití možností stávajících SW nástrojů v oboru hydrotechniky a to jak v oblasti fluidálních tak i napěťových výpočtů včetně jejich vzájemného propojení. Pro konkrétní zadání tvaru klapky je možné i na základě těchto výsledků shrnout následující poznatky. Výsledky fluidální analýzy se velmi dobře shodují s výsledky naměřenými na modelech fyzikálních. Matematický model umožňuje velmi rychlé a přesné zjištění zatížení tvarově složitých hradících konstrukcí od hydrodynamického tlaku pro různé průtoky, resp. přepadové výšky. Jednoduchý přenos výsledků fluidálních výpočtů (tlakového zatížení na návodní hradící konstrukci) do modelu strukturálního je zásadní výhodou pro oblast výpočtů FSI. Strukturální analýza dokázala, že zadaná konstrukce jezové klapky je dostatečně tuhá a deformace jsou s ohledem na provozní podmínky vyhovující. Průběhy napětí v konstrukci ukazují na místa, ve kterých je vhodné konstrukci vyztužit. To se týká zejména střední diafragmy a jejího napojení na spodní ložisko za normálních provozních podmínek a uložení čepu bočního ložiska a hydromotoru za havarijního zatěžovacího stavu.




Uvedený příklad použití FSI velmi dobře dokumentuje flexibilitu a vhodnost tohoto přístupu pro navrhování a optimalizaci hydrotechnických konstrukcí. Tento dvou-oborový problém (hydrodynamika a napjatostní analýza) lze velmi efektivně řešit v jednom komplexním prostředí. Reference BROŽA V., SATRAPA L.: Hydrotechnické stavby 10 – přehrady, pp.78-81, Praha ČVUT, 2000. ČIHÁK F., MEDŘICKÝ V.: Hydrotechnické stavby 20 – navrhování jezů, pp. 96-104, Praha ČVUT, 2001. STROBL T.: Skriptum – Grundkurs “Wasserbau und Wasserwirtschaft“, pp. 127-129, TU München, 2005. ANSYS,Inc., "Release 12.1 Documentation for ANSYS." SYS IP, Inc 2009 Poděkování Tato práce byla podporována a vznikla v rámci projektu Studentské grantové soutěže ČVUT - SGS10/144/OHK1/2T/11 “Modelování proudění a pevnostní analýza jezových uzávěrů“. Kontaktní adresa: Thákurova 7, Praha 6




METODIKA FSI APLIKACE PŘI SIMULACI VYNUCENÉHO KMITÁNÍ LOPATEK PARNÍ TURBÍNY J. JURENKA1, D. GALLOVIČ1, T. MÍŠEK2 1

TechSoft Engineering s.r.o., 2ŠKODA POWER s.r.o.

Abstract: Článek se zabývá problematikou propojení pevnostních analýz (Finite Element Analysis - FEA) se simulacemi proudění tekutin (Computational Fluid Dynamics - CFD) - Fluid Structure Interaction (FSI) - v rámci programového aparátu ANSYS. Je popsána metodika přenosu silových účinků vypočtených pomocí ANSYS CFX/Fluent vznikajících při obtékání aerodynamických profilů lopatek parních turbín do pevnostní, resp. dynamické analýzy lopatek provedené v ANSYS Mechanical APDL. Keywords: FEA, CFD, FSI. 1

Úvod

V inženýrské praxi se poměrně často setkáváme s případy konstrukcí, které jsou namáhány proměnlivými aerodynamickými silami vznikajícími při jejich obtékání okolním médiem. Mezi tyto konstrukce lze zařadit např. lopatky parních a spalovacích turbín, kompresorů nebo větrných elektráren, potrubní trasy, budovy, letecké profily atd. Při pevnostním návrhu nebo kontrole těchto konstrukcí se často vychází z odhadů působících aerodynamických sil, přičemž tyto odhady jsou často spojeny pouze s nominálním provozem daného zařízení. Moderní výpočetní programy dovolují zpřesnit pevnostní analýzy využitím výsledků CFD výpočtů v rámci FEA analýz. Na Obrázku 1 je vidět příklad provázání CFD a FEA analýz.

Obrázek 1 - Provázání CFD a teplotně-pevnostní analýzy parní turbíny s vnitřním chlazením pomocí programů ANSYS CFX a ANSYS Mechanical.

2

Metodika FSI v rámci programového aparátu ANSYS




Využití metodiky FSI je většinou podmíněno jak způsobem implementace (uživatelskou náročností) tohoto výpočetního postupu v daném programovém aparátu, který určuje míru efektivity práce, tak i znalostmi a zkušenostmi uživatele s problematikou CFD a FEA analýz. Prostředí ANSYS Workbench (WB) v sobě integruje jak výpočetní moduly pro provádění CFD výpočtů (ANSYS CFX), tak i moduly pro pevnostní výpočty (ANSYS Mechanical atd.). Díky tomuto provázání je možné poměrně efektivně přenášet vstupy, resp. výstupy mezi jednotlivými typy analýz (CFD <-> FEA). V případě komplexních výpočetních simulací, jako jsou např. nestacionární CFD výpočty na jedné straně a výpočty dynamické odezvy konstrukcí na straně druhé, standardní prostředky ANSYS Workbench nejsou schopny poskytnout dostatečnou efektivitu prováděných simulací a je vhodné je doplnit o vlastní programové skripty. a. Základní FSI simulace V případě stacionárního, nebo kvazistacionárního způsobu namáhání konstrukcí od obtékajícího média je možné sestavit FSI analýzu pomocí standardních nástrojů dostupných v programu ANSYS Workbench, přičemž se může jednat o tzv. 1-way (Obrázek 2) nebo 2-way (Obrázek 3) metodiku výpočtu. V případě 1-way přístupu je nejprve provedena simulace proudění (CFD) a její výsledky (tlaky, teploty) jsou jako vstupy přeneseny do pevnostní analýzy (FEA). Tento přenos probíhá na základě interpolace tlaků z CFD výpočetní sítě na povrchové elementy (SURF154 - v rámci ANSYS WB vytvářeny automaticky) FEA modelu. V případě 2way přístupu se vytváří uzavřená výpočetní smyčka mezi CFD a FEA simulacemi. Výsledky z CFD výpočtu jsou přenášeny do FEA analýzy, kde je vypočtena odezva konstrukce na dané zatížení a výsledné pole posuvů je zpětně přeneseno do CFD analýzy, kde je využito pro aktualizaci tvaru obtékané součásti. Tento postup je možné opakovat, dokud nedojde k saturaci odezvy konstrukce apod. V uvedených základních typech FSI simulací je implicitně zahrnut předpoklad, že přenos zatížení z CFD do FEA simulace probíhá vždy na konstrukcích s topologicky totožnými hranicemi, tzn. tvar konstrukce obtékané v dané CFD simulaci je totožný s tvarem konstrukce ve FEA analýze, na který je následně přenášeno zatížení.

Obrázek 2 - Schéma 1-way přístupu metodiky FSI.




Obrázek 3 - Schéma 2-way přístupu metodiky FSI.

Ukázka metodiky přenosu tlakového zatížení v prostředí ANSYS Workbench pro jednoduché statické zatížení a topologicky totožné CFD a FEA modely je vidět na Obrázku 4. V programu ANSYS Fluent byl proveden stacionární CFD výpočet, jehož výsledkem jsou data uložená v souboru jobname.res. Pokud souřadnice CFD modelu (ANSYS Fluent) a FEA modelu (ANSYS Mechanical) nejsou totožné s ohledem na globální souřadný systém, je nutné provést prostorovou transformaci FEA modelu, tak aby oba modely byly topologicky totožné. Dalším krokem je definice párových obtékaných ploch pro přenos tlaků a to jak v CFD, tak i FEA modelu na základě dříve definovaných množin. Následně lze provést import tlakového zatížení na FEA model aktivováním příkazu "Import Load". Výsledkem je rozložení statického tlaku na dříve definovaných plochách, resp. elementech ve FEA modelu. Je zřejmé, že přenesené tlakové rozložení bude díky interpolačnímu algoritmu záviset na relativní hustotě výpočetních sítí obou modelů. Proto je nutné především v blízkosti náběžné a odtokové hrany FEA modelu patřičně zjemnit výpočetní síť. b. Komplexní FSI simulace V technické praxi se vyskytují případy konstrukcí, kde postupy uvedené v předchozí kapitole nelze přímočaře použít. Důvodem může být jak komplexní nestacionární způsob namáhání konstrukcí při jejich obtékání, tak i např. odlišná topologie modelu použitého pro CFD, resp. FEA analýzu. Příkladem může být oběžná lopatka, jejíž tvar se mění v závislosti na otáčkách rotoru, a která je namáhána proudovým polem za rozváděcími lopatkami statoru. V případech časově proměnlivého namáhání konstrukcí vycházejícího z jejich obtékání daným médiem vede použití metodiky FSI k simulacím vynuceného kmitání těchto konstrukcí, kdy budící účinky odpovídají právě časovým (nestacionárním) průběhům tlakového zatížení. Budící účinky lze v tomto případě získat nestacionární CFD simulací, jejímž výsledkem jsou obecné časové průběhy tlaků v jednotlivých elementech obtékaných povrchů konstrukce. S ohledem na složitost, resp. obecnost získaných časových průběhů je možné volit různé způsoby jejich dalšího zpracování v rámci FEA analýzy: 

V případě zcela obecných průběhů budících účinků lze provést jak statickou, tak i dynamickou implicitní, resp. explicitní analýzu, přičemž budící účinky jsou uvažovány




jako časové průběhy tlaků vypočtené při CDF simulaci a jsou např. pomocí tabulkového zápisu přímo načteny do FEA modelu. 

Pokud lze získané průběhy budících účinků aproximovat např. pomocí Fourierovy řady, je možné provést harmonickou analýzu vynuceného kmitání. Výpočet vynuceného kmitání je v tomto případě nutné provést samostatně pro každý harmonický člen Fourierova rozvoje a výsledná odezva konstrukce je dána superpozicí polí posuvů, resp. napětí a deformací poplatným jednotlivým členům Fourierova rozvoje.

Obrázek 4 - 1 way coupling: Import zatížení v prostředí ANSYS Workbench.

Pro efektivní použití metodiky FSI v případě analýz s časově proměnlivým zatížením a komplexní geometrií je vhodné doplnit standardní nástroje prostředí ANSYS vlastními programovými skripty a to jak pro zpracování nestacionárních CFD výpočtů v rámci např. postprocesoru CFD-Post, tak i pro úpravu a aplikaci budících účinků v rámci FEA analýzy, vyhodnocení získaných výsledků a kontrolu průběhu simulace. Schéma uvedené metodiky je vidět na Obrázku 5.




Obrázek 5 - Schéma metodiky výpočtu - nestacionární CFD + harmonická analýza vynuceného kmitání FEA.

c. Aplikace FSI na lopatku parní turbíny Lopatky nízkotlakého stupně parních turbín jsou typickou součástí, která je namáhaná časově proměnlivým tlakovým polem. Komplexní způsob namáhání a topologická odlišnost CFD a FEA modelu (vycházející z procesu návrhu) vyžadují použití vlastních programových skriptů pro efektivní definování budících účinků a provedení harmonické analýzy. Příklad přípustné topologické odlišnosti CFD a FEA modelu je vidět na Obrázku 6. Z obrázku je zřejmé, že topologické rozdíly mezi modely musí být relativně malé vzhledem k celkovým rozměrům modelů. Dostatečně přesný přenos tlakového zatížení lze potom docílit postupem, který je založen na vhodném rozdělení obtékaných ploch. Následná interpolace tlakového zatížení probíhá pro jednotlivé plochy nezávisle a tím lze eliminovat chyby vycházejí ze špatné identifikace párových ploch. Dalšího zpřesnění budících účinků FEA modelu lze dosáhnout vynucením normálového směru přeneseného tlakového namáhání vzhledem k povrchu FEA modelu. V případě topologicky odlišných modelů je tlakové namáhání přenášeno s orientací odpovídající normálám CFD povrchu, jejichž směr nemusí být ve FEA modelu stejný.




Obrázek 6 - Přípustných topologických rozdílů CFD a FEA modelů (zelené a tmavomodré objemy FEA model, světle modré objemy - CFD model).

Metodika přenosu v případě nestacionárního buzení je poplatná schématu uvedeném na Obrázku 4. Konkrétní realizace uvedené metodiky v programu ANSYS s využitím programu CFD-Post pro zpracování nestacionárního CFD výpočtu a programu MATLAB pro zpracování a přípravu budících účinků FEA modelu je vidět na Obrázku 7. Vzhledem ke způsobu zatěžování uvedené lopatky byl v prvním kroku proveden nestacionární CFD výpočet, jehož výsledkem je časově proměnlivá tlaková distribuce na obtékaných plochách. Výpočetní sít CFD modelu je vidět na Obrázku 8. Během CFD výpočtu v programu ANSYS CFX bylo uloženo celkem přibližně 800 časových okamžiků (každý krok uložen do samostatného souboru typu *.trn), ze kterých bylo následně rekonstruováno buzení FEA modelu. FEA model lopatky byl vytvořen v ANSYS APDL. Byly použity kvadratické elementy SOLID186. Vzhledem k následné harmonické analýze byl uvažován elastický lineární model materiálu. Při tvorbě výpočetní sítě byly zanedbány některé konstrukční detaily především v oblasti stromečkového zámku a oblasti bandáží, které nemají významný vliv na vynucené kmitání lopatky a umožní efektivní výpočet a zpracování dosažených výsledků harmonické analýzy. Na povrchu obtékaných ploch byly definovány povrchové elementy SURF154, které byly následně exportovány do souboru typu *.cdb. V radiálních řezech segmentu disku a oblastech bandáží byly definovány podmínky cyklické symetrie. Výsledný FEA model je vidět na Obrázku 9. Zpracování výsledků nestacionárního CFD výpočtu bylo provedeno pomocí programu CFD-Post. Vzhledem velkému množství zpracovávaných časových okamžiků byl v rámci CFDPost vytvořen vlastní programový skript, pomoci kterého jsou postupně načítány jednotlivé soubory *.trn a výsledky v nich uložené jsou interpolovány na povrchové elementy definované souborem *.cdb. Pro každý časový okamžik je vytvořen soubor typu *.sfe, ve kterém jsou uloženy aktuální hodnoty tlakového namáhání jednotlivých povrchových elementů FEA modelu.




Příklad proměnlivého tlakového namáhání povrchových elementů FEA modelu je vidět na Obrázku 10. Soubory typu *.sfe jsou v následujícím kroku zpracovány pomocí programu MATLAB. Je provedena Fourierova transformace a jednotlivé průběhy jsou aproximovány pomocí Fourierovy řady. Počet uvažovaných harmonických složek je dán složitostí vstupního signálu. Příklad aproximace jedné periody budících tlakových účinků pro tři povrchové elementy FEA model je vidět na Obrázku 11. Jednotlivé harmonické složky výsledného buzení jsou následně uloženy v podobě reálných a imaginárních složek tlakového působení v jednotlivých povrchových elementech FEA modelu do ASCII souborů, které mohou být načteny přímo do ANSYS APDL. Porovnání tlakové distribuce na obtékaných plochách lopatky vypočtené pomocí CFD simulace a distribuce tlakového namáhání aplikovaného ve FEA modelu, které vzniklo výše popsaným způsobem pro daný časový okamžik je vidět na Obrázku 12. Harmonická analýzy vynuceného kmitání lopatky byla provedena v programu ANSYS APDL. Vzhledem ke komplexnosti jak okrajových podmínek a zatížení, tak i následného zpracování dosažených výsledků je celá analýza řízena vlastními APDL skripty. V prvním kroku simulace je proveden statický výpočet předpětí lopatky od rotace a konstantního tlakového namáhání vycházejícím z Fourierova rozvoje. Na statickou analýzu navazuje harmonická analýza s uvažování předpětí, podmínek cyklické symetrie a necyklicky symetrickým zatěžováním jednotlivých segmentů, které odpovídá danému harmonickému indexu a dané (jedné) harmonické složce buzení. Vzhledem k tomu, že ANSYS nepodporuje zpracování výsledků cyklicky symetrických modelů pomocí tzv. loadcase, je nutné harmonickou analýzu opakovat pro každou harmonickou složku buzení. Díky lineární povaze dosažených výsledků je následně možné výsledky jednotlivých harmonických analýz sečíst a tím získat výsledné vynucené kmitání lopatkového kola, resp. reálné a imaginární hodnoty pole posuvů, napětí a deformací. Jednoduchým přepočtem je možné získat časové průběhy (kmitavých složek) uvedených veličin. Superpozicí kmitavých složek s výsledky statické analýzy je možné získat celkové hodnoty uvedených veličin.

Obrázek 7 - Realizace FSI metodiky v programu ANSYS.





Obrázek 9 - MKP výpočetní síť.

Tlak 

Tlak 

Obrázek 8 - CFD výpočetní síť.


Čas 

Obrázek 10 - Příklad časového průběhu tlakového působení na povrchové elementy FEA modelu.


Čas  Obrázek 11 - Příklad aproximace budících účinků pomocí Fourierovy řady pro tři povrchové elementy FEA modelu.


38



Obrázek 12 - Porovnání distribuce tlaku CFD model vs. FEA model po přenosu v daném časovém okamžiku (barevné škály jsou totožné).

Výsledky vynuceného kmitání je možné po provedení všech harmonických analýz zobrazit v závislosti na čase na jedné lopatce nebo na celém lopatkovém kole, případně vykreslit závislost dané veličin na čase pro zvolený uzel výpočetní sítě. Zpracování a vykreslení výsledků harmonických analýz je díky časové náročnosti a komplexnosti řízeno vlastními programovými APDL skripty. Na Obrázku 13 je zobrazena axiální složka pole posuvů osmi lopatek v daném časovém okamžiku. Maximální hodnoty jsou dosaženy na odtokové hraně lopatky jak je naznačeno ve spodní části Obrázku 13. Obdobně lze zobrazit jednotlivé složky pole napětí, resp. deformací, případně HMH napětí apod. 3

Závěr

Metodika FSI patří mezi komplexní analýzy technické praxe. Její zvládnutí vyžaduje jak dostatečný výkon výpočetních prostředků, tak i teoretické znalosti z oblasti CFD, resp. FEA analýz a pokročilé uživatelské znalosti prostředí ANSYS CFX/Fleunt, resp. ANSYS Mechanical/APDL. V prezentovaném článku jsou popsány jak standardní způsoby aplikace metodiky FSI, tak i přístupy použití FSI v případě nestacionárního tlakového zatížení vedoucího na simulaci vynuceného kmitání. Obecně je možné uvedené poznatky shrnout do následujících bodů: 

Použití metodiky FSI je poměrně náročné, a proto by její aplikace na reálné konstrukce měla být podložena relevantními důvody, jako jsou např. významné deformace konstrukce indukované tlakovým zatížením, nebo blízkost frekvence budících účinků a vlastní frekvence konstrukce.



Přístupy při aplikaci metodiky FSI lze rozdělit do dvou skupin vzhledem k stacionární, resp. nestacionární povaze buzení a tedy i způsobu CFD výpočtu. Výsledky stacionárního CFD výpočtu lze použít jako zatížení při statické FEA analýze. Výsledky




nestacionárního CFD výpočtu lze pomocí postupů pospaných v prezentovaném článku aplikovat jako buzení při dynamické transientní nebo harmonické analýze konstrukce.

Obrázek 13 - Axiální složka pole posuvů v daném časovém okamžiku.



V případě stacionárního CFD výpočtu a následné statické FE analýzy lze pro přenos zatížení z CFD do FEA efektivně využít standardní nástroje prostředí ANSYS Workbench. V případě nestacionárního CFD výpočtu a následné např. harmonické FE analýzy je vhodné tyto standardní nástroje doplnit vlastními programovými skripty a tím celý výpočetní postup zefektivnit.






Poměrně robustní interpolační algoritmus implementovaný v prostředí ANSYS selhává, pokud není splněn implicitní předpoklad, že modely CFD a FEA jsou topologicky totožné. Další oblastí, kde není zcela zaručena spolehlivost přenosu, jsou plochy s velkou relativní křivostí (náběžné a odtokové hrany). V případě, že topologická odlišnost CFD a FEA modelů je relativně malá v porovnání s celkovými rozměry je možné dostatečně přesný přenos zatížení (FSI) zajistit vhodným rozdělením obtékaných ploch a případným zjemněním FEA výpočetní sítě v místech s vyšší křivostí.

Literatura Manuál ANSYS R13 Poděkování Tento příspěvek vznikl za laskavého svolení ŠKODA POWER s.r.o. Contact address: TechSoft Engineering s.r.o., Na Pankráci 26, 140 00 Praha 4 ŠKODA POWER s.r.o., Tylova 1/57, 301 00 Plzeň




SOFTWARE ANSYS PRO NÁVRH A OPTIMALIZACI ELEKTRICKÝCH STROJŮ A ZAŘÍZENÍ, MOŽNOSTI MULTIFYZIKÁLNÍCH ANALÝZ JAKUB HROMÁDKA, JINDŘICH KUBÁK Techsoft Engineering spol. s.r.o. Abstrakt: Uvádíme možnosti využití software společnosti ANSYS při návrhu elektromotoru pro elektromobil včetně jeho tepelné analýzy. Zmíněny jsou také možnosti multifyzikálních analýz navazujících na elektromagnetické výpočty. Keywords: ANSYS Maxwell 14.0, ANSYS Mechanical, electromagnetic field, electromagnetic losses, multiphysical analyses, thermal computation. 4

Úvod

V posledních letech se staly v technické praxi zcela běžnými návrhy s pomocí pevnostních výpočtů a analýz proudění až po optimalizaci. S nimi můžeme navrhnout zařízení, která jsou levnější a zároveň dosahují vyšších technických parametrů a umožňují tak získat významnou konkurenční výhodu. Další oblastí, ve které nám mohou technické výpočty výrazně pomoci při konstruování, je kategorie elektrických zařízení jako jsou např. elektromotory, transformátory a celá oblast spotřební elektrotechniky. 5

Software ANSYS Maxwell 14.0 a jeho moduly RMxprt, PExprt

Pro analýzu nejrůznějších elektromagnetických problémů vyvíjí společnost ANSYS, Inc. (USA) program ANSYS Maxwell v nízkofrekvenční oblasti (do 3 kHz) a pro vysokofrekvenční oblast nabízí program ANSYS HFFS. Doplnila tak řadu svých produktů, které nabízejí komplexní řešení problémů v oblasti proudění a pevnostních analýz o analýzy elektromagnetické. ANSYS Maxwell využívá metod MKP (metoda konečných prvků) pro nalezení řešení úplné soustavy Maxwellových rovnic. Konstruktér tak dostává do rukou silný nástroj pro návrh a následnou optimalizaci elektrických strojů a zařízení. a. Popis programu a jeho možnosti Práce s programem ANSYS Maxwell 14.0 je díky grafickému uživatelskému rozhraní velice přátelská a intuitivní. Uživateli je při spuštění nabídnuto 6 typů řešičů, každý přizpůsobený pro efektivní výpočet jiného druhu problému. Je možné provádět analýzy stacionární i nestacionární, analýzy s harmonickými zdroji elektromagnetických polí nebo analýzy frekvenční. Ve výpočtu jsou zahrnuta i elektromagnetická pole indukovaná mechanickým pohybem částí stroje. Program má implementovány inteligentní algoritmy pro návrh a optimalizaci výpočetní sítě (autoadaptive meshing). V několika krocích je síť optimalizována v místech největší chyby řešení; ta je počítána z podmínky nevířivosti magnetického pole. Program dokáže pracovat i s materiály nelineárními s definovanou BH závislostí, s materiály laminovanými složenými s izolovaných plechů, s materiály anizotropními, nehomogenními a teplotně závislými. Proces geometrického návrhu elektrického stroje usnadňují dva expertní moduly: RMxprt - pro točivé stroje, PExprt - pro transformátory. Na základě analytických vztahů jsou provedeny rychlé návrhové výpočty nového stroje, uživatel má během několika sekund základní informace o http://aum.svsfem.cz



svém výrobku. Navržená geometrie se jednoduše exportuje do ANSYS Maxwellu, kde je proveden kompletní FEA výpočet. ANSYS Maxwell také disponuje vlastními optimalizačními nástroji - hledání extrémů cílové funkce, parametrické studie nebo citlivostní studie. Optimalizace tak umožní navrhnout zařízení s vylepšenými technickými parametry při požadavku nízkých výrobních nákladů. b. Elektromagnetismus a multifyzikální analýzy Svět, ve kterém žijeme, je multifyzikální. Fyzikální jevy spolu navzájem souvisí a pro požadovanou přesnost prováděné analýzy není možné je studovat zvlášť, při vývoji nového produktu je potřeba vzít v úvahu všechny fyzikální souvislosti. V případě elektrických strojů se multifyzikální analýza projevuje kombinací elektromagnetických polí, polí teplotních a silových účinků na konstrukci. Praktické využití multifyzikálních analýz nacházíme např. ve velmi aktuální oblasti návrhu elektromotorů pro elektromobily. Takový elektromotor musí být konstruován s důrazem na nízkou hmotnost při dostatečném výkonu a nízkém oteplení. i.

Teplotní analýzy

Průchodem elektrického proudu vodičem vznikají ohmické ztráty, které způsobují zahřívání elektrického stroje. K zahřívání přispívají také hysterezní ztráty vznikající v důsledku přemagnetovávání feromagnetik. Ve vodivých částech elektrických strojů v proměnných magnetických polích jsou generovány i vířivé proudy, které jsou dalším zdrojem ohřevu. Prostředí ANSYS Workbench podporuje vzájemné provázání mezi elektromagnetickými výpočty v ANSYS Maxwellu a teplotními výpočty v ANSYS Mechanical nebo ANSYS CFD (Fluent, CFX). Nejdříve je nutné vytvořit v obou programech geometrický model, který má totožné alespoň díly, u kterých se počítají elektromagnetické ztráty. V dalším kroku se provede výpočet elektromagnetických polí. Výše popsané druhy elektromagnetických ztrát mohou být jednoduchým postupem díky ANSYS Workbench importovány na geometrii v ANSYS Mechanical/CFD. Importovaná data lze využít k deklaraci objemových zdrojů tepla nebo tepelných toků. Specifikováním tepelných materiálových vlastností a okrajových podmínek je možné v ANSYS Mechanical/CFD spustit výpočet teplotních polí. Získáváme rozložení teplot na vytvořeném modelu v ustáleném stavu. V ANSYS Mechanical se řeší prosté vedení tepla v pevné struktuře. V případě výpočtu v ANSYS Fluent je formulován výpočet včetně proudění tekutin v okolí elektromagneticky zahřátých těles (např. ohřev vzduchu v okolí cívky s elektrickým proudem na obrázku) a jejich využití pro chlazení daného zařízení. Elektromagnetické vlastnosti materiálů mohou být teplotně závislé (např. elektrická vodivost). ANSYS Workbench dokáže zpětně přenést vypočtené teplotní pole na model v ANSYS Maxwellu a spustit elektromagnetický výpočet s aktualizovanými okrajovými podmínkami. Výsledkem je nové rozložení ztrát, které může být opět použito k tepelnému výpočtu. V několika iteracích získáváme ustálený stav pro teplotně závislé materiály. ii.

Strukturální výpočty

Elektromagnetická pole mohou být zdrojem vzájemného silového působení mezi jednotlivými díly navrhovaného stroje. Silová pole dokáže ANSYS Maxwell kvantifikovat a díky ANSYS Workbench dokáže hustotu silového působení jednoduše importovat do ANSYS Mechanical. Zde jsou importovaná data užita jako zatížení strukturálních výpočtů. Výsledkem získáváme deformace a napětí v důsledku elektromagnetického silového působení a nebo teplotní roztažnosti materiálu. http://aum.svsfem.cz


iii.


Mechanický pohyb

Samotný program ANSYS Maxwell dokáže zpracovat dva druhy energie. Energii mechanickou (mechanický pohyb komponent elektrického stroje) a energii elektromagnetickou. Elektromagnetická pole generují silová působení na jednotlivé součásti modelu, důsledkem je mechanický pohyb pohyblivých částí. Změna konfigurace zpětně ovlivní rozložení elektromagnetických polí. Navíc mohou být při mechanickém pohybu ve vodivých částech modelu generovány vířivé proudy způsobující ohřev. Všechny tyto interakce dokáže ANSYS Maxwell ve výpočtech zahrnout. 6

Proces návrhu trakčního elektromotoru pro elektromobil

Výše popsaný postup multifyzikálních analýz můžeme prakticky demonstrovat na velmi aktuálním příkladu návrhu trakčního elektromotoru pro elektromobil. Požadavky na technické parametry takového elektromotoru jsou velmi vysoké. Motor musí zajistit dostatečný výkon pro různé jízdní režimy. Projektovaný motor musí mít minimální hmotnost a zároveň musí být schopen vyvinout potřebný kroutící moment. Elektromagnetické ztráty musí být minimalizovány, neboť baterie napájející elektromobil mají omezenou kapacitu a cílem je maximalizovat dojezd elektromobilu. Elektromagnetické ztráty navíc způsobují ohřev motoru. Pro různé jízdní režimy nesmí teplota motoru přesáhnout maximální přípustnou teplotu pro danou tepelnou třídu izolace. Návrh elektromotoru, výpočet elektromagnetických ztát s následnou tepelnou analýzou jsou ukázkovým příkladem využití multifyzikálních možností softwaru ANSYS. Tabulka 1 Požadované parametry motoru

Cílem návrhu je synchronní motor s permanentními magnety. Na začátku projektu jsme měli základní požadavky na navrhovaný motor. Shrnuty jsou v Tabulce 1. V RMxprt jsme na základě těchto parametrů provedli prvotní návrh motoru (rozměry motoru, založení vinutí, materiálové vlastnosti atd.). Vzhledem k časově nenáročné metodě analytických vztahů programu RMxprt lze provádět množství rychlých výpočtů, měnit parametry a pozorovat vliv těchto změn na sledované vlastnosti motoru (účinnost, hmotnost, jmenovitý moment atd.). Tímto postupem jsme navrhli motor, který splnil požadavky uvedené v Tabulce 1. Na Obrázku 1 http://aum.svsfem.cz



vidíme výstup z RMxprt, v Tabulce 2 jsou uvedeny základní hodnoty parametrů navrhnutého motoru.

Obrázek 1 - Výstup z RMxprt Tabulka 2 Výsledné parametry motoru z RMxprt

RMxprt automaticky vytvořil geometrický model navrženého motoru v ANSYS Maxwellu 2D, Obrázek 2, kde byl následně proveden přesnější MKP výpočet. Pro výpočet v ANSYS Maxwellu jsme definovali požadované okrajové podmínky a nastavili vlastnosti mechanického pohybu rotoru s permanentními magnety (počáteční rychlost, počáteční poloha, zátěž, tlumení). http://aum.svsfem.cz



Nejprve jsme simulovali situaci v ustáleném režimu, kdy byl motor zatěžován momentem 5,5 Nm. Tento moment byl ekvivalentní zátěži při námi definovaném jízdním cyklu v běžném městském provozu. Výsledkem výpočtu bylo rozložení elektromagnetických polí, Obrázek 3, a časový průběh ztrát ve vinutí, Obrázek 4. Jejich průměrná hodnota v ustáleném stavu je 230 W.

Obrázek 2 - 2D geometrický model navrženého elektromotoru

Obrázek 3 - Rozložení magnetické indukce v dané pozici rotoru




Obrázek 4 - Časový průběh ztrát ve vinutí v ustáleném stavu

Dalším krokem je export ohmických ztrát výše popsaným způsobem (v prostředí ANSYS Workbench) na 3D geometrii navrhnutého motoru v ANSYS Mechanical. Zde jsou zadány tepelné materiálové vlastnosti jednotlivých částí motoru, počáteční teplota a teplota okolí (40 °C), koeficienty přestupu tepla do okolí. Importované ztráty z ANSYS Maxwellu byly užity k deklaraci objemových zdrojů tepla ve vinutí. Dalším zdrojem tepla jsou vířivé proudy indukující se v permanentních magnetech rotoru, tyto zdroje byly započítány také. Na Obrázku 5 vidíme rozložení teplot v ustáleném stavu na celé geometrii motoru, Obrázek 6 pak ukazuje detail rozložení teplot v tělese statoru. Motor jsme navrhovali s tepelnou třídou izolace F (maximální přípustná teplota 145 °C). Tepelná analýza motoru pro městský jízdní režim ukázala, že navržený motor této třídě izolace vyhoví.

Obrázek 5 - Rozložení teplot na motoru v ustáleném stavu


47



Obrázek 6 - Rozložení teplot na statoru v ustáleném stavu

1

Závěr

Názorný příklad návrhu vhodného elektromotoru pro elektromobil ukazuje, že software ANSYS disponuje vhodnými nástroji pro celý proces vývoje takového zařízení. Od počátečních návrhů, přes přesné elektromagnetické analýzy pro různé jízdní režimy, až po teplotní výpočty výsledného návrhu. Nedílnou součástí jsou také optimalizační nástroje, které dále zefektivňují vyvíjený produkt. Výše uvedený příklad rozebíral možnosti tepelných výpočtů pro ustálený stav. Tepelné analýzy je možné provádět také pro nestacionární jízdní režimy. Další možnosti multifyzikálních analýz nabízí prostředí Simplorer, které umožňuje simulaci komplexních systémů pomocí zjednodušených blokových schémat, elektrických obvodů, vývojových diagramů, rovnic, MKP a CFD modelů apod. Do systému můžeme zahrnout komponenty elektrické, tepelné, elektromechanické, elektromagnetické a hydraulické. Kontaktní adresa: Na Pankráci 26, Praha 4




A steady state approach to calculation of valve pressure rise rate characteristics JÁN ORAVEC Technical centre, Sauer-Danfoss a.s., Považská Bystrica Abstract: Pressure relief valves are widely used in hydrostatic circuits. The pressure rise rate characteristics are one of the most important descriptions of the valves performance. This article proposes design of the basic pressure rise rate characteristic of valve using only a few CFD calculations. These calculations are performed for several poppet strokes and flow values. The method is shown for charge pressure relief valve but can be used similarly for any poppet type valve. Keywords: pressure relief valve, hydrostatics, pressure rise rate characteristic, CFD 1

Introduction

Pressure relief valves are essential components in hydrostatic systems. They protect the systems against excessive pressures or maintain pressures at desired levels. The performances of such valves are represented by their pressure rise rate characteristics. A simplified one is shown in Fig. 1. Note, the valve remains closed while pressure is lower than the cracking pressure. If the pressure exceeds the cracking pressure, pcr, the valve opens to allow escaping of working fluid (hydraulic oil) to a reservoir. A small slope of the characteristic is desired in standard working range of the valve. If the poppet reaches its maximal stroke limit, the valve behaves like a constant orifice for higher flow rates. Because each the point of the characteristic represents different poppet stroke (also flow, pressure…), the standard calculation (as well as real measurement) leads to transient task including remeshing and solution of equation of motion.




Figure - 1 Simplified pressure rise rate characteristic

If we take into account some simplifications: dynamic effects are avoided, valve at fixed poppet stroke behaves like an orifice which has quadratic function behavior for wide range of flow rates, we can design the valve characteristic from several CFD steady state calculations. This approach is shown for charge pressure relief valve.




Figure - 2 Design of charge pressure relief valve

2

Charge pressure relief valve (CPRV)

The function of CPRV is to maintain charge pressure of the closed hydrostatic system at a designated level above case pressure (Sauer-Danfoss, 2011). The CPRV is a direct acting poppet valve (Fig. 2). It opens and discharges fluid to the pump case when pressure exceeds the designated level, pcharge, at a channel called gallery. The valve poppet has a cone which acts against a seat created in the pump end cap (Fig. 2). The poppet guiding is provided by a guide designed in the nut. The 2 springs provide force which acts in direction to close the valve. Pretension of one spring can be adjusted by a bolt to reach desired level of the cracking pressure. The pretension force can be calculated from equation: (1)




(2)

where k is spring rate and s is poppet stroke. 3

CFD model

As mentioned above, the steady state approach is used here. The CFD models (Fig. 3) are prepared for 3 poppet strokes (0.5 mm, 1.5 mm and 3 mm). There is used combination of hexahedral, tetrahedral and wedge mesh types created using ANSYS ICEM CFD. The CFD analysis is performed in ANSYS Fluent. Material properties of hydraulic oil Shell Tellus 46 are used. The fluid is assumed as incompressible. Based on previous investigation, experience and comparison with testing, laminar model is chosen in this specific case. Pressure outlet boundary condition is applied to the pump housing passage. The pressure

is constant for all the calculated cases.

Velocity inlet boundary condition is experienced as suitable and stable for this task. Inlet velocity and flow rate are related to each other by simple relation: (3)

where is inlet flow rate, is inlet velocity magnitude and is inlet area. The velocity magnitude (flow rate respectively) is ramped to several values for each the stroke model. At least, 3 discrete velocity magnitudes are needed for each the stroke to be able approximate results by quadratic function. As a result from each the CFD calculation is: pressure at inlet

,

flow force induced by oil to the poppet (F). For reference, pressure contours and streamlines in region of control passage are shown in Fig. 4 for 1 mm poppet stroke and 70 lmin-1 flow rate case.




Figure - 3 CFD model

-1

Figure - 4 Pressure contours and streamlines in control passage region (1 mm stroke, 70 lmin flow)



4


Post-processing

The post-processing plays big role in this task. We have data from 9 CFD calculations (3 strokes × 3 flow rates per each stroke = 9 calculations) in this specific example. Therefore further data processing in a spreadsheet or other suitable tool is needed. The input and output values can be collected into a table form for further processing. One table proposal is following:

Table 1Proposed table for collection of CFD inputs and outputs

The visualization of data points can be seen in Fig. 5. The triangle points in graph C represent dependency of the pressure drop on flow rate. The different colors relate to different poppet strokes. Similarly, the square points in graph B represent dependency of the flow force on flow rate. The data points can be approximated by quadratic polynomial trendlines (continuous colored lines in graphs B and C) for each the stroke value. The flow force at constant stroke can be then written in form: (4)

where a,b, and c are polynomial coefficients obtained from the spreadsheet trendline properties.It is obvious the flow force has to be in balance with the spring force (2) at each the point of final pressure rise rate characteristic. Therefore we can write: (5)

Solution of this quadratic equation is:

(6)




where the sign plus or minus is chosen in way to obtain physically correct result. So, using (6) we get flow rate at which there is balance between flow force and spring force at constant stroke (compare graphs A and B in Fig. 5 to see balance between spring force and flow force for each the stroke) .

Figure - 5 Visualisation of the pressure rise rate characteristic design

Similarly to the flow forces approximation, we can also approximate the pressure drop data points using quadratic polynomial trendlines in form: (7)

where the polynomial coefficients d, e, and f are obtained again from the spreadsheet trendline properties. If we use the calculated flow rate from (6), we obtain corresponding pressure drop: (8)

This procedure is repeated for all the 3 strokes so we obtain the flow rate and pressure drop pairs: (9)

plus we have the additional point The final pressure rise rate characteristic can be designed from those points – see red line in graph C in Fig. 5.



5


Comparison with measurement

A simplified hydraulic scheme used for measurement of the CPRV is shown in Fig. 6a. An auxiliary pump H provides oil to the testing device. The throttle valve TV is fully open so the pressure indicated by the gauge G1 is small in comparison with the cracking pressure of the tested CPRV. As the TV restricts flow area, the load of CPRV increases and it starts to operate. Pressure before CPRV (G1) and behind it (G2) is recorded together with flow rate (Q). Note the measurement process should be sufficiently slow to avoid dynamic effects of the system. It is possible to perform measurement of the CPRV installed in a simple test block or installed in the pump end cap. The measured valve performance with the one calculated by CFD are shown in Fig. 6b for comparison.

Figure - 6 Measurement hydraulic scheme (a) and comparison CFD with measurement (b)

Conclusion The design of pressure rise rate characteristic of charge pressure relief valve using steady state CFD simulation is shown in the article. The model is prepared for 3 constant strokes. Minimally 3 different flow rate conditions needs to be calculated for each the stroke due to possibility of data approximation by quadratic polynomial functions. Both pressure drop and flow force functions are approximated in relation to flow rate. The final pressure rise rate characteristic of the valve is found using condition of equilibrium of flow force and spring force at each the specific stroke. Because of computing only a few steady state cases, the proposed method saves computing time as well as effort to set up the analysis in comparison to transient CFD computation. Another advantage is the effect of the cracking pressure or spring rate change is recalculated immediately in the spreadsheet without any additional CFD calculation requirement. References “H1 Axial Piston Pumps Single and Tandem”, Basic Information, Sauer-Danfoss, 11062168-Rev BB-Apr 2011. IMRE M., KRIŠŠÁK P., RAHMFELD R., ZAVADINKA P., “Porovnanie simulačných a meraných výsledkov škrtiaceho ventila”, Hydraulika a pneumatika, 1-2/2010, ISSN 1335-5171.



Contact [email protected]





ŘEŠENÍ VNĚJŠÍ AERODYNAMIKY KOLEJOVÝCH VOZIDEL MILAN SCHUSTER Výzkumný a zkušební ústav Plzeň s.r.o. Abstract: This paper deals with numerical simulations of external aerodynamics of rail vehicles, it is focused on aerodynamic interaction between the trains and vicinity infrastructures – mainly train entering to the tunnel. The method to study aerodynamic effects connecting measurements on the real vehicle and simulation on its virtual model is briefly introduced. The main possibilities to solve simulations for this method are discussed. This method can find wider application on research or industry. Keywords: train, aerodynamic load, simulation, measurement, validation. 1

Úvod

Příspěvek obsahuje stručný souhrn možností numerických simulací vnější aerodynamiky kolejových vozidel, sledování aerodynamické interakce jedoucích vozidel s nejbližším okolím tratě a možnosti pro porovnání výsledků simulací s výsledky měření. Numerické simulace umožňují řešit dynamické účinky a aerodynamické vlivy na jedoucí soupravu vlaku za různých podmínek (provozní režimy, rychlost jízdy, uspořádání tratě a okolí, působení bočního větru, úplav a vznos prachu ap.) a dále řešit aerodynamické účinky jedoucí soupravy vlaku na okolí (např. průjezd okolo nástupiště, průjezd okolo svislých a vodorovných konstrukcí nad kolejemi: např. návěstidla, konstrukce trakčního vedení, lávky, protihlukové stěny ap.). Podmínky pro zadání a nastavení numerických simulací a rozsah a způsob vyhodnocení výsledků je velmi často dán doporučujícím nebo závazným předpisem, např. podle TSI pro vysokorychlostní vozidla (TSI = technická specifikace pro interoperabilitu), podle evropských norem EN 14067 a v některých případech podle provozních (nebo i stavebních) předpisů. 2

Výzkum vnější aerodynamiky kolejových vozidel

Jedoucí kolejové vozidlo je za provozu jednak namáháno doplňkovým aerodynamickým odporem a jednak svým průjezdem také ovlivňuje dodatečným silovým zatížením okolní infrastrukturu (tj. stavby a objekty v blízkosti tratě). Znalost vzájemného silového působení mezi jedoucím kolejovým vozidlem a jeho okolím je potřeba pro návrh, konstrukci a zkoušky nového vozidla pro zajištění bezpečného provozu. Experimentální výzkum vnější aerodynamiky kolejových vozidel sebou přináší některé komplikace (měřítko modelu vlaku pro laboratorní experimenty, pronájem tratě pro zkoušky skutečných vozidel, dosažení potřebné rychlosti jízdy, podmínky pro měření a umístění čidel ap.). Alternativou jsou simulace. Metody řešení numerických simulací vnější aerodynamiky kolejových vozidel lze podobně jako měřicí postupy rozdělit do dvou skupin: simulace aerodynamického tunelu, simulace jízdy soupravy vlaku. Úlohy první skupiny jsou simulace obtékání stojícího modelu vlaku a nejbližšího okolí, relativní rychlost proudění vzduchu mezi povrchem vlaku a okolím je dána okrajovou podmínkou. Při řešení simulací v druhé




skupině je potřeba zadat pohyb modelu soupravy, tím řešení s využitím podmínek „pohyblivých sítí“ vede ke složitějším a náročnějším nestacionárním výpočtům. 3

Vybrané řešené situace

Pomocí CFD simulací byly řešeny základní provozní situace spojené s vnější aerodynamikou a nebo aerodynamickou interakcí kolejových vozidel: aerodynamické zatížení povrchu vozidel a nebo zařízení na povrchu vozidel, míjení vlaků, průjezd vlaku tunelem, průjezd vlaku okolo stojící soupravy, průjezd vlaku podél nástupiště, míjení objektů kolejiště v blízkosti jedoucího vlaku, rozptyl spalin z motorového vozu za jízdy atd. V příspěvku jsou uvedeny jen vybrané ilustrativní úlohy simulací, které dokumentují možnosti a výsledky metody simulace skupiny „virtuální aerodynamický tunel“.

Obrázek 1 - Obtékání modelu celého kolejového vozidla



a.


Vnější aerodynamika soupravy vlaku

Simulace vnější aerodynamiky celé soupravy vlaku umožňují určit zejména aerodynamický odpor a zatížení povrchu vozidla, propočítat varianty vnějšího tvaru vozidla a zjistit základní informace o působení na okolí. Na Obr. 1 je výsledek simulací obtékání soupravy elektrické jednotky, jsou zobrazeny proudnice ve vertikální rovině symetrie a v horizontální rovině (v polovině výšky čelního skla). Volba velikosti výpočtové oblasti v okolí vozidla, velikost a zahuštění buněk výpočetní sítě a výběr podstatných detailů povrchu rozhodují o upočítatelnosti simulací tohoto typu. b.

Obtékání zařízení na povrchu vozidel

Cílem numerických simulací obtékání detailů povrchu kolejových vozidel je pro oba směry jízdy soupravy vlaku a pro různé rychlosti jízdy zjistit podmínky obtékání soupravy vlaku a zejména se zaměřit na obtékání a aerodynamické zatížení agregátů umístěných na povrchu soupravy, většinou umístěných na střeše vozidel.

Obrázek 2 - Příklad výsledku simulací tlakového zatížení prvků a zařízení na povrchu kolejového vozidla (části střešní nástavby pro ventilační a chladicí systém vozidla)

Simulace umožňují virtuálním způsobem najít vhodný tvar přídavných prvků a pracovních agregátů umístěných na střeše tak, aby bylo dosaženo zlepšených parametrů funkce těchto zařízení za vybraných provozních situací. Simulace také dovolují zjistit dodatečné aerodynamické zatížení těchto částí na povrchu vozidla. Testovací simulace také dovolují řešit vhodné varianty v uspořádání a v tvarech různých střešních spoilerů a usměrňujících prvků na povrchu vozidel. Účelem je optimalizace proudění v okolí vstupů do klimatizačních a chladicích systémů kolejových vozidel.


60



Obrázek 3 - Příklad výsledku obtékání tvarově složitého povrchu střechy kolejového vozidla

c. Aerodynamika podvozku hnacího vozidla Simulace vnějšího obtékání podvozku kolejového vozidla byly využity ke zjištění tlakových poměrů na vstupu do chladicího systému trakčních motorů. Trakční motory jsou umístěny v podvozku s poměrně hustou zástavbou různými provozními agregáty a vlastní konstrukcí rámu podvozku, uchycení hnacích dvojkolí, odpružení atd. Simulace byly řešeny opět za různých provozních situací: směr a rychlost jízdy, výkon ventilátoru chladicího systému aj. Na Obr. 4 je pohled na výpočtovou oblast a na Obr. 5 jsou ilustrativní výsledky obtékání podvozku. Simulace umožnily také alespoň základní studium možností popisu pohybu drobných částic z okolí jedoucího vozidla. Simulace byly zaměřeny na získání trajektorií fiktivních částic z povrchu kolejiště, které vzniknou zvířením a vznosem při průjezdu vozidla. Výsledkem simulací je popis situace ulpívání nečistot a drobných (poletujících) částic na některých místech kolejových vozidel při jízdě za nejrůznějších provozních stavů, jejichž působení může mít vliv na funkčnost některých zařízení kolejového vozidla.


61



Obrázek 4 - Výpočtová oblast modelu trakčního podvozku soupravy vlaku

Obrázek 5 - Výsledky simulací obtékání modelu podvozku

4

Metoda ověření simulací

Numerické simulace proudění a sdílení tepla se s úspěchem používají v řadě technických úloh spojených s průmyslovými projekty návrhu nových zařízení. Simulace umožňují řešit více variant tvaru zařízení a / nebo více variant provozních parametrů. Výsledky simulací jsou ověřovány měřením, proměřovány by měly být alespoň ty varianty




zařízení, které byly úspěšně simulovány. Porovnávání (validace) výsledků simulací a výsledků měření je vděčné téma mnohého výzkumu s podněty pro obě „strany“. Technická praxe však zná situace, ve kterých měření není možné z důvodů extrémních rozměrů nebo extrémních hodnot parametrů proudících tekutin. V těchto situacích jsou simulace zřejmě jediným zdrojem informací k popisu chování zkoumaných zařízení. Mějme konstrukční nebo provozní situaci zařízení, kdy rozměry nebo parametry umožňují provést měření (stav A). Pro návrh zařízení je však třeba znát informace i o situaci s extrémními parametry (stav B), kdy měření není možné. Za těchto podmínek se nabízí kombinovaná metoda, která spojuje výsledky měření a simulací. Kombinovaná metoda v oblasti aerodynamiky kolejových vozidel nalézá využití při posuzování podmínek nastávajících při průjezdu vysokorychlostní soupravy tunelem. Průjezd vlakové soupravy tunelem sebou přináší přídavné aerodynamické zatížení jednak vlastního kolejového vozidla a jednak okolní stavby tunelu. O velikosti a průběhu těchto silových zatížení rozhodují tvar a délka soupravy a rychlost jízdy soupravy, tvar a rozměry tunelu a ostatní vlivy provozu. Mezinárodní předpis všechny varianty reálných tunelů shrnuje do požadavku na etalonový tunel, tento fiktivní ideální tunel by měl sloužit k posuzování aerodynamických účinků vysokorychlostních kolejových vozidel, předpis určuje podmínky jízdních zkoušek a experimentů potřebných ke schvalování vozidel. V současné síti tratí většiny železničních správ skutečný tunel požadovaných parametrů etalonového tunelu není k dispozici. Kombinovaná metoda může napomoci ke zjištění aerodynamických účinků jedoucí vlakové soupravy v tunelu s tím, že stav A je skutečný tunel a stav B je fiktivní etalonový tunel. Princip metody v několika krocích: 1. Měření (stav A) - průjezd skutečné testovací vlakové soupravy skutečným tunelem 2. Numerická simulace (stav A) - průjezd modelu soupravy modelem skutečného tunelu 3. Validační proces úpravy výpočtových modelů simulací podle výsledků měření 4. Simulace (stav B) - průjezd modelu soupravy modelem etalonového tunelu. Využití navržené kombinované metody vyžaduje splnění a provedení několika podmínek jak v oblasti měření, tak v oblasti simulací: rozdíly mezi testovanou soupravou vlaku pro měření a soupravou modelového vlaku pro simulace, podmínky při měření a jejich zahrnutí v simulacích, postup vyhodnocení výsledků a validace výpočtových modelů. Základní podmínkou je dodržení co nejvíce možného souladu mezi podmínkami měření a možnostmi zahrnutí skutečných podmínek experimentu do zadání numerických simulací. Je snaha minimalizovat dopady případných rozdílů do metodiky simulací. Jako nejvíce riskantní, mimo vlastní proces výpočtů simulací, se zatím jeví otázka vlivu geometrických a tvarových rozdílů detailů povrchu mezi skutečnou a modelovou testovací soupravou, vlivy podmínek výpočetní sítě objemu celé výpočtové oblasti pro simulace a konečně vlivy podmínek vzorkování a zpracování naměřených signálů při porovnání s výsledky získanými metodami pro zpracování simulací. 5

Závěr

V příspěvku jsou stručně zmíněny vybrané úlohy vnější aerodynamiky kolejových vozidel s ilustrativními výsledky numerických simulací nejdůležitějších provozních režimů. Stručně je naznačena možná metoda využití simulací v aplikacích s výskytem extrémních hodnot. Uvedená kombinovaná metoda spojuje postupy a výsledky získané měřením a




simulací skutečných a modelovaných stavů. Jsou zmíněny nejdůležitější podmínky připravované metody pro její využití při zjišťování aerodynamických účinků vlakové soupravy projíždějící tunelem. Postupy a poznatky kombinované metody aplikované na řešení aerodynamiky kolejového vozidla mají také obecnou platnost i pro další úlohy výzkumu proudění. Lze je využít pro ověřování výsledků a validaci výpočtových modelů simulací proudění a přenosu tepla i v dalších oblastech technické výzkumné praxe. Literatura SCHUSTER M., „Simulations of aerodynamic effects of rail vehicles“, 26. konference "Výpočtová mechanika 2010", Nečtiny, ZČU Plzeň, 2010. SCHUSTER, M., NOVÁK, J., KLESA, J., ZÁBRANSKÝ T., „Vnější aerodynamika a aerodynamická interakce kolejových vozidel“, Interní výzkumné zprávy, ŠKODA VÝZKUM, Plzeň 2005 – 2011. ČSN EN 14067-5, Český normalizační institut, 2007. Poděkování V příspěvku jsou uvedeny simulační výpočty s výsledky, které řešili moji kolegové Jan Klesa, Jakub Novák, Tomáš Zábranský v poslední době. Příspěvek vznikl v rámci řešení úkolů Výzkumného centra kolejových vozidel, které je podporováno MŠMT ČR (projekt č. 1M0519). Kontaktní adresa: Výzkumný a zkušební ústav Plzeň s.r.o. Tylova 46 301 00 Plzeň e-mail: [email protected]




SIMULACE PODCHLAZENÉHO VARU VE FLUENTU 13 ING. LADISLAV VYSKOČIL PH.D. Ústav jaderného výzkumu Řež a.s. Abstract: Simulation of subcooled boiling in Fluent 13 The paper presents possibilities for simulation of subcooled boiling under conditions in the core of nuclear pressurized water reactor. All the underlying physical phenomena that have to be modeled are summarized in the first part of the paper. Modeling approach in Fluent 13 CFD code is then described as well as the comparison with capabilities of the older Fluent 6.3 code. Finally, an example simulation of subcooled boiling in heated pipe is presented. The simulation results are compared with experimental data. Keywords: subcooled boiling, CFD, Fluent 13, DEBORA experiments 1

Úvod

Podchlazený var na otápěné stěně při nucené konvekci se využívá v mnoha průmyslových odvětvích. Výhodou podchlazeného varu oproti konvenčním metodám chlazení je velký součinitel přestupu tepla – díky latentnímu teplu a lokálnímu nestacionárnímu omývání stěny při odpoutávání bublin od stěny. Podchlazeného varu se využívá zejména v situacích, kdy je potřeba chladit stěny s vysokými tepelnými toky v řádovém rozmezí 104 – 107 W/m2. Cílem této práce bylo posoudit možnosti modelování podchlazeného varu za podmínek v aktivní zóně tlakovodního jaderného reaktoru. 2

Fyzikální jevy v podchlazeném varu

V nejvíce zatížených kanálech aktivní zóny tlakovodního reaktoru může dojít k podchlazenému varu. Teplota stěny palivového proutku překročí v určité vzdálenosti od vstupu do aktivní zóny teplotu sytosti a na stěně proutku se začnou vytvářet parní bublinky. Teplota kapaliny uprostřed kanálu je přitom v tomto místě stále pod teplotou sytosti. Parní bublinky se odpoutávají od stěny a jsou unášeny doprostřed proudu, kde postupně kondenzují – tj. zmenšují se. Pokud se kapalina uprostřed kanálu zahřeje na teplotu sytosti, dojde k objemovému varu a bublinky uprostřed kanálu se začnou zvětšovat. Další vliv na velikost bublin v proudění má turbulence, díky níž do sebe bubliny narážejí a mohou se spojovat. Na druhé straně turbulence může velké bubliny „trhat“ na menší – dochází k rozpadu bublin. Při modelování varu (a krize varu) v palivové kazetě tlakovodního reaktoru je potřeba zohlednit všechny tyto fyzikální jevy, tj. vznik bublin na stěně a přenos tepla mezi otápěnou stěnou a oběma fázemi; přenos hybnosti, tepla a hmoty mezi bublinami a kapalinou; a interakce bublin. 3

Modelování podchlazeného varu v CFD programu

V této kapitole jsou uvedeny matematické modely jevů, ke kterým dochází v podchlazeném varu. Cílem bylo vybrat co nejpřesnější modely pro dané fyzikální jevy, ale zase http://aum.svsfem.cz



na druhé straně výpočet příliš nekomplikovat, aby daný přístup byl reálně použitelný i pro výpočty ve složitějších geometriích. a.

Eulerovský model

V této práci je pro modelování dvoufázového proudění využit Eulerovský model, protože umožňuje nejpřesněji popsat fyzikální jevy v podchlazeném varu. Uvažují se dvě fáze: primární fáze je kapalina, sekundární fáze jsou parní bubliny. Obě fáze sdílí jedno tlakové pole. Rovnice kontinuity, rovnice zachování hybnosti a rovnice zachování energie se řeší pro každou fázi zvlášť. Proudění je nestlačitelné. Stejný přístup lze ale použít i se stlačitelnou sekundární fází (tj. párou). Fyzikální vlastnosti média jsou funkcí teploty. Turbulence se modeluje pomocí dvourovnicového modelu k-epsilon, upraveného pro dvoufázové proudění. b.

Modelování varu na otápěné stěně

Pokud teplota povrchu stěny překročí teplotu sytosti, mohou se na stěně začít vytvářet bubliny i v případě, že teplota okolní kapaliny je stále pod teplotou sytosti. Místo, kde se na stěně objeví první bublinky, se označuje jako počátek podchlazeného varu. Pro určení počátku podchlazeného varu lze využít například kritérium podle Hsu (1962). Podle tohoto kritéria bude růst bublina ze zárodku v trhlince na stěně, pokud je na vrcholu zárodku teplota kapaliny větší nebo rovná teplotě sytosti při daném tlaku v bublině. Zjednodušený přístup předpokládá, že podchlazený var začne, jakmile teplota stěny překročí teplotu sytosti. Pro stanovení množství páry vznikající na stěně při podchlazeném varu v kanále se používá model dělení tepelného toku. Existuje několik takových modelů. V této práci je použit model dělení tepelného toku podle Kurula a Podowskiho (1990). V uvedeném modelu se tepelný tok ze stěny dělí na tři části: (1) První část je jednofázový přestup tepla (konvektivní tepelný tok): q f  A1 wallfcn Twall  Tl  A1  1  A2

(2) (3)

A1 je podíl plochy povrchu stěny ovlivněný pouze kapalinou, A2 je podíl plochy ovlivněný bublinami vznikajícími na stěně, Tl je teplota kapaliny ve středu buňky přilehlé ke stěně, αwallfcn součinitel přestupu tepla ze stěnové funkce pro teplotu kapaliny. Druhá část tepelného toku qq (tzv. ”quenching”) je teplo předané při nestacionárním omývání stěny v intervalu mezi odpoutáním bubliny od stěny a vznikem další bubliny v dané pozici. Jedná se o velmi rychlý přechodový děj. qq  A2 quench Twall  Tl 


(4)



αquench je součinitel přestupu tepla pro nestacionární omývání stěny (11). Třetí část tepelného toku, qe je použita na vypařování kapaliny: (5) je rychlost vypařování vztažená na jednotku plochy stěny (9), Hlat výparné teplo. Model předpokládá, že průměr oblasti ovlivněné jednou bublinou je dvojnásobkem průměru bubliny odpoutávající se od stěny dw: (6) n je hustota aktivních nukleačních míst. Hustota aktivních nukleačních míst se počítá z přehřátí povrchu stěny – vychází z dat Lemmert a Chawla (1977): (7) Průměr bubliny odpoutávající se od stěny dw se počítá z korelací. Existuje např. poměrně složitá Ünalova korelace (1976) nebo jednoduchá korelace podle Tolubinskiho a Konstančuka (1970). Pro výpočet rychlosti vypařování je nutno spočítat frekvenci odpoutávání bublin od stěny f pomocí následujícího vztahu - Ceumern-Lindenstjerna (1977): (8) Rychlost vypařování je součin hmotnosti bubliny, frekvence odpoutávání bublin a hustoty aktivních nukleačních míst: (9) Druhou část tepelného toku (quenching) analyticky odvodili Del Valle a Kenning (1985). Součinitel přestupu tepla pro nestacionární omývání stěny αquench závisí na délce časového intervalu mezi odpoutáním bubliny od stěny a vznikem nové bubliny. Předpokládá se, že tato „čekací“ doba tw je rovná periodě odpoutávání bublin: (10)

(11) kde al je teplotní vodivost kapaliny (difuzivita) a λl tepelná vodivost kapaliny.



c.


Přenos hybnosti mezi fázemi

Mezifázový přenos hybnosti lze rozdělit na několik částí: odpor, efekt virtuální hmoty, hydrodynamickou vztlakovou sílu (lift) a turbulentní disperzi. Efekt „wall lubrication“ (síla působící na bubliny těsně u stěny), Basetova síla a Marangoniho síla (teplotní gradient podél rozhraní) se v této práci zanedbávají. 3.3.1. Odporová síla Tato síla funguje stejně jako aerodynamický odpor působící na jedoucí auto. Pokud se bubliny pohybují rychleji než kapalina, potom odporová síla brzdí bubliny a urychluje kapalinu. Mezifázová odporová síla se počítá z následující rovnice: (12) db je Sauterův střední průměr bubliny, cD je součinitel odporu. V této práci je použita jedna z variant Ishiiho korelace (1979), která byla vytvořena pro parní bubliny ve vodě. 3.3.2. Hydrodynamická vztlaková síla Tato vztlaková síla nesouvisí s rozdílem hustot fází, ale s rychlostním polem okolo bubliny. Je to stejná síla, jako vztlak křídla letadla. Pokud se bublina pohybuje v neuniformním rychlostním poli, působí na ni hydrodynamická vztlaková síla. Hydrodynamická vztlaková síla se počítá z rovnice: (13) V této práci se pro výpočet součinitele vztlakové síly používá korelace Moraga (1999), která zohledňuje účinky dvou následujících sil působících proti sobě: 1. klasická hydrodynamická vztlaková síla – interakce bubliny a rychlostního pole okolo bubliny – kladný vliv na cL. 2. interakce bubliny s víry v úplavu za bublinou – záporný vliv na cL Pozn.: vztlak daný rozdílem hustot je zohledněn v rovnicích pro zachování hybnosti pomocí členu s gravitačním zrychlením. 3.3.3.

Efekt virtuální hmoty

Pokud sekundární fáze zrychluje relativně vůči primární fázi, setrvačnost primární fáze vyvozuje „virtuální hmotovou sílu“ na sekundární fázi, čímž brání sekundární fázi ve zrychlení. Lze to popsat takto: pokud chceme urychlit bublinu, musíme současně s bublinou urychlit i určitý objem kapaliny okolo bubliny. V literatuře se tento jev někdy také popisuje jako „přidaná hmota“. Rovnice pro virtuální hmotovou sílu vypadá následovně (Drew a Lahey, 1979):




(14) Součinitel virtuální hmotové síly pro zrychlující kulovou bublinu v klidné tekutině je cVM =0.5. Tato hodnota byla odvozena z řešení potenciálního proudění (Batchelor, 1967), stejná hodnota se používá i ve Fluentu. Existuje ovšem řada jiných modelů a korekcí pro tento součinitel v závislosti na objemovém podílu páry, tvaru bubliny nebo vzájemné polohy bublin a stěny. 3.3.4.

Turbulentní disperze

V experimentech s pohybem bublin v turbulentním proudění kapaliny dochází k jevu zvanému turbulentní disperze. Náhodné působení vírů v kapalině má za následek rozptylování bublin z oblasti s vysokým objemovým podílem páry do oblastí s malým objemovým podílem páry. Analogií tohoto jevu je tepelná difúze molekul vzduchu v atmosféře Brownův pohyb. Roli tepelné energie molekul vzduchu zde přebírá turbulentní kinetická energie kapaliny. Turbulentní disperze se počítá z gradientu objemového podílu páry a turbulentní kinetické energie kapaliny kl (Lance, Lopez de Bertodano, 1994):

(15) Součinitel turbulentní disperze cTD se nastaví na konstantní hodnotu, která se odvíjí od použitého modelu turbulence. Ve Fluentu 13 se používá hodnota cTD = 1. Existují ale i sofistikovanější přístupy, kdy se součinitel turbulentní disperze počítá z odporové síly a virtuální hmoty. d. Přenos tepla a hmoty mezi fázemi Model pro přenos tepla mezi fázemi, který je používán v této práci, vychází z korelace Ranz- Marshall (1952). Jedná se o standardní a často používaný model pro přenos tepla na kulových částicích. Korelace je upravena tak, aby teplota mezifázového rozhraní držela teplotu sytosti. V tomto modelu se přenos tepla mezi fázemi dělí na dva problémy: 1. Přenos tepla mezi rozhraním a kapalinou:

(16) Součinitel přestupu tepla hli je:

(17)




db je Sauterův střední průměr bubliny. Objemová hustota mezifázového rozhraní ai je dána vztahem:

(18) Nusseltovo číslo Nu se počítá z korelací, ty jsou různé pro případ vypařování a pro případ kondenzace. 2. Přenos tepla mezi párou a rozhraním: Přenos tepla mezi párou a rozhraním se počítá pomocí metody “návrat do stavu sytosti během časového kroku”. Předpokládá se, že pára v bublinách si udržuje teplotu sytosti pomocí rychlé kondenzace nebo vypařování.

(19) δt je časový krok, cP,v je izobarická tepelná kapacita páry. V literatuře lze nalézt i jiné mnohem složitější modely přenosu tepla mezi fázemi. Přenos hmoty mezi fázemi přímo závisí na přenosu tepla mezi fázemi v objemu a na stěně. Konstantou úměrnosti mezi těmito přenosy je výparné teplo Hlat.

(20) je dáno rovnicí (9), Aface je plocha stěny v buňce přilehlé ke stěně, Vcell je objem buňky přilehlé ke stěně.

e. Modelování velikosti bublin Modely pro přenos hybnosti a tepla mezi fázemi závisí na průměru bublin. Tento průměr bublin je potřeba nějakým způsobem modelovat. V této práci je použita tzv. jednoskupinová metoda, která modeluje přenos mezifázového rozhraní. Využívá se jedna transportní rovnice se zdroji pro nukleaci bublin na otápěné stěně a s modely pro spojování a rozpad bublin. Výsledkem řešení této rovnice je rozložení středního průměru bubliny ve výpočtové oblasti. V každé výpočtové buňce je pak velikost bublin reprezentována jedním středním průměrem, který je ovšem v každé buňce jiný. Pro modelování spojování a rozpadu bublin jsou použity modely Yao a Morel (2004). Existuje samozřejmě řada jiných modelů pro spojování a rozpad bublin. Přenosová rovnice pro hustotu mezifázového rozhraní je z numerických transformována na přenosovou rovnici pro převrácenou hmotnost bubliny X: Převrácená hmotnost bubliny je:


důvodů



(21) db je Sauterův střední průměr bubliny. Transportní rovnice pro skalár X je potom:

(22) Zdrojový člen hustoty počtu bublin kvůli nukleaci na stěně:

(23) „f“ se počítá z rovnice (8) a „n“ z rovnice (7). Aface je plocha povrchu stěny v buňce u stěny,Vcell je objem této buňky. Zánik mezifázového rozhraní kvůli spojování bublin se modeluje takto (Yao a Morel, 2004):

(24)

(25) Kritická hodnota Weberova čísla je Wecr = 1.24. Konstanty jsou Kc1 = 2.86, Kc2= 1.922, Kc3 = 1.017. Vznik mezifázového rozhraní kvůli rozpadu bublin se modeluje takto (Yao a Morel, 2004):

(26) We je definováno rovnicí (25), Wecr = 1.24, konstanty jsou Kb1 = 1.6, Kb2 = 0.42.




Následující rovnice zohledňuje výměnu energie mezi volnou energií rozhraní a turbulentní kinetickou energií kapaliny při spojování a rozpadu bublin:

(27) Sk,l je zdroj v rovnici pro turbulentní kinetickou energii kapaliny, jsou zánik (24) a zdroj (26) mezifázového rozhraní kvůli spojování a rozpadu bublin, σ je povrchové napětí. 4. Možnosti modelování podchlazeného varu ve Fluentu 13 V této kapitole jsou probrány možnosti, jak lze výše uvedené matematické modely zakomponovat do Fluentu 13. 4.1 Var na otápěné stěně Ve Fluentu 13 je zabudován model varu na stěně vycházející z prací Kurula a Podowskiho (1990). Model lze zapnout v uživatelském rozhraní. V tomto modelu se průměr bubliny odpoutávající se od stěny počítá z jednoduché a diskutabilní korelace Tolubinskiho a Konstančuka (1970). V grafickém rozhraní lze jako výsledek zobrazit dílčí tepelné toky qf, qq a qe z rovnice (1). Tyto toky je možné i číst v uživatelem definované funkci (UDF). V grafickém rozhraní nelze zobrazit průměr bubliny odpoutávající se od stěny, lze jej ale dopočítat v UDF ze známých parametrů v rovnicích (5), (7), (8) a (9). V modelu nelze použít jinou korelaci než Tolubinski+Konstančuk. Jako kritérium pro počátek podchlazeného varu na stěně se používá jednoduchá podmínka Tw>Tsat. Pokud chceme použít jinou korelaci pro průměr bubliny odpoutávající se od stěny (např. Ünalovu), nebo jiné kritérium pro začátek podchlazeného varu (např. podle Hsu), potom musíme napsat celý model podchlazeného varu znovu pomocí UDF. Tento přístup je možný, ale velmi pracný. 4.2 Přenos hybnosti mezi fázemi Ve Fluentu 13 lze modelovat všechny 4 základní mechanismy výměny hybnosti (odpor, hydrodynamický vztlak, turbulentní disperzi a virtuální hmotu). Všechny tyto síly lze zapnout v uživatelském rozhraní. 4.2.1 Odporová síla Ve FLUENTu jsou k dispozici tři univerzální korelace pro odporovou sílu: SchillerNaumann, Morsi-Alexander a symetrický model (viz Fluent Theory Guide) a model napsaný přímo pro parní bubliny ve vodě: korelace podle Ishiiho (1979). Poslední model se zdá být nejvhodnější pro modelování proudění s varem. Ve Fluentu 13 nelze číst součinitel odporu cD v uživatelském rozhraní ani v UDF. Je samozřejmě možné jej zpětně dopočítat v UDF. 4.2.2 Hydrodynamická vztlaková síla Ve Fluentu 13 lze v uživatelském rozhraní vybrat korelaci Moraga (1999) nebo korelaci Tomiyama (2002). Obě korelace zohledňují protichůdné účinky (1) interakce bubliny a rychlostního pole okolo bubliny a (2) interakce bubliny s víry v úplavu za bublinou. Je také možné použít konstantní součinitel vztlaku cL. 4.2.3 Efekt virtuální hmoty




Ve Fluentu 13 je pevně nastavený součinitel virtuální hmotové síly cVM = 0.5. Tuto hodnotu nelze měnit. Efekt virtuální hmoty lze pouze vypnout. Bylo by zřejmě možné napsat vlastní virtuální hmotovou sílu pomocí UDF. 4.2.4 Turbulentní disperze Ve Fluentu 13 lze zapnout výpočet turbulentní disperze a to nejjednodušší model Lance, Lopez de Bertodano (1994). Součinitel turbulentní disperze je nastaven na hodnotu cTD = 1. Součinitel cTD lze změnit na jinou konstantní hodnotu z příkazové řádky příkazem (rpsetvar 'mp/boiling/cd-drift' Value), kde „Value“ je nová hodnota cTD. Pokud chceme použít zobecněný model, v němž je cTD funkcí odporové síly, musíme si napsat vlastní UDF. 4.3 Přenos tepla a hmoty mezi fázemi Ve Fluentu 13 lze v uživatelském rozhraní vybrat model přenosu tepla založený na korelaci Ranz-Marshall (1952), který je blíže popsán v oddíle 3.4, a také odpovídající přenos hmoty. Do modelů nelze nijak zasahovat, lze ale napsat vlastní modely pomocí UDF. 4.4 Modelování velikosti bublin Nejjednodušší přístup je použít konstantní střední průměr bubliny – stejný v celé výpočtové oblasti. Tuto hodnotu lze odhadnout buď z hodnot průměrů bublin odpoutávajících se od stěny a nebo z experimentálních měření. Tato metoda může poskytnout slušné výsledky v případě malých objemových podílů páry v oblasti (okolo 0.2 na otápěné stěně). Tento postup ale není vhodný pro případy s vyššími objemovými podíly páry (např. 0.4), kdy lze očekávat velké rozpětí hodnot průměrů bublin (kvůli spojování a rozpadu bublin). V takovém případě už je nutné modelovat velikost bublin pomocí transportní rovnice pro hustotu mezifázového rozhraní s modely pro spojování a rozpad bublin. Ve Fluentu 13 je možné zapnout výpočet přenosu hustoty mezifázového rozhraní spolu s modely pro spojování a rozpad bublin Hibiki+Ishii a Ishii+Kim (Fluent Theory Guide). Z uživatelského rozhraní ale nelze zadat „propojení“ mezi modelem varu na stěně a rovnicí mezifázového rozhraní (vznik bublin na stěně z modelu varu by se měl objevit jako zdroj v rovnici pro mezifázové rozhraní). Tento zdroj lze dopsat pomocí UDF. Autorovi této práce se ale nepodařilo získat tímto způsobem rozumné výsledky, proto musel použít přístup vlastní. Transportní rovnici (22) a modely Yao a Morel (2004) popsané v oddíle 3.5 autor zabudoval do Fluentu pomocí rovnice pro uživatelem definovaný skalár (UDS) se zdrojovými členy napsanými pomocí UDF. 4.5 Nastavení řešiče pro modelování varu Uvažuje se proudění s varem v trubce nebo ve svazku tyčí. Proudění lze modelovat jako nestlačitelné (obě fáze). Použije se tedy „pressure-based“ řešič. Je nutné zohlednit gravitační sílu, protože fáze mají velký rozdíl hustot. Pro modelování turbulence se osvědčily dva různé přístupy: 1. Model k-epsilon pro směs (mixture) V tomto případě Fluent sám přidá do rovnic pro k a epsilon zdroje modelující turbulenci indukovanou bublinami. 2. Složitější model k-epsilon řešený pro obě fáze zvlášť (per phase)




V tomto modelu musí uživatel zdroje popisující turbulenci indukovanou bublinami napsat sám pomocí UDF. V této práci je použit tento složitější přístup. Podle zkušeností autora není velký rozdíl ve výsledcích, pokud zvolíme standardní model k- epsilon nebo model k-epsilon realizable. Fyzikální vlastnosti vody a páry se modelují jako po částech lineární funkce teploty. Nastavení řešiče, které se autorovi osvědčilo, vypadá následovně: Sekvenční řešení polí rychlostí a tlaků: Multiphase coupled Gradienty: Least Square Cell-Based Diskretizační schéma pro konvektivní členy v řešených rovnicích: -

protiproudé schéma 2. řádu pro hybnost, energii, k a epsilon

-

protiproudé schéma 1. řádu pro objemový podíl páry a uživatelem definovaný skalár (pokud se řeší spojování a rozpad bublin v proudu)

Řešič s dvojitou přesností (3ddp) Pozn.: Je sice možné použít schémata 2. řádu i pro objemový podíl a UDS, ale jsou pak velké problémy se stabilitou výpočtu. Ze stejného důvodu je nutno snížit relaxační faktory u těchto dvou veličin na hodnoty zhruba 0.1-0.2. Úloha se řeší jako přechodový proces, s konstatními okrajovými podmínkami v čase, dokud se řešení neustálí. Je možné použít metodu „pseudo-transient“ (implicitní metoda pro řešení stacionárních úloh, v níž se ustalování počítá s různě dlouhými časovými kroky v jednotlivých řešených rovnicích). 5 Možnosti modelování podchlazeného varu ve Fluentu 6 Ve Fluentu 6 je také možné modelovat podchlazený var. Většina potřebných modelů ovšem v programu chybí a uživatel si je musí napsat sám. V programu jsou dostupné pouze následující modely: 

obecné korelace pro odporovou sílu, korelace Ishii pro bubliny chybí



hydrodynamická vztlaková síla (konstantní součinitel cL, korelaci pro cL nutno dopsat v UDF)



virtuální hmotová síla (stejný model jako ve Fluentu 13)

Dále ve Fluentu 6 chybí metoda „Multiphase coupled“ pro sekvenční řešení polí rychlostí a tlaků, je zde jen metoda „Phase Coupled SIMPLE“. Není možné použít metodu „pseudotransient“ pro ustalování. Autor této práce v minulosti napsal a odzkoušel model podchlazeného varu pro programy Fluent 6.1 a Fluent 6.3. Bližší informace lze nalézt např. v publikaci Vyskocil, Macek (2008). Pro zajímavost lze uvést, že uživatelem definovaná funkce (UDF) pro simulaci varu ve Fluentu 6.3 se všemi modely má zhruba 1300 řádek, využívá 20 různých maker DEFINE a 10 uživatelem definovaných pamětí (UDML). Obdobná UDF pro Fluent 13 má 370 řádek, 5 maker DEFINE a 7 UDML. Pokud nepotřebujeme znát parametry jako průměr bubliny odpoutávající se od stěny a součinitel




odporu a pokud nepotřebujeme modelovat přenos hustoty mezifázového rozhraní, lze podchlazený var ve Fluentu 13 simulovat bez jakékoliv UDF. 6. Simulace podchlazeného varu v trubce V této kapitole je ukázka simulace podchlazeného varu v experimentu DEBORA. Experimenty DEBORA byly provedeny v CEA Grenoble, viz ref. Falk et al.(2002), Garnier (1998, 2001). Cílem bylo vytvoření databáze experimentů proudění s varem až do blízkosti krize varu v oblasti parametrů, které jsou v tlakovodních reaktorech. Tato databáze má sloužit k lepšímu porozumění jevům, které probíhají při podchlazeném varu a v blízkosti krize varu1. druhu. Experimentální zařízeníTestovací zařízení je svislá elektricky ohřívaná trubka, kterou proudí freon R12. Trubka se skládá ze vstupní 1 m dlouhé neotápěné sekce, z 3.5 m dlouhé ohřívané části a z 0.5 m dlouhé výstupní neohřívané části. Vnitřní průměr trubky je 19.2 mm. Parametry média R12 jsou zvoleny tak, aby simulovaly proudění vody v aktivní zóně tlakovodního reaktoru (využívá se fyzikální podobnosti mezi prouděním vody a freonu). Výhodou tohoto přístupu je nízký tlak a teplota média. Na vstupu do ohřívané sekce má médium R12 teplotu nižší, než je teplota sytosti. Na ohřívané stěně trubky dochází k podchlazenému varu. Na konci otápěné sekce jsou měřeny po průměru trubky profily objemového podílu páry, hustoty mezifázového rozhraní, velikosti bublin, časově středované osové rychlosti páry a teploty kapaliny. 6.2 Vybrané experimenty V následující tabulce jsou uvedeny parametry experimentů, které byly vybrány pro simulaci v CFD programu Fluent 13. Tabulka 1 Parametry experimentů p G 2 Experiment MPa kg/m /s Case 1 3.008 1005 Case 2 3.007 1004.8 Case 3 3.007 1004.8 Case 4 3.006 994.9

qw 2 W/m 58260 58260 58260 58260

Tinlet °C 67.89 70.14 72.65 73.7

Tsat °C 94.24 94.22 94.22 94.21

Xeq 0.1218 0.1563 0.1953 0.2173

p je tlak, G je hmotový tok, qw je tepelný tok na stěně, Tinlet je teplota na vstupu, Tsat je teplota sytosti při tlaku p a Xeq je rovnovážná suchost na výstupu. Všimněte si, že v této posloupnosti experimentů se zvyšuje teplota na vstupu, tepelný tok na stěně je stejný a hmotový tok je přibližně stejný. Pomocí fyzikální podobnosti freonu R12 s vodou lze spočítat parametry ekvivalentních experimentů ve vodě (viz Coffield, 1969). Například pro experiment „Case 4“ by odpovídající experiment ve vodě měl parametry: tlak 17.3 MPa, průměr trubky 17.5 mm, hmotový tok 1624 kg/m2/s a tepelný tok 1.062 MW/m2.




6.3 Výpočtová síť Případ je osově symetrický. Trubka má vnitřní průměr 19.2 mm. Modeluje se 1 m dlouhá vstupní sekce, 3.5 m dlouhá otápěná sekce a 0.5 m dlouhá výstupní sekce. Síť je ve tvaru 10° klínu ohraničeného rovinami symetrie, viz následující obrázek:

Obrázek 1 - Výpočtová oblast a síť pro simulaci experimentů DEBORA – 14x1000 buněk

Byly vytvořeny 3 sítě s rozlišeními: 10x500, 14x1000 a 19x2000 buněk. Všechny výpočty byly provedeny na síti 14x1000 buněk. Jeden výpočet byl proveden na všech 3 sítích – pro posouzení nezávislosti řešení na síti, viz Obr. 4. 6.4 Okrajové podmínky Vstup: 

stejná rychlost obou fází



teplota kapaliny podle experimentálních dat, teplota páry = teplota sytosti při daném tlaku



objemový podíl páry → 0



intenzita turbulence: 3% (odhad), hydraulický průměr 19.2 mm



převrácená hmotnost bubliny X zadána tak, aby průměr bublin → 0, viz rovnici (21)

Roviny řezu: symetrie Výstup: konstantní statický tlak Stěny: zadaný tepelný tok, nulové rychlosti obou fází, nulová drsnost 6.5 Poznámky k experimentálním datům V experimentálních datech jsou použity dvě charakteristiky průměrné velikosti bublin – Sauterův střední průměr a průměr Dg.




Sauterův střední průměr (Sauter Mean Diameter, SMD) se definuje jako průměr koule, která má stejný poměr objemu a povrchu jako daná částice (tj. zde bublina). Střední průměr Dg je průměr bubliny, který by mělo ekvivalentní dvoufázové proudění se stejnou hustotou počtu bublin, stejnou hustotou mezifázového rozhraní a se všemi bublinami o stejném průměru, viz Garnier (1998). 6.6 Výsledky Na následujících obrázcích jsou znázorněny změřené a vypočtené radiální profily veličin na konci otápěné sekce. V této sérii experimentů nebyla měřena teplota kapaliny.




Obrázek 2 - DEBORA, Case 1 – Tin = 67.89°C, Xeq = 0.1218

Z obrázku je patrný rozdíl mezi změřeným a vypočteným profilem rychlosti páry. Tento problém je pravděpodobně způsoben použitím jednofázových stěnových funkcí ve výpočtu. V literatuře lze nalézt pokusy, jak tento problém odstranit – použitím stěnových funkcí s „umělou drsností“, která se odvíjí od velikosti bublin na stěně. Viz např. práci Končar (2010).




Obrázek 3 - DEBORA, Case 2 – Tin = 70.14°C, Xeq = 0.1563




Všimněte si změřených profilů rozložení středního průměru bubliny. Je z nich patrné, že uprostřed proudu dochází ke spojování bublin. Modely Yao a Morel (2004) dokáží tento efekt předpovědět.

Obrázek 4 - DEBORA, Case 3 – Tin = 72.65°C, Xeq = 0.1953, výpočet na 3 sítích




Z výsledků výpočtu „Case 3“ je vidět, že síť 14x1000 buněk už je dostatečně jemná a další zjemňování sítě nemá vliv na výsledky. Podrobné výsledky pro simulaci experimentu „Case 3“ jsou znázorněny v oddílu 6.6.1.

Obrázek 5 - DEBORA, Case 4 – Tin = 73.7°C, Xeq = 0.2173 (blízko krize varu)




Tento experiment je velmi blízko krize varu (DNB), je to vidět z profilu objemového podílu páry u stěny (prudký nárůst směrem ke stěně). Program velmi dobře zachytil objemový podíl páry v centru trubky, ale nedokázal zachytit „skok“ objemového podílu páry u stěny. Pozn.: Kritérium pro krizi varu podle Weismana (1983) je objemový podíl páry na stěně rovný 0.82.

6.6.1 Podrobné výsledky simulace experimentu „Case 3“




Obrázek 6 - DEBORA, Case 3: Vypočtený objemový podíl páry a střední průměr bubliny

Všimněte si, že v horní části trubky uprostřed proudu dochází ke spojování bublin.

Poznámky k Obr. 6 a Obr. 7:




Střední průměr bubliny je zobrazen pouze v místech, kde je pára. Výpočtová oblast je znázorněna v šikmém pohledu a ortogonálně (bez perspektivy). V bočním pohledu by nebylo nic vidět, protože poloměr trubky je 9.6 mm a délka 5 m (z toho otápěná délka je 3.5 m).

Obrázek 7 - DEBORA, Case 3: Vypočtená teplota a rychlost kapaliny




Obrázek 8 - DEBORA, Case 3: Dělení tepelných toků na stěně, viz rovnici (1)

Otápěná sekce začíná na svislé souřadnici y = 1m a končí na y = 4.5m. Z Obr. 8 je vidět, že na konci otápěné sekce se většina tepelného toku ze stěny využije na vypařování. 7

Závěr

Cílem této práce bylo posoudit schopnosti programu Fluent 13 modelovat podchlazený var v proudění. Nejprve byl proveden rozbor fyzikálních jevů, ke kterým dochází v podchlazeném varu. Dále byly popsány matematické modely uvedených jevů a možnosti jejich modelování ve Fluentu 13. Nakonec byla provedena simulace podchlazeného varu v otápěné trubce. Pro tuto simulaci byly vybrány 4 experimenty DEBORA provedené v CEA Grenoble. Podařilo se dosáhnout dobré shody výsledků s experimentálními daty. Ve Fluentu 13 jsou podstatně lepší možnosti pro modelování podchlazeného varu než ve starší verzi 6.3. Jednoduchou simulaci varu s konstantním průměrem bubliny lze ve Fluentu 13 provést i bez jakékoliv UDF. Některé detaily ale stále ve Fluentu 13 chybějí. Chceme-li například modelovat spojování a rozpad bublin nebo znázorňovat parametry jako průměr bublin odpoutávajících se od stěny nebo součinitel odporu, musíme použít UDF. Literatura ANSYS Fluent Theory Guide, Release 13.0, November 2010 ANSYS Inc, Southpointe, 275 Technology Drive, Canonsburg,PA 15317 BATCHELOR, G.K.: An Introduction to Fluid Dynamics. Cambridge University Press 1967, 1973, 2000 CEUMERN- LINDENSTJERNA, W.C.: Bubble Departure Diameter and Release Frequencies During Nucleate Pool Boiling of Water and Aqueous NaCl Solutions. Heat Transfer in Boiling, editors E. Hahne and U. Grigull, Academic Press and Hemisphere 1977 COFFIELD R.D., ROHRER W.M., TONG L.S.: A Subcooled DNB Investigation of Freon-113 and Its Similarity to Subcooled Water DNB Data, Nuclear Engineering and Design 11, pp 143153, 1969




DAVIDSON, M. R.: Numerical Calculations of Two-Phase Flow in a Liquid Bath with Bottom Gas Injection: The Central Plume. Appl. Math. Modelling, 14, 67-76, 1990. DEL VALLE, V.H., KENNING, D.B.R.: Subcooled Flow Boiling at High Heat Flux, International Journal of Heat and Mass Transfer, vol. 28, no.10, pp. 1907-1920, 1985 DREW, D. A., LAHEY, R. T. JR.: Application of General Constitutive Principles to the Derivation of Multidimensional Two-Phase Flow Equations, Int. J. Multiphase Flow, Vol.7, pp.243-264, 1979 FALK F., CUBIZOLLES G., HUGONNARD R.: Rapport d’essais DEBORA-FNU – Campagne d’essais USINU 3.6 m et 3.17 m SETEX/LTAC/02-155. CEA 2002 GARNIER J.: DEBORA topologie d’écoulement Mise en forme et description de la banque finalisée. SETEX/LTDF/98-65, CEA 1998 GARNIER J., MANON E., CUBIZOLLES G.:Local measurements on flow boiling of refrigerant 12 in a vertical tube Multiphase Science and Technology, Vol. 13, 2001 HSU, Y.Y.: On The Size Range Of Active Nucleation Cavities On A Heating Surface, J. Heat Transfer 84, pp.207-216, 1962 ISHII, M., ZUBER, N.:Drag coefficient and relative velocity in bubbly, droplet or particulate flows, AIChE Journal Vol.25, No.5, pp. 843-855 (1979) KONČAR B., BORUT M.:Wall function approach for boiling two-phase flows, Nuclear Engineering and Design 240, pp 3910–3918, 2010 KURUL, N., PODOWSKI, M.Z.:Multidimensional Effects in Forced Convection Subcooled Boiling,9th International Heat Transfer Conference, Jerusalem, Israel, 21-26 August 1990 LANCE M., LOPEZ DE BERTODANO M.:Phase distribution phenomena and wall effects in bubbly two-phase flows, Multiphase Science and Technology 8, 1994 LEMMERT, M., CHAWLA, J.M.:Influence of Flow Velocity on Surface Boiling Heat Transfer Coefficient, v Hahne E., Grigull: Heat Transfer in Boiling, Academic Press and Hemisphere, 1977 MORAGA, F.J., BONETTO, F.J., LAHEY, R.T.:Lateral forces on spheres in turbulent uniform shear flow, Int. J. Multiphase Flow 25, pp.1321-1372, 1999 RANZ, W. E., MARSHALL, W. R. JR.:Evaporation from Drops, Part I+II. Chem. Eng. Prog., 48(3): 141{146, March 1952, Chem. Eng. Prog., 48(4): 173{180, April 1952. Kontaktaktní adresa: Husinec-Řež 130, 250 68 Řež, [email protected]




VLIV ČÁSTEČNÉ SMÁČIVOSTI NA PROUDĚNÍ A KAVITAČNÍ OBLAST V LAVALOVĚ DÝZE B. FRODLOVÁ1, P. RUDOLF2, L. ZAVADIL3, M. KOZUBKOVÁ4, J. RAUTOVÁ5 1,2,3,4,5

Vysoká škola báňská – Technická univerzita Ostrava, 2VUT v Brně

Abstract:In the introduction, the contribution shows the term of cavitation and partial wetting of the walls. Next step is how the partial wetting of the walls influences the flow in Laval nozzle and how the cavitation field which is forming behind the nozzle narrow is influenced. Dimensions of the nozzle and boundary conditions for numerical experiments are used from experimental measurement which was done at VSB-TUO. Cavitation model was used for numerical modelling of flow with boundary condition of partial wetting on the wall. This model is also compatible with Mixture model and Euler´s model for multiphase flow. Keywords: Cavitation field, partial wetting,Laval nozzle narrowing. 1

Úvod Kavitace

Kavitace je složitý jev charakterizovaný vznikem a zánikem dutin v proudící kapalině. Vzniká při snížení tlaku kapaliny na tlak nasycených par, odpovídající teplotě kapaliny. Za těchto termodynamických podmínek se začne kapalina odpařovat a tvoří se velmi malé bublinky. Tyto bubliny jsou unášeny proudící kapalinou, a jakmile se dostanou do oblasti vyššího tlaku, pára v bublinách kondenzuje a vznikají kavitační dutiny. Do těchto dutin vniká okolní kapalina velkou rychlostí a po zaplnění dutiny dochází k velkému rázu. Kavitační bubliny, které vyplňují část proudu kapaliny, tvoří kavitační oblast. Při malém poklesu tlaku pod tlak nasycených par vzniká počáteční kavitace, jejíž oblast se periodicky zvětšuje a zmenšuje, jedná se tedy o děj dynamický. Pokud se tlak bude dále snižovat, kavitační oblast se zvětší a ustálí. Při popisu kavitační oblasti se zkoumá její tvar, místo výskytu a stálost. Účinky kavitace na materiál se nazývají kavitační napadení nebo rozrušení. Částečná smáčivost Právě materiál a jeho povrch ovlivňují proudění kapalin, ať už svými fyzikálními vlastnostmi nebo dodatečnou povrchovou úpravou. Částečná smáčivost je vlastnost, kdy kapalina špatně smáčí povrch druhé, tomto případě pevné, fáze. Kapalina tedy neulpívá na pevném povrchu. Rozdíl při proudění kapalin kolem smáčivého a částečně smáčivého povrchu je zřejmý. Při dobrém smáčení kapalina ulpívá na povrchu tělesa a relativní rychlost proudící kapaliny je na stěně nulová. Ovšem při částečné smáčivosti stěn, kdy platí podmínka prokluzování kapaliny na povrchu tělesa, není relativní rychlost kapaliny na stěně nulová. Pro obecně zakřivené plochy ve styku s kapalinou odvodil prof. Pochylý vztah (1), kdy se předpokládá, že vektor smykového adhezního napětí leží v rovině určené vektorem vnější normály k povrchu a vektorem rychlosti (Pochylý, 2006), viz obrázek 1 .




(1)

kde k je adhezní součinitel [Pa∙s∙m-1]. Na základě této představy lze předpokládat, že vektor smykového napětí na částečně smáčivém povrchu je úměrný rychlosti kapaliny.

Obrázek 1 - Smykové napětí na obecně zakřiveném povrchu

Aby bylo možné tuto vlastnost zkoumat na složité geometrii Lavalovy dýzy, byl nejprve proveden numerický experiment. Při něm byl zkoumán vliv částečné smáčivosti povrchu na proudové pole v potrubí kruhového průřezu při laminárním proudění, poté byly výsledky porovnány s teorií. 1.1 Rychlostní profil v potrubí kruhového průřezu při laminárním proudění s uvažováním částečné smáčivosti stěn Na částečně smáčivém povrchu je předpokládána nenulová relativní rychlost kapaliny. Rychlostní profil v potrubí kruhového průřezu při laminárním proudění s částečně smáčivou stěnou tedy bude odlišný od rychlostního profilu v potrubí se smáčivou stěnou, viz obrázek 2.

Obrázek 2 - Rychlostní profil pro smáčivý (vlevo) a částečně smáčivý povrch (vpravo)

Vztah pro výpočet rychlostního profilu v potrubí kruhového průřezu v cylindrických souřadnicích, rovnice (2) vlevo, lze odvodit z Hagen-Poiseuilleho zákona. Maximální rychlosti dosáhne kapalina v ose potrubí. Vzhledem k tomu, že profil má tvar paraboly, lze jednoduše




odvodit i střední rychlost kapaliny v potrubí, rovnice (2) vpravo, jelikož ta je rovna polovině rychlosti maximální:

(2)

Podstatným rozdílem u částečně smáčivých stěn je nenulová rychlost na stěně potrubí, což se projeví na snížení tlakové ztráty v potrubí. Vztah pro výpočet rychlostního profilu je pak dán vzorcem (3) vlevo. Tak jako v předchozím případě má kapalina maximální rychlost v ose potrubí. Střední rychlost nelze stanovit z úvahy, že je rovna polovině rychlosti maximální, tak jako tomu je u plně smáčivé stěny. Z obrázku 2 lze vyčíst, že rychlostní profil není tvořen pouze paraboloidem, musí být rozdělen na válcovou část a paraboloid. Střední rychlost je pak úměrná celé válcové části a polovině paraboloidu (Pochylý, 2006):

(3)

Úloha částečně smáčivých stěn v potrubí kruhového průřezu měla za úkol ověřit možnost použití okrajové podmínky pro částečně smáčivé povrchy v softwaru Fluent. Protože se podmínka částečné smáčivosti v softwaru nevyskytuje, je třeba ji definovat do výpočtu pomocí UDF funkce, která může být využita jak ve 2D, tak i ve 3D úlohách. Touto UDF se napětí na stěně při výpočtu změní tím, že se rychlost v první buňce od stěny vynásobí adhezním součinitelem k. Změnou parametru adhezního součinitele v UDF lze měnit míru částečné smáčivosti stěn, přičemž platí, že čím menší hodnotu adhezní součinitel má, tím je povrch méně smáčivý a kapalina na něm více prokluzuje, obrázek 3.

Obrázek 3 - Srovnání rychlostních profilů z výsledků z Fluentu a empirických rovnic




Obrázek 3 zobrazuje podobu rychlostních profilů z 3D simulací jak pro smáčivý povrch, tak pro částečně smáčivý povrch při různých hodnotách adhezního součinitele k. Lze vidět, jak skluzová délka se snižující se hodnotou adhezního součinitele roste, povrch se tedy stává méně smáčivý. Při laminárním režimu proudění má rychlostní profil přesně daný parabolický tvar. Z grafu je zřejmé, že maximální rychlost na stěně, kterou lze při použití podmínky částečné smáčivosti stěn a následném snižování adhezního součinitele dosáhnout, je střední rychlost proudění. Naopak při zvyšování hodnoty adhezního součinitele je možné se maximálně přiblížit rychlostnímu profilu proudění kolem smáčivé stěny. Protože při turbulentním proudění v dýze nelze přesně určit tvar rychlostního profilu, musíme vycházet z těchto předpokladů a po výpočtu s určitou hodnotou adhezního součinitele vždy kontrolovat, zda se hodnota rychlosti na stěně pohybuje ve vytýčeném rozmezí relevantních hodnot, tedy od nulové hodnoty rychlosti na stěně po hodnotu střední rychlosti proudu, což by odpovídalo nulovému napětí na stěně. Pro zvyšující se hodnoty adhezního součinitele k, které se mohou blížit nekonečnu, se rychlostní profily shodují. Analyticky jsou tyto úlohy řešitelné, numericky pro větší hodnoty k výpočet diverguje. Divergenci je možné eliminovat zjemněním výpočtové sítě u stěny, avšak toto zjemňování sítě nelze provádět do nekonečna, zvláště pro složitější geometrie. Proto je třeba najít limitní hodnotu adhezního součinitele, při které je rychlostní profil téměř totožný s rychlostním profilem proudící kapaliny kolem smáčivých stěn pro danou hustotu sítě. 2

Experimentální zařízení s Lavalovou dýzou

Kavitační oblast byla pozorována v zúžení Lavalovy dýzy, která je vyrobena z průhledného materiálu Tecanat. Geometrie potřebná pro numerické výpočty vycházela zcela z geometrie dýzy použité při fyzikálním experimentu, viz obrázek 4.

Obrázek 4 - Geometrie Lavalovy dýzy

Kavitační oblast se jevila jako nestálá, v průběhu fyzikálního i numerického experimentu měnila svoji velikost i tvar a bylo možné sledovat periodické kmitání této oblasti. Tvar a délka kavitační oblasti je na obrázku 5.


90



Obrázek 5 - Kavitační oblast v Lavalově dýze (Babejová, 2010)

3

Matematický model pro řešení kavitace v dýze

Kavitace vzniká při proudění v oblastech vyšších Reynoldsových čísel, proto proudění v dýze s Re = 75000 považujeme za turbulentní. Úlohy byly řešeny metodou konečných objemů s použitím dvourovnicového turbulentního modelu k-ε RNG. Modelovaná úloha byla řešena jako časově závislá. Kavitace, díky jednotlivě modelovaným fázím vody a páry, je definována tzv. vícefázovým modelem proudění, kdy byl použit viskózní model směsi - Mixture. Jednotlivé fáze byly voda a pára, popřípadě další modelovanou fází může být vzduch. Software Ansys Fluent nabízí několik modelů kavitace, které jsou odlišné v přístupu řešení a zadávání vstupních parametrů – kavitační model Singhal, Schneer-Sauer a Zwart-Gerber-Belamri. U výpočtů tedy bylo proudící médium uvažováno jako směs vody a páry. Na vstupu do Lavalovy dýzy byla nastavena průtoková podmínka (mass-flow-inlet) odpovídající průtoku 3 kg·s-1. Na výstupu z trysky byla definována tlaková podmínka (pressure outlet) 105000 Pa. Ostatní hranice byly definovány jako stěna (wall), viz obázek 6.

Obrázek 6 - Nastavení okrajových podmínek na výpočetní geometrii dýzy

4

Výsledky vlivu částečné smáčivosti stěny dýzy na proudění a kavitační oblast v dýze

Zajímavým výsledkem numerických experimentů je srovnání tvaru a délky kavitační oblasti u smáčivého a částečně smáčivého povrchu. Z obrázku 7 a obrázku 8 je zřejmé, že čím větší byla hodnota adhezního součinitele, tím více se délka kavitační oblasti blížila k podmínce „no slip“, která odpovídá plně smáčivému povrchu. Naopak čím menší byla hodnota adhezního součinitele, tím byla kavitační oblast delší.


91



Obrázek 7 - Délka kavitační oblasti z výsledků z Fluentu

Obrázek 8 - Kavitační oblast u různých povrchů stěny dýzy – plně nesmáčivá stěna při τ = 0, částečně smáčivé stěny s k = 2, 5, 15 a 20 a plně smáčivá stěna s podmínkou „no slip“

Na následujícím obrázku 9 je možno sledovat rychlostní profily v příčném řezu v blízkosti vstupu do dýzy. Na částečně smáčivé stěně není rychlost nulová, což odpovídá okrajové podmínce. Maximum rychlosti nepřevyšuje maximální rychlost pro smáčivou stěnu. Pro nesmáčivou stěnu je rychlost konstantní v celém průřezu.




Obrázek 9 - Srovnání rychlostních profilů před zúžením dýzy

5

Závěr

Cílem práce bylo využít teorie částečné smáčivosti povrchů obtékáných kapalinou při numerickém modelování. Částečná smáčivost povrchu byla definována UDF funkcí a otestována na případě laminárního proudění v trubce. Druhým úkolem bylo zkoumat vliv částečné smáčivosti stěn na velikost kavitační oblasti, přitom bylo využito geometricky jednoduchého zařízení, tj. Lavalovy dýzy, kde je kavitace běžným jevem. Byla potvrzena domněnka, že oblast kavitace se zvětšuje při snižující se smáčivosti, tj. snižujícím se součiniteli smáčivosti. Literatura BABEJOVÁ, S. „Experimentální zařízení pro vizualizaci vzniku a vývoje kavitace”. Diplomová práce. 69 s. Ostrava 2010. BIRD, R. B., STEWART, W. E., LIGHTFOOT, E. N. „Transport Phenomena”. 2nd Edition. New York: John Wiley & Sons, Inc. 2002. 914 p. ISBN 0-471-41077-2. BOJKO, M., KOZUBKOVÁ, M., BABEJOVÁ, S., BÍLÝ, D. „Experimentální a numerické řešení vzniku a vývoje kavitace v Lavalově dýze”. In sborník mezinárodní konference „XXIX. setkání kateder mechaniky tekutin a termomechaniky“. Rožnov pod Radhoštěm: VŠB-TU Ostrava, 23.26. června 2010. s. 9-12. ISBN 978-80-248-2244-0 BRENNEN, C. E. „Fundamentals of Multiphase Flow”. 1st Edition. USA: Cambridge University Press, 2005.345 p. ISBN-13 978-0-521-84804-6. KOZUBKOVÁ, M., „Modelování proudění tekutin FLUENT, CFX”. Ostrava: VŠB-TUO, 2008. 153 s. [online]. Datum poslední revize 12. 12. 2008. Dostupné z http://www.338.vsb.cz/PDF/Kozubkova-Fluent.pdf. KOZUBKOVÁ, M. „Matematické modely kavitace a hydraulického rázu“. 1. vyd. Ostrava : VŠBTUO, 2009. 130 s. ISBN 978-80-248-2043-9. POCHYLÝ, F., HABÁN, V., „Smáčivost kapalin vůči pevným povrchům”. Výzkumná zpráva VUTEU13303- QR-14-06, Brno: VUT v Brně, FSI, 2006, 36 s. POCHYLÝ, F., FIALOVÁ, S., KOZUBKOVÁ, M., ZAVADIL, L. „Assessment of cavitation creation depending on the surface wettability”. Proceedings of the 25th IAHR Symposium on Hydraulic




Machinery and Systems. Volume 12 of the IOP Conference Series: Earth and Environmental Science, 012106 doi:10.1088/1755-1315/12/1/012106




PROBLEMATIKA NUMERICKÉ SIMULACE FENOMÉNŮ SPOJENÝCH S ŽÁROVÝMI NÁSTŘIKY ING. TOMÁŠ ZÁBRANSKÝ, ING. JAN HAMZA Výzkumný a zkušební ústav Plzeň s. r. o.

Abstract: Paper presents ways of numerical simulations and processes occurring in flame spray technology mainly in the combustion area. Problems regarding to chemical reactions and heat transfer under high velocities appear by modeling these processes with ANSYS FLUENT Keywords: CFD modeling, flame spray, heat transfer. 1

Úvod

Technologie žárového nástřiku je technologie umožňující nanášení povlaků o tloušťce obvykle větší než 50µm. Nanášený materiál, obvykle se jedná o prášek, je nataven a urychlen směrem k povrchu v plameni zařízení.Po dopadu na povrch se nanášený materiál rychle ochlazuje a tuhne. Na povrch tedy dopadají celá natavené zrnka materiálu (kapky). Hlavními parametry technologii žárových nástřiků, ovlivňující technologické parametry povrchu jsou: 

Natavení částice a s tím související doba setrvání v plameni a teplota plamene



Rychlost dopadu částice na povrch

Žárové nástřiky vyhrávají nad ostatními technologiemi nanášení povlaků jednoduchostí, cenou a širokou škálou nanášených materiálů. Hlavním problémem je menší dopadová rychlost částic a pro nanášení keramiky malá teplota plamene a z toho vyplývající malé natavení dopadajících částic. To vede k horším parametrům povlaků než u jiných technologií. Z těchto důvodů se jeví jako perspektivní vyvinout zařízení zlepšující parametry žárových nástřiků za současného zachování jednoduchosti a nízkých provozních nákladů.




Obrázek 1 - Schéma původního zařízení

Obrázek 2 - Schéma nového zařízení 2UZ nástřiku plamenem

Použitím tohoto zařízení se předpokládá zlepšení depozičních parametrů a tím i parametrů nanášených povlaků. 2

Popis úlohy

Nástřikové zařízení používané v současnosti ve firmě VZÚ Plzeň je typu METCO 6P-II a výrazně se neliší od jiných zařízení pro nástřiky plamenem. Schéma zařízení je na obrázku 3 (obrázek použit z METCO online brožury zařízení). Důležitou vlastností zařízení je možnost plynule regulovat průtoky plynů hořících v plameni i plynu dopravujícího do plamene nanášené částice. V tomto případě jde o hoření acetylenu s kyslíkem, jejichž předmíchaná směs vystupuje z trysky vnějšími otvory a dusíku jako nosného plynu viz obrázek 5

Obrázek 3 - Schéma zařízení pro žárové nástřiky.




Obrázek 4 - Fotografie zařízení

Obrázek 5 - Koncová část trysky

3

Modelování a výsledky V oblasti plamene dochází současně k několika fyzikálním jevům. Tyto jevy obsahují například: 

Proces hoření



Proudění směsi plynů



Přestup tepla



Vícefázové proudění obsahující letící částice






Natavování letících částic a jejich ohřev



Interakce dopadajících částic s povrchem, na který jsou nanášeny



Interakce letících částic s vnitřním povrchem navrhovaného zařízení (obrušování, zanášení)

Celý problém byl řešen jak trojrozměrně, tak i dvojrozměrně jako osově souměrná úloha. Z důvodů časové náročnosti byly vývojové výpočty prováděny jako 2D axisymetrický model. a. Proces hořeníProces hoření v plameni byl modelován jako směs plynů s uvažováním chemické reakce (volumetric reaction) acetylenu s kyslíkem o následující rovnici: 2 C2H2 + 5 O2 → 4 CO2 + 2 H2O Tato reakce postihuje stav, kdy je plamen nasycený a v oblasti hoření je dostatek kyslíku. Tento stav nastává v procesech žárových nástřiků jen pro nanášení některých povrchů. FLUENT má při použití této možnosti problém s tím, že reakce se sama nezpustí. To se v praxi obchází imaginární jiskrou zadanou jako „patch“ teploty před výpočtem v malé oblasti. Nastavovaná teplota musí být vyšší než teplota hoření acetylenu. b. Model proudění Hlavním úkolem vyvíjeného nasazeného urychlovacího a usměrňovacího zařízení je zvýšení rychlosti za výstupem ze zařízení, a pokud to bude možné i předvídatelnost rychlostního profilu v různých vzdálenostech. Problém byl modelován jako stlačitelné turbulentní proudění. Rychlosti se v modelu blíží 1 machu pro pokojové teploty. Model je výrazně turbulentní nejen vlivem vysokých rychlostí ale také chemické reakce, která zvyšuje intenzitu turbulence v modelu. Vypočtený profil je na obrázku 6.




Obrázek 6 - Rychlosti v oblasti za tryskou s nasazeným zařízením

c. Modelování přestupu tepla Jedna z důležitých věcí pro úspěšné fungování vyvíjeného zařízení je účinné chlazení jeho povrchu, což se týká i vhodného výběru materiálu pro jeho výrobu. V našem případě byl vybrán materiál s dobrou tepelnou vodivostí a jako chladící medium byl použit vzduch. Pro modelování byla zvolena jednoduchá válcová stěna a obecný průřez pro chladící medium. Oblastí ohraničené válcovou stěnou protéká horký vzduch o teplotě 3000 K (oblast hoření) a vnější stěna je obtékána chladícím médiem o teplotě 300 K (oblast chlazení). Jako horní mez pro teplotu zařízení nebyla brána teplota tání zvoleného materiálu, ale teplota nižší, nad níž již materiál ztrácí své mechanické vlastnosti. Byl použit k-epsilon model. Proměnná v jednotlivých simulacích byla rychlost chladícího media, kdy se začalo s nulovou rychlostí (zde musela být zapnuta gravitace, jelikož se jedná o volnou konvekci), která se postupně zvyšovala až do hodnoty, kdy teplota zařízení klesla pod požadovanou mezní hodnotu. Jak je patrno z obrázku 7, pro oblast hoření a oblast chlazení byla použita inflační metoda síťování směrem od rozhraní dané oblasti s povrchem zařízení. Pro tvorbu výpočetní sítě samotné válcové stěny inflační metoda použita nebyla, jelikož z důvodu dobré tepelné vodivosti zvoleného materiálu je teplotní pole téměř konstantní. Po odladění simulace na válcové stěně byla provedena simulace na modelu reálného zařízení se zanedbáním složitostí tvaru (obrázek 8).




Obrázek 7 - Teplotní pole příčného řezu výpočtovou oblastí se zobrazením sítě

Obrázek 8 - Teplotní pole v podélném řezu výpočtovou oblastí

4

Závěr




Popsané simulace budou využity pro vývoj přídavného zařízení aparatury pro žárové nástřiky, na kterém budou provedeny experimenty pro porovnání výsledků jednotlivých simulací. Další práce budou zaměřeny na vícefázové simulace, konkrétně na započtení diskrétních částic, jejich natavování a depozice na substrát. V oblasti modelování přestupu tepla se další práce budou soustředit na simulaci činnosti reálného zařízení s uvažováním složitosti tvaru. Literatura FIALA J., MENTL V., ŠUTTA P., „Struktura a vlastnosti materiálů,” Academia, Praha 2003. SAZIMA M., KMONÍČEK V., SCHENELLER A. a kol., „Teplo,“ SNTL, Praha 1989. ANSYS FLUENT 12.0 User’s Guide (2009) ANSYS FLUENT 12.0 Theory Guide (2009) Poděkování Příspěvek vznikl v rámci řešení výzkumného projektu se státní podporou MPO ČR č. FRTI2/561(Inovace technologie nástřiku materiálu plamenem). Kontaktní adresa: Tylova 46, 301 00 Plzeň, e-mail: [email protected], [email protected]




CFD CODE AS A HELPER FOR THE DESIGN OF THE WATER JET CUTTING DEVICE ZDENĚK ŘÍHA Ústav geoniky AV ČR, v.v.i, Abstract: The efficiency of the amplification of high pressure pulsations in the geometric transition between larger and smaller pipes is analysed in this paper. Three types of geometric transitions are compared in the paper: step change, conical and radius change of pipe diameters. In other words, water flow in a device for water jet cleaning and cutting is evaluated. There are ultrasonic pulsations of high pressure in the given domain. Keywords: water jet cutting, ultrasonic pulsations of high pressure, geometric transition between pipes. 1

Introduction

Our work is focused on high pressure compressible water flowing in the domain of a device for water jet cutting. The inner space of the device filled by water is shown in Figure 1. We can say that such the device is mostly a system of cylinders and cones placed consecutively in line, except for the input hole, see Figure 1. It is possible to find front wall of a vibrating body termed waveguide, which brings into the fluid flow high pressure pulsations through its periodic linear length changes (mainly in the direction of rotation axis). Then, pressure pulsations are transferred into the output nozzle, where the energy of pressure transforms to kinetic (velocity) energy and dissipation energy. Pulsations of the energy create pulsations of the water jet at the exit from the nozzle. Hence the water jet divides into clusters of droplets at certain standoff distances from the nozzle exit, see (Foldyna, 2004). This type of jet is called a pulsating water jet. The advantage of such a pulsating jet over a continuous one is based on the fact that the initial impact of clusters of droplets of the pulsating jet on the target surface generates impact pressure pi that is several times higher than the stagnation pressure ps generated by the action of a continuous jet under the same working conditions. In addition, the action of a pulsating jet induces fatigue stress in the target material due to the cyclic loading of th e target surface. This further improves the efficiency of the pulsating liquid jet in comparison with the continuous one (Foldyna, 2004). The shape of the transition region between large and small pipes could strongly influence the value of pressure pulsation amplification. The influence of selected possible shapes of the transition area between larger and smaller pipes on the amplification of pressure pulsations is investigated. The efficiency of the high pressure pulsation amplification is studied with regard to the shape of the transitional region.




Figure 1 – The device and its inner space for the water jet cutting and given transition regions

2

Configurations

We can choose some realizable interval of transition types, because it is very hard to find the best solution from all possibilities available. In the first step, we can define three basic transitions with regard to their geometry and manufacturing (see Figure 1): 1. simple step change of diameters; 2. conical change of diameters; and 3. radius change of diameters. We would like to find a certain configuration of the geometry, which will be used in standard practice with a high efficiency of pressure pulsations transfer. Therefore, the requirement of simple manufacturing represents one of the main parameters to be taken into account. The above-mentioned transitional regions create mutually different variants of the geometry and they also offer the possibility of simple manufacturing. The device for water jet cutting consists of several parts, such as main body, waveguide, nozzle holder, etc, see Figure 1. Body and its inner space represent the main part of the tool assembly. The front wall of the waveguide vibrates axially and, due to this change of the geometry, in time creates pulsations of pressure in the flowing water. The inner space of the nozzle produces the final water jet. The total length of a large diameter cylinder can be adjusted to tune a given domain. The place where water flows from the inner space of body to




the nozzle holder, represents our crucial point of interest. The above-mentioned transitional regions have been placed at that point.

3

CFD model

The water flow in the tool for water jet cutting is a complicated physical problem. Turbulent fluid flow was physically defined using the Navier-Stokes equation adjusted by the RANS method. A realizable k- turbulence model was used in the calculations. The model should be able to describe with sufficient accuracy the given type of regions with low and high velocities. The following system of equations was solved:

(1)

The first expression (1) shows the shape of the continuity equation for compressible flow, where: 

is density,

t

is time,

xj

is component of coordinate,

ui

is component of velocity.

(2)

The second equation (momentum equation) shows the special shape of the NavierStokes equation for turbulent fluid flow, where: p

is pressure,



is viscosity,

t

is turbulent viscosity,

k

is turbulent kinetic energy.

Equations (1) and (2) also have been called Reynolds-averaged equations. The terms in equation (2) placed on the second line define Reynolds stresses. There are other



unknown

parameters,

that

can

be

described


by

equations

(3),

(4),

and

(5):

(3)

(4)

(5)

where: ε

is turbulent dissipation energy,

for other terms 13.0, 2010).

see

(ANSYS:

ANSYS

FLUENT

Theory

Guide,

Release

Equations (1)-(5) create a closed system for one phase compressible turbulent isothermal fluid flow. More exactly, the above-mentioned system of equations needs to define compressibility ((p), a(p)) and equations describing turbulent flow need to be complemented with some equations describing fluid flow at the wall. This is necessary because the ratio of molecular viscosity and turbulent viscosity is completely different at the wall and in the middle stream. In this case, equations describing turbulent flow have been completed with nonequilibrium wall function, see (ANSYS: ANSYS FLUENT Theory Guide, Release 13.0, 2010). This approach substitutes exact solution of the fluid flow in the given sub-layers of boundary layer with the help of empirical formulas. Equations (6) and (7) complete the above-mentioned system of equations with the possibility to solve multi phase flow,

(6)

(7)

where: q

is volume fraction of the qth phase,

is density of the qth phase, for other terms, see (ANSYS: ANSYS FLUENT Theory Guide, Release 13.0, 2010). q

Due to the necessity to calculate water flow as compressible, it is obligatory to complement the above-mentioned equations with relations for density and sound speed as a function of pressure. It is possible to use equation (8) as a base to derive equation (9).




(8)

(9)

where: K

is bulk modulus of water.

For the description of the sound speed, a simple linear equation (10) will be used, see [3]. abpd

(10)

Where: a

is sound speed,

b, d

are constants and p is absolute pressure.

In the CFD model, the part termed waveguide was represented only by the front circular area, see Figure 1. The space or distance between the large pipe diameter and the diameter of the waveguide arm was neglected. The movement of the head circular waveguide wall has been defined by the following equation:

v pulsatorA cos t  ,

(11)

where: vpulsator is velocity moving of the waveguide wall in direction of its rotational axis, A is amplitude of oscillations,  is angular velocity, T is time.

Equations (9), (10), (11) have been implemented in CFD code as user-defined functions. These equations create a closed system and allow the solution of 3D two-phase turbulent compressible flow of water in the tool and out of the tool in incompressible air. The effect of cavitation was neglected in this case. The above-mentioned partial differential equations are solved with the help of the numerical method known as method of control volumes. Using these procedures, a suitable model for studying high pressure pulsations propagation in the given geometry was created.




Figure 2 - Inner space of the nozzle

4

Calculated results

The nozzle, or its cylinder with the smallest diameter, is the place where pressure energy or pressure pulsations are transformed to kinetic and dissipation energy, see Figure 2. Pulsating kinetic energy of water in the nozzle cylinder with the smallest diameter creates the final shape of the water jet at the output of the tool or of the nozzle. It is valid that higher amplitude of pressure pulsations in front of the output nozzle produces higher amplitude of velocity in the smallest cylinder of the output nozzle. Therefore, the evaluated parameter was velocity in the smallest cylinder of the nozzle, more exactly, maximum value of velocity on the depicted cross- sectional area, see Figure 2. Figure 3 illustrates evaluated velocities in time. It is possible to see that the variant with radius transition between big and small diameter (see Figure 1) produces maximum amplitude of normalized velocity vz. Then, the variant with cone transition follows, and the variant with diameter jumping produces minimum amplitude of normalized velocity vz, see Figure 3. The size of the normalized velocity amplitude defines the efficiency of the high pressure pulsations transfer in the given geometry of the tool for water jet cutting, and also defines efficiency of cutting. Higher amplitudes of velocity at the output from the output nozzle will produce pulsating water jet with well-developed pulses at certain standoff distances from the nozzle exit. Subsequently, such a pulsating water jet will produce higher impact pressures on the surface area of the cut materials. Contrarily, low amplitude of normalized velocities will produce poorly developed pulsating water jet, with very low impact pressures on the surface




area of the cut material, resulting in low cutting efficiency, and its effects will roughly correspond to the effects of a continuous water jet.

Figure 3 - Maximum values of the normalized velocity vz in direction of the output nozzle rotational axis (cross-sectional area in the smallest cylinder of the output nozzle)

5

Conclusion

This study shows that it is also necessary to choose the correct shape of the transition region between large and small pipe diameter, if high efficiency of cutting is requested. The design can be supported by numerical modeling of flow with the help of CFD code ANSYS. The objective of the calculations was to determine the amplitude of velocity in the smallest cylinder of the nozzle as a response to high pressure pulsation, and to compare given geometrical variants with regard to the efficiency of high pressure pulsation transfer in the tool for water jet cutting. The highest amplitude was calculated for variants with the radius transition. It is possible to expect that the variant with radius transition will produce the best pulsating water jet for given boundary conditions (input pressure, amplitude of waveguide vibrations, nozzle diameter, etc.). The other two variants of the transition region produce lower amplitudes of velocity at the nozzle exit, that will result most likely in poorly developed pulsating water jets. In this case, a higher amplitude of vibrations of the waveguide would be required to obtain the same amplitude of velocity as in the case of radius transition. This would require more energy for vibrations of waveguide and better material of the waveguide, due to the presence of higher stresses. Three fixed geometrical variants have been resolved here. Generally, it is possible to find cases with the highest pulsations of pressure and velocity if the length of the large cylinder




is changed by nozzle holder moving with respect to the body. In other words, global extremes of pulsations should be found for given variants, and then such results could be compared with each other. This approach will be used in future studies of the problem. References ANSYS: ANSYS FLUENT Theory Guide, Release 13.0, November 2010, www.ansys.com. FOLDYNA, J., SITEK, L., ŠVEHLA, B., ŠVEHLA, Š.: Utilization of ultrasound to enhance highspeed water jet effects. Ultrasonics Sonochemistry, Vol 11/3-4, pp 131-137, 2004. Acknowledgements The research reported here has been carried out in connection with the Institute project on clean technologies for mining and utilization of raw materials for energy use, reg. no. CZ.1.05/2.1.00/03.0082 supported by Research and Development for Innovations Operational Programme, financed by the Structural Funds of the European Union, and from the state budget of the Czech Republic. This work was also supported by the Academy of Sciences of the Czech Republic, project No. AV0Z30860518. The authors are thankful for the various sources of support for their research. Contact address: Drobného 28 Brno 60200, [email protected]




NUMERICKÁ SIMULACE PROUDĚNÍ STUPNĚM S VYROVNÁVACÍMI ŠTĚRBINAMI BARTOLOMĚJ RUDAS, ZDENĚK ŠIMKA, PETR MILČÁK, LADISLAV TAJČ, MICHAL HOZNEDL ŠKODA POWER, A Doosan Copany Abstract: The numerical simulation of flow through the turbine stage including stator and rotor seals and equalizing gaps is presented. The pressure field in particular parts of the stage is described. The main focus is attended to the parts with maximum increasing of entropy production. The distribution of velocities at the inlet and outlet from the rotor blades is observed as well as the flow through equalizing gaps. From these data the final thermodynamical efficiency is determined. Keywords: steam turbine, efficiency, numerical simulation, equalizing gaps 1

Úvod

Koncepce bubnového uspořádání rotoru vede k prodloužení lopatek a ke zvětšení jejich štíhlosti. Výsledkem úpravy je snížení podílu okrajových ztrát a zlepšení termodynamické účinnosti stupně. Aby se omezil negativní dopad úniku páry přes hřídelovou ucpávku na účinnost, používají se pod patou oběžných lopatek odlehčovací štěrbiny, které umožní odvést veškerou páru z hřídelové ucpávky za stupeň. Průtok páry štěrbinami je silně závislý na přetlaku na oběžné lopatce, resp. na patní reakci stupně. Čím větší je přetlak, tím více páry může unikat přes vyrovnávací štěrbiny. Je žádoucí, aby velikost štěrbiny umožnila průtok páry jen z hřídelové ucpávky a aby odsávání páry z lopatkové části stupně bylo nulové. Cílem provedeného výpočtu je prověřit, zda se tento záměr povedl. Numerická simulace může rovněž ukázat na rychlostní a tlakové poměry ve všech částech stupně, může odhalit i určité rezervy v návrhu lopatkové části i vyrovnávacích štěrbin. 2

Výpočtový model turbínového stupně

Výpočtová oblast je zachycena na obrázku 1. Představuje reálný návrh turbínového stupně včetně ucpávek a vyrovnávacích štěrbin. Jedná se o mírně přetlakový stupeň s 3D tvarováním rozváděcích lopatek kombinujícím řízený průtok (controlled flow) se složeným náklonem (compound lean). Základní charakteristické údaje stupně se nacházejí v tabulce 1.




Tabulka 1 Základní charakteristické údaje stupně

Rozváděcí lopatka

Oběžná lopatka

29

20

Tětiva

b [mm]

Štíhlost

l/b [-]

1,44

2,21

Poměrná rozteč

t/b [-]

0,68

0,68

Tabulka 2 Provozní parametry stupně

Obrázek 1 - Výpočtová oblast

Provozní parametry uvažované při výpočtu zachycuje tabulka 2. Hřídelová ucpávka byla modelována dvěma posledními břity, proto parametry páry na vstupu do ní byly nastaveny tak, aby průtok odpovídal skutečnosti. K výpočtu byl použit program ANSYS – FLUENT, verze 13. Numerická síť byla generována pomocí softwarů GridPro a Gambit. Celkový počet buněk dosáhl ~8.106. K výpočtům byl použit pressure based coupled řešič s druhými řády přesnosti. Proudění bylo uvažováno jako vazké a turbulentní. K výpočtů byl použit k-ω SST turbulentní model. Proudícím médiem byla uvažována pára chovající se jako ideální plyn. Okrajové podmínky byly definovány pomocí celkového tlaku a celkové teploty na vstupu do statorové mříže a statickým tlakem na výstupu z rotorové mříže. Protékající množství páry přes hřídelovou ucpávku bylo zadáno pomocí konstantního hmotnostního toku a směrového vektoru v radiální mezeře předposledního břitu.


111


3


Poznatky z numerické simulace

Rozložení celkového tlaku na středním řezu lopatkovými řadami se nachází na obrázku 2. Ke změnám celkového tlaku dochází především v oběžné lopatkové mříži. V rozváděcí mříži se pokles celkového tlaku týká pouze mezní vrstvy na povrchu profilů a úplavů za lopatkami. Na rozhraní mezi pevnou rozváděcí mříží a oběžnou mříží se úplavy stírají. Je to dáno volbou interface “mixing plane”, kde se pro proudění v oběžné mříži nastaví střední hodnoty z výstupu z rozváděcí mříže. Pro sledování pohybu úplavu přes oběžnou lopatkovou mříž by se musel uskutečnit nestacionární výpočet s rozhraním “sliding mesh”. Pro účely této studie je však použitý typ výpočtu dostačující.

Obrázek 2 - Rozložení celkového tlaku ve stupni

Obrázek 3 - Rozložení statického tlaku ve stupni

V oběžné lopatkové mříži lze vystopovat plynulý pokles celkového tlaku. Změna statického tlaku je zachycena na obrázku 3. V tomto případě nastávají rozhodující změny v rozváděcí mříži. V oběžné lopatkové mříži nejsou změny díky nízké reakci dostatečně průkazné. Návrhový stupeň reakce na patě lopatkování je 18%. Tomu odpovídá i rozložení tlaků pod patou lopatek a především proudění přes vyrovnávací štěrbiny.

Obrázek 4 - Rozložení entropie ve vyrovnávacích štěrbinách

Obrázek 5 - Rozložení statického tlaku ve vyrovnávacích štěrbinách

Na obrázku 4 je zachycen průběh entropie na středním řezu přes vyrovnávací štěrbiny a v mezeře mezi břitem hřídelové ucpávky a oběžným kolem. Bezprostředně za břitem je úsek s nárůstem entropie. Břit na rotoru napomáhá vzniku vírové oblasti s pochopitelným nárůstem ztrát a tudíž i entropie. Proud páry vytékající ze štěrbiny nad břitem se ohýbá směrem dolů a vstupuje do jednotlivých štěrbin. Vstupní úhel neodpovídá geometrii vyrovnávacích štěrbin. Z tohoto důvodů dochází na přetlakové


112



straně štěrbiny k odtržení proudu a k nárůstu entropie. Další oblastí vyšších ztrát je pasivní úsek za oběžným kolem. I zde vzniká charakteristická zavířená oblast, kde se pára postupně dostává do hlavního proudu z lopatkové části stupně.

Obrázek 6 - Relativní rychlosti na vstupu do vyrovnávacích štěrbin

Obrázek 7 - Rozložení statického tlaku ve vyrovnávacích štěrbinách

Rozložení statického tlaku před štěrbinou a ve štěrbině ukazuje obrázek 5. Zřetelný je přetlak působící na stěny oběžného kola. V samotném kanálu vyrovnávacích štěrbin je patrný pokles tlaku v místě utržení proudu ve vstupní partii štěrbiny a nárůst tlaku na protější straně. Obdobný efekt v rozložení tlaku nastává v místě zlomu štěrbiny. V tomto místě však rozložení tlaku na povrchu stěny kanálu napomáhá k přenosu energie na rotující hřídel. Ve vstupní části štěrbiny se projevují brzdný účinek toku páry. Detailní pohled na vektory rychlosti na vstupu do štěrbiny ukazuje obrázek 6. Vedle odtržení proudu od stěny je dobře vidět i rozdílnost vstupního úhlu proudu a geometrie štěrbiny. Změny rychlosti proudu páry v místě zlomu štěrbiny ukazuje obrázek 7. Na podtlakové straně proudí pára rychleji, což se projevuje poklesem statického tlaku a tažným účinkem na rotor. Rozložení statického tlaku na jedné a druhé straně vyrovnávacích štěrbin ukazují obrázky 8 a 9. Rozdíly v radiálním směru nejsou podstatné a existují především ve vstupní partii štěrbin. Tok proudu páry ze štěrbiny nad břitem hřídelové ucpávky dopadá především do spodního úseku vyrovnávacích štěrbin.

Obrázek 8 - Statický tlak na podtlakové straně štěrbin

Obrázek 9 - Statický tlak na přetlakové straně štěrbin


113



Rozložení tlaku v jednotlivých řezech podél štěrbiny ukazuje obrázek 10. Rozdíly tlaků v radiálním směru v ústí a v místě ohybu štěrbiny pak ukazuje obrázek 11. Ukazuje se, že by bylo vhodné upravit vstupní geometrii štěrbiny. Jednalo by se zejména o změnu vstupního úhlu štěrbiny, který by lépe odpovídala úhlu proudu z hřídelové ucpávky.

Obrázek 10 - Statický tlak v jednotlivých řezech vyrovnávací štěrbiny

Obrázek 11 - Tlaky v ústí a v ohybu štěrbiny

Z rozboru obrazu entropie na vstupu do oběžného kola, který se nachází na obrázku 12, je vidět, že úseky pod patou a nad špičkou lopatek jsou nejvýraznějším zdrojem ztrát. Jsou to místa s výskytem vírů i třecích ztrát na neaktivních plochách. Rozměry nadbandážové ucpávky jsou dané a nelze je výrazněji zmenšovat. Rozměry pod patou jsou závislé na velikosti vyrovnávacích štěrbin. V daném případě je možné rozsah patního úseku zmenšit.

Obrázek 12 - Rozložení entropie na vstupu do oběžného kola

Obrázek 13 - Rozložení axiální složky rychlosti na vstupu do oběžného kola

Rozložení axiální složky rychlosti na vstupu do oběžného kola se nachází na obrázku 13. Největší rychlosti se vyskytují v místě výtoku páry z hřídelové ucpávky. Tento proud je nasměrován přímo na vyrovnávací štěrbiny. V úseku pod patou se vyskytují místa se zápornou hodnotou axiální složky rychlosti. Svědčí to o výskytu vírů v daném místě. Rovněž v úseku nadbandážové ucpávky jsou ve sledovaném řezu minimální axiální rychlosti. Pro patní i špičkovou oblast jsou charakteristické periodické změny rychlosti, které odpovídají lopatkové frekvenci. V lopatkovém úseku jsou rychlosti vyrovnané. http://aum.svsfem.cz

114


Obrázek 14 - Radiální složka rychlosti před oběžným kolem


Obrázek 15 - Tangenciální složka rychlosti před oběžným kolem

Jak to vypadá s radiální a tangenciální složkou rychlosti ve stejném řezu před oběžnou lopatkovou mříží ukazují obrázky 14 a 15. Pod patou existují úseky, kde pára proudí směrem dolů i nahoru. Výrazný je pohyb páry z hřídelové ucpávky směrem dolu k vyrovnávacím štěrbinám. V nadbandážovém úseku proudí pára k okrajům ucpávkových břitů. V lopatkové části je radiální směr proudění zanedbatelný.

Obrázek 16 - Vstupní a výstupní úhel relativní rychlosti proudu páry na oběžném kole Obrázek 17 - Rozložení hmotnostního toku na oběžném kole


115



Tangenciální složka rychlosti je však v tomto místě maximální. K rotaci proudu dochází i v místech mimo lopatkovou část. Tangenciální složka rychlosti však má menší intenzitu než v lopatkovém úseku. Jaký je vstupní a výstupní úhel relativní rychlosti proudu na oběžném kole ukazuje obrázek 16. Úhel je orientován od osového směru. V místě vyrovnávacích štěrbin je odklon od konstrukčního úhlu 30 a více stupňů. V místech nad 90° se vyskytuje zpětné proudění. Rozdílné charaktery proudění v lopatkové a mimolopatkové části jsou zřetelné. O rozložení hmotnostního toku vypovídá obrázek 17. Zachycuje proudění přes vyrovnávací štěrbiny, lopatkovou částí i únik páry přes nadbandážovou ucpávku. Úniky páry jsou minimální. V celkové bilanci využití energie se však nedají zanedbat.

Obrázek 18 - Průběh reakce na stupni

Obrázek 19 - Účinnost turbínového stupně

Uvažovaný stupeň je navržen s mírně zvětšenou reakcí na patě. Její průběh je uveden na obrázku 18. Není to tedy typický rovnotlaký ale ani typický přetlakový stupeň. Reakce v rozsahu 22-30% je vyhovující z hlediska snížení ztrát v oběžné lopatkové mříži. Jaká je výsledná termodynamická účinnost ukazuje obrázek 19. Ta je vyhodnocena v řezu za stupněm. Zachycuje tudíž páru i z nadbandážové i z hřídelové ucpávky. V převážné části průtočné plochy dosahuje hodnoty 95%. Vliv okrajových ztrát je potlačen vhodným tvarováním lopatkové části i zvětšením štíhlosti při aplikaci bubnového uspořádání rotoru. 4

Závěr

Numerická simulace proudění komplexním provedením turbínového stupně představuje vhodný nástroj umožňující optimalizaci návrhu stupně i kontrolu kvality finálního provedení. Napomáhá k odhalení i případných rezerv v účinnosti. Výpočet prokázal kvalitní provedení lopatkové části i jisté rezervy v návrhu tvaru vyrovnávacích štěrbin.Největší nárůst ztrát je spojen s úseky pod patou a nad špičkou lopatkové části. Je žádoucí minimalizovat rozměry těchto částí. Výsledná termodynamická účinnost splňuje požadavky zadání. Poděkování Příspěvek vznikl za finanční podpory MPO České republiky při řešení grantového úkolu TIP FRTI3/432. Autoři příspěvku děkují MPO za snahu napomoci při vývoji nové generace turbínových stupňů s vyšší termodynamickou účinností.

http://aum.svsfem.cz1.25






A STEADY STATE APPROACH TO CALCULATION OF VALVE PRESSURE RISE RATE CHARACTERISTICS JÁN ORAVEC Technical centre, Sauer-Danfoss a.s., Považská Bystrica Abstract: Pressure relief valves are widely used in hydrostatic circuits. The pressure rise rate characteristics are one of the most important descriptions of the valves performance. This article proposes design of the basic pressure rise rate characteristic of valve using only a few CFD calculations. These calculations are performed for several poppet strokes and flow values. The method is shown for charge pressure relief valve but can be used similarly for any poppet type valve. Keywords: pressure relief valve, hydrostatics, pressure rise rate characteristic, CFD 1

Introduction

Pressure relief valves are essential components in hydrostatic systems. They protect the systems against excessive pressures or maintain pressures at desired levels. The performances of such valves are represented by their pressure rise rate characteristics. A simplified one is shown in Figure 1. Note, the valve remains closed while pressure is lower than the cracking pressure. If the pressure exceeds the cracking pressure, pcr, the valve opens to allow escaping of working fluid (hydraulic oil) to a reservoir. A small slope of the characteristic is desired in standard working range of the valve. If the poppet reaches its maximal stroke limit, the valve behaves like a constant orifice for higher flow rates. Because each the point of the characteristic represents different poppet stroke (also flow, pressure…), the standard calculation (as well as real measurement) leads to transient task including remeshing and solution of equation of motion.

Figure 1 - Simplified pressure rise rate characteristic




If we take into account some simplifications: 

dynamic effects are avoided,



valve at fixed poppet stroke behaves like an orifice which has quadratic function behavior for wide range of flow rates,

we can design the valve characteristic from several CFD steady state calculations. This approach is shown for charge pressure relief valve.

Figure 2 - Design of charge pressure relief valve

2

Charge pressure relief valve (CPRV)

The function of CPRV is to maintain charge pressure of the closed hydrostatic system at a designated level above case pressure (Sauer-Danfoss, 2011). The CPRV is a direct acting




poppet valve (Figure 2). It opens and discharges fluid to the pump case when pressure exceeds the designated level, pcharge, at a channel called gallery.

The valve poppet has a cone which acts against a seat created in the pump end cap (Figure 2). The poppet guiding is provided by a guide designed in the nut. The 2 springs provide force which acts in direction to close the valve. Pretension of one spring can be adjusted by a bolt to reach desired level of the cracking pressure. The pretension force can be calculated from equation: (1)

are evident from Figure 2 As the poppet moves, the spring force corresponds to equation: (2)

where k is spring rate and s is poppet stroke. 3

CFD model

As mentioned above, the steady state approach is used here. The CFD models (Figure 3) are prepared for 3 poppet strokes (0.5 mm, 1.5 mm and 3 mm). There is used combination of hexahedral, tetrahedral and wedge mesh types created using ANSYS ICEM CFD. The CFD analysis is performed in ANSYS Fluent. Material properties of hydraulic oil Shell Tellus 46 are used. The fluid is assumed as incompressible. Based on previous investigation, experience and comparison with testing, laminar model is chosen in this specific case. Pressure outlet boundary condition is applied to the pump housing passage. The pressure pout = pcase is constant for all the calculated cases. Velocity inlet boundary condition is experienced as suitable and stable for this task. Inlet velocity and flow rate are related to each other by simple relation:

(3)

where q in is inlet flow rate, vin is inlet velocity magnitude and A in is inlet area. The velocity magnitude (flow rate respectively) is ramped to several values for each the stroke model. At least, 3 discrete velocity magnitudes are needed for each the stroke to be able approximate results by quadratic function. As a result from each the CFD calculation is: 

pressure at inlet (pin),



flow force induced by oil to the poppet (F).

For reference, pressure contours and streamlines in region of control passage are shown in Figure 4 for 1 mm poppet stroke and 70 l min-1 flow rate case.




Figure 3 - CFD model

-1

Figure 4 - Pressure contours and streamlines in control passage region (1 mm stroke, 70 l min flow)



4


Post-processing

The post-processing plays big role in this task. We have data from 9 CFD calculations (3 strokes × 3 flow rates per each stroke = 9 calculations) in this specific example. Therefore further data processing in a spreadsheet or other suitable tool is needed. The input and output values can be collected into a table form for further processing. One table proposal is following: Table 1 Proposed table for collection of CFD inputs and outputs

The visualization of data points can be seen in Figure 5.The triangle points in graph C represent dependency of the pressure drop on flow rate. The different colors relate to different poppet strokes. Similarly, the square points in graph B represent dependency of the flow force on flow rate. The data points can be approximated by quadratic polynomial trendlines (continuous colored lines in graphs B and C) for each the stroke value. The flow force at constant stroke can be then written in form: (4)

where a,b and c are polynomial coefficients obtained from the spreadsheet trendline properties. It is obvious the flow force has to be in balance with the spring force (2) at each the point of final pressure rise rate characteristic. Therefore we can write:

(5)

Solution of this quadratic equation is:

(6)




where the sign plus or minus is chosen in way to obtain physically correct result. So, using (6) we get flow rate at which there is balance between flow force and spring force at constant stroke (compare graphs A and B in Figure 5 to see balance between spring force and flow force for each the stroke).

Figure 5 - Visualisation of the pressure rise rate characteristic design

Similarly to the flow forces approximation, we can also approximate the pressure drop data points using quadratic polynomial trendlines in form:

(7)

where the polynomial coefficients d,e, and f are obtained again from the spreadsheet trendline properties. If we use the calculated flow rate from (6), we obtain corresponding pressure drop:

(8)

This procedure is repeated for all the 3 strokes so we obtain the flow rate and pressure drop pairs: (9)

plus we have the additional point [qs0 = 0,pcr]. The final pressure rise rate characteristic can be designed from those points – see red line in graph C in Figure 5 . 5

Comparison with measurement




A simplified hydraulic scheme used for measurement of the CPRV is shown in Figure 6a. An auxiliary pump H provides oil to the testing device. The throttle valve TV is fully open so the pressure indicated by the gauge G1 is small in comparison with the cracking pressure of the tested CPRV. As the TV restricts flow area, the load of CPRV increases and it starts to operate. Pressure before CPRV (G1) and behind it (G2) is recorded together with flow rate (Q). Note the measurement process should be sufficiently slow to avoid dynamic effects of the system. It is possible to perform measurement of the CPRV installed in a simple test block or installed in the pump end cap. The measured valve performance with the one calculated by CFD are shown in Figure 6b for comparison.

Figure 6 - Measurement hydraulic scheme (a) and comparison CFD with measurement (b)

6

Conclusion

The design of pressure rise rate characteristic of charge pressure relief valve using steady state CFD simulation is shown in the article. The model is prepared for 3 constant strokes. Minimally 3 different flow rate conditions needs to be calculated for each the stroke due to possibility of data approximation by quadratic polynomial functions. Both pressure drop and flow force functions are approximated in relation to flow rate. The final pressure rise rate characteristic of the valve is found using condition of equilibrium of flow force and spring force at each the specific stroke. Because of computing only a few steady state cases, the proposed method saves computing time as well as effort to set up the analysis in comparison to transient CFD computation. Another advantage is the effect of the cracking pressure or spring rate change is recalculated immediately in the spreadsheet without any additional CFD calculation requirement. References “H1 Axial Piston Pumps Single and Tandem”, Basic Information, Sauer-Danfoss, 11062168-Rev BB-Apr 2011. IMRE M., KRIŠŠÁK P., RAHMFELD R., ZAVADINKA P., “Porovnanie simulačných a meraných výsledkov škrtiaceho ventila”, Hydraulika a pneumatika, 1-2/2010, ISSN 1335-5171. Contact address:



[email protected]





MODELOVANIE SPAĽOVANIA HNEDÉHO UHLIA V KOTLY A ĎALŠIE PRÍKLADY ING. RASTISLAV BENIAK, ING. CSABA LATIKA SES. a.s. Tlmače, Slovenská Republika Abstract: Main section use Ansys- Fluent in our company is modelling combustion lignite and hard coal with focusing on Low NOx by primary arrangements. In this report is mention combustion lignite. In second part are in short show two typical cases from our work: Reviewing influence natural convection air on offtake heat with insulation parts boiler and unsteady solution flow water in flow piping boiler with respect to temperature course on exit from piping. Keywords: Combustion, lignite, natural convection 1

Úvod

Spaľovanie uhlia v súčasnosti stále zohráva významnú úlohu pri výrobe elektrickej energie a tepla aj napriek rozvoju obnoviteľných zdrojov energie a výraznej plynofikácii. Tomto príspevku je popísaný postup pri riešení spaľovania hnedého uhlia v granulačnom ohnisku kotla výkonu 325 t/h prehriatej pary. Na dopravu zomletého prášku do kotla nám slúžia ventilátorové mlyny. Kotol je vybavený štyrmi ventilátorovými mlynmi, z ktorých každý zásobuje palivom jeden práškový horák umiestnený v rohoch spaľovacej komory. Horák je po výške rozdelený na viacej častí. Horáky sú tangenciálne nasmerované na fiktívnu kružnicu v strede spaľovacej komory . 2

Materiál a metodika

Samotný výpočet bol prevedený vo viacerých variantách, kde bol vybratý variant s najnižšími emisiami NOx. Palivo spaľované v kotle je hnedé uhlie s výhrevnosťou 18485 kJ/ kg . Sušenie paliva je zabezpečené spalinami, ktoré sú odoberané v hornej časti spaľovacej komory a teplota primárnej zmesi na výstupe z triediča je vo výpočte uvažovaná 140°C. Úloha bola riešená v 3D a model je prevedený skutočnej veľkosti, model zahŕňal celú spaľovaciu komoru . Geometria spaľovacej komory, povrchová a objemová sieť bola vytvorená programom GAMBIT. Výpočtová sieť je vytvorená zo štvorstenov a obsahuje 655215 elementov. Vykonané výpočty: Variant 0: Pôvodný variant Variant 1: Upravený variant Zadanie počiatočných, okrajových podmienok a fyzikálnych vlastností pre spaľovanie: Výpočet prebiehal pri týchto podmienkach : Stacionárna úloha, turbulentný model kepsilon, radiácia P1, operačný tlak p0=101325, nestlačiteľný plyn. Uhlie je modelované pomocou diskrétnej fázy- DPM model. Prietok uhlia je uvažovaný rovnomerne cez primárne hubice horákov. Na vstupe do kolta bolo uvažované so strednou veľkosťou zrno 90µm.




Elementárne zloženie paliva v pôvodnom stave

V Tabulke 1 sú vypísané hodnoty emisií NOx a hodnoty emisií CO a čo daný variant obsahoval. Na obrázku 2 a 3 môžeme vidieť kontúrové zobrazenie teplôt. Z porovnania jednotlivých variant nám emisie klesli z hodnoty 348 na 277 mg/Nm3, čo predstavuje 20 % pokles.

Tabulka 1 Hodnoty emisií NOx a emisií CO v jednotlivých variantách 0, 1




Obrázok 1 - Geometria spaľovacej komory

Obrázok 2 - Pôvodný Variant




Obrázok 3 - Upravený Variant

1.2 Modelovanie voľnej konvekcie. Úlohou bola namodelovať teplotné pole vonkajšej plechovej steny a bandáže chladiča fluidnej vrstvy v oblasti, kde sú vyvedené hady do komôr, a celá oblasť je izolovaná vonkajšou izoláciou. Na základe teplotného poľa potom posúdiť vhodnosť riešenia resp. materiálu steny a bandáže. Vnútorný priestor môžeme považovať za uzavretú oblasť, kde môže dôjsť voľnej konvekcii. V chladiči fluidnej vrstvy je dokonale zmiešaný popol, ktorý má 900°C. Vo fluidnej vrstve je vysoký koeficient prestupu tepla, tak môžeme použiť podmienku rovno teplotnú, t.j. teplota povrchu steny chladiča bude 900°C Stena chladiča sa skladá z troch vrstiev zámurovky a vonkajšieho plechu. Všetky sú modelované aj s príslušnými koeficientmi vedenia tepla. Tam kde nie sú vyvedené hady, model končí druhou okrajovou podmienkou teplota voľného prúdu 30°C a koeficient prestupu tepla 10W/m2K. Tam ,kde sú vyvedené hady, je oblasť vzdušniny, ktorá je ohraničená stenou chladiča na ľavej strane, v pravých rohoch komorami a vo zvyšku izolačnou vrstvou Nobasil 80. Izolačná vrstva má podobné okrajové podmienky zvonka ako voľný povrch chladiča, t.j. teplota voľného prúdu 30°C a koeficient prestupu tepla 10W/m2K. Vo vnútri oblasti vzdušniny sú ešte rúry a predkomôrky, všetky majú médium rovnakej teploty vo vnútri, a keď si uvedomíme charakter modelu nízkymi tepelnými tokmi, tak môžeme danú podmienku umiestniť aj na vonkajší povrch rúr. Danú úlohu modelovali zjednodušene v 2D. Spracovali sme dve varianty výpočtu, a na základe výsledkov vidieť malá citlivosť na detailné riešenie vnútra pri daných okrajových podmienkach, a tým aj zbytočnosť riešiť danú úlohu v 3D detailne. 1.3 Varianty výpočtu Boli spracované 2 varianty výpočtu. Variant 1 – súčasný stav. 2D rez chladičom krajné rúry sú modelované ako čiary o teplote 276°C, a povrch komôr má tiež 276°C. V chladiči je 900°C a na vonkajšom povrchu je konvekcia pri teplote 30°C a prestupe 10 W/m2K. Variant 2 – súčasný stav. 2D rez chladičom. Všetky rúry sú zanedbané len povrch komôr má 276°C. V chladiči je 900°C a na vonkajšom povrchu je konvekcia pri teplote 30°C a prestupe 10 W/m2K.




1.4 Nestacionárne riešenie prúdenia vody v napájacom potrubí kotla s ohľadom na priebeh teplôt na výstupe z potrubia Cieľom bolo zistiť teplotné gradienty vo vstupnej komore počas jeho nábehu zo studeného stavu. Jedná sa o nestacionárny dej. Bolo modelovaná celá dĺžka napájacieho potrubia. V danom potrubí je voda ktorá má teplotu 55°C a takú teplotu má aj potrubie vrátane izolácie. V čase t=0s začne prúdiť na začiatku potrubia voda o teplote 188,2°C. Potrubie má rozmery

a má dĺžku 158m.

Materiál potrubia je 15NiCuMoNb5-6-4. Izolácia je Nobasil 60mm s hustotou 100kg/m3.

Modelovaná je osmina prierezu potrubia s izoláciou v celej dĺžke v skutočnej veľkosti ale bez kolien. Na konci napájacieho potrubia je prechodový kus a začiatok komory, aby bolo možné zistiť priebeh ohrevu vnútorného povrchu steny komory pri nábehu kotla. Výpočtová sieť je štruktúrovaná hexahedral s 311 000 bunkami. Pozdĺžne roviny rezu sú typu symetry. Začiatok aj koniec steny rúry aj izolácie sú typu wall. Začiatok potrubia je typu velocity inlet a koniec potrubia je typu pressure outlet. Vonkajší plášť izolácie je typu wall s podmienkou prestup tepla pri teplote 15°C s αvz 10 w/m2.K. Vodný objem potrubia je typu fluid priradený materiál water liquid z databáze Fluentu. Plášť potrubia je typu solid, materiál nadefinovaný steel 15NiCuMoNb5 s nasledujúcimi vlastnosťami: hustota 7850kg/m 3 koeficientom vedenia tepla závislej na teplote (podľa zadaných bodov teplota [°C], koeficient [w/m.K]: 20, 38, 100, 41, 200, 42) a tepelnou kapacitou závislej na teplote (podľa zadaných bodov teplota [°C], Cp [J/kg.K]: 20, 430, 100, 500, 200, 540).




Izolácia je typu solid, materiál nadefinovaný Nobasil basalt s nasledujúcimi vlastnosťami: hustota 100 kg/m3 koeficientom vedenia tepla 0,05 w/m.K a tepelnou kapacitou 871 J/kg.K. Modelovaná časť komory má rozmery a má dĺžku 500mm a prechodový kus je plynulý prechod medzi rúrou a komorou s materiálom komory na dĺžke 500mm. Materiál komory je X10CrMoVNb9-1. Plášť komory je typu solid, materiál nadefinovaný steel X10CrMoVNb9-1 s nasledujúcimi vlastnosťami: hustota 7770kg/m3 koeficientom vedenia tepla závislej na teplote (podľa zadaných bodov teplota [°C], koeficient [w/m.K]: 20, 22.6, 50, 23.2, 100, 24, 150, 24.8, 200, 25.4) a tepelnou kapacitou závislej na teplote (podľa zadaných bodov teplota [°C], Cp [J/kg.K]: 20, 456, 50, 472, 100, 493, 150, 515, 200, 531). Výpočet prebiehal sekundovým krokovaním s 20 iteráciami na krok po 1200 s čiže vypočítaný bol prvých 20 minút. V prvej variante je daná rýchlosť vody (zodpovedá prietoku 13,88 t/h. Ďalšie výpočty sú pre prietok 75 t/h (v=2,967 m/s) a 100 t/h (v=3,955 m/s). 3

Výsledky

1.5 Stacionární řešení Ako vidieť z výsledkov v obrázkovej prílohe, v oboch variantoch je dominantná teplota média v rúrkach. Je to spôsobené tým, že vyšetrovaný priestor je dobre odizolovaný. Teplotné rozdiely v skúmanej oblasti sú pomerne malé a teplota v ľavej časti je vyššia ako v pravej. Ak by sme prúdenie úplne zastavili, tak rozdiely by boli vyššie, ale keď si uvedomíme, že celá oblasť je vyplnená výparníkovými rúrami vrátane prechodu rúr cez zámurovku, tak je jasné, že výrazne iná teplota nie je možná nikde v danom priestore ako tých 276°C. Dané konštrukčné riešenie vyhovuje, lebo konštrukčný materiál pri danom zaťažený vydrží dlhodobo nad 290°C.

Obrázok 4 - Vektorové pole vo variante 1




Obrázok 5 - Teplotné pole vo variante 1

Obrázok 6 - Vektorové pole vo variante 2




Obrázok 7 - Teplotné pole vo variante 2

1.6 Nestacionární řešení Teplotu vody sledujeme na piatich miestach rovnomerne rozmiestnením po dĺžke rúry. Miesta sú značené podľa vzdialenosti od vstupu: 31,6 m, 63,2 m 94,8 m 126,4 m a 158 m je značené ako výstup. V prvom grafe vidieť priebehy teplôt vody vo všetkých sledovaných prierezoch a vnútornom povrchu steny komory v čase od 0 s po 1200 s. V druhom grafe je vybraná najdôležitejšia časť t.j. priebeh teploty vody na výstupe a vnútornom povrchu steny komory od 240 s po 840 s. Porovnávacia tabuľka obsahuje maximálne hodnoty rýchlosti ohrevu vnútorného povrchu steny komory v jednotlivých variantoch. Tabulka 2 Maximálne hodnoty rýchlosti ohrevu v [ m/s ] Q [ t/h ]

tOmax [ s ]

Variant 1

0,549

13,88

0,47

285-315

Variant 2

2,967

75

3,1

50-60

Variant 3

3,955

100

4,0

36-46

Model s osminou potrubia neumožňuje modelovať kolená na potrubí, ktoré miestne zintenzívnia prestup tepla. Podľa tabuľky vyplýva, že je možné použiť aj najrýchlejší nábeh, čiže variant 3.




Obrázok 8 - Geometria modelu.

Obrázok 9 - Časový priebeh teplôt vody na konci potrubia a na vstupe do komory od 0s po 600s vo variante 1,2 a 3

Literatura V.ČERNÝ, B.JANEBA, J.TEYSSLER .: Parní kotle Praha 1983 DR.INŽ.JOZEF ČERMÁK .: Spalovací zarízení a kotle Praha 1964 MIROSLAV RÉDR, MIROSLAV PŘÍHODA.: Základy tepelné techniky Praha 1991 Users Guide Help Copyright 2010 ANSYS Inc.




PARNÍ STROJ - VÝPOČTY S POHYBLIVÝMI SÍTĚMI P. ZÁCHA ČVUT v Praze, Fakulta strojní, Ústav energetiky Abstract: This paper briefly analyse problems arising from the steam flow in the piston machine. There are discussed difficulties linked with moving and dynamic meshing and with wet steam formation. At the end of the paper some results of performed simulations are presented. Keywords: piston machine, sleeve valve, wet steam. 4

Úvod

Parní stroj už vynalezli naši předci před čtvrt tisíciletím, ale i v dnešní době stále nacházíme nevyřešené otázky, které se pokoušíme zodpovědět. Firma PolyComp s.r.o. již několik let vyvíjí pístový parní stroj, který by dokázal převést byť jen malou část jinak zmařené energie na energii elektrickou a tím ušetřit nemalé finance. Kde se takové ztráty vyskytují? Představte si, že vlastníte firmu, která pro své technologické procesy vyrábí ostrou páru pro několik rozdílných tlakových hladin. Nabízí se dvě cesty, jak si takovou páru zajistit: a) vyrábět ji pro každou pracovní hladinu samostatně, b) vyrábět veškerou páru pro nejvyšší tlakově-teplotní hladinu a pro ostatní procesy pak parametry páry upravit škrcením. První řešení znamená výrazně vyšší investice, to druhé zase zvyšuje oproti prvnímu náklady na provoz. Často se však jedná již o výrobní procesy v provozu a způsob získání páry byl vyřešen cestou „b)“. Hledání úspor se tedy řeší dodatečně. Jak víme, škrcení je proces, při kterém veškerou přebytečnou energii v páře ztrácíme. Proto se nabízí otázka, zda lze sestrojit takový stroj, který bude natolik investičně nenáročný a bude mít takovou účinnost, že si zajistí svou brzkou ekonomickou návratnost. Pokud je výkon natolik malý, že investice do malé parní turbíny je příliš vysoká, stávají se zde vhodným kandidátem upravené sériově vyráběné pístové motory (např. Liaz), kde úlohu vaček přebírají rotační šoupátka. Novátorské řešení má však svá úskalí, jako jsou např. těsnění šoupátek, teplotní namáhání stroje při startu nebo pevnostní meze klikového mechanismu motoru. 5

Cíl práce

Tento příspěvek shrnuje základní poznatky ze simulace proudění páry v pístovém jednočinném stroji pro dané podmínky.



6


Popis parního pístového stroje

Na Obr.1 je řez pístovým strojem zahrnující válec s pohyblivým pístem, vstupní a výstupní komoru s rotačními šoupátky a kanály spojující komory s válcem. Ve skutečnosti vstupní (admisní) a výstupní (emisní) komora zahrnuje 6 paralelně zapojených šoupátek zajišťujících spojení do 6-ti válců. Vždy dva válce mají shodně natočená šoupátka, další dvojice je pak pootočena vždy o 120°. Tím je zajištěna plynulejší distribuce páry mezi jednotlivé válce a přenos momentové síly na klikovou hřídel. Pístový stroj takto pracuje jako dvoutakt, přičemž na vstupu je pára o vyšší tlakové hladině, na výstupu pak logicky na nižším tlaku. Při otevření vstupního šoupátka pára expanduje do válce (píst se pohybuje směrem dolů) a působením na píst koná práci. Při otevření výstupního šoupátka se pak píst pohybuje směrem nahoru a dochází k vypouštění páry do výstupní komory. Vše se děje při 450 ot./min.

Obrázek 1 - Definice okrajových podmínek a výpočtových oblastí pro upravenou geometrii

1

Materiál a metody

Při bližším pohledu na řešený problém lze konstatovat, že se jedná o nestacionární úlohu s pohyblivou sítí, ve které dochází ke kondenzaci páry. Takto komplikovanou úlohu jsme na ústavu v minulosti ještě neřešili. Vzhledem k výpočetní náročnosti je úloha řešena v dvojrozměrném prostoru. K dosažení správného popisu proudění páry v pístovém motoru pomocí CFD kódu je dále nutné souběžně vyřešit několik dílčích problematik: pohyblivé části - řešená oblast zahrnuje píst a rotační šoupátka; dvoufázové proudění – při prudké expanzi přehřáté páry může dojít k jejímu podchlazení až pod teplotu sytosti a je tedy nutné použít model dvoufázového proudění, který dokáže postihnout vznikající kondenzaci; nestacionární děj – úloha je v čase výrazně proměnná, což vyžaduje při výpočtu použít velmi krátký časový krok.


136



a. Problém č. 1 – pohyblivá síť a štěrbiny Aby při zavřeném stavu docházelo k co nejmenšímu úniku páry, je snaha štěrbinu mezi šoupátkem a bronzovou vložkou (viz Obr.2 a Obr.3) zmenšit na minimum. Výpočty byly prováděny se štěrbinou o tloušťce 0,05 a 0,1 mm. To je o 3 řády menší rozměr, než jsou rozměry komor. Pro oblast šoupátka bylo vhodné použít pohyblivou síť pomocí metody pohyblivé soustavy souřadnic. Rotační část je zde od nepohyblivé sítě oddělena pomocí podmínky Interface. Součástí této podmínky je nutnost vytvořit v jejím okolí několik řad buněk, jenž umožní správný popis proudění v bezprostřední blízkosti. Na Obr.2 je zobrazena výpočetní oblast v okolí šoupátka. Podmínka Interface je dobře patrná na detailu v Obr.3 (červená čára). Z každé strany této podmínky jsou vytvořeny 2 až 4 řady buněk. Výsledkem této první části studie bylo stanovení množství uniklé páry ze vstupní komory v zavřeném stavu.

Obrázek 2 - Výpočtová síť řešené oblasti

Obrázek 3 - Detail výpočtové sítě v oblasti štěrbiny mezi rotačním šoupátkem a bronzovou vložkou

Při zachování podmínky doporučeného maximálního poměru stran buněk (1:5) celou síť, tvořenou řadami buněk v okolí podmínky Interface, výrazně komplikuje. Buňky zde dosahují řádově 0,01x0,05mm. Z těchto důvodů byl při generování výpočtové sítě kladen důraz na minimalizaci počtu buněk při současném zachování kvality popisu proudění v oblasti štěrbiny.


137



b. Problém č. 2 – dynamická síť Dynamická síť se využívá v případech, kdy se mění tvar výpočtové sítě pohybem jedné nebo více stěn, a tedy pro pohyb pístu ideální volba. Je zde využita metoda dynamického vrstvení (Dynamic layering). Ta je založena na spojování resp. dělení nejbližších dvou vrstev buněk u pohyblivé stěny, viz Obr.4. Vrstva j je zde rozdělena nebo spojena s vrstvou i v závislosti na hodnotě výšky h vrstvy j. Funkce dělení resp. spojování se pak již odvíjí pouze od definovaných proměnných (minimální výšky a dělícího parametru). Program k popisu dynamické sítě současně nabízí i definici pohybu pístu na základě specifikace jeho pohybu (tj. rychlost otáčení klikové hřídele, výšku zdvihu a časový krok v podobě úhlu pootočení hřídele).

Obrázek 4 - Metoda dynamického vrstvení

c. Problém č. 3 – nestacionární úloha a čas Celý výpočet byl zatížen vysokou náročností výpočetního času. Při nastavení objemových okrajových podmínek na počátku každého výpočtu dochází v celém modelu k předdefinování uniformní hodnoty všech parametrů (tlak, teplota, relativní vlhkost atd.). Tato nastavení způsobují, že počáteční 2 až 3 otáčky se ustalují hodnoty časově proměnných veličin. Teprve poté bylo možné zahájit sledování konvergenčních kritérií. Kritérium konvergence bylo považováno za splněné, mají-li průběhy obou sledovaných veličin (teplota, tlak) ve dvou po sobě následujících otáčkách stejnou charakteristiku a současně nedochází k výrazným relativním odchylkám. Pro jednotlivé výpočty nepřevýšily maximální relativní odchylky tlaků v jedné poloze u dvou po sobě následujících otáčkách hodnotu 2%, stření relativní odchylky pak byly do 0,3%. U teplot dosahovaly relativní odchylky poněkud vyšších hodnot, a to v místech a časech prudké změny tlakových podmínek. Přesto střední relativní odchylka teplot zde nepřesáhla u žádného výpočtu 5% hranici. Příklad zkonvergovaného průběhu teploty pro oblast válce je znázorněn na Obr.5. Z výše uvedeného vyplývá, že bylo třeba pro každý výpočet vykonat 5-6 otáček. Protože z důvodů stability výpočtu podél pohyblivé sítě bylo nutné udržovat časový krok s pootočením šoupátka o 0,1° (při 450 ot./min tedy 3,704.10-5 s), dostáváme se k hranici 20.000 časových kroků, přičemž každý časový krok zahrnoval 20 iterací. Jeden výpočet na 16-ti procesorech Intel Nehalem znamenal 80 až 100 hodin čistého výpočetního času. Během výpočtů bylo navíc nutné provádět optimalizace konvergentních mechanizmů tak, aby nedošlo k numerickým oscilacím a následné divergenci řešení. Náchylnost k divergenci vznikala zejména v oblasti škodlivého prostoru (je velmi úzký a pohyb poblíž horní úvrati bylo nutné řídit „ručním“ nastavováním podrelaxačních faktorů) a v momentech prudké změny vlhkosti páry (při otevírání šoupátek) - a to byl 4. problém k řešení…


138



Obrázek 5 - Příklad průběhu průměrné teploty ve válci v čase (výpočet č. 4)

d. Problém č. 4 – kondenzace páry Tlakově-teplotní podmínky admisní páry jsou u technologických procesů často blízké teplotě sytosti. Tak tomu bylo i v modelových případech. Při rychlé expanzi, kdy dochází k poklesu teploty pod bod varu, může pára krátkodobě kondenzovat. Proces expanze způsobuje nejprve podchlazení přehřáté páry a následně tvorbu kapiček a tedy přechod do stavu mokré páry. Z tohoto důvodu byl použit dvoufázový model mokré páry. Modelování mokré páry je důležité zejména při analýzách a konstrukci nízkotlakých dílů parních turbín. Pro popis tvorby mokré páry používá program Fluent model Euler-Euler. Míšení páry a vody je modelováno Navier- Stokesovými rovnicemi stlačitelného proudění doplněných dvěma transportními rovnicemi pro stanovení hmotnostního podílu vodní fáze () a počtu vodních kapiček v jednotce objemu (). Model změny fáze zahrnuje tvorbu vodních kapiček v homogenním nerovnovážném procesu kondenzace, který je založen na klasické neizotermní teorii nukleace. Omezení použitelnosti modelu mokré páry je následující: 

je dostupný pouze při použití numerického schématu založeném na metodě korekce hustoty;



nastavení materiálových vlastností vodní páry je omezeno, neboť hodnota vlastností je funkcí použitého modelu;



rozdíl v rychlosti obou fází lze zanedbat (skluz ~ 1), interakce mezi kapičkami je zanedbatelná a hmotnostní podíl vznikající vodní fáze je nízký ( < 0,2);



předpokládá se, že velikost kapiček je velmi malá (v průměru mezi 0,1 až 100 m), tj., objemový podíl vodní fáze je zanedbatelný.

Výše uvedená omezení jsou v dobré shodě s podmínkami řešené úlohy (hodnota hmotnostního podílu  je překračována jen lokálně), a proto lze tento model použít.


139


2


Nastavení výpočtu

Okrajové a provozní podmínky jsou shrnuty v Tab.1. Dále je definováno, že úhel pootočení 0° značí píst v horní úvrati. Nastavení základních parametrů výpočtu jsou uvedena v Tab.2. Uvedené výpočty se lišily buď v tlaku na výstupu resp. nastavením velikosti předstihu. Tabulka 1 Okrajové a provozní podmínky č. výpočtu 2, 4, 5 a 6

Parametr 1 Operační tlak na vstupu [MPa] Vstup

0,1 Tlak na vstupu

typ okrajové podmínky přetlak (gauge pressure) [MPa]

2,0

teplota [°C]

213

Rotační šoupátka

rychlost otáčení [otáčky/min]

450

Píst

rychlost otáčení [otáčky/min]

450

zdvih [mm]

150

délka ojnice [mm] Výstup

typ okrajové podmínky přetlak (gauge pressure) [MPa]

Časový krok [s] / úhel pootočení kliky [°]

3

250 Tlak na výstupu 0,5

1,0

1,5

-5

3,704.10 / 0,1

Tabulka 2 Nastavení základních parametrů pro výpočty Parametr Numerické schéma Časová závislost Model turbulence Stěnová funkce Diskretizační schémata

3

Nastavení Metoda korekce hustoty Přechodový děj Standard k-ε model Standardní stěnová funkce „Upwind“ schémata druhého řádu

Výsledky výpočtů Z provedených výpočtů bylo pro původní zprávu zpracováno: 

stanovení množství uniklé páry ze vstupní komory v zavřeném stavu – byl nalezen vztah mezi velikostí štěrbiny a množstvím unikající páry v zavřeném

stavu; 

indikátorové diagramy - slouží pro stanovení účinnosti celého cyklu při různých variantách (rozdílné předstihy, rozdílné tlakové hladiny) a pro stanovení velikosti sil působících na píst;



průběhy teplot na stěnách - umožňují lépe definovat teplotní namáhání v konstrukčních materiálech motoru (pevnostní výpočty);



velikost hm. podílu vodní fáze – kontrola, zda vlhkost páry nepřesahuje mezní hodnoty ( < 0,2);



videosekvence sledovaných veličin – vizualizace pro lepší pochopení celého procesu.

Významná část výsledků podléhá obchodnímu tajemství, některé výsledky se do příspěvku prostě nevejdou. Proto je dále uveden pouze příklad indikátorových diagramů pro dva různé předstihy při shodné tlakové hladině (Obr.6) a pro jednotlivé sledované veličiny je uveden jeden vybraný stav (Obr.7 až Obr.9).


140



a) b) Obrázek 6 - Indikátorové diagramy a) výpočet č. 2 (pin = 2 MPa; pout = 1 MPa; předstih 7°); b) výpočet č. 4 (pin = 2 MPa; pout = 1 MPa; předstih 5°)

Obrázek 7 - Průběh tlaku v řešené oblasti pro vybraný úhel pootočení


141



Obrázek 9 - Průběh teploty v řešené oblasti pro vybraný úhel pootočení

Obrázek 8 - Průběh hmotnostního podílu vody v řešené oblasti pro vybraný úhel pootočení

Literatura Výkresová dokumentace pístového stroje. GAMBIT 2.4 User’s Guide. Fluent Inc. May 2007. FLUENT 6.3 User’s Guide, Volume 1, 2, 3 and 4. Fluent Inc., September 2006. ZÁCHA P., GREGOR K., ŽELEZNÝ V., „Simulace proudění páry v pístovém stroji“ Výzkumná zpráva č. 12107- JE/2009/01, Praha : ČVUT v Praze, 2009. ZÁCHA P., GREGOR K., „Simulace proudění páry v pístovém stroji – 2. část“ Výzkumná zpráva č. 12107- JE/2010/01, Praha : ČVUT v Praze, 2010. Poděkování Příspěvek vznikl za podpory projektu MPO TIP FR-TI1/521 Kontaktní adresa: Ing. Pavel Zácha, Ph.D. Ústav mechaniky tekutin a energetiky Fakulta strojní



České vysoké učení technické v Praze Technická 4, 166 04 Praha 6 E-mail: [email protected]





CFD SIMULACE TEPLOTNĚ-HYDRAULICKÉ CHARAKTERISTIKY NA MODELU PALIVOVÉ TYČI V OBLASTI DISTANČNÍ MŘÍŽKY D. LÁVIČKA Západočeská univerzita v Plzni, Katedra energetických strojů a zařízení Abstract: This paper looks into detailed investigation of thermal-hydraulic characteristics of the flow field in a fuel rod model, especially near the spacer. The area investigate represents a source of information on the velocity flow field, vortex, and on the amount of heat transfer into the coolant all of which are critical for the design and improvement of the fuel rod in nuclear power plants. The flow field investigation uses commercial CFD software ANSYS/FLUENT. Keywords: fuel rod model, thermal-hydraulic characteristics, spacer, ANSYS/FLUENT, heat transfer. 4

Úvod

Palivový článek je součástí složitého systému, který představuje reaktor jako celek. Na palivový článek je vždy zaměřena největší pozornost, protože patří k nejexponovanějším a technicky nejvýznamnějším částem jaderného reaktoru. Distanční mřížky jsou jedna z důležitých komponent používaných na svazku palivových tyčí. Jedním z hlavních úkolů je držení palivových tyčí a vytváření mezery pro chladící médium v aktivní zóně jaderného reaktoru. Dalším z důležitých úkolů distanční mřížky je plnit funkci turbulizátoru. V místě distanční mřížky dochází k jejímu obtékání, rozdělení proudu podle tvaru distanční mřížky a také k zavíření proudu. Velikost a intenzita zavířené oblasti přímo ovlivňuje přestup tepla ve směru proudění mezi palivovou tyčí a chladícím médiem v aktivní zóně. K vyšetření nejen tohoto úkolu budou použity CFD simulace pomocí komerčního softwaru ANSYS FLUENT. 5

Okrajové podmínky a výpočetní model a. Výpočetní model

Výpočetní geometrie modelu palivové tyče vychází z rozměrů skutečného palivového proutku v jaderném reaktoru. Numerické simulace lze ověřovat na experimentálním zařízení vybudovaného na Západočeské univerzitě v Plzni (Lávička, 2010). Toto experimentální zařízení pro studium přestupu tepla při dvoufázovém proudění obsahuje jednu palivovou tyč, která je tvořena nerezovou trubkou o vnějším průměru 9.1mm. Nerezová trubka je umístěna do skleněné trubice o vnitřní průměru 14.5mm. Takto vzniklý mezikruhový prostor představuje shodný prostor pro proudící chladící médium jako ve svazku palivových tyčí umístěných v reaktoru. Popis výpočetního modelu s umístěním distanční mřížky a základními rozměry průtočného kanálu s nerezovou trubkou je popsán na obrázku 1. V programu GAMBIT 2.4.6 byly vytvořeny 3 výpočetní sítě pro 3 tvary distančních mřížek. Jednotlivé tvary distančních mřížek jsou uvedeny na obrázku 1. Základem modelu je plocha s tvarem distanční mřížky promítnutá kolmo na vstupní plochu. Na této ploše je vytvořena povrchová síť pomocí quad elementů. Počet povrchových elementů se na ploše pohybuje okolo 25000. Povrchová síť distanční mřížky je zobrazen na obrázku 1. http://aum.svsfem.cz



Tyto povrchové elementy jsou vytaženy pomocí funkce „cooper“ (obdoba funkce sweep) v programu GAMBIT do objemové sítě.. Počet elementů na výšku je rozdělen nerovnoměrně s větším zahuštěním buněk v místě distanční mřížky a jejím bezprostředním okolím. Objemová síť je tvořena cca 1,3 miliónem buněk.

Obrázek 1 - Model palivové tyče a navržené tvary distančních mřížek s povrchovou sítí

b. Okrajové podmínky a nastavení řešiče Okrajové podmínky pro numerickou simulaci jsou zobrazeny na obrázku 1. Ve spodní části průtočného kanálu je definována vstupní okrajová podmínka mass-flow-inlet pro zadání průtočného množství. Průtočné množství bylo nastaveno na hodnotu 0,0007kg/s (zhruba 2,5l/hodinu). Tato hodnota průtoku chladícího média je používána pro zaplavování palivových proutků v havarijních stavech např. pokud dojde k úbytku chladící kapaliny-vody v reaktoru v primárním okruhu. Na výstupní ploše byla nastavena okrajová podmínka pressure-outlet. Zdroj tepla byl nastaven na stěny nerezové trubky a distanční mřížky ve velikosti 5000 W/m 2. Na rozhraní okolí a skleněné trubice byl nastaven koeficient přestupu tepla 10W/m2K a teplota okolí 300K. Numerická simulace byla prováděna s nastaveným modelem turbulence „RSM – LowReynolds Stress-Omega”. a v nestacionárním režimu v druhém řádu přesnosti o 4000 časových krocích s časovým krokem 0.05s a 0.1s. Vyhodnocování získaného proudového a teplotního pole probíhalo po posledních 500-ti časových krocích s aktivovanou funkcí datasampling. V každém časovém kroku probíhalo 20 iterací, které stačily pro ustálení konvergenčních kritérií.




Výpočet a vyhodnocení tepelně-hydraulických charakteristik pro jednotlivé distanční mřížky byl prováděn v komerčním softwaru ANSYS FLUENT 13.0. 6

Výsledky

Z numerické simulace byly do vyhodnocení dopočteny další veličiny pro lepší hodnocení distančních mřížek. Některé dopočtené výsledky jsou uvedeny v tabulce 1. Důležitým údajem je tlaková ztráta distanční mřížky, která byla zjišťována odečtením tlaku před mřížkou a za mřížkou. Dalším výsledkem je ustředěná průměrná teplota na ploše uvnitř mezikruhové kanálu v oblasti s distanční mřížkou. V tabulce 1 je z uvedených hodnot vypočten součinitel místní ztráty ξ [-] pro jednotlivé typy distančních mřížek. Hodnota součinitele místní ztráty byla vypočtena podle vztahu (1) a Reynoldsovo číslo Re [-] na vstupu do kanálu podle vztahu (2).

(1)

(2) Další vztah (3) definuje výpočet hydraulického průměru, DH, který je použit pro výpočet Reynoldsova čísla. Reynoldsovo číslo vyjadřuje vliv vnitřního tření v důsledku viskozity dané kapaliny při proudění a také slouží pro určení, zda je proudění laminární nebo turbulentní. Hydraulický průměr je definován vztahem (3) jako poměrem čtyřnásobku vnitřního průřezu daného profilu potrubí „A“ a jeho omočeného obvodu „P“.

(3) , kde D je vnitřní průměr vnější skleněné (plexisklové) trubice a d je vnější průměr nerezové tyče uvnitř modelu palivové článku. Tabulka 1 Charakteristické údaje a výsledky distančních mřížek VAR 01 VAR 02

VAR 03

Tlaková ztráta [Pa]

0.58

0.64

0.65

Rychlost vY* [m/s]

0.00738

0.00756

0.00747

313.2

316.5

313.2

21.3

22.4

23.3

Teplota [K] Součinitel hydraulické (místního odporu) ξ* [-]


ztráty


Reynoldsovo číslo Re [-]

45


45

45

Z uvedených výsledků je patrné, že nejmenší odpor vykazuje distanční mřížka označená jako VAR 01 a naopak distanční mřížka VAR 02 má největší odpor. Tato mřížka také vykazuje největší prohřátí chladící kapaliny na průřezu mezikruhového kanálu a teplota chladícího média je o cca 3K vyšší než u ostatních variant. Ostatní varianty mají shodnou průměrnou teplotu na řezu mezikruhovým kanálem.

Obrázek 2 - Rychlostní a teplotní pole pro distanční mřížky VAR 01 a VAR 02

Hodnota rychlosti z tabulky 1 ukazuje u geometrické varianty VAR 02 zvětšenou rychlost ve směru y oproti variantě VAR 01 a VAR 03. Proud chladícího média u varianty VAR 03 je nejvíce dělen na jednotlivé menší proudy, a proto tato mřížka vykazuje nejvyšší tlakovou ztrátu. Varianta VAR 02 se vyznačuje excentrickým mezikruhovým kanálem. Tato excentricita způsobuje vyšší odpor distanční mřížky, ale také vyšší průměrnou rychlost v ose y na řezu mezikruhovým kanálem. Tato charakteristika je detailněji zachycena na obrázku 2 ukazující rychlostní a teplotní proudové pole v jednotlivých řezech pro jednotlivé geometrické tvary distančních mřížek. Grafy na obrázku 3 ukazují rozložení teploty a rychlosti na středové kružnici v mezikruhovém kanále. Rozložení teplot bylo provedeno na rovině y=0.242 mm. Z grafu je patrné, že varianta VAR 01 a VAR 03 dosahují přibližně stejných maximálních a minimálních teplot. Excentricky umístěná distanční mřížka u varianty VAR 02 ovlivňuje rozložení teploty na řezu mezikruhovým kanálem. Varianta VAR 02 dosahuje vyšší maximální teploty zhruba o 15K a minimální teplota je nižší zhruba o 5K. Toto teplotní rozložení na mezikruhovém kanále se nese celým proudovým polem za distanční mřížkou. Rozložení rychlosti ve směru osy y na obrázku 3 je vytvořeno na stejné středové kružnici jako na předchozím grafu a opět na rovině y=0.242 mm. Tento graf




ukazuje podobné rozložení rychlosti u varianty VAR 01 a VAR 03. U varianty VAR 02 lze pozorovat výrazné nerovnoměrné rozložení rychlostí v mezikruhovém kanále.

Obrázek 3 - Rychlostní a teplotní průběh pro jednotlivé distanční mřížky

Poslední graf na obrázku 4 ukazuje průběh součinitele přestupu tepla podél nerezové trubky. Tento průběh byl získán odečtením teplot v jednotlivých bodech podél nerezové trubky na povrchu stěny a v mezikruhovém kanále. Na základě velikosti definovaného zdroje tepla a odečtených teplot byl vypočten součinitel přestupu tepla. Nejnižší součinitel přestupu tepla v oblasti distanční mřížky je u varianty VAR 03. Součinitel přestupu tepla je menší zhruba o cca 5%. Naopak tato varianta se vyznačuje nejvyšším součinitelem přestupu tepla za distanční mřížkou. Varianta VAR 01 a VAR 02 mají velmi podobný průběh a velikost součinitele přestupu tepla. Nejlepší variantou z hlediska přestupu tepla lze označit variantu VAR 01.

Obrázek 4 - Součinitel přestupu tepla pro jednotlivé distanční mřížky

1

Závěr




Na základě těchto zkušeností lze provádět optimalizace tvaru mřížek, které povedou k lepšímu rozložení teplotního profilu po mezikruhové ploše kanálu a ke zlepšení přestupu tepla mezi palivovou tyčí a chladícím médiem. Nejlepšími výsledky teplotněhydraulické charakteristiky se vyznačuje varianta VAR 01. Tato varianta VAR 01 má nejnižší tlakovou ztrátu než ostatní varianty a tlaková úspora se výrazně zvětší při vynásobení rozdílu tlakové ztráty počtem palivových prutů a počtem distančních mřížek na jeden palivový prout. Tato varianta se také vyznačuje nejvhodnějším průběhem součinitele přestupu tepla v mezikruhovém kanále podél nerezové trubky. Tyto geometrické modely budou dále využity k rozšíření numerických simulací o tzv. var, který zahrnuje vícefázové proudění jako např. pára, vzduch. U těchto numerických simulací se budou vyhodnocovat podobná kritéria jako v tomto článku, ale také další parametry týkající se vícefázového proudění. Takto získané výsledky lze validovat a porovnávat s měřením z experimentálního zařízení pro studium dvoufázového proudění. Literatura LÁVIČKA D., „CFD simulation the thermal-hydraulic characteristic within fuel rod bundle near grid”, „International Conference on Computational Fluid Dynamics 2011”, WASET - World academy of science, engineering and technology, Issue 0079: 2011, Paris, France, 2011, pp. 168-173. ANSYS, Inc., 2011, „ANSYS FLUENT Theory Guide”, Release 13.0, Nov. 2010. INCOPERA, DEWITT, BERGMAN, LAVINE, Fundamentals of Heat and Mass Transfer. 6th edition, October 2006, USA. ISBN 978-0-471-45728-2. KOLEV N.I. Multiphase Flow Dynamics – 4 Nuclear Thermal Hydraulics.1ST edition, Springer. ISBN 978-3-540-92917-8. LÁVIČKA, D., „Popis řešení numerických simulací v mezikruhovém průtočném kanále okolo nerezové trubky”. Stretnutie katedier mechaniky tekutin a termomechaniky, Jasna, Demanovska dolina, Slovensko, červen 2009. ISSN 1335-2938. LÁVIČKA D., „Rychlostní a teplotní proudové pole v mezikruhovém průtočném kanále okolo vyhřívané tyče”. XVII. Aplikácia experimentálnych a numerických metód v mechanike tekutín a v energetike, Bojnice, Slovenská republika, duben 2010. ISBN 978-80-554-0189-8. LÁVIČKA D., „Od jaderného reaktoru k experimentálnímu modelu chlazení palivového proutku”. Modelování a měření v energetice - Tepelné cykly, jaderně energetická zařízení, Hrad Nečtiny, Česká republika, květen 2010. ISBN 978-80-02-02239-8. LÁVIČKA D., „A Comparison of Flow Field Characteristics from PIV Experiment Measurement to Numerical Simulation behind a Spacer in a Vertical Pipe”, Applied and Computational Mechanics, vol. 4, no°1, pp.57-66, Czech republic, July 2010. Poděkování Tato práce vznikla za finančního přispění Grantové agentury ČR v rámci postdoktorského projektu GAČR 101/09/P056. Kontaktní adresa: Západočeská univerzita v Plzni Katedra energetických strojů a zařízení Univerzitní 8 306 14 Plzeň email: [email protected]




PREDICTION OF FLOW THROUGH THE BURNER IN THE INDUSTRIAL-SIZE COMBUSTION CHAMBER JIŘÍ VONDÁL Institute of Process and Environmental Engineering, Faculty of Mechanical Engineering, Brno University of Technology Abstract: It is a well-recognized fact that reliable predictions of turbulent swirling non-premixed flames are very difficult, especially for practical cases where application of LES methods is not feasible. It is also clear that detailed measured data of inlet velocity profiles for swirling combustion air are unavailable in practical applications. Therefore is a need for validation of industry-standard codes for the prediction of flow through swirl generators. In this work are predictions validated by published experimental data for a swirler with guide vanes, which is similar to a typical flame holder in a staged-gas low-NOx burner. Validation of predicted local wall heat fluxes follows with measured data from industrial-scale burner testing facility. Computations are done in ANSYS Fluent v12 code using a range of frequently used momentclosure turbulence models. Keywords: swirl flow, CFD, combustion, heat flux 1

Introduction

Study of flame structure is the subject of long-lasting interest of the combustion modelling community. Detailed in-flame measurements of temperature, velocity and species concentrations have served for the validation of all existing combustion models. Unlike the in-flame properties, wall heat fluxes have been used for model validation only rarely. Heat flux measurements reported in the literature are either spot measurements or global heat transfer rates. Spot measurements however mostly provide just the thermal irradiation flux, not the actual radiative or total heat transfer rate (e.g. in the study of industrial furnaces and boilers (Kobayashi et al. 2002; Ströhle 2004; Hayes et al. 2001)). Likewise, global heat transfer rates calculated from the total hot water (steam) production are insufficient for the validation of detailed predictions. In contrast to that, the interest of engineering community focuses primarily on local heat fluxes and pollutant emissions. Emissions are studied namely to ensure compliance with legislative regulations (e.g. directive (Anon 2001)), while heat fluxes are required to check proper furnace design and to ensure safe operation. It is thus apparent that the correct prediction of local heat fluxes on heat transfer surfaces is one of the most important aspects of practical combustion simulations that should receive adequate attention. Swirl-stabilized non-premixed flames are frequently used in industrial burners, but they present one of the most difficult problems to predict computationally. Only with the advances in large eddy simulations (LES), successful predictions of in-flame properties were reported (Fureby et al. 2007; Sadiki et al. 2006; James et al. 2007). The LES approach is unfortunately still too computationally expensive for the simulation of large-scale fired heaters due to their huge dimensions (on the order of 10 m) and the need to resolve fine features like gas nozzles with diameters on the order of 1 mm. The only viable alternative for practical




predictions in the present as well as for a number of years to come thus consists of models based on first or second-order turbulence closures.




The key question in predicting swirling diffusion flames however is, whether the prediction of swirl using geometry of swirl generator is dependable. In the literature, only scarce instances may be found of measurements suitable for the validation of such swirl generation predictions (Mak and Balabani, 2007; Fernandes et al., 2005). In most cases of advanced predictions of swirling flows including those mentioned above, boundary conditions on the inlet are typically specified using measured velocities and velocity fluctuations. Predictions validated by experimental data are almost nonexistent in peerreviewed journals. Occasionally, swirl is even specified by geometric swirl number, i.e. by inclination of swirl generator vanes (helixes) (So et al., 1985). Neither of these approaches is however suitable for most cases of practical predictions of swirlstabilised gas and liquid fuel burners, due to the large variety of swirl generator designs used by burner vendors and due to the unavailability of detailed measurements. This work therefore includes a simulation of swirl generator validated by an experiment recently published in a peer-reviewed journal. 2

Swirl prediction and validation

For the validation of swirl prediction capability was used a recently published work (Mak and Balabani, 2007). The geometry of experimental setup was further clarified in personal communications with one of the authors (Balabani, 2010). Here we include only a brief summary of the experiment and its results.

Figure 1 - Experimental expansion chamber with guide vane swirl generator

The measurements were performed for a vane swirl generator by optical method (particle image velocimetry, PIV). Geometry of the computational domain including the swirl generator is displayed in Figure 1. Inclination of the guide vanes in the present case is 45°. The experimental work was focused on analyzing flow features in a sudden expansion and its deeper analysis by proper orthogonal decomposition (POD), but they


152



measured also velocity components just above the expansion (x/D = −0.44) in order to determine accurately the amount of swirl in the expanding flow. These velocity measurements above the expansion were used in the present work to validate computational predictions. Working medium was water. Simulations are in this work performed with moment closure turbulence models and a transient formulation (URANS). The tested turbulence models include mainly variants of the popular k-ε model previously validated for different swirl generators e.g. in (Escue & Cui 2010; Fudihara et al. 2007; Widmann et al. 2000). The RNG k-ε model has in those works displayed some superiority over the standard k-ε model since it predicted more accurately recirculation zones in the combustors. Another model applied in this work is the shear stress transport (SST) k-ω model (Menter 1994), which removes the need for wall functions and was previously applied for isothermal flow in a vortex combustor simulation in (Ridluan et al. 2007). Reynolds stress model (RSM) was also included in the analysis, which is more computationally intensive but handles anisotropy of turbulence and has been reported to outperform two-equation eddyviscosity models in highly swirling flows. Finally, the realizable variant of the k-ε model (Shih et al. 1995) was tested as well. All simulations were done with discretization of second order for momentum. All the simulations were performed using a commercial code ANSYS Fluent v12 (ANSYS 2009). Reynolds number at the inlet (before swirl generator) was 10,000 and swirl number calculated from the measured velocity data was 0.64. Two meshes were created for the simulation. One used simplified geometry, taking advantage of the rotational periodicity of the flow domain, which covered 90° section of the whole domain. This simplified geometry was meshed with a high-quality all- hexahedral mesh with 1,138,000 cells. The second model included the whole flow domain and it was meshed by a more coarse mesh containing 621,000 polyhedral cells. Results from the simulation are displayed in Figure 2. The plots show averaged velocity over sufficient time period normalized with the average velocity magnitude. Position of the reporting cross-section lies downstream from the swirl generator but still upstream from the sudden expansion (see Figure 1). The results of simulations display consistent local disagreements with measured data. While simulation suggests no reverse flow, the measurement contradicts it. The key problem is in the central jet prediction (short circuit through the swirler). The most consistent result among two grid sizes gives the RSM model. The biggest grid type sensitivity shows RNG k-ε. Generally the full model predicts lower axial velocity at the r = 0 m than the quarter model.


153



Figure 2 - Comparison of measured and calculated velocity profiles for several turbulence models. Lines with symbols are for full model and lines without symbols for the quarter model with the periodic conditions. 3

Burner testing facility

The construction of the semi-industrial experimental facility for burners up to 2 MW was aimed at providing variable length of combustion chamber and accurate heat flux and emission measurements. The main feature distinguishing the test facility at Brno University of Technology from others is the ability to measure local heat transfer rates to the cooled walls, which is enabled by segmental design of the combustion chamber. There are up to seven water cooled segments of the combustion chamber, see Figure 3. All internal segments have the same flame-facing area of 1.57 m2, whereas the first and seventh section has 1.26 m2 and 3.14 m2, respectively. The last three segments are removable which allows adjustment of the combustion chamber length. The burner was a low-NOx design with staged gas supply (Figure 3) and axial swirl generator (Figure 4), fired by natural gas. Swirl generator with diameter 240 mm consisting of 8 vanes acts as a flame holder. Gas inlet consists of twelve primary nozzles and eight secondary nozzles. Eight of the primary nozzles have diameter 2.6 mm and the other four 3.0 mm. All the primary nozzles are drilled in a nozzle head located on the burner axis. Secondary gas injection is performed by four additional nozzle heads located in regular intervals around annular air channel which surrounds the primary nozzle head. Each of the four secondary nozzle heads has two nozzles with a diameter of 3.3 mm. http://aum.svsfem.cz

154



Flame ignition and stabilization is performed by a small (25 kW) premixed naturaldraft pilot burner. Its thermal duty was included in the total thermal duty.

Figure 3 - Combustion chamber and main parts of the data acquisition system

Figure 5 - Swirl burner with two gas stages

Figure 4 - Swirl generator (flame holder)


155


4


Heat flux measurements

The burner performance was set to 745 kW (Case 1) and 1120 kW (Case 2) with excess air ratio 1.1. Several measurements previously performed at the same testing facility, although with different objectives than in this work, were described in (Bělohradský et al. 2008; Kermes et al. 2008b, Kermes et al. 2008a). Each of the two studied regimes (Case 1 and 2) was measured twice to test repeatability Measured heat fluxes (extracted heat rates in each of the seven sections along the combustion chamber) are shown in Figure 6. It can be seen that reproducibility is good even though conditions were slightly changed e.g. air temperature 19.3 °C vs. 4.3 °C in Case 1. The deviation in the total extracted heat flux into the water was 0.16 % for the Case 1 and 0.88 % in the Case 2. The maximum load difference was in 3rd section 2.7 % for the Case 2 and 2.4 % in 5 th section for the Case 1. The total extracted heat from flame (Qtotal) was evaluated from natural gas flow rate, temperature and pressure. The data displayed in the graphs are averaged values over 5 minutes (Case 1) and 15 minutes (Case 2). Averaging was applied to remove random fluctuations in the measured values. The fuel distribution among primary and secondary nozzles was determined from numerical simulation of gas distributor based on total fuel flow rate since only the gas main is equipped with flow meter. Fuel flow rate in Case 1 was 3.84 10-3 kg/s for primary and 1.1 10-2 kg/s for secondary nozzles. In Case 2 the distribution was 5.79 10-3 kg/s to the primary and 1.65 10-2 kg/s to the secondary nozzles. 5

Heat flux prediction and validation

Grid was created with great care to ensure high quality. Most of the volume of computational domain was meshed by hexahedral cells and only in the vicinity of nozzles were used tetrahedral elements. Total number of grid cells was approximately 1,300,000. Commercial code ANSYS Fluent was utilized for simulation. Turbulence was modeled by second-order turbulence closure that dominates current engineering applications of CFD in the area of reactive flows. These include namely the realizable k−ε model (Shih 1995) and shear stress transport k-ω model (Menter 1994). For the description of chemical reactions and interaction with turbulence was used eddydissipation (ED) model (Magnussen & Hjertager 1977) with a one-step global reaction. Radiation was included using the discrete ordinates method and radiative gas properties were calculated by the weighted-sum-of-grey-gases method with data based on (Coppalle & Vervisch 1983) and (Smith et al. 1982). Emissivity of the walls was set at 0.9.Boundary conditions were utilized according the real measured data. Walls other than water cooled are assumed adiabatic. The solutions were performed using a pressure-based solver option of FLUENT 12.1, with the SIMPLEC solution algorithm and generalized QUICK thirdorder scheme (ANSYS Academic Research, Release 12.1 2009) for momentum and density, and pressure staggering option. The remaining model equations, namely the turbulence model equations, were discretized using first order differencing. Steady solution was sought, but the solution did not converge even in most stabilizing numerical settings with first order upwind differencing and strong relaxation. Therefore a more rigorous treatment was adopted by adding time into the problem formulation. The resulting unsteady-RANS (URANS) simulations display smaller fluctuations and complete convergence is achieved in every time step with step size 0.002 s. In order to obtain results independent of the initial conditions at least 3 s of physical time were always simulated, as the residence time in the combustion chamber is about 2 s and only afterwards it was possible to start collecting data for statistics.


156



Wall heat fluxes predicted by the simulation are shown in the Figure 6 for 745 kW and Figure 7 for 1120 kW. There are deviation in all sections – from 3 % to 16 %, but the last section where the deviation reaches 31 % and 20 % for the realizable k-ε model and SST k-ω model respectively. In the case of 1120 kW the last section has differed significantly from measured data (39 %).

Figure 7 - Heat flux measurements and predictions for the 745 kW (Case 1)

6

Figure 6 - Heat flux measurement and prediction for the 1120 kW (Case 2)

Conclusions

Computed flow profiles of the axial guide vane swirler display consistent local disagreements with measured data. While simulation suggests no reverse flow, the measurement contradicts it. The most consistent result among two grid sizes gives the RSM model. The biggest grid type sensitivity shows RNG k-ε. Generally the full model predicts lower axial velocity at the r = 0 m than the quarter model. The results might be influenced by short-circuit coarse grid at the inlet to the swirler. Consequently it affects the downstream flow predictions. Even though there are uncertainties in swirl flow modeling the complex simulation of heat fluxes was performed. Local wall heat fluxes were investigated in a cylindrical, water-cooled large experimental combustion chamber for non-premixed swirling natural gas flames. Measured data are provided for two firing rates (745 kW and 1120 kW) together with a schema of the fluid flow domain including air duct, staged-gas burner and combustion chamber. The measurements are used to validate heat flux predictions. The results show ability of the utilized models to predict heat fluxes with acceptable accuracy. Unsteady RANS computations are required to obtain well-converged results due to internal flow instabilities. The prediction generally differs from measured local wall heat fluxes in all segments by 7 % but the last water cooled section which differs by up to 39 %.


157



Reference ANSYS Academic Research, Release 12.1, 2009. Help System, Users’s Guide. ANON, 2001. Directive 2001/80/EC of the European Parliament and of the Council on the limitation of emissions of certain pollutants into the air from large combustion plants, Available at: http://eur- lex.europa.eu/LexUriServ/LexUriServ.do?uri=CELEX:32001L0080:EN:HTML. BĚLOHRADSKÝ, P., KERMES, V. & STEHLÍK, P., 2008. Design and analysis of experiment for low-NOx burners design for process industries. V 8th European Conference on Industrial Furnaces and Boilers. 8th European Conference on Industrial Furnaces and Boilers. Vilamouria, Portugal. COPPALLE, A. & VERVISCH, P., 1983. The total emissivities of high-temperature flames. Combustion and Flame, 49(1-3), s.101-108. ESCUE, A. & CUI, J., 2010. Comparison of turbulence models in simulating swirling pipe flows. Applied Mathematical Modelling, 34(10), s.2840-2849. FERNANDES, E.C., HEITOR, M.V. & SHTORK, S.I., 2006. An analysis of unsteady highly turbulent swirling flow in a model vortex combustor. Experiments in Fluids, 40(2), s.177-187. FUDIHARA, T., GOLDSTEIN, L. & MORI, M., 2007. A numerical investigation of the aerodynamics of a furnace with a movable block burner. BRAZILIAN JOURNAL OF CHEMICAL ENGINEERING, 24(2), s.233-248. FUREBY, C. ET AL., 2007. An experimental and computational study of a multi-swirl gas turbine combustor.Proceedings of the Combustion Institute, 31(2), s.3107-3114. HAYES, R.R. ET AL., 2001. Crown incident radiant heat flux measurements in an industrial, regenerative, gas- fired, flat-glass furnace. Experimental Thermal and Fluid Science, 24(1-2), s.35-46. JAMES, S., ZHU, J. & ANAND, M.S., 2007. Large eddy simulations of turbulent flames using the filtered density function model. Proceedings of the Combustion Institute, 31(2), s.1737-1745. KERMES, V., BĚLOHRADSKÝ, P. & STEHLÍK, P., 2008. Influence of burner geometry with staged gas supply on the formation of nitrogen oxides. V 8th European Conference on Industrial Furnaces and Boilers. 8th European Conference on Industrial Furnaces and Boilers. Vilamouria, Portugal. KERMES, V., BĚLOHRADSKÝ, P., ORAL, J., ET AL., 2008. Testing of gas and liquid fuel burners for power and process industries. Energy, 33(10), s.1551-1561. KOBAYASHI, H. ET AL., 2002. Performance of High Temperature Air Combustion Boiler with Low NOx Emission. JSME International Journal Series B Fluids and Thermal Engineering, 45(3), s.481-486. MAGNUSSEN, B.F. & HJERTAGER, B.H., 1977. On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion. Symposium (International) on Combustion, 16(1), s.719-729. MAK, H. & BALABANI, S., 2007. Near field characteristics of swirling flow past a sudden expansion. Chemical Engineering Science, 62(23), s.6726-6746. MENTER, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32, s.1598-1605. RIDLUAN, A., EIAMSA-ARD, S. & PROMVONGE, P., 2007. Numerical simulation of 3D turbulent isothermal flow in a vortex combustor. International Communications in Heat and Mass Transfer, 34(7), s.860-869. SADIKI, A. ET AL., 2006. Unsteady methods (URANS and LES) for simulation of combustion systems. International Journal of Thermal Sciences, 45(8), s.760-773 SHIH, T.-H. ET AL., 1995. A new k-[epsilon] eddy viscosity model for high reynolds number turbulent flows. Computers & Fluids, 24(3), s.227-238. SMITH, T.F., SHEN, Z.F. & FRIEDMAN, J.N., 1982. Evaluation of Coefficients for the Weighted Sum of Gray Gases Model. Journal of Heat Transfer, 104(4), s.602-608.




SO, R.M.C., AHMED, S.A. & MONGIA, H.C., 1985. Jet characteristics in confined swirling flow. Experiments in Fluids, 3(4), s.221-230. STRÖHLE, J., 2004. Spectral Modelling of Radiative Heat Transfer in Industrial Furnaces, Aachen, Germany: Shaker Verlag GmbH, Germany.

WIDMANN, J.F., RAO CHARAGUNDLA, S. & PRESSER, C., 2000. Aerodynamic study of a vane-cascade swirl generator. Chemical Engineering Science, 55(22), s.5311-5320. Acknowledgement The authors gratefully acknowledge financial support of the Ministry of Education, Youth and Sports of the Czech Republic within the framework of project CZ.1.07/2.3.00/20.0020 and research plan MSM 0021630502. Contact Address: Jiří Vondál Institute of Process and Environmental Engineering Faculty of Mechanical Engineering Brno University of Technology Czech Republic [email protected]




VYUŽITÍ 3D NUMERICKÉHO MODELOVÁNÍ PRO STANOVENÍ EFEKTU PŘEKÁŽEK V KANALIZAČNÍM POTRUBÍ KAREL KŘÍŽ*, JAROSLAV POLLERT ML.*, MARTIN KANTOR** * České vysoké učení technické v Praze, Fakulta stavební, Katedra zdravotního a ekologického inženýrství, ** České vysoké učení technické v Praze, Fakulta stavební, Katedra hydrauliky a hydrologie Abstract: The operation of the sewer systems requires high operational costs. Therefore, an improvement and an optimization of the support decision software tool is a very important issue to help spend the budget more effectively. The 3D CFD numerical models are used to determine the effect of failures in the pipes on their hydraulic capacity. In the pipe of gravity sewer network, the flow regime can be changed between pressurized and free surface flow regime. Therefore, those two flow regimes were simulated. This paper describes the modelling procedure and presents the partial results in comparison with the hydraulic models. Keywords: sewer systems, 3D CFD modelling, hydraulic capacity, free surface flow, pressurized flow 1

Úvod

Vysoké náklady na provoz stokových sítí nutí provozovatele rozhodovat se maximálně efektivně, k čemuž často využívají podpůrné softwarové nástroje určující prioritní úseky pro sanaci. Úroveň těchto podpůrných nástrojů je velmi různá. Liší se v detailnosti výstupu a tím pádem samozřejmě i náročností na vstupní data. Jedním z mnoha důležitých aspektů posouzení stavu stokové sítě je její kapacita. Ta je dnes běžně posuzována pomocí 1D srážko-odtokových modelů. Kalibrace srážko-odtokovým modelů klasických (nepoškozených) stokových sítí byly podrobně testovány a popsány v řadě publikací a existují přesné metodiky a doporučení pro výběr a nastavení modelů povrchového odtoku (hydrologický model) i samotného proudění trubním systémem (hydraulický model) (EPA, 2002). Problémem jsou však stokové sítě poškozené. Strukturální poruchy vnáší do modelů řadu nejistot (Kříž, 2007). Přitom mohou vyvolávat poměrně výrazné energetické ztráty proudění a snižovat hydraulickou kapacitu systému. Problematika zahrnující vliv poruch stokových sítí zjištěných při kamerové inspekci do 1D hydrodynamického modelu byla řešena již v evropském projektu (ukončen v roce 2004), jehož výsledkem byl jeden z nástrojů manažeru CARE-S (Computer aided Rehabilitation for water and storm water networks) (Saegrov, 2006). Tento modul nazvaný Degradation Tool přepočítává efekt poruch zjištěných při kamerové inspekci pomocí zvýšené hodnoty Manningova drsnostního součinitele n. Při přepočtu však nejsou uvažovány vzájemné vzdálenosti mezi překážkami, tj. poruchy jsou uvažovány v izolaci i v pozicích, ve kterých se vzájemně ovlivňují. Dalším zjednodušením je použití rovnic zjištěných pro tlakové proudění i pro proudění s volnou hladinou (Pollert, 2004) (NTNU, 2004). S ohledem na popsané nedostatky jsou prováděny práce na zpřesnění tohoto nástroje.



2


Metodika

Vliv poruch stokových sítí na jejich hydraulickou kapacitu je pro získání potřebného množství rovnic vyjadřující vliv velikosti, polohy a tvaru překážky na energetické ztráty pro tlakové proudění i proudění s volnou hladinou podrobně zkoumán pomocí 3D CFD numerického modelování (Fluent verze 12.0.16). Výsledky jsou porovnávány s hodnotami změřenými nebo vypočtenými z modelování hydraulického. Tvorba výpočetních mřížek byla provedena s programu Gambit (verze 2.4.6) s měřítkem 1:1 k oběma hydraulickým modelům. a.

Fyzikální modelování překážek v tlakovém proudění

Pro stavbu fyzikálního modelu bylo využito prostor vodohospodářské haly Fakulty stavební. V rámci projektu byla vybudována trubní linka. Potrubí z průhledného PVC D200 (DN194) v celkové délce 12,18m bylo osazeno na ocelovou konstrukci umožňující změnu sklonu. Z obou stran experimentální části jsou přírubově osazeny redukce 200/150, na které je namontována přítoková a odtoková část D150. Pro zajištění tlakového režimu je odtoková část vedena obloukem nad záklenek potrubí experimentální části linky (viz. Obr. 1). V současné fázi projektu byly simulovány překážky cihel o různých rozměrech (kvádrové plochy) a přípojek zaústěných do průtočného profilu (válcové plochy). Pro měření tlakových ztrát byla linka s odtokovým obloukem vyrovnána do nulového spádu a pro stanovení rozdílu tlaku mezi dvěma profily je používán diferenciální tlakoměr Siemens Sitrans P. Průtok, který je měněn elektrickým uzávěrem v rozmezí hodnot 0 – 30 l/s (odpovídá průřezovým rychlostem 0 – 1,2m/s), je měřen magneticko - indukčním průtokoměrem Krohne. Hodnoty průtoku a rozdílu tlaků byly zaznamenávány pomocí přístroje Micromec MultiSens a softwaru MMGrafix7. Měření probíhalo vždy minimálně 30 sekund a z naměřených hodnot tlakového rozdílu a průtoku byly spočteny průměry. Během všech měření byla měřena a zaznamenávána teplota vody pro stanovení kinematické viskozity. V první fázi experimentů byla stanovena tlaková ztráta potrubí bez překážek (ztráta třením). Tato hodnota je referenční hodnotou experimentálního úseku, který je následně odečítán od tlakových ztrát zjištěných po vložení překážek tak, aby byl získán samotný efekt překážek.

Obrázek 1 – Schéma trubní linky pro tlakové proudění



b.


Fyzikální modelování překážek v proudění o volné hladině

Pro stanovení hydraulických ztrát v potrubí při proudění s volnou hladinou bylo využito fyzikálního modelu trubní linky v hydraulické laboratoři university v Bradfordu (Velká Británie), kde je vyvíjena metoda stanovení hydraulických ztrát akustickou metodou (Romanova, 2010). Použitá trubní linka z plexi potrubí D300 (DN290) má délku 20,1m. Experimenty byly prováděny při sklonech 1% a 0,5%. Nátok do trubní linky z nadzemní nádrže je veden přes uklidňovací objekt, výtok je do volna. Ustálené rovnoměrné proudění bylo zajišťováno pomocí stavítka tak, aby sklon hladiny odpovídal sklonu potrubí. Voda, přepadající přes stavítko, protéká přes měrný objekt průtoku do hlavní podzemní nádrže, odkud je čerpána zpět do nádrže nadzemní. Ve dně potrubí je po délce linky umístěno celkem 30 manometrů, dle kterých byla zaznamenávána hladina a pomocí stavítka na konci potrubí byla vyrovnávána hladina rovnoměrného proudění (sklon hladiny = sklon potrubí) pro každý řešený průtok. Pro stanovení hydraulických ztrát v popisované trubní lince byla zvolena „tracerová metoda“ (Kříž, 2009), (Reichert, 1994), (Rieckermann, ). S ohledem na délku potrubí a řešené průtoky (Q = 3,0; 5,0 a 8,0l/s) byl používaný dávkovaný roztok o koncentraci 30mg/l. Odezva v potrubí byla zaznamenávána v 5ti příčných profilech ve staničení od 7,4 o 17,65m. Vodivost byla zaznamenávána s časovým krokem 0,1s. Hodnoty Manningovy drsnosti vypočtené pro různé kombinace velikostí a vzdáleností překážek byly vztaženy k referenční hodnotě drsnosti pro potrubí bez překážek. Pro kalibraci numerického modelu byl čtením manometrů zaznamenáván průběh hladin. c. Numerické modelování překážek v tlakovém proudění Tvorba výpočetní mřížky V první fázi bylo vytvořeno potrubí odpovídající reálné lince v hydraulické laboratoři (shodná délka, DN, geometrie nátoku a výtoku). Tato síť byla určena pro kalibraci matematického modelu na základě výsledků z hydraulického modelování. Ve druhé fázi byla připravena výpočetní mřížka modelu se stejnou průřezovou plochou potrubí jako u modelu hydraulického, ale byla změněna délka experimentálního úseku modelu. Pomocí této mřížky byly simulovány situace se dvěma různými překážkami v různých vzdálenostech. Délka experimentální části numerického modelu je 50 m. Výpočetní mřížka byla v obou případech generována pomocí „velikostní funkce” (Create Size Function), kde zdrojem byly plochy vytvořené překážky (překážek), které byly vztaženy k objemu celého potrubí. Velikost mřížky na plochách překážek byla 0,002m a zvětšovala se poměrem (growth rate) 1.1na maximální velikost buňky 0,02 m. V experimentální části potrubí byl použit trojúhelníkový typ mřížky (Tet/Hybrid), částech nátoku a odtoku z modelu byly použita čtvercová síť ve velikost (Hex/Wedge) ve velikosti výpočetní buňky 0,03m (ANSYS, 2009). Nastavení modelu Výpočetní mřížka byla importována do softwaru Fluent, kde byl výpočet řešen jako ustálené proudění (typ solveru „pressure-based”). Na základě testování turbulentních modelů byl zvolen model k-omega Standard, s využitím funkce Shear Flow Correction. Okrajovými podmínkami byl vtok řešený jako velocity-inlet. Simulace byly prováděny pro vtokovou rychlost v rozmezí 0,3 – 2,0 m/s. druhou okrajovou podmínkou byl výtok do volna řešený jako pressure-outlet. Parametry stěny (stationary wall) byly Roughness height = 0m a




drsnostní konstanta Roughness Constant = 0,5. Byly použity defaultní referenční hodnoty (ANSYS, 2009). Konvergenční kriteria pro všechny hodnoty residuí (kontinuita, složky rychlostí x, y, z, a turbulentních parametrů k a omega) byla nastavena hodnotou 1.10-4. Vyhodnocení ztrát Po dokončení simulace byly vyhodnocovány rozdíly statického tlaku v příčných řezech odpovídajících umístění tlakoměru na fyzikálním modelu. Pro všechny simulované rychlosti byl nejprve vyhodnocen rozdíl tlaků prázdného potrubí (odpovídá ztrátě třením) a následně byly vyhodnoceny modely s překážkami. Po odečtení ztráty třením byla získána tlaková ztráta vyvolaná simulovanými překážkami. d. Numerické modelování překážek v proudění o volné hladině Výpočetní síť pro numerický model pro proudění s volnou hladinou odpovídá geometrii reálné trubní linky v hydraulické laboratoři (University of Bradford). Tvorba výpočetní mřížky Obdobně jako u výpočetní mřížky pro tlakové proudění bylo pro sestrojení modelu využito definování velikosti mřížky pomocí velikostní funkce (Size Function). Zdrojem byly plochy objektu stavítka (start size 0,0035 m, growth rate 1,1 a max. size 0,04m). Pro výpočetní mřížku potrubí s překážkami byly navíc definovány další velikostní funkce pro každou z překážek. Parametry této funkce byly totožné s funkcí pro hradítko (start size 0,0035 m, growth rate 1,1 a max. size 0,04m). S ohledem na použití velikostních funkcí byl opět použit trojúhelníkový typ mřížky (Tet/Hybrid) (ANSYS, 2009). Nastavení modelu Model byl pro všechny sklony geometricky řešen s nulových spádem. Sklon byl do výpočtu zahrnut pomocí složek gravitačního zrychlení. Stejně jako pro simulace tlakového proudění byl použit turbulentní model k- omega standard s „Shear flow correction”. Do dvoufázového proudění byl nadefinován materiál vzduchu a vody. Okrajovými podmínkami byl „mass flow inlet” definovaný jako „open channel”. Samotný průtok byl pro obě fáze definován hodnotami v kilogramech za sekundu. Druhou okrajovou podmínkou byl opět výtok do volna „pressure outlet”. Stejně jako v modelech tlakového proudění byly uvažovány defaultní hodnoty stacionární stěny. Konvergenční kriteria pro všechny hodnoty residuí (kontinuita, složky rychlostí x, y, z, turbulentních parametrů k a omega a také vf phase vody) byla nastavena hodnotou 1.10 -3. Pro nadefinování hladiny byla pro Region Adaption nastavena hodnota hladiny na vtoku získaná s fyzikálního modelu (ANSYS, 2009). 3

Výsledky

Vzhledem ke stále probíhajícím pracím na projektu jsou v článku uvedeny pouze dílčí výsledky. V této kapitole budou uvedeny výsledky pro tlakové proudění potrubím.




Obrázek 2 - Kalibrační graf 3D matematického modelu tlakového proudění

Z obrázku 2 je patrný průběh tlakové ztráty čistého potrubí (ztráty třením) v závislosti na rychlosti proudění, kde první světlá křivka znázorňuje tlakovou ztrátu získanou matematickým modelováním, tmavá křivka byla získána z modelu fyzikálního a čárkovaná tmavá křivka představuje kontrolní teoretickou vypočtenou hodnotu dle Altšula pro hydraulicky hladké potrubí. Na obrázku 3 jsou znázorněny hodnoty součinitele místní ztráty k v závislosti na velikosti překážky zjištěné z matematického modelování. Osa x znázorňuje poměr ploch příčného řezu překážky Sb a průřezové plochy potrubí Sp (viz. Obr. 2 a Obr. 3A).

Obrázek 3 – Součinitel místní ztráty pro různé velikosti překážek (platí pro Re 50000 až 130000)




Energetická ztráta místní vyvolaná překážkou definovanou poměrem Sb/Sp je pak vypočtena dle základní hydraulické rovnice:

(1) Uvažování poměru Sb/Sp je možné pouze u překážek, které svým tvarem v příčném profilu kruhového potrubí tvoří kruhovou výseč (viz. Obr 4A). Jedná se nejčastěji o sediment, či materiál napadaný do stokové sítě z povrchu (kód poruchy BBC) (ČSN EN 13508-2). V případě, kdy je v příčném profilu potrubí překážka typu B, C, nebo D (Obr. 4 B, C, a D), je nutné řešit nejen, poměr Sb/Sp, ale také poměr h (výšky překážky) a w (šířky překážky). Tento typ překážky v reálu odpovídá např. zaústěné kanalizační přípojce do průtočného profilu (kód poruchy BAG), kusy trubního materiálu (kód poruchy BBE) atd. (ČSN EN 13508-2).

Obrázek 4 – Typy příčných profilů překážek o stejné průřezové ploše Sb

Na obrázku 5 A je znázorněn vývoj součinitele tlakové ztráty k, pokud uvažujeme pouze poměr Sb/Sp, na obrázku 5 B pak závislost k na tvarovém součiniteli ps (zjištěno ze CFD simulací)




Obrázek 5 – Závislost součinitele místní ztráty k na tvaru a velikosti překážky

Součinitel ps je počítán dle následující experimentálně zjištěné rovnice:

(2)

Obrázek 6 – Vývoj tlakové ztráty v závislosti na Re pro různé vzdálenosti mezi překážkami

Parametry tvaru překážky použité v rovnici (2) jsou znázorněny na obrázku 4. Platnost rovnice byla ověřována pro rozmezí Re 50000 – 200000.




Z obrázku 6 je patrný vývoj tlakových ztrát v potrubí s tlakovým prouděním pro 2 překážky zjištěný pomocí matematických modelů (typ překážek viz obr. 4 A), jejichž výšky odpovídají 20% vnitřního průměru potrubí. Aby bylo možné obecné využití výsledků, jsou grafy zpracovány s bezrozměrnými osovými hodnotami λ (bezrozměrný součinitel tření) a Reynoldsovo číslo. Stejně tak i vzdálenosti mezi překážkami jsou obecně řešeny jako n*D, kde n je násobek a D vnitřní průměr potrubí. 4

Závěr

V článku byla popsána metodika stanovování vlivu poruch potrubí na jeho hydraulickou funkci pro účely zpřesnění nástroje Degradation Tool, který přepočítává vliv nalezených poruch při kamerové inspekci na parametr zahrnující zvýšení hydraulických ztrát do 1D hydrodynamického modelu. Na základě matematických simulací ověřených na fyzikálním modelu byly nalezeny závislosti velikostí překážek pro rozmezí Re od 2.104 do 2,1.105 na vývoj tlakových ztrát. Dle uvedených dílčích výsledků pro tlakové proudění je zřejmé, že při stanovování velikosti energetických ztrát pro určité typy poruch je nepřesné uvažovat pouze poměr velikostí ploch příčného profilu potrubí a překážky.Je vhodné uvažovat vliv tvaru překážky (poměr h a w). Na základě výsledků simulací jsou získány také rovnice popisující vliv vzájemné polohy mezi překážkami. Veškeré poruchy je možné více či méně přesně matematicky popsat a na základě záznamů inspekce potrubí pak uvažovat energetickou ztrátu v reálném potrubí. Toto přiblížení realitě opět sníží tolik nežádoucí nejistoty pro modelování pomocí srážkoodtokových modelů porušených stokových sítí. Zatím nedořešenou otázkou zůstává minimální nutná přesnost a detailnost hydraulického popisu poruch ve vztahu k maximální možné schematičnosti pro využití popisované metodiky v běžné praxi. Tato analýza bude součástí závěrů popisovaného projektu. Literatura ANSYS. (n.d.). ANSYS FLUENT 12: User's Guide. ANSYS. (n.d.). GAMBIT 2.4 User's Guide . ČSN EN 13508-2 Posuzování stavu venkovnich sítí a kanalizačních přípojek – Část 2: Kódovací systém pro vizuální kontrolu. Praha: Český normalizační institut, 2004. 116 s. KŘÍŽ, K., BAREŠ, V., POLLERT, J., STRÁNKÝ, D., Vliv poddolování na hydraulickou funkci stokové sítě. Wastewater 2007, Brno, ISBN 978-80-239-9618-0. Page 97-104. 2007 KŘÍŽ, K., BAREŠ, V., POLLERT, J., STRÁNKÝ, D., Risk analysis of sewer system operational failures caused by unstable subsoil, Risk management of water supply and sanitation systems. Dodrecht: Springer, 2009, p.25-75. ISBN 978-90-481-2364-3 KŘÍŽ, K., BAREŠ, V., POLLERT, J., Tracer method for experimental determination of the hydraulic roughness. In: 1 st Eastern European Regional Young Water Professionals Conference. Oxford: IWA, 2009, vol. 1, p. 374-380. ISBN 978-985-525-145-4. NTNU, SINTEF. Manual for Degradation Tool. SINTEF. 2004 POLLERT, J., UGARELLI, R., SAEGROV, S., SCHILLING, W., DI ZEDERICO, 2004. Change of the Hydraulic Capacity of Sewer Systems due to Deterioration, UDM – Dresden, PP 719 – 727 REICHERT, P., 1994. Aquasim—A tool for simulation and data analysis of aquatic systems. Water Sci. and Technol., 30(2), 21–30.




RIECKERMANN, J., NEUMANN, M., ORT, C., HUISMAN, J.L., GUJER W., 2005. Dispersion coefficients of sewers from tracer experiments Water Sci. and Technol., 52(5), 123 - 133. ROMANOVA, A., 2010 Hydraulic head loss monitoring and determination of roughness by acoustic instrumentation in partially full sewer pipes, 6th International conference on sewer processes and networks, Australia SAEGROV, S. CARE-S Computer aided rehabilitation for water and storm water networks. London: IWA Publishing, 2006, Počet stran 140. ISBN 1843391155. Wastewater planning users group. Code of practice for hydraulic modelling of sewer systems, Wastewater planning users group. 2002 Poděkování Tento článek byl vypracován za podpory výzkumných projektů MSM 6840770002, SGS10/147/OHK1/2T/11 a projektu GAČR 103/08/0134. Kontaktní adresa: Thákurova 7, 166 29 Praha 6, tel.: +420 224 354 344, email:[email protected]




MODELOVÁNÍ PŘEPADU VODY PŘES POHYBLIVOU KLAPKOVOU KONSTRUKCI V. JIRSÁK, M. KANTOR, P. SKLENÁŘ České vysoké učení v Praze, Fakulta stavební Abstract: The project deals with the issue of detailed CFD (computational fluid dynamics) analysis of flow over movable fish-belly flap gate. The concept of project is based on execution of set of physical and correspondent numerical experiments, where the fundamental measured and evaluated variables are compared with the previous hydraulic research. Determination of hydraulic load acting on water-side surface of the movable gate using numerical simulations is fully justified for flap gates with atypical shape or for structures with complicated weir bottom. Keywords: CFD, Multi-phase ﬂows, physical modelling, flap gate. 1

Úvod

Klapkové uzávěry jsou dnes často používanými typy pohyblivých hradících konstrukcí na jezech (příklad vyobrazen v Obr. 1). Svou konstrukcí jsou tvořeny dvěmi zaoblenými plechy na návodní a vzdušné straně a příčně vyztuženy diafragmami přes které je klapka pomocí ložiska připojena ke spodní stavbě. Klapky udržují vzdutou hladinu na požadované úrovni v jezové zdrži vhodným sklopením. Pohyb a ovládání bývá zajištěno pomocí závěsných nebo podpíraných hydromotorů.

Obrázek 1 - Příklad klapkového uzávěru na jezu v Berouně

Pro vlastní dimenzování hradících plechu a diafragmy je velice důležité znát průběh a rozložení hydrodynamických sil po vlastní konstrukci. Pro návrh ložisek a hydromotoru je důležité znát celkové sklopné momenty, výslednice sil a jejich ramena. 2

Koncept řešení

Základní koncepce spočívala ve fyzikálním a numerickém modelování proudění vody o volné hladině přes klapku. Model duté klapky pro fyzikální a numerický experiment vycházel z klapky o parametrech: hrazená výška H = 0,192 m, sklon tečny k zaoblené hradící stěně ω = 68,1° a poloměr návodní válcové plochy R = 2,25·H (schématicky zobrazeno v Obr. 2).




Obrázek 2 - Schématický řez klapkovým uzávěrem

Obrázek 3 - Pohled na experimentální zařízení

a.

Fyzikální modelování

Fyzikální model této klapky byl umístěn ve hydraulickém žlabu vodohospodářských laboratoří FSv (vyobrazeno na Obr.3). Na fyzikálním modelu byl měřen přepadový průtok a tvar přepadového paprsku pro jednotlivé polohy sklopení klapky (vyjádřeno poměrem přepadové výšky ku hrazené výšce - h/H). b.

Matematické modelování

Matematické modelování proudění přes geometricky shodnou klapku bylo provedeno CFD programem ANSYS CFX za využití modelu pro vícefázové proudění vody a vzduchu. Geometrie je zjednodušena na 2D výsekový model o tloušťce 10 mm. CFD analýza je provedena na 2D výpočetní síti tvořené hexaedry (detail sítě na Obr. 4): 

Mesh 1 = hrubá síť, velikost elementů od 10 mm po 5 mm, počet všech elementů je 13 tisíc;



Mesh 2 = jemná síť, velikost elementů od 2,5 mm po 1 mm, počet všech elementů je 362 tisíc.

Okrajové a počáteční podmínky numerického výpočtu vycházely z obdobných předpokladů jako u fyzikálního experimentu (schématický pohled na Obr. 4): 

Vstup vody do modelu definován rovnoměrným hydrostatickým rozložením tlaku po výšce;



Výstup vody z modelu definován tlakovou okrajovou podmínkou – výtok do volna;



Strop modelu definován otevřenou tlakovou podmínkou s 0 Pa přetlak/podtlak;


rychlostním

polem

s





Zavzdušnění prostoru pod klapkou definovanou otevřenou tlakovou podmínkou s 0 Pa přetlak/podtlak;



Dno a povrch klapky definován jako stěna;



Boční stěny výsekového modelu definovány jako podmínka symetrie.

V rámci vícefázového proudění vody a vzduchu byly testovány různé modely Interphase Transfer s defaultním nastavením: 

Homogeneous Model



Free Surface Model



Mixture Model

 Základní nastavení výpočetního modelu: 

Řešič – CFX;



Model turbulence SST k-ω;



Residua – RMS 1e-05;



Stacionární řešení s Pfysical Timescale – dle stability výpočtu 0,05s – 0,25s;



Řád přesnosti – High Resolution;



Celkový počet iterací dle ustálení průtokových podmínek a momentového zatížení klapky (imbalance do 1% z průměrné hodnoty);

Obrázek 4 – Schématické vyobrazení parametrů CFD modelu




Vyhodnocením numerického experiment byl zjištěn přepadový průtok, tvar přepadového paprsku, rozložení hydrodynamického tlaku po návodní straně, sklopný moment a hydrodynamická síla působící na klapku pro jednotlivé polohy sklopení klapky (h/H). 3

Vyhodnocení experimentu

Vyhodnocení projektu spočívalo v základním porovnání dat získaných z fyzikálního a numerického experimentu s daty z předchozího hydraulického výzkumu (Čihák, 2001), (Strobl, 2005). a. Porovnání tvaru přepadového paprsku Nejlepší shody mezi fyzikálním experimentem a CFD modelováním ve tvaru přepadového paprsku (vyobrazeno v Obr. 5) je dosaženo pro vícefázové modely Free Surface Model a Mixture Model na jemné síti – Mesh 2. Model vícefázového proudění Homogeneous Model vycházející ze společného rychlostního pole pro obě fáze napočítává velké strhávání vzduchu paprskem do prostoru pod klapkou, a proto dochází k mírnému oddálení paprsku od klapky v místě dopadu paprsku na dno.

Obrázek 5 – Porovnání tvaru přepadového paprsku pro h/H = 0,3




Porovnání průběhu hladiny mezi fyzikálním experimentem a CFD simulací pro různé úrovně sklopení klapky je vyobrazen v Obr. 6, které demonstruje dobrou komparaci Free Surface Modelu s fyzikálním experimentem pro různé úrovně sklopení klapky.




Obrázek 6 – Tvar přepadového paprsku pro různé polohy sklopení klapky


174



b. Vyhodnocení součinitele přepadu Při přepadu vody přes klapkovou pohyblivou konstrukci jezu se mění přepadová výška a tvar přelivné plochy, a proto se mění i součinitel přepadu – µv. Součinitel přepadu – µv závisí na poloze sklopení a na geometrickém tvaru vlastní klapky a má přibližně hodnotu µp = 0,55 – 0,75 (Čihák, 2001). Na hodnotu přepadového součinitele má vliv úroveň vody v podjezí (vliv zatopení přepadu) a způsob zavzdušnění prostoru pod klapkou. Při experimentálním stanovení součinitele přepadu µp pro různé polohy sklopení klapky se vychází z konzumpční křivky přepadu vody (Obr. 7) a rovnice přepadu:

(1)

Obrázek 7 – Konzumpční křivka přepadu vody přes pohyblivou klapku pro různé polohy sklopení

V Obr. 8 jsou porovnány průběhy součinitele přepadu experimentálně zjištěné z fyzikálního a CFD modelování s daty z literatury (Čihák, 2001), (Strobl, 2005).




Obrázek 8 – Průběh součinitele přepadu v závislosti na poloze sklopení klapky




c. Vyhodnocení zatížení pohyblivé klapky hydrodynamickým tlakem Údaje o hydrodynamickém zatížení klapky zjištěné CFD simulací jsou porovnány s údaji z literatury (Čihák, 2001) v Obr. 9 a je zde patrná dobrá shoda. Největší zatížení hradící konstrukce klapky nastává ve střední poloze klapky při jejím sklopení přibližně o 0,55 H od hladiny horní vody. V této poloze se mění zatížení od vodního tlaku trojúhelníkového tvaru na zatížení křivočarého lichoběžníku, který se již blíží obdélníkovému tvaru.

Obrázek 9 – Průběh hydrodynamického zatížení po návodní hradící stěně klapky

Krouticí momenty potřebné pro dimenzování průřezu klapky a pohybovacích mechanismů se stanoví z hodnot sklopných momentů (Čihák, 2001) nebo mohou být výsledkem CFD simulace. 4

Závěr

Práce se zabývala fyzikálním a numerickým modelováním proudění vody o volné hladině přes zmenšený model pohyblivého klapkového uzávěru. Vyhodnocení experimentu spočívalo v základním porovnání dat získaných z fyzikálního a numerického modelování s daty z předchozího hydraulického výzkumu, kde byla prokázána velice dobrá shoda jak v nasimulování tvaru přepadového paprsku, velikosti přepadového průtoku a proudem generovaného hydrodynamického zatížení na těleso klapky. CFD modelování je v tomto směru velice silným nástrojem, který umožňuje zjištění hydrodynamického zatížení geometricky atypických jezových uzávěrů a to i při extrémních




průtokových stavech. Výsledky a závěry z této práce na zmenšené modelové klapce jsou využity pro FSI analýzu skutečných klapkových uzávěrů.




Literatura Čihák F.,Medřický V.: Hydrotechnické stavby 20 – navrhování jezů, pp. 96-104, Praha ČVUT, 2001. Strobl T.: Skriptum – Grundkurs “Wasserbau und Wasserwirtschaft“, pp. 127-129, TU München, 2005. Poděkování Tato práce byla podporována aktivitami výzkumného centra CIDEAS v rámci projektu 1M0579 MŠMT ČR a projektu SGS10/034/OHK1/1T/11 “Numerická simulace přepadových a výtokových paprsků u pohyblivých hradících konstrukcí“. Kontaktní adresa Fakulta stavební, Thákurova 7, 166 29 Praha 6




FSI ANALÝZA BRZDOVÉHO KOTOUČE TRAMVAJE MICHAL MOŠTĚK TechSoft Engineering, s.r.o. Abstrakt: Tento příspěvek vznikl ze vzorového příkladu pro tepelný výpočet brzdových kotoučů tramvaje, jehož vypracování bylo jednou z podmínek výběrového řízení na dodávku CFD software pro firmu DAKO-CZ, a.s. Třemošnice. Řeší aktuální potřeby ve společnosti pro návrh brzdového kotouče tramvaje ve dvou režimech, kde je počítáno oteplení kotouče tramvaje, u které je udržována 1) konstantní rychlost na svahu mechanickými brzdami a při příjezdu do cílové stanice je provedeno zabrzdění maximální brzdnou silou a 2) opakované brzdění na rovině. Na závěr CFD výpočtu jsou výsledky importovány do FEA řešiče pro pevnostní analýzu. Klíčová slova: brzdění, FSI analýza, oteplení 1

Úvod

Výpočet obsahuje nestacionární teplotní analýzu kotouče tramvaje v různých režimech. Do výpočtu je zahrnut vliv rychlosti nabíhajícího vzduchu, otáčení brzdového kotouče a jeho součástí při jízdě, pohyb okolí. Výsledky výpočtu zobrazují teplotní pole na povrchu kola a rychlostní a teplotní pole v řezu okolo proudícího vzduchu. Analýza byla provedena pomocí CFD programu ANSYS FLUENT v13.0. Následně byla spočtena FEA analýza v programu ANSYS Mechanical. Podle zadání byly provedeny dva CFD výpočty proudění pro Úlohu 1 a Úlohu 2. 2

Zadání úlohy č. 1 – udržování konstantní rychlosti na svahu – určení výkonů pro tepelný výpočet kotouče

Cílem řešení byl výpočet oteplení kotouče tramvaje, u které byla udržována konstantní rychlost na svahu (11,1 m/s) mechanickými brzdami, a při příjezdu do cílové stanice bylo provedeno zabrzdění maximální brzdnou silou. Brzdění bylo prováděno normálnou silou působící na jednu brzdovou desku N = 18 499 N a dobrzdění bylo prováděno maximální normálnou silou působící na jednu brzdovou desku N = 40 000 N. Doba jízdy tramvaje ze stanice V bokách do stanice Hlubočepy byla vypočtena na 128,6 s, dobrzdění pak na 9,07 s. Cílem řešení byl výpočet oteplení kotouče tramvaje, u které byla udržována konstantní rychlost na svahu (11,1 m/s) mechanickými brzdami, a při příjezdu do cílové stanice bylo provedeno zabrzdění maximální brzdnou silou. Brzdění bylo prováděno normálnou silou působící na jednu brzdovou desku N = 18 499 N a dobrzdění bylo prováděno maximální normálnou silou působící na jednu brzdovou desku N = 40 000 N. Doba jízdy tramvaje ze stanice V bokách do stanice Hlubočepy byla vypočtena na 128,6 s, dobrzdění pak na 9,07 s.




Tepelný tok jednou brzdící deskou byl vypočten podle vzorce:

(1)

Při dobrzdění, které bylo uvažováno po projetí tohoto úseku od času 128,6 s, byl tepelný tok jednou brzdící deskou vypočten v každém okamžiku podle vzorce:

(2)

Rychlost okolních stěn byla počítána v každém časovém okamžiku podle vzorce:

(3)

3

Zadání úlohy č. 2 – opakované brzdění na rovině – určení výkonů pro tepelný výpočet kotouče

Oproti předchozí úloze, kde byla téměř po celou dobu konstantní rychlost soupravy, je u této úlohy řešen tento brzdný režim, kde byly jednotlivé veličiny počítány podle následujících vzorců:




(4)

(5)

(6) (7)

(8) (9)

(10) (11) (12) (13)

(14) (15)




(16) (17)

(18)

(19) (20) (21)

4

Geometrie

Zadavatelem byla předána data 3D geometrie kotouče a jeho uložení na nápravě ve formátu STEP, viz Obr. 1 a Obr. 2. Pro simulaci byly použity tyto součásti: brzdový kotouč, dvě brzdové desky, náboj, distanční podložka a část nápravy. Kromě těchto součástí bylo vytvořeno okolí respektující uložení brzdného kotouče na nápravě.

Obrázek 1 - Geometrie kotouče a jeho uložení na nápravě – dodaný 3D model




Obrázek 2 - Detail geometrie kotouče a jeho uložení na nápravě – dodaný 3D model

Schéma modelované oblasti je uvedeno na Obr. 3. Geometrie byla importována z předaného 3D CAD modelu ve formátu STEP, upravena dle zvyklostí pro tento typ výpočtu, a kolem byla vytvořena oblast reprezentující okolí kola. Oblast okolo kola má tvar kvádru s čelní vstupní plochou a zadní výstupní plochou. Kolo je umístěno na výšku a hloubku uprostřed modelované oblasti, na šířku pak podle umístění na nápravě, kde okraje jsou 0,5 m od středu brzdného kotouče.

Obrázek 3 - Simulovaná výpočetní oblast




Pro vstup vzduchu je použita okrajová podmínka mass-flow inlet. Pro výstup vzduchu je použita podmínka pressure-outlet. Pro zbylé plochy "obalující" oblast okolo kola je použita podmínka stěny z důvodu použití sdílení tepla radiací a pohybu země pod vozem. 5

CFD model

Na základě upraveného geometrického modelu byl vytvořen CFD model pevných součástí vozidla a okolí. Výpočetní síť modelu obsahuje 6.000.000 buněk. Výpočetní oblast je tvořena převážně šestistěny. Síť byla vhodně zahuštěna v kritických oblastech, kde se předpokládají vysoké gradienty teploty a rychlostí. Součásti navazující na brzdové desky nebyly modelovány, což ovlivní nejenom proudění v této oblasti ale především přenos tepla těmito částmi. Ukázka použitá povrchové výpočetní sítě na brzdovém kotouči je na obrázku 4.

Obrázek 4 - Ukázka povrchové výpočetní sítě na brzdovém kotouči

a.

Okrajové podmínky

Proudění bylo modelováno jako nestacionární, turbulentní se standardním k-eps modelem, s konvektivním, konduktivním sdílení tepla a sáláním pomocí modelu P1. Rotace kola je řešena pomocí Multiple Reference Frame (MRF) modelu. Podstata MRF modelu spočívá v rozdělení celé výpočetní oblasti na domény, z nichž může být kterákoliv rotační nebo posuvná ve vztahu k další navazující doméně. Na rotující oblasti byl použit MRF model se zadanou rychlostí rotace, která byla vypočtena z rychlosti tramvaje. Ta byla konstantní pro Úlohu č. 1 (s výjimkou dobrzdění ve stanici) a časově závislá podle režimu tramvaje a) až e) pro Úlohu č. 2. Stejně tak stěny ohraničující výpočetní oblast byly modelovány jako pohyblivé v závislosti na rychlosti tramvaje. Hmotnostní průtok na vstupu do oblasti byl modelován jako vzduch při teplotě 20 °C a byl upravován v závislosti na rychlosti tramvaje.




Z vypočtené normálové síly na jednu brzdovou destičku byl v každém okamžiku vypočten tepelný tok generovaný brzděním, který byl dosazován na stěnu mezi brzdovou destičkou a brzdovým kotoučem jako okrajová podmínka. Ten byl opět závislý na rychlosti tramvaje. Závislost tepelného toku a rychlosti na čase pro Úlohu č. 2 pro režim od a) do e) je znázorněn na obrázku 5 a 6.

Obrázek 5 - Závislost tepelného toku na čase pro Úlohu č. 2

Obrázek 6 - Závislost rychlosti soupravy na čase pro Úlohu č. 2

6

Výsledky úlohy č. 1

Na následujících Obr. č. 7 - 15 je zobrazeno rozložení teploty na plochách mezi vzduchem a stěnami zařízení s časovým krokem 5 s až do okamžiku 135 s, včetně zobrazení v čase 128,6 s (na konci konstantního brzdění) a 137,68 s (po samotném dobrzdění tramvaje). Na levé části obrázků jsou zobrazeny mapy teplot v intervalu od minimální do maximální hodnoty teploty na povrchu brzdy, na pravé části obrázků pak hodnoty od minimální hodnoty teploty na povrchu brzdy po maximální hodnotu teploty brzdy v čase 128,6 s.




Obrázek 7 - Rozložení teplot [K], čas = 0 s













Obrázek 14 - Rozložení teplot [K], čas = 128,6 s (konec konstantního brzdění)

Obrázek 15 - Rozložení teplot [K], čas = 137,68 s (po dobrzdění)

V průběhu simulace byla ukládána průměrná hodnota teploty na vnější ploše levého a pravého brzdného kotouče. Na Grafu 1. je patrný průběh teploty až do úplného zastavení tramvaje.

Graf 1 - Úloha 1: Teplotní profil na stěnách brzdového kotouče - průměrná hodnota

Na Obr. 16 je zobrazena mapa rychlostního pole v řezu procházejícím středem brzdového kotouče pro čas 128,6 s, tedy na konci brzdění s konstantní silou. Na Obr. 17 jsou pak zobrazeny vektory rychlosti na stejném řezu ve stejném čase. Na Obr. 18 je doplněno odpovídající teplotní pole na stejné ploše.




Obrázek 16 - Mapa rychlostního pole v čase 128,6 s

Obrázek 17 - Vektory rychlosti v čase 128,6 s




Obrázek 18 - Mapa teplotního pole v čase 128,6 s

7

Výsledky úlohy č. 2

Na Obr. 19 je znázorněn průběh rychlosti soupravy na čase a vyznačeny body časových úseků, kde jsou na následujících Obr. 20 - 26 zobrazeny mapy teploty mezi vzduchem a stěnami zařízení od 0 s, s časovým krokem 20 s do 100 s, resp. do ukončení cyklu stání e) v čase 113,05 s. Na levé části obrázků jsou zobrazeny mapy teplot v intervalu od minimální do maximální hodnoty teploty na povrchu brzdy, na pravé části obrázků pak hodnoty od minimální hodnoty teploty na povrchu brzdy po maximální hodnotu teploty brzdy v čase 10 s.

Obrázek 19 - Rychlost soupravy s vyznačenými časy pro zobrazené obrázky 20 – 26













Obrázek 26 - Rozložení teplot [K], čas = 113,05 s

Stejně jako u Úlohy 1 byla v průběhu simulace ukládána průměrná hodnota teploty na vnější ploše levého a pravého brzdného kotouče. Na Grafu 2 je patrný průběh teploty až do času 113,05 s.




Graf 2 - Úloha 2: Teplotní profil na stěnách brzdového kotouče - průměrná hodnota

Pro CFD výpočet simulace periodického opakování brzdného cyklu b) až e) do ustálení teploty kotouče je vhodné užití stacionární simulace s nastavením průměrných hodnot veličin vypočtených z režimu b) až e). Jedná se o velice efektivní využití CFD programu, které v porovnání s nestacionární úlohou trvá jen zlomek výpočetního času. Pro tento ustálený stav byly vypočteny tyto hodnoty: 

tepelný tok jednou deskou, Qavg = 2 471 W



průměrná rychlost tramvaje, vavg = 5,35 m/s



rotační rychlost brzdy, ωavg = 17,55 rad/s

Obrázek 27 - Rozložení teplot [K], průměr

Pro porovnání hodnot vypočtených stacionární simulací se simulací nestacionární je na následujícím obrázku (Obrázek 28) zobrazeno rozložení teplot na jedné ploše kotouče. Vlevo je obrázek stacionárního ustáleného stavu po dosažení ustálení teploty na kotouči, vpravo pak rozložení teploty v čase 80 s nestacionární simulace.




Obrázek 28 - Rozložení teplot na 1 ploše brzdového kotouče [K], porovnání stacionární a nestacionární simulace (80 s)

8

MKP výpočet teplotních deformací

Deformace od teplotních polí spočítaných CFD analýzou v programu ANSYS FLUENT jsou řešeny v programu ANSYS Mechanical. Přenos zatížení z CFD výpočtu do MKP modelu je proveden automaticky prostřednictvím prostředí ANSYS Workbench. během přenosu dochází k interpolaci dat z CFD sítě na MKP síť pro každou vybranou součást. V ukázce je použit závěrečný čas brzdění po svahu (t = 128,6 s). Výpočet je proveden na sestavě kotouče s částí nápravy. Ve výpočtu není řešena kontaktní úloha, všechny díly jsou na stykových plochách pevně spojeny.

Obrázek 29 - Schéma výpočtu v prostředí ANSYS Workbench




Obrázek 30 - Okrajové podmínky a zatížení, interpolované teploty na pravou část kotouče

Náprava je na jednom konci uchycena ve směru osy otáčení a v obvodovém směru, tedy může se volně roztahovat ve směru radiálním. Zatížení je od teplotního pole a od rotace s rychlostí 36,39 rad/s.

Obrázek 31 - Okrajové podmínky a zatížení

Obrázek 32 - MKP model

Kotouč se vlivem rozdílných teplot na pravé a levé straně „naklápí“. To je spolu s vysokým gradientem teplot v uchycení kotouče do náboje zdrojem vysoké napjatosti v uchycení kotouče.




Obrázek 33 - Celková posunutí v čase 128,6 s

Obrázek 34 - Redukované napětí HMH v čase 128,6 s, řez modelem

Je potřeba zdůraznit, že je nutná pečlivá verifikace CFD výsledků, odladění okrajových podmínek a analýza vlivu velikosti výpočetní domény pro CFD výpočet na odvod tepla. Vzhledem k nedostatku času na výpočet a vzhledem k tomu, že se jedná především o výpočet ukázkový, který má prokázat principiální schopnost řešit problematiku ohřevu a chlazení brzdových kotoučů při různých režimech jízdy včetně navazujících pevnostních analýz, nebyly verifikace výsledků a odladění modelu provedeny. Výsledky jsou předvedeny tak, jak byly získány „na první pokus“, a nelze je považovat za výsledky finální. 9

Závěr

Provedená analýza prokazuje, že vzorový výpočet podle zadávací dokumentace k úlohám k výběrovému řízení je úspěšně řešitelný při použití konfigurace ANSYS FLUENT a ANSYS Mechanical. Výsledky jsou ovlivněny zadáním úlohy. Pro přesnější výpočty je potřeba především uvažovat veličiny pevných částí závislé na teplotě, znát vlastnosti pevných částí z pohledu radiačního modelu a uvažovat zahrnutí součástek kolem kola ovlivňující proudění a přenos tepla (zejména z brzdových desek). Provedená analýza nebyla srovnána s experimentálním měřením. Má za cíl ukázat možnosti dnešních moderních CFD programů řešit problematiku ohřívání a chlazení brzdných kotoučů.




THERMODYNAMIC DESIGN OF THE ROOF ING.ARCH. ING. MILAN PALKO, PHD. Slovak University of Technology in Bratislava , Faculty of Civil Enginering Abstract: The paper deals with optimal thermal design of lightweight roofs. The main supporting structure is a whip element. The carrier is compared with the thermal insulation material far worse thermal properties. Such a roof construction creates thermal bridges. Properly optimized design can reduce the loss of roof shell up to 20% compared to conventional roof construction. Application of appropriate design is particularly important when the roofs of highenergy standards, such as low energy, passive houses and zero houses. Keywords: roof, thermal insulation, wood 1

Introduction

Fulfillment of requirements of heat-moisture and proper choice of constructional materials solutions are essential for a quality flat roof construction. Global trends in energy leads to a reduction in their consumption. This fact is reflected by increasing the thermal properties of all of the heat-exchange building envelope. When the total energy concepts and optimization of thermal-moisture plays a roof structure and the segmental part of the essential role. Reflection is that prescriptive requirements are stringent to the roof. 2

Material basis for solving thermal-moisture problems

In terms of thermal-moisture materials can be incorporated according to their most important characteristics into three basic groups: -

thermal insulation,

-

vapor barrier (vapor brakes)

-

premium waterproofing.

The primary ability of the insulation to minimize heat loss and create a thermal insulation barrier between the interior and exterior. As thermal insulation in pitched roofs are used: -

mineral wool

-

extruded polystyrene,

-

polyurethane foam,

-

multi-layer reflective insulation

-

fibreboard,

-

cellulose

-

straw

-

hemp



-

cork

-

wool,

-

and others.


Mineral wool is the most commonly used thermal insulation material pitched roofs and is used in the form of sheets or webs. Extruded polystyrene is produced in the form of pressed board covered with aluminum foil with integrated galvanized beam for fixing tiles. Polyurethane foam is a highly effective thermal insulation material used in the form of plates of hard polyurethane foam. Multilayer reflective insulation are based on the principle of multiple layers separated by air reflective sheeting. Fibreboard is made from coniferous wood garnetted. Cellulose thermal insulation materials made from recycled newsprint, the basic raw material is wood, the very beginning. The insulation is applied by blowing. Straw is one of the oldest and long-used insulating materials. It is used in the form of packages or pressed board. Its advantage is rapid recoverability. Hemp is one of the frequently used technical plants. Its biggest strength is like the straw fast recoverability. The fibers of this plant to produce insulation in the form of sheets or webs. Cork is obtained from the bark of cork pine. It is used in the form of granules or cork tiles. Sheep's wool is produced in the form of mats for thick wool insulation as a gland.That thermal energy-saving materials in use at the time of construction is obvious, but it is not the only energy relationship. For thermal insulation materials can evaluate the energy required for their production and placing on the site. This evaluation involves extracting raw materials, transportation of basic materials, manufacturing, transport to the construction and installation. Primary energy requirement for production of selected insulation materials in MJ / kg of finished product is for example: -

mineral wool 17.5 MJ / kg (30-150 kg/m3),

-

extruded polystyrene 101 MJ / kg (37 kg/m3),

-

polyurethane (hard foam) 100 MJ / kg (30 kg/m3),

-

cellulose fibers 4.2 MJ / kg (30-70 kg/m3),

-

insulating panels of wood wool 4.1 MJ / kg (400-530 kg/m3),

Steam barrier (steam brake) is a layer, you must include any roof. Has a high diffusion resistance. Proper placement of vapor barriers respectively steam brakes are very important (Fig. 1). Its position has to be physical but also structurally correct. From the physical point of view, this layer should be placed as far as possible into the interior. From a structural point of view should respect the fact that this layer is preferably at least perforated anchors and technological devices (eg, electrical switch)




Figure 1 - Location steam brakes, vapor barriers in the construction of roof (A - inadequate design, Badequate design)

Insurance waterproofing is usually located under the ventilated air cavities. It is characterized by low diffusion resistance and the ability to divert water to the upper cheek. A secondary function is to prevent the cooling of thermal insulation composed of fibers due to cold air flow. Insurance hydro-known non-contact and contact. Under the non-contact must be kept air layer and not touch or thermal insulation layer (if prejudice is not waterproof). In case of contact insurance waterproofing is important to avoid the bulge in the direction of the air cavity (print may cause thermal insulation). The locations of perforations should be applied selfadhesive sealing tape. Appropriate handling of insurance waterproofing layer is shown in Fig. 2.A. The case is airtight and prevent cooling thermal insulation layer is important prelepovanie connections between bands (Fig. 2.B). Special importance in terms of vapor and insured waterproofing details to manage the changing flights through the roof structures such as casing pipes, a chimney, but also construction roof windows and climbed out.




Figure 2 - Insurance waterproofing. A - Use self-adhesive sealing tape prism, B - Mutual insurance adhesive waterproofing relabelling belts parallel to the eaves

3

Examples of solutions of thermal bridges in the rig

Roofs for low energy and passive standard can be addressed as well as inclined flat material, the thickness of the insulation, which in these cases is 260-400 mm (depending on the type of insulation). In terms of choice of the structure difficult to distinguish (eg, reinforced concrete structure) and light construction (eg construction with lightweight truss). Construction of pitched roofs in this standard is substantially similar to the construction of exterior walls in wooden house. As structural elements are used perpendicular beams or planks in the form of gratings, I-beams and other shapes (a combination of prism and OSB). Sense to optimize the design rafters in the standard of construction is very important. The reason for such a solution is that the thermal conductivity of wood is 0.18 W / (mK) and thermal insulation is 0.04 W / (mK), which is more than fourfold difference. Heat flow through this optimized element rafters can be reduced by up to 20%.

Figure 3 - Examples of design beam (A-inappropriate, B, C, D, a suitable solution)

4

Computer simulation Models for different types of structure are the following figures.




Figure 4 - Examples of design beam (A) - isotherm and total heat flux

Figure 5 - Examples of design beam (B)- isotherm and total heat flux




Figure 6 - Examples of design beam (C) - isotherm and total heat flux




Figure 7 - Examples of design beam (D) - isotherm and total heat flux

5

Conclusion

Correct application of all knowledge gained through long years of development roofs using technology and using today's materials engineering principles with predestination materials may be a roof protecting the ecology of environment. Protects maintain optimal indoor climate change with regard to energy conservation. Reviewed: Ing. Michal Frimmer, PhD. Reference STN 73 0540 – 2, Tepelnotechnické vlastnosti stavebných konštrukcií a budov. Tepelná ochrana budov.Časť 2: Funkčné požiadavky.




TYWONIAK,J.: Nízkoenergetické domy princípy a príklady. Vydala Grada Publishing,a.s. Praha 2005.




WINDOW FOR PASSIVE HOUSES AND ITS LOCATION IN THE WALL ING.ARCH. ING. MILAN PALKO, PHD. Slovak University of Technology in Bratislava, Faculty of Civil Enginering Abstract: The paper deals with the optimal position of the window in the wall of a passive house. Windows are the weakest in terms of heat loss cladding elements. Thermal transmittance of the window is five to eight times worse than other parts of the heat transfer surface of the building. In type buildings with the need to minimize heat loss is necessary to optimize the window geometry, material and construction with regard to the surrounding structures and environment. Keywords: window, passive houses 1

Introduction

Window as transparent building element is a significant design and architecture, resistance to forced entry, the effect of sun, heat, wind, cold, rain, or the ability to resist mechanical loads, fire, acoustic and so on. Currently the window, as part of the heat-exchange jacket design packaging important element of the energy concept of the building (especially for passive house). When designing the location and size of glazed surfaces is important to look beyond the avoidance of heat and solar gain. It would be misleading to believe that the transparent parts of the building envelope interaction internal environment - outdoor climate affects only the size of the glazed area. An important feature of quality window is the correct solution of structural parts of the window and its incorporation into the envelope. 2

Passive house

Passive house is building, to ensure that thermal comfort does not need conventional power supplying system. To achieve this, we must have this house heat energy demand than conventional house reduced by 90%. This density heat demand (kWh / (m2.year)) is the normal modern house about 100, in an old apartment building in an average of 200, low energy house in under 50, in a passive house below 15 kWh/m2.rok. Specific heat demand below 15 is an essential characteristic of the passive house and shows her calculation. 3

Thermal properties and requirements for a passive house window

Thermal properties of windows are characterized by two variables. Heat transfer coefficient (W / m2. K) and the inner surface temperature (° C). These values can be obtained by measurement or calculation. Measurement of the window construction (frieze frame and sash, glass system, glass-frame contact system) determine the declared value of the coefficient of heat transfer around windows, the boundary conditions of the internal air temperature, outside air temperature and relative humidity. For the measuring procedure is EN ISO 12567-1 - Thermal performance of windows and doors. Determination of thermal transmittance hot box method. By means of guarded hot box temperature creates a steady state to be on the basis of equation (1) characterized by a transparent product specific dimensions. http://aum.svsfem.cz



(1)

The new approach towards the quantification of thermal characteristics of opening constructions used in European standards included in the Slovak technical standards system (e.g. STN EN ISO 10077-1 Thermal characteristics of windows, doors and shutters. Calculation of heat transfer coefficient) enable designers to show the value of heat transfer coefficient for each window construction based on the values declared by an accredited company using the following equation:

(2)

where:

For windows it is relevant what kind of thermal characteristics and geometry the material used in joint glazing – casement has, as the joint influences the additional thermal flow, which has significant impact on indoor surface temperature (condensation and mould is possible), heat losses through window. U g value (W/(m2 . K)) is determined for a concrete glazing




system – without taking into consideration the edge influence (Ug,NED  1,1 W/(m2 . K)). Similarly, also the Uf value (W/(m2 . K)) is determined for frame without glazing. Linear thermal transmittance g (W/(m . K)) expresses the impact of joint glazing – casement on frame thermal flow, that is the deformation of temperature field by additional thermal flow. g is a function of the used glazing system, frame material and spacer system. It is necessary to consider the linear thermal transmittance is not a material constant factor, but is dependant on specification of joint frame – glazing. From the above stated it follows that the frame as well as glazing systems have constant declared calculation values unlike "characteristics" of glazing edge. The same material of spacer will have different additional thermal flow in wooden, PVC or aluminium frame depending on combination frame – glazing, and, consequently, it will have an impact on Uw value and, thus, also the energy balance. Using the glazing with improved spacer system – warm edge opposite to aluminium one, the effect of direct dependence of additional thermal flow formation is significantly pressed decrease of g value. STN EN ISO 10077-2 Thermal performance of windows, doors and shutters. Calculation of thermal transmittance. Part 2: Numerical method for frames, describes the way of determination of linear thermal transmittance also in joint of glazing as follows:

(3)

where:

4

Effect of stocking location of windows for heat loss

As already mentioned in properly designed passive house is necessary to place emphasis on segmental parts (such as stocking location window in the outer wall). Cladding passive houses may have variants (light wooden structure and thermal insulator, heavy and monolithic thermal insulator). Of these variants also receives the correct type of window recess in the outer wall. In this paper elaborates on the shoulder in a simplified box (homogenized) outer wall for two types of frame friezes. a. Geometric parameters of the computer model Geometry window is determined by five window locations in the outer wall in relation interior and exterior. Models are two alternatives EUROPROFIL 78 mm and 92 mm (Fig. 1).




Figure 1 - EUROPROFIL 78 mm and 92 mm

Glazing has been used as thermal triple the same building and physical quantification of both systems. Position of installation of the windows are shown in Fig.2

Figure 2 - Position of the window recess




b. Results of computer simulation

Figure 3 - The temperatures for the models' Euro 78 with a recess 78-1,78-2,78-3,78-4,78-5




Figure 4 - The temperatures for the models' Euro 92 with a recess 92-1, 92-2, 92-3, 92-4, 92-5

5

Conclusion

Increased thermo-technical quantification fragment walls and window elements eliminates the problems of surface temperatures. For this reason, is limited to a stocking location in terms of minimizing heat loss. The most appropriate location of the installation of windows in terms of minimizing losses is a model 78-3 and 92-3. Given the characteristics of segments of the outer wall are very small differences, to be negligible. Can essentially rule the argument that in a simplified homogeneous wall can be fitted with a window in terms of heat loss in any place.




Figure 5 Example embedding windows in external wall

Reviewed: Ing. Michal Frimmer, PhD. Reference CHMÚRNY, I.: TEPELNÁ ochrana budov. Vydavateľstvo Jaga group - Bratislava 2003




HEAT AND MOISTURE ANALYSIS OF FUNCTIONAL GAPS FOR WOODEN WINDOWS ING.ARCH. ING. MILAN PALKO, PHD. SLOVAK UNIVERSITY OF TECHNOLOGY IN BRATISLAVA , Faculty of Civil Enginering Abstract: The paper deals with thermal-moisture performance of a wooden window in standard external and indoor climate conditions in a gap between casement and window frame. Condensation or ice coating in the gap is the main problem of wooden windows. Thus, destructive processes on the wood surface finish causing the moisture penetration into the wood occur. The surface condensation is sufficiently covered by the requirements in currently valid standards. The condensate generation in the gap is not described neither limited by the standards nor regulations. The frequent problem occurrence in practice is to be dealt with and solved. The experimental and computing investigation was carried out using four types of MIRADOR 682, 682T, 783, 923 euro-profiles. For computing investigation the ANSYS programme was applied. Keywords: wooden window, water vapour, condensate 1

Introduction

Thermal-moisture performance of a window in external and indoor climate conditions is a subject of our investigation. Water vapour condensation on surface of window structure is not desirable phenomena. The condensate generation in functional gap is not described neither limited by the standards nor regulations. Seriousness of the problem varies depending on the material base. Considering the wooden window structures the degradation of surface finish occur and its thermal characteristics change. Generally, there are five condensation types occurring in different intensities: 

indoor glazing condensation,



glazing joint condensation,



functional gap condensation,



external glazing condensation,



supplementary elements condensation (shutters).

The moisture transport through gap between casement and window frame at condensation temperature is the main cause of condensation of water vapour. The transport causes:

2



Air pressure differences (taking water vapour) – infiltration a exfiltration,



Diffusion of water vapour (different saturation of water vapour in exterior and indoor air).

Infiltration and exfiltration

Infiltration and exfiltration are basically permeability of window structure towards indoor or outdoor. The airflow through window structure is a function of pressure differences




and window structure design and its shoulder into building structure. The phenomena causing the infiltration or exfiltration can be divided into the following groups: 

wind impact (building exterior aerodynamics),



air temperature impact (due to gravity),



forced ventilation (building services),



combination of the above given phenomena.

Wind impact – if there is an obstacle on the airflow path, the airflow on windward is stopped, on both sides is accelerated and there is turbulence on leeward. The pressure on windward is, therefore, positive and on all other sides the pressure is negative (suction). The building modifies the airflow in its environment not only due to its shape, but also by surface structures, geometry of elements overhanging the façade as well as physical characteristics of its envelope and partitions (air permeability). The infiltration or exfiltration occurs due to inequality of pressure distribution on building envelope and building tightness (including indoor partitions). Air temperature impact – At different pressure due to various air temperatures the air flows using the natural convection principle. The pressure difference values due to natural convection of warm air are from 0 to 18 Pa, depending on the thermal difference and the distance, on which the pressure difference occurs. The distance is the height between input and output apertures. Considering a multi-floor building the principle is, to a large extent, applied by staircase hall or elevator shaft. The case makes worse the air quality chiefly on upper floors, where mainly the exfiltration occur. The polluted infiltrated air is coming from lower floors. Forced ventilation – Building services (pipes and ducts) are applied at force ventilation. Certainly, central air conditioning is not considered (air offtake and intake), but only local venting of polluted or used air by forced ventilation (fan). The air is carried away, for example, from bathrooms, toilets and kitchens. In such cases there is a vacuum at exhaustion, which is compensated by infiltration through openings. The air intake for combustion in fireplaces is the other partly similar way (if the air intake is not provided by a pipe connected to the exterior). Combination – all the above given reasons are relevant at infiltration or exfiltration impact. The particular factors play different roles and they can be summed up or balanced out. Condensation of water vapour – phenomena for balancing water vapour partial pressure by interaction of molecules. The diffusing water vapours are moving from places having higher pressure to places with lower pressure. The indoor water vapour partial pressure (20°C and air humidity 50%) is calculated by: (1)

where: is indoor water vapour partial pressure in Pa, is relative humidity of indoor air in %, is saturated water vapour partial pressure of indoor air in Pa.




The outdoor water vapour partial pressure ((-11°C and air humidity 83%) is calculated as follows:

(2)

where: is outdoor water vapour partial pressure in Pa, is is relative humidity of outdoor air in %, is saturated water vapour partial pressure of outdoor air in Pa. It follows from the above given that the pressure difference between interior and exterior is 971,8 Pa.

Figure 1 - Schematic representation of heat transmission, diffusion of water vapour and airflow through window structure

The warmer air can contain more water vapour as the cooler one, e.g. the air with 20 °C can contain up to 17,25 g/m2 of vapour, but the air with -11 °C can have only 1,96 g/m2 of moisture. Interactive accumulation of infiltration or exfiltration and diffusion of water vapour leads to condensate or ice coating in winter period. 3

Experimental measurements and computer model

The experimental measurements were carried out in heat engineering laboratories at the Faculty of Civil Engineering. Measurement models represent real window structures. The big climatic chamber illustrated in Figure 2 was used at measurements. The chamber A represents exterior climate (outdoor temperature - 11°C, pressure difference between exterior and interior from 0 to 2000 Pa, heat transfer coefficient 25 W/m2K, relative humidity 50%). The chamber B represents balancing chamber for HOT-BOX and simulate the indoor climate (indoor temperature + 20°C, relative humidity 50%). The HOT-BOX is used for measurement of heat transmission using measured element simulating the indoor conditions (indoor temperature +




20°C, relative humidity 50%, heat transfer coefficient 7,7 W/m2K). Part D is a masking panel, in which the measured window structure is imbedded. In case of our measurement the HOT-BOX was not used as the conditions with infiltration or exfiltration were to be investigated.

Figure 2 - Scheme of laboratory equipment of big climatic chamber

Thermal and humidity parameters for ambient conditions are constant for all models. The indoor air temperature is 20°C and the relative air humidity is 50 %. The outdoor air temperature is -11°C and the relative air humidity is 83%. The pressure differences are 200, 100, 75, 50, 25, 10, 5, and 0 Pa for infiltration and the similar pressure differences are also for exfiltration. The measurements were performed on completed panel having real windows. Such panel having four windows with dimensions 540x695 mm, is illustrated in Figure 3. The MIRADOR 682 with no external sealing is window No. 1, MIRADOR 682 with external sealing is window No. MIRADOR 783 is window No. 3 and MIRADOR 923 is window No. 4.




Figure 3 - Geometric parameters of window set ups in masking panel

Measurement sensors are divided into two groups. The first group measures the surface temperatures (PT 100). The second one measures the temperatures and relative air humidity (SHT 75). The sensors distribution is presented in Figure 4. For windows No. 1, 2, 3 the temperatures and relative air humidity are measured (SHT 75). For window No. 4 the surface temperatures (PT 100), temperatures and relative air humidity are measured (SHT 75). In Figure 4 for MIRADOR 923 profile the placement of SHT 75 sensors on the left side and PT 100 on the right side are shown. Computer simulation was carried out for MIRADOR 923 profile that is used mostly at low-energy and passive building construction. The geometry of computer model is identical with the measured window. Boundary conditions and material characteristics are also compatible with the experimental model in such a way that the comparison can be done. For the calculation the ANSYS programme was used. Geometry and netting of MIRADOR 923 window computer model is given in Figure 5.




Figure 4 - Geometric parameters of window profiles and placement of measurement sensors




Figure 5 - Geometry and netting of MIRADOR 923 window computer model

4

Results of experimental measurements and computer model

The measurement and computer models results are classified into two groups. At experimental measurements in big climatic chamber the occurrence and amount of water vapour and ice coating condensate was found out after each completed measurement phase. The internal surface temperatures and air temperature including relative air humidity in gap measurement points (see Figure 4) belonged into the second valuation data group. The measured values are presented in Graphs 1 and 2. Considering the fact that it is not possible to distribute the sensors in ideal positions in gap the computer model was established. The computer simulation results show more detailed temperature distribution and enable more thorough investigation of condensation possibility in required positions. The simulation results are shown in Figure 6 and Figure 7.




Graph 1 - Course of surface temperature in time for the pressure difference (MIRADOR 923)

Graph 2 - Temperatures and relative humidity in time for the pressure difference (MIRADOR 923)




Figure 6 - Surface temperature fields and surface temperatures in investigated points (°C) of MIRADOR 923 window (upper part)

Figure 7 - Surface temperature fields and surface temperatures in investigated points (°C) of MIRADOR 923 window (lower part)




Table 1 Comparison of measured and calculated surface temperature values for MIRADOR 923

Table 1 presents the comparison of the results gained from experimental measurement and computer simulation. The surface temperature results in investigated gap in identical points for simulation and experiment are evaluated. It can be seen from the numerical difference in the last column that the values for both methods are comparable.




Figure 8 - Demonstration of ice coating and condensation in gap between casement and window frame at exfiltration and 200Pa pressure difference (lower part)

Figure 9 - Demonstration of condensation in gap between casement and window frame at exfiltration and 200Pa pressure difference (upper part)

1

Conclusion

At an initial investigation stage it was known that condensation originated only in equilibrium state and at different pressure impact causing the exfiltration. However, it was found out during the measurements that the condensation also originates due to the infiltration. The condensation boundary for water vapour was being finding after the correct set up of circuit forging and, thus, also the sealing compression. After the measurements the value of differences of pressures between cool and warm chamber at which the condensate formation stops in the range 60 to 75 Pa. The detailed computer simulation of the investigated problem is enormously demanding and the further investigation is needed. Although the mathematic algorithms for moisture transport are relatively well handled in ANSYS programme, it is not possible to continue in simulation without knowing the characteristics of used sealing materials. In order to succeed in




this area it is inevitable to know the air permeability of sealing and window structure connection, diffusion constant of sealing and window frame connection. The air permeability has been already surveyed, but the diffusion constant of sealing for windows has not been investigated so far. The thermal model, which was compared with the experimental measurements in the previous chapter, has been authentically calibrated in the ANSYS programme. The total temperature distribution on the window frame was obtained by the simulation. Reviewed: Ing. Michal Frimmer, PhD. Reference HAUSER, G. – KEMPKES, CH.: Der Einfluss von windinduzierten Druckschwankungen auf das thermisch- hygrische Verhalten von durchstromten Leckagen. Kassel, 2005 HENS, H.: Building Physics – Heat, Air and Moisture. Fundamentals and Engineering Methods with Examples and Exercises. Berlin:Ernst &Sohn A Wiley Company, 2007.ISBN 978-3-43301841-5




ALUMINUM WINDOW FOR THE PASSIVE HOUSE ING.ARCH. ING. MILAN PALKO,PHD. Slovak University of Technology in Bratislava , Faculty of Civil Enginering Abstract: The article deals with creation of residential buildings in passive standard. Lack of energy concentrated in fossil fuels and environmental pollution, leads us to create buildings in passive standard. The planning phase for this type of building is very important. It is necessary to create good energy concept. Transparent elements and their connection to other parts of the envelope must be optimized. A wide range of construction materials is forcing them to consider the choice and interoperability. Incorporation of materials should be properly coordinated. It is important to build high-quality building materials. Keywords: aluminum window, passive house 1

Introduction

Passive House is the object to ensure thermal comfort does not need conventional power supplying system. To achieve this, we must have this house needs heat to heat the house than conventional reduced by 90%. This density heat demand (kWh / (m2.rok)) is normally around 100 modern house in an old apartment building in an average of 200, lowenergy house below 50, the passive house under 15 (kWh / (m2.rok)). Specific heat demand by 15 (kWh / (m2.rok)) is an essential characteristic of the passive house and shows her calculation. As a transparent window of the building is a significant element of design and architecture, resistance to forced entry, the effect of sun, heat, wind, cold, rain, or the ability to resist mechanical loads, fire, acoustic and so on. Currently the window as part of the heatexchange structures of container housing an important element of the energy concept of the building (especially for passive house). When designing the location and size of the glazed surfaces is important to find the boundary between heat loss and solar gain. It would be misleading to believe that the transparent part of the building envelope interaction indoor - outdoor climate affects only the size of the glass-surface. An important feature of quality window is the correct solution of individual structural parts of the window and its incorporation into the envelope. 2

Requirements for transparent structures

To date, standard used windows for passive buildings are inadequate. For these houses should be applied to windows Uw ≤ 0.8 W / (m2. K). The window as a whole has reached the value of the coefficient of heat transfer must friezes window frame and sash to achieve value Uf ≤ 1.0 W / (m2. K). The friezes on the wood-based material it is a special prism thickness> 70 mm, where the median strip replaced by a highly efficient insulator. Appropriate to consider whether we eat, combined wood-aluminum base, where just part of the wooden slats with a highly efficient insulator can provide the desired value of the coefficient of heat transfer window and add an architecturally and functionally. Decisive influence on the thermal properties of the transparent window is part of the glazing, which accounts for 70% to 80% of the window. The window as a whole achieved the necessary heat transfer coefficient




value of the glazing is achieved (depending on the size of the window) Ug ≤ 0.6 to 0.35 W / (m2. K). That value can be achieved with selective triple insulating layers, filled with krypton. 3

Computer model and implementation

As already mentioned in properly designed passive house is necessary to place emphasis on segmental parts (such as position in the window recess siding). Cladding passive houses may have variants (light wooden structure and thermal insulator, heavy and monolithic thermal insulator). Of these variants also receives the right shoulder construction window cladding. In this paper we shoulder the heavy construction window envelope for non-traditional frame material basis for the EPD and the aluminum with intermittent heat bridge. On Figure 1 shows the various options in the window recess siding. The composition of the interior walls are plaster, concrete walls, thermal insulation of polystyrene, thin-film rendering system.

Figure 1 - Examined the position of the window recess in heavy cladding

Figure 2 - Surface temperature fields - model 1.1.1 (-11 ° C outdoor, indoor 20 ° C)




Figure 3 - Surface temperature fields - model 2.1.3 (-11 ° C outdoor, indoor 20 ° C)

Figure 4 - Example of implementation of the optimal embedding window detail (model 2.1.2)

4

Conclusion

Aluminum construction is non-standard windows for passive energy buildings, but also usable. Given that aluminum frames are of frame construction with hidden windows can be insulated by selecting the appropriate detail and his shoulder to achieve interesting results. The computer models, it is clear that the best model in terms of minimizing the heat flows and thus heat loss model is 2.1.3. When deciding on the appropriate model for the realization of entering other factors (eg, aesthetic solution, the possibility of additional applications mosquito nets ...). Taking into account all the factors in detail Replacing windows siding implemented the model 2.1.2 (Fig. 3). Reviewed: Ing. Michal Frimmer, PhD. Reference Chmúrny, Ivan. Tepelná ochrana budov. Vydavateľstvo Jaga group, v.o.s., Bratislava 2003 PUŠKÁR, A., SZOMOLÁNYIOVÁ, K., FUČILA, J., VAVROVIČ, B. Okná, zasklené steny, dvere, brány. Bratislava: Jaga Group, 2008. 266 s. ISBN 978-80-8076-062-5.




Title: 19th SVSFEM ANSYS Users' Group Meeting and Conference 2011 Date and place of the event: 19th – 20th of October 2011, Congress Center in Prague, Czech Republic Publisher: SVS FEM spol. s r.o., SVS FEM s.r.o., Škrochova 3886/42, 615 00 Brno-Židenice & TechSoft Engineering, spol. s r.o., Na Pankráci 322/26, 140 00 Praha 4-Nusle 255, 775 000 Editor: Ing. Miroslav Stárek, Ing. Antonín Tuček, CSc. Pages: 228 Edition: 1. Year: 2011, Prague, Czech Republic

ISBN: 978-80-905040-0-4


19th SVSFEM ANSYS Users' Group Meeting and Conference 2011

Recommend Documents