Vysoká škola chemicko-technologická v Praze Ústav plynárenství, koksochemie a ochrany ovzduší Technická 5, 166 28 Praha 6
Únik zemního plynu z potrubí a jeho následky při havárii na plynovodu Semestrální projekt
Vypracoval: Školitel:
Zdeňka Rohanová Doc. Ing. Václav Koza CSc. Praha, květen 2007
Souhrn Cílem tohoto projektu je matematické modelování výtoku plynu z vysokotlakého potrubí. Výpočtem rychlosti unikajícího plynu v kritické a podkritické oblasti a výpočtem hmotnostního toku, který je založen na numerickém řešení diferenciálních rovnic, lze následně zjistit hodnotu výtokového součinitele. Jeho hodnota vyjadřuje korekci na chování reálných tekutin oproti ideálnímu stavu při výtoku z otvoru. Měření bylo uskutečněno na několika modelových trhlinách různých tvarů. Pro omezení rizika při práci byl jako plyn použit vzduch. Data pro výpočet byla získána velice perspektivní metodikou CFD (computational fluid dynamics) a samotný výpočet byl pak proveden pomocí několika závislostí v programu MS Excel.
Obsah 1 2
Úvod ...................................................................................................................................2 Teoretická část....................................................................................................................3 2.1 Havárie .......................................................................................................................3 2.1.1 Příčiny havárií ....................................................................................................3 2.1.1.1 Vnější zásah....................................................................................................3 2.1.1.2 Koroze ............................................................................................................4 2.1.1.3 Ostatní příčiny ................................................................................................5 2.2 Zemní plyn..................................................................................................................5 2.2.1 Historie ...............................................................................................................5 2.2.2 Fyzikální vlastnosti.............................................................................................6 2.2.2.1 Univerzální plynová konstanta R ...................................................................6 2.2.2.2 Tlak p..............................................................................................................6 2.2.2.3 Teplota T ........................................................................................................7 2.2.2.4 Hmotnost neboli zádrž plynu m .....................................................................7 2.2.2.5 Hustota ρ.........................................................................................................7 2.2.2.6 Izoentropický koeficient κ ..............................................................................7 2.2.2.7 Střední molární hmotnost M..........................................................................8 2.3 Výtok plynu z otvoru v nádobě ..................................................................................8 2.3.1 Otvor v nádobě ...................................................................................................8 2.3.2 Kritický poměr tlaků...........................................................................................9 2.3.3 Rychlost výtoku plynu v podkritické oblasti......................................................9 2.3.4 Rychlost výtoku plynu v kritické oblasti..........................................................10 2.4 Hmotnostní tok .........................................................................................................11 •
2.4.1 Hmotnostní tok měřený m měě. ...........................................................................12 2.4.1.1 Ze vzorce pro výpočet derivace....................................................................12 2.4.1.2 Polynomická .................................................................................................12 2.4.1.3 Exponenciální ...............................................................................................13 •
2.4.2 Hmotnostní tok teoretický m teor . .......................................................................14 2.4.2.1 Hmotnostní tok (model 1) ............................................................................14 2.4.2.2 Hmotnostní tok (model 2) ............................................................................15 2.5 Hustota hmotnostního toku G a výtokový koeficient α............................................15 3 Cíle práce..........................................................................................................................16 4 Experimentální část ..........................................................................................................17 4.1 Experimentální zařízení............................................................................................17 ..................................................................................................................................................19 4.2 Princip měření ..........................................................................................................20 4.3 Příklad měření ..........................................................................................................20 4.4 Výsledky a diskuze...................................................................................................23 5 Závěr.................................................................................................................................27 Seznam symbolu.......................................................................................................................29 Seznam obrázku........................................................................................................................31 Seznam grafů ............................................................................................................................31 Seznam tabulek.........................................................................................................................31
1
1 Úvod Havárií na vysokotlakém plynovodu se rozumí vznik mimořádné situace, kdy může docházet ke vzniku trhlin v matriálu nebo protržení plynovodu ve velkém rozsahu a následnému úniku zemního plynu do zeminy obklopující potrubí nebo do atmosféry, kde může dojít ke vzniku požáru. Dochází tak mnohamilionovým škodám, zraněním či úmrtím osob a v neposlední řadě jsou tyto havárie nežádoucí také z hlediska ekologického. Významným ukazatelem přispívajícím k předcházení těchto havárií je matematické modelování průběhu a následků havarií. Klíčovým a zárověň problematickým faktorem se jeví, při vyhodnocování úniku tekutiny z plynovodů, zjištění vlastností otvoru, jeho geometrie, průřez otvoru (kruh, čtverec, nepravidlný tvar,…), ale i profil ve směru toku (ostrohranný otvor, zaoblení hran, atp.). Proto je důležité zjištění výtokového součinitele, který vyjadřuje korekci na chování reálné tekutiny oproti ideálnímu stavu při výtoku z otvoru (clona, otvor v nádobě). Ztráty mechanické energie jsou zapříčiněny tvarem otvoru, vlastnostmi tekutin a výtokovou rychlostí. Hodnota výtokového koeficientu se pohybuje od 0,3-0,8. V literatuře jsou hodnoty výtokového koeficientu uvedeny pouze pro pravidelné tvary (clona, dýza) a pro nepravidelné tvary, kterými jsou trhliny v potrubí, chybí. Zjištění výtokového koeficientu je velmi důležité pro přesný výpočet množství unikajícího plynu z narušeného plynovodního potrubí.
2
2 Teoretická část
2.1 Havárie Havárie je ve zprávách EGIG (European Gas pipeline Incident data Group) definována jako neúmyslný únik zemního plynu z pozemního přepravního potrubí o jmenovitém tlaku vyšším než 15 bar. V úvahu se berou jen úniky ze samotného potrubí, nikoli tedy z přidružených prvků jako armatur ani uvnitř oplocených instalací jako trasových uzávěrů, předávacích, kompresních nebo regulačních stanic. Havárie se dále ve zprávách EGIG dělí na třídy podle velikosti vzniklého otvoru: dírka/prasklina (pinhole/crack) - rozměr menší než 2 cm díra (hole)- od 2 cm do průměru potrubí roztržení (rupture) - rozměr větší než průměr potrubí. [1]
2.1.1 Příčiny havárií Příčiny vztažené na havárie v letech 1970-2001 jsou uvedeny v tabulce 1.
Tabulka 1. EGIG – Příčiny havárií Příčina
%
Vnější zásah
50
Konstrukční defekty/vada materiálu
17
Koroze (z toho 79 % způsobeno vnější korozí)
15
Pohyb zeminy
7
Potrubí omylem navrtané pod tlakem
5
Jiné
6
2.1.1.1 Vnější zásah Vnější zásah, který je způsoben zemními pracemi (kopání), těžkou technikou (pluh, buldozer, kotva). Závisí na průměru potrubí a druhu závady. Menší potrubí je nejsnadněji
3
poškozeno z důvodu větší aktivity pozemních prací, slabších stěn a nižší pevností třídou materiálu. U potrubí, které má tloušťku stěny větší než 15 mm nebyly havárie zjištěny.
2.1.1.2 Koroze Koroze se podle místa i příčin rozlišuje na vnitřní a vnější. Normálně končí velikostí dírky. U potrubí s tloušťkou stěny nad 15 mm nebyly pozorovány žádné korozivní havárie. [2]
2.1.1.2.1 Vnitřní koroze Vnitřní
koroze
vznikají
vlivem
korozních
složek
plynu
(
především
H 2 S , NH 3 , SO2 , HCN , CO 2 , O2 ). V plynné fázi nejsou tyto složky nebezpečné, ale jakmile v potrubí dojde ke kondenzaci vodní páry, dochází k jejich rozpuštění a vytváří tak vrstvovité nánosy v potrubí, až několik centimetrů vysoké ve formě usazenin nebo tvoří prach . Korozní zplodiny se skládají převážně z oxidu železitého,vody a sloučenin síry .
2.1.1.2.2 Vnější koroze Vnější koroze může být způsobena bludnými proudy nebo elektrochemickým vlivem okolní zeminy, někdy i za spoluúčasti některých druhů mikroorganismů (půdní koroze). Bludnými proudy označujeme elektrickou energii, která uniká ze svého el. obvodu (z kolejí,zpětných kabelů atp.) do okolní půdy. Tyto proudy jsou sbírány dobře vodivými konstrukcemi a jimi vedeny do míst nejpříhodnějších k návratu zpět do původního el. obvodu. Na potrubí rozlišujeme dvě základní charakteristické oblasti. Oblast katodická, což je oblast vstupu bludných proudů, a oblast anodická, kde bludné proudy z potrubí vystupují. Korozně nebezpečná je anodická oblast, která je často koncentrována na relativně velmi malé ploše potrubí (v místech poškození izolace). Půdní korozí rozumíme narušování kovových konstrukcí vlivem okolního prostředí. V podstatě jde o elektrochemickou reakci, při které je napadená část kovu anodou. Z korozního hlediska jsou nebezpečné půdy charakterizovány -
špatným provzdušněním, tzn., že jednotlivé části půdy jsou tak v těsné blízkosti, že zamezují přístup vzduchu do půdy,
-
vodí dobře elektrický proud,
-
obsahují velké množství vody nebo rozpustných látek, které dobře disociují, 4
-
střídáním různorodých půd podél trasy potrubí.
2.1.1.3 Ostatní příčiny Jsou například stavební závady, vady materiálu nebo pohyb zeminy (eroze, povodeň, sesuv půdy, dolování…) Jiné a neznámé příčiny, které lze rozdělit na řadu dalších předem definovaných příčin.
2.2 Zemní plyn Zemní plyn je směs plynných uhlovodíků a nehořlavých složek (zejména dusíku a oxidu uhličitého). Jeho charakteristickým znakem je vysoký obsah metanu. Zemní plyny typu H, které jsou využívány ve většině evropských zemích obsahují zpravidla více než 90% metanu a méně než 5% nehořlavých látek.
2.2.1 Historie První zmínka o zemním plynu se pravděpodobně týká “věčných ohňů“ na území dnešního Iráku, o kterých se zmiňuje Plutarchos ve 2. stol. n. l. Na své využití však musel čekat 1 500 let. Dříve, než ho člověk začal využívat jako prakticky hotový zdroj energie, používal svítiplyn, který se vyráběl z uhlí tzv. karbonizací. Svítiplyn se poprvé začal používat k veřejnému osvětlení ulic a domů v Anglii v 18. stol., protože neexistoval rozvod plynu do domácností. Roku 1821 v USA ve státě New York byl vyvrtán první vrt s úmyslem získat zemní plyn. Vrt byl hluboký pouze 3 m. Koncem 19. stol. se začala využívat elektřina a plynové osvětlení začalo být nahrazováno. V tomto století byl vynalezen kahan, jehož principem bylo mísení vzduchu a plynu před jeho hořením. Tento princip se stal základem všech plynových spotřebičů pro topení i vaření a podnítil rozvod a těžbu zemního plynu. První dálkový plynovod byl vybudován v roce 1891. Byl dlouhý75 km a vedl z plynových polí ve střední Indianě do Chicaga. Výroba dálkových plynovodů však byla náročná a tak do 40 let 20. stol. jich bylo vybudováno jen několik. Až pokrok ve svařování, výrobě oceli a trub v 50. a 60. letech 20. stol. umožnil mohutný rozvoj plynovodů. V České republice byla výroba svítiplynu ukončena v roce 1996 a od té doby je zde dodáván pouze zemní plyn. První vrty u nás byly provedeny v roce 1899 až 649 m, přesto se 5
však těžby nedosáhlo. Až v roce 1908 vrtem nedaleko Slavkova bylo nalezeno ložisko plynu. Česká republika však nemá velké zásoby zemního plynu, vlastní těžba nedosahuje ani 1% z celkové spotřeby. Proto se musí dovážet především z Ruska a od roku 1977 také z Norska. Než se plyn dostane ke spotřebiteli,musí projít systémem sběrných plynovodů do úpraven, kde se zbavuje nežádoucích látek a odtud se přepravuje prostřednictvím sítí do jednotlivých systémů. [5]
2.2.2 Fyzikální vlastnosti Při vyjadřování s dopravou, kompresí, regulací a měřením množství plynu se opíráme o stanovení jejich fyzikálního stavu, především tlaku, teploty, měrné hmotnosti a o závislosti jejich změn daných stavovými rovnicemi ideálního a reálného plynu.
2.2.2.1 Univerzální plynová konstanta R Rovnici spojeného zákona Boyle-Gay-Lussacova je možno také psát ve tvaru: p n ⋅ Vmn = konst. = R = 8,314 To
p n KK je normální tlak 101 325 Pa Vmn KK = 22,4 dm^3/mol, normální molový objem 1 molu plynu při T0 KK = 273.15 K Základní jednotkou univerzální konstanty J/mol.K
2.2.2.2 Tlak p Vzniká nárazy molekul na stěnu nádoby a vyjadřuje poměr velikosti síly F, působící kolmo na rovinnou plochu a rovnoměrně spojitě rozloženou po této ploše a obsahu této plochy S, tedy:
p=
F S
Hlavní jednotkou v SI soustavě je pascal (Pa) a je to tlak, který vyvolává sílu 1 newtonu působící kolmo na plochu o obsahu 1 m2 .Tlak plynu v základních vztazích se udává převážně 6
v absolutní hodnotě, která se rovná součtu (rozdílu) tlaku atmosférického a změřeného přetlaku (podtlaku).
2.2.2.3 Teplota T Je to veličina vyjadřující míru tepelného stavu daného tělesa. V termodynamice se převážně užívá absolutních teplot udávaných v Kelvinově stupnici. 0 stupňů v Celsiově stupnici je rovna 273,15 K.
2.2.2.4 Hmotnost neboli zádrž plynu m Zádrž plynu je vyjádřena ze stavové rovnice ideálního plynu jako:
m=
p ⋅V ⋅ M R ⋅T
Hlavní jednotkou hmotnosti je kg.
2.2.2.5 Hustota ρ Je veličina, která vyjadřuje hmotnost objemové jednotky látky. Základní jednotkou je kg/m^3 a lze vypočítatat ze vztahu :
ρ=
m V
ρ=
p⋅M R.T
, nebo se dá také vypočítat ze stavové rovnice id.plynu:
2.2.2.6 Izoentropický koeficient κ Také nazývaný Poissonovou konstantou definovanou jako:
κ=
cp
cV
7
Je to veličina bezrozměrná a je závislá na teplotě a tlaku. Pro náš výpočet jsme si to zjednodušili a brali jsme tento koeficient konstantní, tedy nezávislý na tlaku a teplotě. κ (zemní plyn) = 1,3 κ (vzduch) = 1,4
2.2.2.7 Střední molární hmotnost M M vz = 28 ⋅ 10 −3 K[kg mol ] M zp = 16,04 ⋅ 10 −3 K[kg mol ]
2.3 Výtok plynu z otvoru v nádobě
2.3.1 Otvor v nádobě Nádoba je třírozměrný útvar, kde rozměr v žádném směru nepřevládá a plyn vytékající z nádoby se k otvoru jakoby sbíhá ze všech stran. Otvor v nádobě je charakterizován svým průřezem – velikostí plochy otvoru a součinitelem zúžení proudu α, který závisí na tvaru otvoru. Schématické znázornění výtoku plynu ukazuje obr. 1.
Obr. 1 Otvor v nádobě
8
2.3.2 Kritický poměr tlaků Je to poměr tlaku, který vytéká z nádoby do okolí, k atmosférickému tlaku. Výtoková rychlost je závislá na tomto poměru a roste pouze tak dlouho, dokud tento poměr nedosáhne tzv. kritické hodnoty.Tento poměr je definován veličinou β: κ
p κ + 1 κ −1 β = 0 = pa 2 Po dosazení κ = 1,3 získáme pro zemní plyn hodnotu tohoto poměru 1,8357.
2.3.3 Rychlost výtoku plynu v podkritické oblasti Při výtoku plynu z nádoby musí platit zákon zachování energie, který bývá vyjádřen Bernoulliho rovnicí,říkající, že v každém okamžiku musí být součet kinetické a potenciální energie konstantní. Potenciální energie je dána součtem tlakové a gravitační energie a platí [3]: p
v2 dp + ∫ + g ⋅ h = konst. 2 p0 ρ - g.h je gravitační potenciální energie,která převažuje v případě kapaliny, avšak je zcela zanedbatelná pro plyn, jehož potenciální energie je skryta především v jeho stlačeném objemu. Předpokládáme-li, že plyn vytéká z otvoru z místa, v němž je v klidu, pak můžeme konstantní hodnotu na pravé straně rovnice položit rovno nule a dostáváme: p
dp v2 = −∫ 2 p0 ρ Dále předpokládáme-li, že expanze plynu probíhá adiabaticky,tak z rovnice adiabaty pro změnu stavu plynu 9
ρ0
ρ =
p0
1κ
p
=
T0
1 dp dp = ∫p 0 ρ ∫ ρ 0 ⋅ ( p p )κ 0 p
1 (κ −1)
T
− κ p0 = ⋅ κ −1 ρ 0
κ −1 p κ ⋅ 1 − p 0
a dosazením za hustotu
ρ0 =
p0 ⋅ M R ⋅ T0
dostaneme pro podkritickou rychlost plynu
v podkrit .
κ −1 2 ⋅ κ R ⋅ T0 p κ = ⋅ ⋅ 1− κ − 1 M p 0
2.3.4 Rychlost výtoku plynu v kritické oblasti Rychlost výtoku plynu v kritické oblasti nemůže být větší než lokální rychlost zvuku (lokální, jelikož je závislá na teplotě). Přitom pro rychlost zvuku platí vztah [3]:
c=
dp dρ
ze kterého opět dosazením hustoty z rovnice adiabaty obdržíme přes
dp =
ρ p 0 ⋅ κ ⋅ ρ0
κ −1
⋅ dρ
ρ0
10
a dosazením za dp, získáme pro rychlost zvuku, která je rovna kritické rychlosti výtoku vztah:
c2 = v2 =
p p 0 ⋅ κ ⋅ p 0
κ −1 κ
ρ0
kde dále dosazením kritického poměru tlaků a vyjádřením hustoty ze stavové rovnice získáváme vzorec pro výpočet kritické rychlosti
v krit . = 2 ⋅
R ⋅ T0 κ +1 M
κ
⋅
2.4 Hmotnostní tok Hmotnostní tok definujeme jako úbytek plynu z nádoby v závislosti na čase. Tok se postupně zmenšuje, protože dochází k úbytku tlaku i množství zádrže. Lze ho vypočíst dle vztahu [1]: •
m = v ⋅ S ⋅ ρ ⋅α
Jelikož neznáme velikost součinitele zúžení proudu α , který je předmětem našeho zkoumání, zjišťujeme hmotnostní tok jako derivace hmotnosti za čas •
m=
dm dτ
a to několika způsoby.
11
•
2.4.1 Hmotnostní tok měřený m měě. Určuje se ze závislosti hmotnosti na čase proložením dat vhodnou funkcí a následným odvozením vztahu pro výpočet této derivace.
2.4.1.1 Ze vzorce pro výpočet derivace Jako nejjednodušší cesta k určení derivace se může jevit rovnice:
•
m1mëreno =
m j − m j +k
j = 0, 1, 2, …
τ j − τ j +k
k = j+1, j+2, j+3, …
2.4.1.2 Polynomická Naměřená data byla aproximována polynomickou rovnicí pátého stupně m = a ⋅τ 5 + b ⋅τ 4 + c ⋅τ 3 + d ⋅τ 2 + e ⋅τ + f
kde a, b, c, d, e, f -
parametry aproximační funkce
Derivace hmotnosti za čas odvozená z této rovnice je pak: •
m 2 mereno = 5 ⋅ a ⋅τ 4 + 4 ⋅ b ⋅τ 3 + 3 ⋅ c ⋅τ 2 + 2 ⋅ d ⋅τ + e
V následujícím grafu je pro ukázku zobrazeno proložení naměřených dat polynomickou funkcí.
12
Graf 1 Polynomická závislost hmotnosti na čase Závislost hmotnosti na čase Závislost hmotnosti na čase Polynomický (Závislost hmotnosti na čase)
0,40 y = -8E-06x 5 + 0,0004x 4 - 0,006x 3 + 0,0537x 2 - 0,2899x + 0,9711 R2 = 0,9998
0,35
m (kg)
0,30 0,25 0,20 0,15 0,10 0,05 0,00 3
5
7
9
11
13
t (s)
2.4.1.3 Exponenciální Naměřená data byla po celé délce rozdělena do několika úseků a následně aproximována sérií exponenciál ve tvaru:
m = a ⋅ e −b⋅τ
kde a, b
-
jsou parametry aproximační funkce
Derivace hmotnosti za čas odvozená z této rovnice je pak: •
m 3mereno = −b ⋅ a ⋅ e −b⋅τ
V následujícím grafu je pro ukázku zobrazeno proložení naměřených dat serií exponenciál.
13
Graf 2 Exponenciální závislost hmotnosti na čase Závislost hmotnosti na čase
y = 0,2591e-0,0519x R2 = 0,9996
hmotnost (kg)
0,370
kriticka 1 y = 0,4012e-0,1785x R2 = 0,9989
podkriticka 1
y = 0,3878e-0,1451x R2 = 0,9992
podkriticka 3
0,320
y = 0,309e-0,0828x R2 = 0,9963
y = 0,3771e-0,1327x R2 = 0,998
0,270
y = 0,2307e R2 = 0,9997
y = 0,3602e R2 = 0,9983
-0,0347x
-0,1171x
podkriticka 2 podkriticka 4 podkriticka 5 podkriticka 6 podkriticka 7 podkriticka 8 podkriticka 9 Exponenciální (podkriticka 2)
y = 0,3434e-0,1048x R2 = 0,9991
0,220
Exponenciální (podkriticka 3) Exponenciální (podkriticka 4) Exponenciální (podkriticka 5) Exponenciální (kriticka 1)
0,170 0,0
2,0
4,0
6,0
8,0
čas (s)
Exponenciální (podkriticka 6) 10,0 Exponenciální (podkriticka 7) Exponenciální (podkriticka 1)
•
2.4.2 Hmotnostní tok teoretický m teor .
2.4.2.1 Hmotnostní tok (model 1) Tento vztah pro hmotnostní průtok vychází z rovnice kontinuity a kritického průtoku [1]:
•
m1teor .
κ +1
2 κ −1 = α ⋅ S ⋅ Ψ0 ⋅ pi ⋅ ρ i ⋅ κ ⋅ κ +1
Koeficient Ψ0 je roven 1 při kritickém toku. Při podkritickém proudění, kdy dochází k poklesu výtokové rychlosti z maximální hodnoty pak nabývá tvaru: 14
2
κ +1
2 κ + 1 κ −1 p a κ Ψ0 = ⋅ ⋅ κ −1 2 pi
κ −1 κ p a ⋅ 1 − pi
2.4.2.2 Hmotnostní tok (model 2) Nejprve získáme časově závislou proměnnou τr z naměřených hodnot tlaku [1]: 2⋅κ
1 κ −1 pi = p 0 ⋅ 1 + ⋅ (κ − 1) ⋅ τ r 2
Zde je průtok vyjadřující časovou závislost:
1 = α ⋅ S ⋅ β ⋅ pi ⋅ ρ i ⋅ 1 + ⋅ (κ − 1) ⋅τ r 2
•
m 2teor
−
κ +1 κ −1
Při podkritickém toku je potřeba do výpočtů zahrnout Ψ0 podle předešlého vztahu.
2.5 Hustota hmotnostního toku G a výtokový koeficient α Hustotu hmotnostního toku získáme poměrem hmotnostního toku a průřezem otvoru •
m G= S
[kg / s ⋅ m ] 2
Koeficient α získáme poměrem Gmereno a Gteor . •
α=
Gmereno = Gteor .
m mereno
.
•
m teor .
15
3 Cíle práce Cílem projektu bylo sestavení experimentálního zařízení pro zjištění výtokového koeficientu. Hodnota výtokového koeficientu je důležitá pro zjištění výtokového koeficientu při havárii plynovodního potrubí. Výtokový koeficient vyjadřuje korekci na chování reálné tekutiny oproti ideálnímu stavu při výtoku z otvoru (clona, otvor v nádobě). Ztráty mechanické energie jsou zapříčiněny tvarem otvoru, vlastnostmi tekutin a výtokovou rychlostí. Hodnota výtokového koeficientu se pohybuje od 0,3-0,8. V literatuře jsou hodnoty výtokového koeficientu pouze pro pravidelné tvary (clona, dýza) a pro nepravidelné tvary, kterými jsou trhliny v potrubí, chybí. Zjištění výtokového koeficientu je velmi důležité pro přesný výpočet množství unikajícího plynu z narušeného plynovodního potrubí. Trhliny v potrubí mohou být různě velké: od velikosti špendlíkové hlavičky až po destrukci velké části potrubí například rozevřením šroubovicového svaru.
16
4 Experimentální část 4.1 Experimentální zařízení Schéma experimentálního zařízení vidíme na obrázku 2, fotografie zařízení na obrázku 3 a na obrázku 5 jsou ukázány modelové trhliny.
Obr. 2 Schéma aparatury
Zařízení sestává z následujích součástí: 1- Bezolejový kompresor Gude. Přípojka motoru 230 V, výkon 1.1 kW, sací výkon max. 220 l/min, skutečné dodávané množství 125 l/min. Max tlak 8 barů. Kompresor je hadicí s rychlouzávěrem připojen na plastové potrubí 1/2´´. 2- Kulový kohout, připojení 1/2´´. 3- Zaslepený T-kus pro případné napojení dalších komponent 4- Kulový kohout, připojení 1/2´´. 5- Plnozdvižný pojistný ventil nastavený na tlak 10 barů. 6- Tlaková nádoba.
17
Jako tlakovou nádobu jsme použili standardní expanzní nádobu Aquamat vyráběnou v Dukle Trutnov. AQUAMAT P je svařená ocelová tlaková nádoba, stojatá s připojovacími nátrubky. Celá nádoba je žárově zinkována (vně i uvnitř). Tlaková nádoba je konstruována na přetlak 10 bar, objem nádoby je 150 l. Nádoba má 4 otvory G 1´´ - vstup, výstup a dva otvory na snímače tlaku a teploty. 7 – Snímač teploty. Snímač Pt 100, hlavicový typ. Sériově vyráběný snímač byl výrobcem upraven – stonek z nerez oceli byl otevřen kvůli rychlejší odezvě při měření. Snímač teploty má teplotní rozsah od -40 - 100°C. Na snímač teploty je napojen proudový převodník PTC/I, který převádí změny odporu Pt 100 na proudový signál 4 ÷ 20 mA. K napájení snímače byl použit napájecí zdroj ZS 24 ZPA Ekoreg. 8- Pro snímání tlaku byl použit průmyslový snímač tlaku DMP 331. Převádí tlak plynů na elektrický signál na proudový signál 4 ÷ 20 A. 9- Manometr. 10- Příruba. Polyethylenová příruba G 1´´, která je zajištěna čtyřmi šrouby. Obě části příruby je možné od sebe snadno oddálit pomocí kladkostroje viz obr. č 3. Mezi obě části příruby jsou vkládány vzorky – ocelové plechy s otvory simulujícími narušení plynovodního potrubí. Na následujících obrázcích je vidět několik vzorků. 11 – 2/2 cestný magnetický ventil s membránovým uzávěrem, nepřímo řízený. Ventil je možné otevřít naplno za velmi krátký časový okamžik. Elektrické výstupy ze snímačů teploty a tlaku jsou vedeny do svorkovnice a dále do měřící karty v řídící jednotce. [4]
18
Obr. 3 Aparatura a detail aparatury
Obr. 2 Fotografie modelových trhlin
Vzorek č. 3
Vzorek č. 4
19
4.2 Princip měření Do příruby 11 je vsunut vzorek měřeného otvoru a magnetický ventil 10 nastaven na zavřeno. Nádoba o objemu V je natlakována pomocí kompresoru na max. 8 barů (maximální možný tlak kompresoru). Kulový kohout 4 je uzavřen a tím přerušeno tlakování nádoby. Pomocí manometru je zaznamenáván tlak a pomocí teplotního čidla měřena aktuální teplota se vzorkovací periodou až 10 ms (tj. 100x za sekundu). Po natlakování nádoby se magnetický ventil 11 otevře a z nádoby začne unikat plyn přes vzorek s otvorem. Doba vyprazdňování závisí na výchozím tlaku v nádobě a na velikosti a tvaru otvoru, který je charakterizován výtokovým koeficientem.
4.3 Příklad měření K výpočtům byl použit Microsoft Excel. V následující tabulce č. 2 jsou počáteční podmínky pro měření.
Tabulka 2 Počáteční podmínky pro měření přetlak na poč.
pi0
448 515
[Pa]
absolutní tlak
Pi0
549 840
[Pa]
průřez otvoru
S
0,00020
[m^2]
izoentropický koeficient
κ
1,4
[-]
molární hmotnost
M
0,02880
[kg/mol]
objem nádoby
V
0,15
[m^3]
Z naměřených dat vytvoříme závislost absolutního tlaku Pabsolutní na čase τ
20
Graf 3 Závislost tlaku na čase
Závislost tlaku na čase zavislost tlaku na čase
550 000 500 000 450 000 tlak (Pa)
400 000 350 000 300 000 250 000 200 000 150 000 100 000 0
10
čas (s)
20
30
40
50
Po spočtení zádrže mi trhlinou o průřezu S vytvoříme graf závislosti zádrže na čase a aproximujeme vhodnou rovnicí.
21
Graf 4 Závislost hmotnosti na čase
Závislost hmotnosti na čase Závislost hmotnosti na čase
0,90 0,80 0,70 m (kg)
0,60 0,50 0,40 0,30 0,20 0,10 0,00 2
7
t (s)
12
17
Získáme tak hmotnostní průtoky : •
m num ,mer - hmotnostní průtok získaný ze základní rovnice pro derivaci (viz. 2.4.1.1) •
m pol ,mer - hmotnostní průtok získaný z derivace polynomu (viz. 2.4.1.2) •
m exp,mer - hmotnostní průtok získaný z derivace exponenciál (viz. 2.4.1.3)
Teoretické hmotnostní průtoky získáme ze vzorců z Yellow book – model 1 a model 2 •
•
m1teor . a m 2teor . .
Hustoty hmotnostního toku získáme podělením hmotnostního průtoku průřezem trhliny: •
Gi ,mereno
m i ,mereno = S •
Gi ,teor .
m = i ,teor S
22
Výtokové součinitele pak získáme ze vzorce: Gi , mereno Gi ,teor
α= kde:
α num,1=
Gnum ,mer
α num, 2 =
G1teor
Gnum ,mer
α pol ,1 =
G pol ,mer
α pol , 2 =
G pol ,mer
α exp,1 =
Gexp,mer
α exp, 2 =
Gexp,mer
G2teor .
G1teor G2teor G1teor . G2teor .
4.4 Výsledky a diskuze V této části jsou shrnuty výsledky výpočtů, které byly získány z naměřených dat a dále zpracovány. Měření bylo prováděno při několika různých přetlacích ( až 5 barů) a přes různé typy trhlin. Pro ukázku jsem vybrala měření při přetlaku 5 barů a pro vzorek č.3.
23
Graf 5 Závislost teoretické hustoty hmotnostního toku na tlaku.
Závislost teoretické hustoty hmotnostního toku na tlaku G 1teor. G 2teor 1400,0
G (kg/s.m^2)
1200,0 1000,0 800,0 600,0 400,0 200,0 0,0 0
100 000
200 000
300 000
400 000
500 000
600 000
P (Pa)
Graf 6 Závislost měřené hustoty hmotnostního toku na tlaku.
Závislost měřené hustoty hmotnostního toku na tlaku 1400,0 1200,0 1000,0
G (kg/sm^2)
800,0 G 1mereno 600,0
G 2mereno
400,0
G 3mereno
200,0 0,0 -200,0
0
100 000
200 000
300 000
400 000
500 000
600 000
-400,0 P (Pa)
Vypočtené hodnoty výtokového součinitele jsou zobrazeny v následujících grafech. Je z nich patrné, že u výpočtu numerické alfy je zapotřebí vyhlazení dat a u polynomické alfy je 24
zřejmé, že v krajních hodnotách prudce vzrůstá. Také u těchto dvou typů výpočtů je hodnota alfy větší než jedna, což by být nemělo.
Graf 7 Závislost
α num,1 a α num, 2
na tlaku.
Závislost alfa num. na tlaku alfa num1 alfa num2 4,00 3,00
alfa [-]
2,00 1,00 0,00 -1,00
0
100 000
200 000
300 000
400 000
500 000
600 000
-2,00 -3,00 tlak [Pa]
25
Graf 8 Závislost
α pol ,1 a α pol , 2
na tlaku.
Závislost alfa pol. na tlaku alfa pol,1 alfa pol,2 4 3,5
alfa [-]
3 2,5 2 1,5 1 0,5 0 0
100 000
200 000
300 000
400 000
500 000
600 000
tlak [Pa
Graf 9 Závislost
α exp,1 a α exp, 2
na tlaku
Závislost alfa exp. na tlaku alfa 31 alfa 32 0,90 0,80 0,70 alfa [-]
0,60 0,50 0,40 0,30 0,20 0,10 0,00 0
100 000
200 000
300 000
400 000
500 000
600 000
tlak [Pa]
26
5 Závěr Zabývali jsme se modelováním výtoku plynu z tlakové nádoby. Získali jsme sérii časové závislých hodnot tlaku plynu unikajícího z různých typů trhlin. Tato data byla následně použita pro výpočet výtokových součinitelů modelového plynu. Jejich hodnoty lze v praxi použít pro výpočet hmotnosti unikajícího plynu při havárii vysokotlakých zařízení. Zvolené matematické zpracování naměřených dat úspěšně popisuje většinu studovaného děje, avšak při nízkých hodnotách tlaku v potrubí (konec děje) dochází k větším odchylkám od teorie. Proto je třeba problém nadále studovat a upřesnit použitý matematický model.
27
Seznam použité literatury 1. Pasková I., „Únik zemního plynu a jeho následky při havárii na plynovodu “,VŠCHT, Praha, (2005). 2. Potužák K., „Dálková doprava a rozvod plynu “, SNTL, Praha, (1981) 3. Kalčík J., „ Technická termodynamika “, NČAV, Praha, (1963) 4. Koza V., Pasková I., „ Kolokvium 2007”, Praha, (2007) 5. http://www.quido.cz/objevy/zemniplyn.htm
28
Seznam symbolu Proměnné
hustota hmotnostního toku
kg ⋅ m 2 ⋅ s −1
m
hmotnostní tok
kg ⋅ s −1
m
hmotnost
kg
M
molární hmotnost
kg ⋅ mol −1
P
absolutní tlak
Pa
p
přetlak
Pa
R
univerzální plynová konstanta
J ⋅ mol −1 ⋅ K −1
S
průřez otvoru
m2
T
absolutní teplota
K
V
objem
m3
α
výtokový koeficient
-
β
kritický tlakový poměr
-
κ
izoentropický koeficient
-
τ
čas
s
v
rychlost
m ⋅ s −1
ρ
hustota
kg ⋅ m −3
F
síla
N
EIGIG
European gas popelíne incident data group
G •
29
Indexy
0
počátek v čase
i
nátok do únikového otvoru
a
atmosféra
krit.
kritická
podkrit.
podkritická
mereno
měřená hodnota
teor.
Teoretická hodnota
num
numerická derivace
exp
derivace zexponenciál
pol
derivace polynomu
30
Seznam obrázku Obr. 1 Otvor v nádobě ................................................................................................................8 Obr. 2 Schéma aparatury ..........................................................................................................18 Obr. 3 Detail aparatury .............................................................................................................20 Obr. 4 Fotografie modelových trhlin........................................................................................20
Seznam grafů Graf 1 Polynomická závislost hmotnosti na čase .....................................................................13 Graf 2 Exponenciální závislost hmotnosti na čase ...................................................................14 Graf 4 Závislost tlaku na čase ..................................................................................................21 Graf 5 Závislost hmotnosti na čase ..........................................................................................22 Graf 6 Závislost teoretické hustoty hmotnostního toku na tlaku..............................................24 Graf 7 Závislost měřené hustoty hmotnostního toku na tlaku..................................................24 Graf 8 Závislost α num,1 a α num, 2 na tlaku. .................................................................................25 Graf 9 Závislost α pol ,1 a α pol , 2 na tlaku....................................................................................26 Graf 10 Závislost α exp,1 a α exp, 2 na tlaku ..................................................................................26
Seznam tabulek Tabulka 1. EGIG – Příčiny havárií.............................................................................................3 Tabulka 2 Počáteční podmínky pro měření..............................................................................20
31
