MATEMATICKÉ MODELOVÁNÍ VÝBUCHU METANU V RODINNÉM DOMKU V KAMENNÉ POMOCÍ SW FLUENT

PŘÍSPĚVKY

THE SCIENCE FOR POPULATION PROTECTION 2/2012

MATEMATICKÉ MODELOVÁNÍ VÝBUCHU METANU V RODINNÉM DOMKU V KAMENNÉ POMOCÍ SW FLUENT MATHEMATICAL MODELLING EXPLOSION METHAN IN FAMILY HOUSE IN KAMENNA USING FLUENT SOFTWARE Milada KOZUBKOVÁ, Jaroslav KRUTIL, Marian BOJKO, Otto DVOŘÁK [email protected], [email protected], [email protected], [email protected] Došlo 5. 5. 2012, přijato 6. 6. 2012. Dostupné na http://www.population-protection.eu/ attachments/041_vol4n2_kozubkova_krutil_bojko_dvorak.pdf.

Abstract This paper describes the issues of risks conditions associated with the explosion of gaseous mixtures. Method of mathematical modeling by ANSYS FLUENT software for solving is used. The work contains two alternatives of solution. Variant A not considering damage room, and variant B considering destruction room, in which the blast was initiated. Problem solution caused by the generation of explosion pressure waves is not easy. It should be mentioned that it is very important to orient oneself on physical knowledge, related to fluid flow, technical drawings, thermophysical properties and material knowledge. It is important to note that a major impact on solving of similar tasks depends on correct determination of chemical reaction constants as activation energy and pre-exponential factor. These constants have significantly influence on chemical reaction of gaseous fuel with an oxidizer. In the final summary the comparison of calculated data with the results of both experimental measurements and also with other problem oriented numerical software (FLACS) is evaluated. The essence of such works is verification of mathematical models for the fire technical expertise. This work should also contribute to better understanding of burning behavior of gaseous fuel mixtures in confined spaces and thereby significantly reduce the risk of such situations or prevent them. Keywords ANSYS FLUENT, CFD, explosion, methane, numerical simulation. ÚVOD Tento článek podává informace o tvorbě matematického modelu a následné simulaci výbuchu rodinného domku. Matematický model se snaží co nejvěrohodněji napodobit skutečnou velkorozměrovou požární zkoušku rodinného 1


PŘÍSPĚVKY

domu v Kamenné u Milína, která byla provedena dne 8. 10. 2010 společně institucemi MV – GŘ HZS ČR, TÚPO a VŠB – TUO, FBI [1]. Úkolem této práce je tvorba matematického modelu v programu ANSYS FLUENT pro stanovení vývinu/šíření tlakových vln a jejich ničivých účinků. Cílem této práce je ověření těchto modelů s výsledky získanými prostřednictvím experimentu a výsledky simulace provedenými v programu FLACS. Jádrem samotného matematického modelování je řešení hoření směsi plynného paliva a z toho vyplývající generace tlakové vlny, přestupu tepla, proudění plynů a zkoumání teplotních polí v zasažené oblasti. Problematika modelování výbuchu je velmi složitá a v programu ANSYS FLUENT existuje několik možných přístupů k její realizaci (akustický model, řešení pomocí přetlakového signálu, model využívající chemických reakcí). V tomto případu, kdy dochází ke generaci tlakové vlny v důsledku hoření plynné směsi, přichází v úvahu pouze možnost s využitím obecného modelu proudění plynů s chemickou reakcí (species transport and chemical reaction model). Avšak tato varianta byla rozpracována do dvou dílčích bodů. V první části byla úloha řešena tak, že nebylo uvažováno poškození střechy, oken a dveří v důsledku vzniku tlakové vlny (což bylo pozorováno v experimentu). V druhé variantě bylo toto poškození uvažováno pomocí tzv. porézní vrstvy. V této části řešení je také popsána teorie porézní vrstvy a vysvětleno její použití. 1

TEORIE

1.1

Matematický model proudění s přestupem tepla

Při výpočtech turbulentního proudění se u praktických inženýrských úloh využívá časově středovaných veličin, jež jsou v následných odstavcích a rovnicích označeny pruhem nad danou fyzikální veličinou. Je to způsobeno faktem, že při vysokých hodnotách Reynoldsova čísla nelze s ohledem na nutný počet buněk sítě a možnosti výpočetní techniky použít metodu DNS [2]. 1.2

Rovnice kontinuity pro proudění stlačitelné tekutiny

Rovnice vyjadřující zákon zachování hmotnosti se nazývá „rovnicí kontinuity“. Pro neustálené, tedy časově závislé proudění stlačitelných tekutin ji lze v diferenciálním tvaru vyjádřit takto:

∂ρ ∂ (ρ u j ) + =0 ∂t ∂x j

(1) ,

kde u i je „časově středovaná složka rychlosti“ proudění. 2

PŘÍSPĚVKY

1.3


Pohybové rovnice pro proudění stlačitelné tekutiny

Rovnice vyjadřující zákon zachování hybnosti se nazývají „NavierStokesovy rovnice“. Pro výpočty turbulentního proudění je však potřeba použít časově středovaných veličin. Po dosazení časově středovaných veličin do NavierStokesových rovnic nabývají tyto rovnice tvaru tzv. „Reynoldsových rovnic“. Rovnice pro přenos hybnosti u stlačitelných tekutin mají tedy tvar:

