NUMERICKÝ MODEL NESTACIONÁRNÍHO PŘENOSU TEPLA V PALIVOVÉ TYČI JADERNÉHO REAKTORU VVER 1000 SVOČ – FST 2014 Miroslav Kabát, Západočeská univerzita v Plzni, Univerzitní 8, 306 14 Plzeň Česká republika ABSTRAKT Práce pojednává o návrhu, fungování a programování vlastního numerického modelu sloužícího k simulacím nestacionárního přestupu tepla v dvou-dimenzionálních rotačně symetrických geometriích s proměnlivými veličinami v radiálním a axiálním směru s uvažováním proudění tenkých vrstev tekutin. Numerický model vychází z diskretizace energetické rovnice na rovnici Fourier-Kirchhoffovu, která je upravena na explicitní bilanční výpočet ve výpočetní síti aplikované na geometrii palivové tyče jaderného reaktoru VVER 1000. KLÍČOVÁ SLOVA numerický model, přenos tepla, simulace, jaderný reaktor, palivová tyč ÚVOD V tomto příspěvku naleznete informace o principech fungování jednotlivých částí velmi kompaktního zdrojového kódu programu numerického modelu pro SW Matlab sloužícího k simulacím nestacionárního přestupu tepla v dvou-dimenzionálních rotačně symetrických geometriích s proměnlivými veličinami v radiálním a axiálním směru s uvažováním proudění tenkých vrstev tekutin. Numerický model vychází z diskretizace energetické rovnice na rovnici Fourier-Kirchhoffovu, která je upravena na explicitní bilanční výpočet ve výpočetní síti aplikované na geometrii palivové tyče jaderného reaktoru VVER 1000. Dle charakteru tohoto numerického modelu jej lze použít univerzálně na obdobné geometrie jako je například izolované potrubí či trubkové výměníky tepla. ZADÁNÍ Tento projekt vznikl jako semestrální práce k předmětu Termodynamika jaderného reaktoru v ZS 2013/14. Zadání: "Vytvořte numerický model pro simulaci přestupu tepla palivovým článkem reaktoru VVER1000." ŘEŠENÍ Pro vyřešení zadané úlohy byl využit SW MATLAB od společnosti MathWorks. Numerický model byl vytvořen nestacionární dvojdimenzionální explicitní bilanční metodou v rovině dané v axiálním směru výškou a v radiálním směru poloměrem. Samotný zdrojový kód algoritmu byl rozdělen do dílčích úseků:
Geometrie Síťování Výpočet o Počáteční podmínky o Vlastní výpočet Zobrazení výsledků
Celý výpočet funguje na principu vytvoření trojrozměrné matice "Z", kde dva směry udávají polohu výpočetních buněk (radiální a axiální směr) a třetí směr slouží k ukládání informací o buňce. Tyto informace jsou následně použity k naplnění série rovnic pro řešení přenosu tepla mezi jednotlivými buňkami. Na obr. 1 a obr.2 je velmi dobře patrný axiální směr (rovnoběžný se sloupci) a radiální směr (rovnoběžný s řádky). Informace o buňkách jsou pak v jednotlivých vrstvách nad vrstvou výpočetní sítě.
GEOMETRIE [2] Geometrie představuje složený plný válec modelovaný od osy symetrie po střední hodnotu poloměru vnější vrstvy tvořené chladivem po celé výšce palivového článku. Geometrie má 4 výpočtové oblasti: FUEL – palivová tableta 𝑈𝑂2 - 𝑟𝜖〈0; 3,775〉 𝑚𝑚 GAP – inertní vrstva 𝐻𝑒 - 𝑟𝜖〈3,775; 3,9〉 𝑚𝑚 CLAD – krycí vrstva 𝑍𝑟 − 𝑁𝑏 - 𝑟𝜖〈3,9; 4,55〉 𝑚𝑚 OUT – chladivo 𝐻2 𝑂 - 𝑟𝜖〈4,55; 6,735〉 𝑚𝑚 Každá z těchto částí je vysoká - ℎ𝜖〈0; 2960〉 𝑚 Na každou palivovou tyč připadá vrstva chladiva ve tvaru šestiúhelníku. Jelikož je výpočet provádět dvojdimenzionálně, byla geometrie zjednodušena středním poloměrem tohoto šestiúhelníku a považována za válec. SÍŤOVÁNÍ Všechny výpočtové oblasti geometrie byly rozděleny na 5 buněk v axiálním směru a na rozdílný počet buněk v radiálním směru (FUEL 10 buněk, GAP 1 buňka; CLAD 2 buňky; OUT 1 buňka) (viz. obr. 1) Jednotlivý počet buněk pro každou výpočtovou část byl zvolen z důvodu co nejvyšší tepelné kapacity každé buňky s ohledem na tepelný tok, který každou buňkou prochází. Velmi jemná síť vyžaduje velmi malý časový krok, který neúměrně prodlužuje čas výpočtu. K takto zvolené výpočetní síti byly dodány buňky navíc z každé její strany. Tyto buňky reprezentují okrajovou podmínku a zjednodušují tak výpočetní algoritmus. Každou výpočetní část i okrajovou podmínku reprezentuje index buňky, který je uložen do první vrstvy třírozměrné matice 𝑍. Index 1 = FUEL Index 2 = GAP Index 3 = CLAD Index 4 = OUT
Index 6 = okrajová podmínka „dno“ Index 7 = okrajová podmínka „poklop“ Index 8 = okrajová podmínka „povrch“ Index 9 = okrajová podmínka „osa“
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
9 9 9 9
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
2 2 2 2
3 3 3 3
3 3 3 3
4 4 4 4
8 8 8 8
9
1
1
1
1
1
1
1
1
1
1
2
3
3
4
8
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
Obrázek 1: Indexy výpočtových oblastí ve výpočetní síti (1. vrstva matice Z)
VÝPOČET Počáteční podmínky Do vrstev 2 a 4-12 matice 𝑍 byly vepsány počáteční podmínky výpočtu. Vrstva 2 reprezentuje počáteční teplotu každé buňky [°𝐶]. Další vrstvy reprezentují vlastnosti potřebné pro další výpočty: Vrstva 2 – počáteční teplota buňky [°𝐶] Vrstva 4 – tloušťka buňky [𝑚] Vrstva 5 – střední poloměr buňky [𝑚] Vrstva 6 – minimální poloměr buňky [𝑚] Vrstva 7 – maximální poloměr buňky [𝑚] Vrstva 12 – objem buňky [𝑚3 ] Vlastní výpočet [1] Výpočet vyhodnocuje nestacionární vlastnosti buněk v simulovaném čase po časových krocích. Simulovaný čas byl zvolen 600 [𝑠] . Délka časového kroku pak 0,0001 [𝑠]. Celkem tedy proběhne 6 000 000 iterací výpočtu, což trvá přibližně 15 minut. Každá iterace přepočítá fyzikální vlastnosti každé buňky závislé především na teplotě a tlaku. Tlak je považován za konstantní (15,7 𝑀𝑃𝑎). Hodnoty fyzikálních vlastností buněk jsou zapisovány do vrstev trojrozměrné matice. Vrstva 8 – hustota [𝑘𝑔/𝑚3 ] Vrstva 9 – měrná tepelná kapacita [𝐽/𝑘𝑔. 𝐾] Vrstva 10 – součinitel tepelné vodivost [𝑊/𝑚. 𝐾] a součinitel přestupu tepla [𝑊/𝑚2 . 𝐾] Vrstva 11 – vnitřní zdroj tepla [𝑊/𝑚3 ] Podstata explicitní bilanční metody spočívá ve vytvoření tepelné bilance každé buňky z vlastností, které o nejbližších sousedních buňkách a samotné počítané buňce známe. Tepelná bilance je vyhodnocena z 6 zdrojů tepla a informativně zapsána zbylých vrstev matice 𝑍: N - tepelný tok v axiálním směru shora dolů [𝑊] (vrstva 13) S – tepelný tok v axiálním směru zdola nahoru [𝑊] (vrstva 15) E – tepelný tok v radiálním směru zprava doleva [𝑊] (vrstva 14) W – tepelný tok v radiálním směru zleva doprava [𝑊] (vrstva 16) Q – vnitřní zdroj tepla [𝑊] (vrstva 17) M – tepelný tok míšením [𝑊] (vrstva 18) Pro tepelnou bilanci N a S byla upravena rovnice pro prostup tepla složenou rovinnou stěnou. 𝑆∙∆𝑇
(𝑁; 𝑆) = 𝑘 ∙ 𝑆 ∙ ∆𝑇 = 𝛿1
𝛿 + 2 𝜆1 𝜆2
[𝑊]
(1)
Pro tepelnou bilanci E a W byla upravena rovnice pro prostup tepla složenou válcovou stěnou. (𝐸; 𝑊) = (𝐸) =
2∙𝜋∙𝑙∙∆𝑇 1 𝑟 1 𝑟 ln 2 + ln 3 𝜆1 𝑟1 𝜆2 𝑟2
2∙𝜋∙𝑙∙∆𝑇 1 𝑟 1 ln 2 + 𝜆1 𝑟1 𝛼∙𝑟2
[𝑊]
[𝑊]
(2) (3)
Pro tepelnou bilanci Q byla využita hodnota vnitřního zdroje. 𝑄 = 𝑄 [𝑊]
(4)
Pro tepelnou bilanci míšením bylo využito rovnice míšení dvou rozdílných látek. (5)
𝑀 = 𝑆 ∙ 𝑤 ∙ (𝜌1 ∙ 𝑐𝑝1 ∙ 𝑇1 − 𝜌2 ∙ 𝑐𝑝2 ∙ 𝑇2 ) [𝑊] Z tohoto pak plyne celková tepelná bilance.
(6)
𝑄𝑏𝑖𝑙 = (𝑁 + 𝐸 + 𝑆 + 𝑊 + 𝑄 + 𝑀) [𝑊] A následný výpočet teploty v buňce po uplynutí časového kroku. 𝑇𝑛𝑜𝑣á = 𝑇𝑧𝑚ě𝑛𝑎 + 𝑇𝑝ů𝑣𝑜𝑑𝑛í =
(𝑁+𝐸+𝑆+𝑊+𝑄+𝑀)∙Δ𝑡 𝑉∙𝜌∙𝑐𝑝
(7)
+ 𝑇𝑝ů𝑣𝑜𝑑𝑛í [°𝐶]
Nová teplota se na konci iterace stává původní a proces se opakuje do dosažení výpočtového času. ZOBRAZENÍ VÝSLEDKŮ Hlavním výsledkem je rozložení teplot po uplynutí výpočtového času. 𝑟𝑠 [𝑚𝑚]
0,37
0,75
1,13
1,51
1,88
2,26
2,64
3,02
3,39
3,77
3,90
4,22
4,55
6,73
Typ
FUEL
FUEL
FUEL
FUEL
FUEL
FUEL
FUEL
FUEL
FUEL
FUEL
GAP
CLAD
CLAD
OUT
2,97 m 523,29 521,85 519,17 515,21 509,97 503,48 495,75 486,83 476,73 465,51 422,95 385,76 382,80 328,63 528,14
526,60
523,73 519,51
513,92
506,99
498,76
489,25
478,50
466,55
421,27
381,69
378,53 321,04
532,43
530,79
527,75 523,25
517,32
509,96
501,21
491,12
479,72
467,06
419,01
377,01
373,66 312,97
536,31
534,57
531,35 526,58
520,30
512,51
503,26
492,58
480,53
467,16
416,33
371,87
368,32 304,43
0 m 539,85 538,01 534,60 529,57 522,94 514,72 504,96 493,70 481,01 466,94 413,29 366,35 362,60 295,42
Obrázek 2: Graf rozložení teplot ve výpočetní síti [°C]
550,00
Teplota [°C]
500,00 450,00 Víko
400,00
Dno
350,00 300,00 250,00 FUEL FUEL FUEL FUEL FUEL FUEL FUEL FUEL FUEL FUEL GAP CLAD CLAD OUT Obrázek 3: Graf rozložení teplot ve středu buňky u "dna" a "víka" reaktoru
Teplota [°C]
550 500
500-550
450
450-500
400
400-450
350
350-400
300
300-350
250
250-300
Obrázek 4: 3D Graf rozložení teplot ZÁVĚR A DOPORUČENÍ Podařilo se vytvořit numerický model věrohodně odpovídající skutečnému stavu. Tento model lze použít nejen pro palivové tyče jaderných elektráren, ale pro všechny nestacionární, dvojdimenzionální, rotačně symetrické úlohy, jako je například trubkový výměník tepla či ohodnocení nestacionárních tepelných ztrát v izolovaném potrubí a podobně. LITERATURA Knižní publikace: [1] Michejev, M. A.: Základy sdílení tepla. Praha: Průmyslové vydavatelství, 1952. [2] Bečvář, J.: Jaderné elektrárny. Praha: Nakladatelství technické literatury, 1978.