∂ (ρ u i ) ∂ (ρ u i u j ) ∂p ∂ + =− + ∂t ∂x j ∂xi ∂x j

⎛ ∂u i ⎜ μt ⎜ ∂x j ⎝

⎞ ⎟ + ρ δ i 3 g + ρ f c ε ij 3 u j + ρ f j ⎟ ⎠ ,

(2)

což odpovídá diferenciálnímu tvaru „rovnice pro přenos hybnosti“, kde g = − 9,81 m ⋅ s −1 je „gravitační zrychlení“ v případě účasti vztlakových sil. Rovnicemi pro vyjádření turbulentních veličin jsou myšleny rovnice pro „turbulentní kinetickou energii“ k a „rychlost disipace“ ε . Rovnici pro k lze odvodit z Navier-Stokesových rovnic a má tvar: ∂ ul ∂ u l′ ∂ u l′ ∂k ∂ujk p ′ ⎞⎤ ∂ ⎡ ⎛ u l′u l′ ∂ 2k ⎟⎟⎥ + ν 2 − u l′u ′j + =− + δ jl −ν ⎢u ′j ⎜⎜ ∂t ∂ xj ∂ x j ⎣⎢ ⎝ 2 ρ ⎠⎥⎦ ∂ x ∂ xj ∂ xj ∂x j j

(3)

„Turbulentní kinetická energie“ k je uvedená v rovnici (3) a je definována jako:

k=

(

)

1 2 1 u1′ + u 2′ 2 + u 3′ 2 = u ′j 2 2 2

(4)

Rovnici pro ε lze opět odvodit z Navier-Stokesových rovnic a má tvar:

∂ ε ∂ u jε ∂ + = ∂t ∂ xj ∂ xj

⎛ ∂u ⎛ νt ∂ε ⎞ ⎜ . ⎟ + C1εν t ⎜ j + ∂u l ⎜ ∂x ⎜σ ∂ x ⎟ j ⎠ ⎝ l ∂x j ⎝ ε

⎞ ∂u l ε2 ⎟ − C 2ε ⎟ ∂x k ⎠ j

(5)

Vztah pro „turbulentní viskozitu“ ν t je pak definován takto:

ν t = Cν

k2

ε

(6)

1.3.1

Rovnice energie Rovnice energie vyjadřuje „zákon zachování energie“, podle kterého je „celková změna energie“ E tekutiny v určitém „objemu“ V dána změnou 3


PŘÍSPĚVKY

„vnitřní energie“ a „kinetické energie“ a tokem obou energií „plochou“ S omezující objem V . Výsledná rovnice má tvar:

∂ ( pu j ) ∂ (τ j l u j ) ∂ q j ∂ [ρE ] + ∂ [ρu j E ] = ρu j f j − + − + Sh ∂t ∂x j ∂ xj ∂ xl ∂ xj

(7)

1.3.2

Transportní rovnice pro přenos příměsí FLUENT počítá v modelu „časově středované hodnoty lokálních hmotnostních zlomků příměsí“ Yi′ , které jsou popsány podobnou bilanční rovnicí, jako je tomu u rovnice energie (7) zahrnující řešení konvektivní a difúzní složky přenosu. Je využíváno vztahu, který má v konzervativní formě tvar:

∂ (ρYi′ ) + ∂ ⋅ (ρ u j Yi′ ) = − ∂ J j ,i′ + Ri′ + S i′ , ∂t ∂x j ∂xi

(8)

kde ui je „časově středovaná složka rychlosti“ proudění a na pravé straně je Ri′ „rychlost produkce příměsí i ′ “ vlivem chemické reakce a S i′ „rychlost tvorby

přírůstku z distribuované příměsi“. Výše uvedená rovnice platí pro N − 1 příměsí, kde N je „úplný počet komponent“ prezentovaných v matematickém modelu. Distribuce příměsí může být realizována za různých podmínek, obecně lze rozlišovat distribuci za laminárního nebo turbulentního proudění. J j ,i′ představuje „difúzní tok i ′ -té komponenty směsi“. Při turbulentním proudění FLUENT pro vyjádření difúzního toku i ′ -té složky uplatňuje vztah:

⎛ μ ⎞ ∂Y J i′ = −⎜⎜ t ⎟⎟ i′ ⎝ Sct ⎠ ∂x j ,

(9)

kde Sc t je „Schmidtovo turbulentní číslo“ (přednastaveno na hodnotu 0,7). 1.4

Matematický model řešení chemických reakcí

FLUENT používá pro řešení rychlosti produkce příměsí i′ vlivem chemické reakce několika modelů: laminární model (Laminar finite-rate model), turbulentní model (Eddy-Dissipation model), kombinovaný model (Finiterate/Eddy-Dissipation model) a EDC turbulentní model (Eddy-DissipationConcept). Každý z těchto modelů je vhodný pro určité podmínky průběhu chemické reakce [3], [10]. 4

PŘÍSPĚVKY


1.4.1

Laminar finite-rate model (laminární model) Model počítá chemické zdrojové členy užitím Arrheniova vyjádření, kdy turbulentní fluktuace jsou zanedbatelné. Laminární model je dostatečně přesný pro spalování s relativně pomalou dobou vlastní chemické reakce a zanedbatelnými turbulentními fluktuacemi jako je nadzvukové hoření. Zdrojový člen Ri′ z důvodu chemické reakce v rovnici pro příměs i′ je počítán jako součet N R reakčních zdrojových členů příměsí, které se na reakci podílejí

⎞ ⎛ E NR N − k N η′ η ′′ Ri′ = M i ′ ∑ (ν i′′ ,′k − ν i′′ ,k ) ⎜ Ak T βk e RT ∏ C j′ j ′ ,k − kb,k ∏ C j′ j ′ ,k ⎟ ⎟ ⎜ 14243 j′ =1 k =1 j ′ =1 k f ,k ⎠ ⎝

[ ]

[ ]

(10) ,

kde N je počet chemických příměsí, ν i′′, k stechiometrický koeficient pro reaktant

i′ v k -té reakci, ν i′′′, k stechiometrický koeficient pro produkt i′ v k -té reakci,

M i ′ je molární hmotnost příměsi i′ , k f ,k rychlostní konstanta pro k -tou přímou (dopřednou) reakci, k b ,k rychlostní konstanta pro k -tou zpětnou reakci, C j ′ látková koncentrace každého reaktantu a produktu příměsi j ′ v k -té reakci, η ´j´,k rychlostní exponent (řád reakce) pro reaktant a produkt j ′ v

η ´´j´,k

k -té přímé reakci,

rychlostní exponent pro reaktant a produkt j ′ v k -té zpětné reakci, Ak pre-

edexponenciání faktor Arrheniova výrazu, β k je teplotní exponent, E k aktivační energie reakce, R univerzální plynová konstanta a T teplota. Reakce může probíhat v homogenní fázi, mezi fázemi jednotlivých příměsí, nebo na povrchu, jejíž výsledkem je usazování nebo vznik fáze [4], [10]. 1.4.2

Eddy-Dissipation model (turbulentní model) Probíhá-li chemická reakce rychle, tak celková rychlost reakce je řízená turbulentním směšováním. Rozlišujeme dva typy reakcí, s nepromíchanými a promíchanými reaktanty. FLUENT poskytuje model chemické turbulentní interakce založeného na Magnussen a Hjertager (nazvaný eddy-dissipation model) [4]. Střední rychlost chemické reakce tvorby produktu i′ -té příměsí v k -té reakci je dána menší hodnotou ze dvou vyjádření

⎧ ⎪ ⎛ Y ⎞ ε ε ⎪ Ri ′ = M i ′ ∑ min ⎨ν i´´,k Aρ min⎜ ´ R ⎟ , ν i´´,k ABρ R ν k M k ⎝ R ,k i´,R ⎠ k =1 ⎪ ⎪⎩ NR

⎫ ⎪ ⎪ P ⎬ N ∑j´ ν ´´j´,k M j´ ⎪⎪ ⎭,

∑Y

P

(11)

5


PŘÍSPĚVKY

kde YP je hmotnostní zlomek jednotlivých produktů příměsí ( P ), YR je hmotnostní zlomek konkrétních reaktantů ( R ), A je empirická konstanta (rovna 4) a B je empirická konstanta (rovna 0,5). ρ je měrná hmotnost i′ -té příměsi. Rychlost chemické reakce je řízená časovým měřítkem k / ε směšování velkých víru na základě Spaldingova modelu eddy-breakup (rozpad víru). Proces chemické reakce probíhá, jestliže je proudění turbulentní ( ε / k < 0 ) [4], [10]. 1.4.3

Finite-rate/Eddy-Dissipation model (kombinovaný model) Dále je možno uvažovat kombinovaný finite-rate/eddy-dissipation model, kdy se rychlost reakce určí jak podle Arrhenia (10), tak podle Eddy-dissipation rovnice (11). Lokální rychlost reakce je udávána jako minimální hodnota z těchto dvou rovnic [4]. Přestože FLUENT dovoluje několika stupňové reakční mechanismy pro eddy-dissipation a finite-rate/eddy-dissipation model, lze u reakčních mechanizmů vyšších řádů očekávat ne příliš přesné řešení. Příčinou je, že několika stupňové reakční mechanizmy jsou založeny na Arrheniových rychlostech, které jsou rozdílné pro každé reakce. V eddy-dissipation modelu mají všechny reakce stejnou rychlost, proto by měl být tento model použit pouze pro jednokrokové (reaktant → produkt) nebo dvoukrokové (reaktant → přechodný produkt → produkt) obecné rovnice. Model nemůže předpokládat kinetickou kontrolu příměsí, jako jsou radikály [4], [10]. 1.4.4

Eddy-Dissipation-Concept (EDC) model (turbulentní model) V tomto modelu je zahrnuta kinetika několika krokového chemického mechanismu v turbulentním proudění. Předpokládá chemické reakce, jejichž děj probíhá v malých turbulentních strukturách, zvaných fine scaled. Zdrojový člen Ri′ vlivem chemické reakce pro příměs i′ je počítán podle (7), kde Yi′ je hmotnostní zlomek příměsi i′ , Yi′ * hmotnostní zlomek příměsi i′ pro fine scaled [4], Cξ je konstanta objemového zlomku (2,1377), C r je konstanta časového měřítka (0,4082), ν je kinematická viskozita [10].

⎛ ⎞ 1,5 ⎜ ⎟ 2 ⎛ν ε ⎞ ρ Cξ ⎜ 2 ⎟ ⎜ ⎟ ⎝k ⎠ * ⎟ Ri ′ = Yi ′ − Yi ′ ⎜ 2.25 ⎜ ν ⎡ ⎛ 3 ⎛ν ε ⎞ ⎞ ⎤ ⎟ ⎢1 − ⎜ C ⎜ ⎟ ⎟ ⎥ ⎟ ⎜ Cr ε ⎢⎣ ⎜⎝ ξ ⎝ k 2 ⎠ ⎟⎠ ⎥⎦ ⎠ ⎝

(

)

(12)

Pokud jsou modelovány laminární reakční systémy užitím laminárního finite-rate modelu, budeme pravděpodobně využívat řešič s názvem coupled, 6

PŘÍSPĚVKY


pokud bude/nebude výpočet konvergovat. V EDC modelu je nutné nastavit: pravděpodobnost limitní rychlosti pro teplotu (formuluje změnu mezi teplotou vypočtenou a teplotou z předchozího kroku výpočtu, implicitně nastaveno na 0,2), časový krok teplotního redukčního faktoru (limit pro lokální CFL číslo, pro případ, kdy je změna teploty příliš výrazná, přednastaveno na 0,25) a přípustné maximum poměru časové měřítko/chemické časové měřítko (limitní lokální hodnota CFL čísla, přednastaveno na 0,9). Implicitně nastavené hodnoty jsou vhodné pro velkou škálu modelů [4], [10]. V matematickém modelu bylo využito rovnice pro dokonalé spalování metanu:

CH 4 + 2O2 → CO2 + 2H 2 O

(13)

Matematické modely řešící chemické reakce plynů jsou založeny na řešení transportních rovnic pro hmotnostní zlomky příměsi s definovaným reakčním mechanismem chemické reakce. Rychlosti reakce, která se objeví jako zdrojové členy v rovnicích pro přenos příměsi, jsou počítány v případech laminárního proudění z Arrheniových výrazů pro rychlost, v případech turbulentního proudění z modelu turbulentní (eddy) disipace dle Magnussena Hjertagera nebo z EDC (Eddy-dissipation-concept) modelu [5]. Proto mají zásadní vliv na správnou realizaci výpočtu konstanty aktivační energie a pre-exponenciálního faktoru. V odborné literatuře je mnoho variací těchto konstant např. Zambon Chelliah (ZC), Puri-Seshadri (PS), WD (Andersen-et-al), CM2 (Bibrzycki-Poinsot) atd… V našem případě se nejvíce osvědčily konstanty dle Zambon Chelliah, které nabývají těchto hodnot: Pre-exponencial factor: 500000000000 J.kmol-1 Activation energy: 47800 cal.mol-1 Řešení matematických modelů je provedeno v programu ANSYS FLUENT (CFD - Computational Fluid Dynamics). Tento program je založen na metodě konečných objemů a řešení základních rovnic umožňující komplexní řešení úloh z oblasti turbulence, přenosu tepla atd. 1.5

Fyzikální vlastnosti

U hustoty byl nastaven parametr ideal-gas. To znamená, že hustota plynu je počítána podle stavové rovnice.

p ⋅V = konst T ⋅ ri′

(14)

Ostatní fyzikální veličiny se definují v závislosti na teplotě experimentálně zjištěnými závislostmi, jako polynom, tabulka atd. Podle kinetické energie [6] ideálního plynu mohou být definovány následující fyzikální vlastnosti jednotlivých plynů a parametry [7]: 7


PŘÍSPĚVKY

• viskozita, • tepelná vodivost, • měrná tepelná kapacita, • koeficienty difúze hmoty (pro multi-speciální druhy směsi). Definice dynamické viskozity μ i′ při užití kinetické teorie je následující:

μi′ = 2.67e − 6

MT σ Ω μ , i′

(15)

2

,

kde

( )

Ω μ , i′ = Ω μ T *

a T = *

T (ε / k B ) .

(16)

Funkce Ω μ , i′ je experimentálně určenou závislostí na bezrozměrné teplotě, např. [8]:

Ω μ ,i′ =

1.16145

(T )

* 0.14874

+

0.52487 2.16178 + * exp(0.77320T ) exp(2.43787T * )

(17)

Vzorec pro měrnou tepelnou kapacitu plynů c p , i′ je vyjádřen jako funkce teploty polynomem n-tého řádu:

c p ,i , (T ) = A1 + A2 ⋅ T + A3 ⋅ T 2 + A4 ⋅ T 3 + A5 ⋅ T 4

(18)

Tepelná vodivost λ i′ při užití kinetické teorie je vyjádřena takto:

λi′ =

15 R ⎡ 4 c p ,i′ M 1 ⎤ μ i′ + ⎥ 4 M ⎢⎣15 R 3⎦

(19)

Fyzikální vlastnosti směsi plynů se pak určí podle směšovacích zákonů. 2

NUMERICKÉ MODELOVÁNÍ

2.1

Definice geometrie a výpočetní sítě oblasti

Geometrie a výpočetní síť vychází z přiložené výkresové dokumentace, avšak není zcela totožná. Jelikož je z výsledků experimentu patrné, že v obývacím pokoji dochází jen k velmi nepatrným změnám tlaku v závislosti na čase, je tato 8

PŘÍSPĚVKY


místnost vynechána z matematického modelu. Výkresová dokumentace také podává detailní informace o umístění tlakových senzorů. Umístění tlakových senzorů v matematickém modelu je shodné s požární zkouškou. Z důvodu lepší orientace a čitelnosti při srovnávání výsledků jsou tyto měřící body na modelu označeny stejně jako tlakové senzory při experimentu. Výpočetní síť domu byla vytvořena v programu Workbench 13.0 s celkovým počtem buněk 75110. Jedná se o nestrukturovanou síť vytvořenou pomocí prvků mnohostěnu různorodých tvarů. Geometrie modelu byla zjednodušena tím, že se neuvažuje nábytek, příslušenství a vybavení uvnitř místností. Naproti tomu byly pro variantu uvažující poničení části místnosti vytvořeny specifické oblasti, které budou definovány jako porézní zóny. Tyto oblasti jsou umístěny podél oken, dveří a pod střechou a budou v modelu simulovat vysklení oken, vyražení dveří a nadzdvižení spolu s poškozením střechy komory. Geometrie matematického modelu domku je vyobrazena na obr. 1. Kde jsou zeleně vyznačeny porézní zóny. Varianta, u které neuvažujeme poničení, pochopitelně tyto oblasti neobsahuje. U této varianty jsou okna, dveře a část střechy uvažovány jako tlakové výstupy (pressure outlet).

Obr. 1 Geometrie a síť matematického modelu (vlevo) Detail řezu sítě a zobrazení nejčetnějších tvarů buněk (vpravo) 2.2

Okrajové podmínky

V podstatě byly veškeré důležité okrajové podmínky pro obě varianty totožné a byly kompletně převzaty z dostupných podkladů o experimentu (tedy z přiložené výkresové dokumentace). Okrajové podmínky pro všechny tři výstupy (okno, dveře venkovní, dveře do ložnice) jsou definovány jako tlakový výstup (pressure outlet) a jsou označeny červeně, viz obr. 1. Všechny další stěny jsou 9


PŘÍSPĚVKY

definovány jako Wall tedy pevné stěny. Teplotní podmínky na výstupech a uvnitř řešené oblasti byly definovány z dostupných informací o povětrnostních podmínkách jako T = 10°C. Proudícím médiem celé oblasti je plynná směs vzduch s metanem a je definována tímto složením N2 = 72,68 %, O2 = 19,32 % a CH4 = 8 %. K nastavení plynné směsi do řešené oblasti bylo využito funkce inicializace. 2.3

Definice porézní vrstvy

Jak už bylo uvedeno výše, u varianty uvažující i simulace vysklení oken, vyražení dveří a poničení střechy byla řešena za pomoci modelu proudění média přes porézní vrstvu. Tato funkce v programu ANSYS FLUENT zahrnuje dvě možnosti definování oblasti, přes které se výpočet řeší. A jsou to tyto: plocha - je možné využít plochy dané geometrie a těmto plochám nastavit porézní koeficienty; objem - využije se v případech 3D těles, kde je nezbytné nastavit porézní vrstvu o daných koeficientech v celém objemu. V našem případě bylo využito objemového zadání porézní oblasti. V prostředí ANSYS FLUENT se 3D porézní zóna řeší v menu Cells Zone, kde se vybraná oblast nastaví jako porézní (Porous Zone) [9]. Teorie porézního prostředí vychází z tzv. Darcyho zápisu porézního prostředí. V rovnici průtoků pro porézní prostředí je přidaný člen S i složený ze dvou částí laminární a turbulentní ztráty. V případě homogenního porézního prostředí má tento tvar [4]:

1 μ S i = −( v i + C 2 ⋅ ⋅ ρ v v i ) α 2 v

velikost rychlosti

α

propustnost součinitel odporu

C2

(20)

Laminární ztráty v porézním prostředí (Viscous resistence) V případě laminárního proudění skrz porézní vrstvu je ztráta tlaku úměrná rychlosti a konstanta C2 je nulová. Model porézního prostředí se poté redukuje na Darcyho zákon [4], [9]:

∇p = −

μ ⋅v α

Ztráta tlaku pro směry x, y, z má pak následný tvar: 10

(21)

PŘÍSPĚVKY


μ ⋅ v j Δn x j =1 α xj 3 μ Δp y = ∑ ⋅ v j Δn y j =1 α yj 3 μ Δp z = ∑ ⋅ v j Δn z j =1 α zj 3

Δp x = ∑

(22)

1 / α ij

položka v matici D

vj

složka rychlosti x, y, z

Δn x, y , z

tloušťky média v osách x, y, z

Turbulentní ztráty v porézním prostředí (Inertial resistence) U případu např. děrovaného plátu je možné eliminovat propustnost alfa a použít pouze inertní ztráty pro směry x, y, z [4], [9]: 3 1 Δp x ≈ ∑ C 2 xj Δn x ⋅ ( ⋅ ρv j v ) 2 j =1 3 1 Δp y ≈ ∑ C 2 yj Δn y ⋅ ( ⋅ ρv j v ) 2 j =1 3 1 Δp z ≈ ∑ C 2 zj Δn z ⋅ ( ⋅ ρv j v ) 2 j =1

(23)

V našem případě byly koeficienty viscous a inertial resistance faktor orientačně odvozeny na základě naměřených hodnot (měření bylo převzato z [9]) rychlosti a ztráty tlaku pomocí následujícího postupu: A) Z naměřených dat rychlosti a tlaku proudícího vzduchu se vytvoří graf a proloží se rovnicí regrese. Tabulka 1 Tabulka naměřených hodnot z experimentu v (m.s-1) 0 1,57 2,37 3,55 4,64 5,61

p (Pa) 0 125 254 455 780 941 11


PŘÍSPĚVKY

Graf 1 Závislost rychlosti na tlakové ztrátě B)

Z následujících vzorců se vypočítají koeficienty 1 / α a C2 :

μ μ ⋅t ⋅t ⇒ α = = 3964180 m −1 α 61,643 1 20,006 ⋅ 2 20,006 = C 2 ρ ⋅ t ⇒ C 2 = = 32,66286 m −2 ρ ⋅t 2

61,643 =

2.4

Definice materiálových vlastností

V oblasti matematického modelu jsou uvažovány dva základní materiály. Jsou to materiály proudícího média (fluid) a materiál stěn (solid). Proudícím médiem je uvažována směs, která odpovídá fyzikálním parametrům o tomto složení N2 = 72,68 %, O2 = 19,32 % a CH4 = 8 %. Tato směs je definována parametry, které jsou naznačeny v tab. 2. Jelikož se směs skládá z jednotlivých plynů, je také nezbytné nadefinovat fyzikální parametry u samotných prvků směsi. Definování těchto parametrů je u všech prvků směsi totožné a proto jsou naznačeny materiálové vlastnosti pouze pro jeden plyn (viz tab. 3). 12

PŘÍSPĚVKY


Pevný materiál (solid), použitý v modelu, je cihla. Tento materiál je nastaven pro všechny pevné a nepohyblivé části v domě (strop, podlaha a jednotlivé stěny uvnitř domu). Fyzikální parametry jsou detailně popsány v tab. 4. Tabulka 2 Fyzikální parametry zadávané do programu FLUENT pro proudící směs Složení směsi (vzduchu) na vstupu

objemová %

dusík

%

72,68 %

kyslík

%

19,32 %

metan

%

8%

%

100 %

Suma Hustota směsi (ρ)

kg.m-3N -1

ideal-gas

-1

Specifické teplo (Cp)

J.kg .K

Tepelná vodivost (λ)

W.m-1.K-1

ideal-gas-mixing-law

Viskozita (η)

kg.m-1.s-1

ideal-gas-mixing-law

Hmotnostní rozptyl (D) Koeficient teplotní difúze (DI)

2 -1

m .s

mixing-law

kinetic -theory

-1 -1

kg.m .s

kinetic-theory

Tabulka 3 Fyzikální parametry popisující jeden plyn směsi (O2 – kyslík) Specifické teplo (Cp)

J.kg-1.K-1

piecewise-polynomial

Tepelná vodivost (λ)

-1

-1

kinetic-theory

-1 -1

W.m .K

Viskozita (η)

kg.m .s

kinetic-theory

Molekulová hmotnost (Mr)

kg.kgmol-1

31.9988

j.kgmol-1

0

Entalpie (H) Entropie (S)

-1

j.kgmol .k

-1

205026,9

K

298,15

L - J charakteristická délka

Angstrom

3,458

L - J energetický parametr

K

107,4

Referenční teplota (Tref)

13


PŘÍSPĚVKY

Tabulka 4 Fyzikální parametry popisující pevné materiály v modelu Hustota (ρ) Specifické teplo (Cp) Tepelná vodivost (λ)

3

kg.m-3N J.kg-1.K-1 W.m-1.K-1

750 1040 0,86

VYHODNOCENÍ MATEMATICKÉ SIMULACE VÝBUCHU – ŘEŠENÍ ČASOVĚ ZÁVISLÉ ÚLOHY PROUDĚNÍ PLYNŮ S CHEMICKOU REAKCÍ

Grafické zhodnocení je založeno na časové závislosti změny teplotních, rychlostních, tlakových atd. polí. Předmětem zájmu bylo především sledování změn tlakového pole. Byly provedeny dva základní způsoby vyhodnocení dosažených výsledků. V první fázi byly vyhodnoceny formou vyplněných kontur výsledky výpočtu, již výše zmiňovaných polí. Znázornění těchto tlakových polí je pro lepší viditelnost provedeno pomocí vhodně umístěných řezů přes modelovanou oblast. Byly přidány další veličiny, které označují úbytek hmotnostního zlomku metanu vlivem hoření a teplo vznikající z chemické reakce směsi (místo, kde dochází k hoření metanu). Pro variantu A je navíc provedeno srovnání průběhu výbuchového tlaku v čase v programech ANSYS FLUENT a programu FLACS. A ve variantě B je provedeno srovnání průběhu výbuchového tlaku v programu ANSYS FLUENT s experimentem. 3.1

Varianta A neuvažující poškození místnosti

Na obrázcích 2 až 6 jsou graficky zobrazeny výsledky varianty, která neuvažuje poničení oblasti. Porovnání je provedeno v časovém kroku 0,125 s a 0,135 s, kdy je nejlépe vidět vznik tlakové vlny.

Obr. 2 Tlakové pole v čase 0,125 a 0,135 s [Pa] 14

PŘÍSPĚVKY


Obr. 3 Rychlostní pole v čase 0,125 a 0,135 s [m.s-1]

Obr. 4 Teplotní pole v čase 0,125 a 0,135 s [K]

15


PŘÍSPĚVKY

Obr. 5 Úbytek hmotnostního zlomku metanu vlivem hoření plynné směsi v čase 0,125 a 0,135 s

Obr. 6 Místo hoření plynné směsi v čase 0,125 a 0,135 s (Heat of reaction [W])

16

PŘÍSPĚVKY


Graf 2 Srovnání průběhu výbuchového tlaku (nahoře ANSYS FLUENT dole FLACS)

17


3.2

PŘÍSPĚVKY

Varianta B uvažující poškození místnosti

Na obrázcích 7 až 11 jsou znázorněny výsledky matematického modelu, který zahrnuje poškození místnosti. Z obrázku 7 je patrný vznik a šíření tlakové vlny v čase 0,3 s (vlevo) a v čase 0,4 s (obrázek vpravo). Tento průběh potvrzují také rychlostní pole, která jsou naznačena na obr 8. Rychlostní profily jsou opět naznačeny ve stejném časovém intervalu. Je také patrné, že se tlaková vlna šíří všemi směry stejně. Na snímcích 9, které popisují teplotní pole, je velmi dobře vidět vznik plamene v řešené oblasti. Šíření teploty postupuje s mírným zpožděním, než postupuje tlaková vlna. Další obrázky 10 a 11 naznačují, jak dochází k úbytku metanu v závislosti na hoření plynné směsi, respektive identifikují místo, kde dochází k hoření. Tyto poslední obrázky jsou zaznamenány v časovém kroku 0,4 s a 0,5 s.

Obr. 7 Tlakové pole v čase 0,3 a 0,4 s [Pa]

18

PŘÍSPĚVKY


Obr. 8 Rychlostní pole v čase 0,3 a 0,4 s [m.s-1]

Obr. 9 Teplotní pole v čase 0,3 a 0,4 s [K]

19


PŘÍSPĚVKY

Obr. 10 Úbytek hmotnostního zlomku metanu vlivem hoření plynné směsi v čase 0,4 a 0,5 s

Obr. 11 Místo hoření plynné směsi v čase 0,4 a 0,5 s (Heat of reaction [W]) V této části byl vyhodnocen průběh změny tlaku v závislosti na čase v místech umístění tlakových senzorů. Jak už bylo uvedeno výše, umístění tlakových senzorů v experimentu se shoduje s umístěním v numerickém modelu. Pro lepší orientaci ve výsledcích je označení tlakových senzorů v matematickém modelu totožné jako u experimentu. Z grafu 3 je patrné, že bylo dosaženo velmi dobré shody experimentu s numerickým modelem, a to nejen u maximálních hodnot tlaku, ale také u průběhu 20

PŘÍSPĚVKY


tlakové vlny v závislosti na čase. Tlakový snímač s označením PS 2 (zelená křivka) dosahuje v matematickém modelu ve srovnání s experimentem velmi dobrého maximálního výbuchového tlaku. V experimentu dosahuje maximálních hodnot 1800 Pa a v numerickém modelu je dosaženo maximální hodnoty 1834 Pa. Průběh tlakové vlny v závislosti na čase lze považovat také za zdařilý, neboť maximálních hodnot je v experimentu dosaženo v čase 35,2 s a v modelu byl tento čas 35,01 s. Velmi dobře byl zachycen nejen začátek tlakové vlny, ale i její prudký pokles vlivem vysklení oken, vyražení dveří a poničení střechy v závěru vlny. Je potřeba také zmínit, že vznikly mírné nesrovnalosti, které byly zjištěny na špičce vlny, kdy v numerickém modelu nastal nejprve mírný pokles tlaku a až posléze bylo dosaženo maximální hodnoty výbuchového tlaku. Posledního snímače PS 3 (světle fialová křivka) nebylo v modelu uvažováno z důvodu, že se tlak v závislosti na čase mění jen velmi málo. Závěrem lze říci, že dosažené výsledky lze považovat za velmi uspokojivé.

Graf 3 Srovnání tlakového průběhu PS 2 (vlevo numerická simulace vpravo experiment) 4

ZÁVĚR

Článek řeší výbuch směsi metanu a vzduchu a šíření takto vzniklé tlakové vlny v rodinném domku s využitím numerické simulace. Matematický model domku byl vytvořen na základě experimentu týkajícího se výbuchu ve skutečném domě v Kamenné u Milína. Bylo také využito technické dokumentace a meteorologických dat k definici okrajových podmínek. 21


PŘÍSPĚVKY

V první části této práce je rozpracován podrobný matematický model celé úlohy, kdy je popsán nejprve matematický model proudění s přestupem tepla a ten je následně rozšířen o model řešící chemické reakce v plynech. V této kapitole jsou popsány základní rovnice, které jsou využity k řešení takového problému a popis fyzikálních vlastností objektu, vzduchu a paliva. Následuje oddíl zabývající se geometrií a tvorbou sítě v programu ANSYS FLUENT 13.0. Dále jsou specifikovány okrajové podmínky. Následně jsou určeny potřebné koeficienty k definici porézní vrstvy získané měřením. Závěr této sekce je zaměřen na definici materiálových vlastností. Při řešení této úlohy jsou uvažovány dvě varianty. Varianta A, která neuvažuje poničení domu a varianta B, která uvažuje s částečným poničením domu. Pro variantu A bylo provedeno vzájemné srovnání dvou softwarů (SW ANSYS FLUENT a SW FLACS). Bylo zjištěno, že oba tyto softwary dosahují velmi výrazné shody výbuchového tlaku v závislosti na čase. Ovšem je nutno podotknout, že se takto řešená úloha výrazně liší od experimentu. Maximální hodnoty tlaků v obou simulacích jsou shodné s experimentem, ale jejich průběh je výrazně rychlejší, než jak tomu bylo ve skutečnosti (asi 4x rychlejší). Naproti tomu varianta s porézní vrstvou dosahuje velmi dobré shody v průběhu tlakové vlny a v maximálních hodnotách tlaku v porovnání s experimentem. Je potřeba říct, že mírné odchylky výsledků mezi numerickými a experimentálními daty byly zjištěny. Jsou důsledkem špatného odhadu konstant porézní vrstvy, což má zásadní vliv na správnou realizaci numerické simulace. Hlavním přínosem práce je získání podrobných informací a dat o vývinu/šíření tlakových vln a jejich ničivých účinků v obytných objektech. Na základě porovnání vypočtených a naměřených hodnot sledovaných veličin v definovaných pozicích a čase lze takto získané informace dále využít a zobecnit pro potřebu technických expertiz, k záchraně lidských životů a zabránění škod na hmotném majetku. Résumé Article solves explosion of methane and air mixture and spread the resulting pressure waves in a house using numerical simulations. A mathematical model of house was created based on experiment concerning to the explosion in a real house in Kamenná in Milína. It was also used technical documentation and meteorological data to define the boundary conditions. In the first part of this thesis the detailed mathematical model of the entire job was elaborated. At first the mathematical model of flow with heat transfer was described and then was extended to model solving chemical reactions in gases. This chapter the basic equations are introduce used for solving the problem and description of the physical properties of an object, air and fuel. The following section deals with the geometry and mesh generation in ANSYS FLUENT 13.0. Further boundary conditions are specified. Subsequently, the coefficients are needed to define the porous layer and are obtained from 22

PŘÍSPĚVKY


measurements. The conclusion of this section is focused on the definition of material properties. In solving this task, two variants are considered. Variant A not considering damage of the house and variant B, which is considering a partial damage of house. For variant A the comparison between two software (ANSYS FLUENT software and FLACS) was done. It was found that both of these softwares reach very strong consensus of explosion pressure versus time. However, it should be noted that this solved task significantly differs from the experiment. Maximum values of pressures in both simulations are consistent with experiment, but their progress is much faster than it was in fact (about 4x faster). In contrast, a variant with a porous layer has a very good agreement in the course of pressure waves and maximum pressure values in comparison with experiment. It is necessary to say that a slight deviation between the numerical results and experimental data has been identified. They result from miscalculation of constants porous layer, which fundamentally affects the correct implementation of numerical simulation. The main contribution of this work is to obtain detailed information and data on the evolution / spread of pressure waves and their destructive effects in residential buildings. Based on the comparison of calculated and measured monitored variables depending on defined position and time thus obtained information can be used for technical expertise, to save lives and prevent damage to property. Poděkování Tento článek vznikl v rámci výzkumného projektu č. VD20062010A07 "Výzkum moderních metod pro ZPP a hodnocení nebezpečných účinků požárů na osoby, majetek a životní prostředí". Literatura [1] DVOŘÁK, O., DUDÁČEK, A. Zpráva o výsledcích požární zkoušky v rodinném domku v Kamenné u Milína dne 8. 10. 2010. Praha, Ostrava: TÚPO MV – GŘ HZS ČR a FBI VŠB – TU Ostrava, 2010. 30 s. [2] KOZUBKOVÁ, M. Numerické modelování proudění FLUENT I. [online]. Ostrava: VŠB – TU Ostrava, c2003.116 s. Poslední revize 6. 1. 2005. Dostupné z: . [3] PLATOŠ, P. Aplikace modelu v oblasti ekologie. [Disertační práce]. Ostrava: Katedra hydromechaniky a hydraulických zařízení, Fakulta strojní VŠB – Technická univerzita Ostrava, 2005. 118 s. [4] Ansys, Inc. ANSYS FLUENT 13.0 - Theory Guide. 2010. [5] BOJKO, M., KOZUBKOVÁ, M., MICHALEC, Z. Mathematical Model of the Low-Temperature Oxidation of Coal in Coal Stockpiles and Dumps. In Twenty – Seventh Annual International Pittsburgh Coal Conference. Istanbul (Turkey), 2010. 23


PŘÍSPĚVKY

[6] STULL, R. B. An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, 1994. s. 251-289. [7] KOZUBKOVÁ, M., BLEJCHAŘ, T., BOJKO, M. Modelování přenosu tepla a hybnosti. Ostrava: VŠB – TU Ostrava, 2011. 173 s. [8] BIRD, R. B., STEWART, W. E., LIGHTFOOT, N. N. Transport Phenomena. 2. vyd. John Wiley & Sons, 2002. 914 s. ISBN 0-471-41077-2. [9] RASZKA, J. Matematické modelování proudění přes porézní vrstvu v ANSYS FLUENT 13.0. Ostrava: VŠB – TU Ostrava, 2011. 22 s. [10] KOZUBKOVÁ, M., BOJKO, M., ZAVILA, O. Zpráva řešení modelování požáru daného tepelným výkonem a chemickou reakcí. Ostrava, 2009. 45 s.

24

MATEMATICKÉ MODELOVÁNÍ VÝBUCHU METANU V RODINNÉM DOMKU V KAMENNÉ POMOCÍ SW FLUENT

Recommend Documents