CITLIVOSTNÍ ANALÝZA JAKO NÁSTROJ PRO VERIFIKACI CFD MODELU A OPTIMALIZACI KONKRÉTNÍHO PRVKU

Simulace budov a techniky prostředí 2006 4. konference IBPSA-CZ Praha, 7. listopadu 2006

CITLIVOSTNÍ ANALÝZA JAKO NÁSTROJ PRO VERIFIKACI CFD MODELU A OPTIMALIZACI KONKRÉTNÍHO PRVKU Ondřej Hojer1, Jiří Bašta1, Jan Hensen2 1

Ústav techniky prostředí, Fakulta strojní, ČVUT v Praze e-mail: [email protected], [email protected] 2 Building Physics & Systems, Technische Universiteit Eindhoven www.bwk.tue.nl/bps/hensen

ABSTRAKT

SIMLAB

V článku je ukázáno využití citlivostní analýzy na příkladu nejprve jednoduchého modelu přenosu tepla radiací a následně složitějšího případu, tvořeného kombinací programů Gambit – Fluent. K této analýze byl použit software Simlab, který plně vyhovuje všem požadavkům a je doporučován v celé řadě prací zabývajících se touto tématikou. Na obou příkladech je detailně popsáno nastavení Simlabu a jsou uvedena základní doporučení pro praktické použití. V závěru je také jeden odstavec věnovaný použití citlivostní analýzy jako nástroje k verifikaci konkrétního CFD modelu a jsou zde diskutovány další možnosti využití.

Jedná se o program, který je navržený pro zjišťování nejistot a citlivostí modelu metodou Monte Carlo. Celý postup můžeme rozdělit na tři postupné kroky – přípravná fáze (pre-processing), vlastní model systému a analýza získaných výsledků (post-processing) – viz obr. 1.

ÚVOD Citlivostní analýza je nástroj, který určuje závislost výstupních parametrů na parametrech vstupních. Jedná se o velmi cennou pomůcku hlavně pro optimalizační postupy, neboť umožňuje mimo jiné určit pořadí důležitosti jednotlivých vstupních veličin a tedy dále se pak při optimalizaci soustředit pouze na ty nejdůležitější veličiny. Velmi často se může stát, že jsme při práci postaveni před problém optimalizace prvku, výpočetního postupu nebo metody. To si většinou žádá důkladnou analýzu celého systému, zjištění vzájemných vazeb a závislostí. V neposlední řadě i zhodnocení jak je výsledek citlivý na vstupy. A právě k tomuto účelu slouží citlivostní analýza. Citlivostní analýza je obecně určena ke zjišťování působení změn vstupních parametrů modelu na výsledné hodnoty (výsledky) tohoto modelu. Používá se ze tří základních důvodů: 1. 2. 3.

porozumění spolehlivosti výsledků; rozpoznání nejlepší regulovatelnosti systému; detailnější zaměření budoucích experimentů.

Citlivostní analýza má celou řadu velmi zajímavých a hlavně užitečných využití. Dále v článku se však budeme zabývat pouze její aplikací v oboru techniky prostředí, konkrétně jejím využitím jako pomůcky pro verifikaci a optimalizaci.

Obr. 1 – Postup výpočtu v programu Simlab V přípravné fázi vkládáme do Simlabu vstupní parametry, jejich střední hodnoty a standardní odchylky. Zadáváme jejich statistická rozdělení (normální, exponenciální, gama nebo jiná další) a volíme tzv. vzorkovací metodu. Vzorkovací metodou je myšleno vytvoření matice MxN, kde pro každou z M vstupních proměnných vygeneruje Simlab N náhodných hodnot podle předchozího nastavení. Pak je zapíše do externího souboru. V druhé fázi je třeba zadat, zda je model interní (je-li ho možné popsat jednoduchou rovnicí pomocí dříve zadaných parametrů a konstant) a nebo externí, kdy navzorkovaná vstupní data budou vložena do nějakého jiného programu, který provede výpočty. Výhodou je, že tento externí program (model) může být naprosto libovolný (např. Matlab, Excel, IES, Energy+ nebo Fluent). Uvnitř zvoleného modelu tedy proběhne N potřebných výpočtů (cyklů) a výstupy musí být poté opět zpětně uloženy do souboru se specifickou strukturou, aby byl Simlab schopen tyto hodnoty znovu použít. Slovem cyklus budeme tedy dále označovat jeden výpočet modelu, kde jsou jako vstupní veličiny zadány hodnoty příslušného řádku matice vstupních veličin.

A konečně poslední, třetí fáze provádí vlastní vyhodnocení a to analýzu nejistoty (určující nejistotu v předpokladech, které model poskytuje) a analýzu citlivostní, hodnotící vztah mezi nejistotami nezávislých vstupních parametrů a nejistotami závislých výstupních parametrů. Obě obsahují celou řadu výstupů a to jak grafických tak i tabelovaných.

PŘÍKLAD POUŽITÍ Použití výše uvedeného postupu bude ukázáno na jednoduchém příkladu - přestupu tepla sáláním mezi dvěma přesně definovanými povrchy. Situace je zobrazena na obr. 2. Destička na podlaze je umístěna 0,1 m od počátku jak ve směru x tak y. Další technický popis je následující: rozměry prostoru (1 x 1 x 1) m sálající ploška S1 = 0,0138 m2 t1 = 350 °C e1 = 0,90 ploška na podlaze (0,1 x 0,1) m t2 = 10 °C e2 = 0,85

Obr. 2 – Situace pro příklad přenosu tepla radiací s využitím citlivostní analýzy v Simlabu Na popis tohoto problému bude stačit jednoduchá rovnice pro přestup tepla sáláním mezi dvěma povrchy

[

Q12 = σ ⋅ ε 1 ⋅ ε 2 ⋅ S1 ⋅ ϕ12 ⋅ T14 − T24

]

(1)

Hodnotu všech veličin až na poměr osálání již známe, nicméně ten je možné snadno zjistit nebo vypočítat. V dalším kroku určíme nejistotu všech parametrů. Pro tento příklad bude zvolena pro všechny hodnoty relativní nejistota 5 %. Všechny parametry i s příslušnými standardními odchylkami jsou shrnuté v tabulce 1. Pro všechny tyto parametry byla zvolena vzorkovací metoda Random [3] a pro první přiblížení 100 cyklů. Manuál Simlabu udává, že počet cyklů by rozhodně neměl být menší než 1,5 násobek počtu proměnných (1,5 x M). Doporučuje hodnotu 10 x počet proměnných. Jak se však dále ukáže, ani tato hodnota není k dosažení správných

hodnot dostačující. Statistické rozdělení bylo u všech hodnot zvoleno normální. Tab. 1. – Hodnoty parametrů a jejich nejistot pro příklad přestupu tepla radiací Parametry e1 [-] e2 [-] S1 [m2] T1 [K] T2 [K] f12 [-] s [W.m-2.K-4]

Střední Standardní hodnota odchylka 0,90 ± 0,05 0,85 ± 0,04 0,0138 ± 0,0007 623 ± 31 283 ± 14 0,0029 ± 0,0002 5,67·10-8

Když je popis systému jednoduchý jako v tomto případě, je výhodné využít interního modelu Simlabu a rovnici (1) vytvořit pomocí kalkulačky uvnitř Simlabu. Neobjevují se žádné komplikace s výstupními ani vstupními soubory a celý výpočet je tak velmi zrychlený a snadno opakovatelný pro libovolný počet cyklů. Jak je patrné z tabulky 2, analýza byla současně provedena několika technikami a tedy hodnocena několika různými citlivostními indexy (PEAR, SPEA, PCC, PRCC, SRC, SRRC a Smirnov). Rozdíl mezi jednotlivými citlivostními indexy spočívá ve způsobu jejich vyhodnocení. PEAR a SPEA jsou založeny na vyhodnocení bodových diagramů (scatterplot) závislé na každé jednotlivé nezávislé proměnné, kde PEAR je pro hodnocení lineárních modelů, zatímco SPEA se používá pro hodnocení nelineárních. Další hodnotící index využívá k řešení regresní analýzu. Měřítkem citlivosti jsou standardizované regresní koeficienty SRC. Používá se metoda lineární regrese. PCC neboli částečné korelační koeficienty jsou založeny na základech korelace a částečné korelace. Udávají váhu jednotlivých vstupů na výstupu bez uvažování jakékoli vazby s jiným vstupem. Regresní analýza často nepracuje spolehlivě, pokud vztah mezi vstupními proměnnými není lineární. V tomto případě je vhodné použít transformované indexy PRCC nebo SRRC. Transformace spočívá v nahrazení vlastních dat jejich pořadím (nejnižší hodnotě přiřadit číslo 1, atd.). Regresní analýza je pak provedena na těchto transformovaných hodnotách. Podle různých indexů lze tedy dojít k různým závěrům a je hlavně na autorovi analýzy, aby byl schopen zvolit vhodnou techniku, podle které případ vyhodnotí. Na základě hodnocení výhod a nevýhod jednotlivých indexů, jak jsou rozepsány v manuálu Simlabu [3] bylo zvoleno pro tento případ hodnocení indexem SRC. Podle absolutní velikosti hodnot v jednotlivých řádcích vidíme, že výsledek Q12 je nejvíce citlivý na T1, to znamená teplotu sálající plošky 1. Takový výsledek se vzhledem k charakteru rovnice (1) dal předpokládat.

Tab. 2 – Výsledky první analýzy přestupu tepla radiací pro 100 cyklů

Tab. 3 – Výsledky první analýzy přestupu tepla radiací pro 65000 cyklů

Další pořadí ovšem na první pohled nemusí být tak jasné a zde nám právě mohou pomoci výsledky zpracovávané analýzy. Z výsledků vidíme, že následuje poměr osálání f12, sálající plocha S1, emisivita této plochy e1, emisivita přijímající plochy e2 a nakonec teplota přijímající plochy T2. Rozdílná hodnota indexů SRC u T1 a T2 je způsobena jejich rozdílnou absolutní hodnotou a nejistotou. Záporná hodnota vyjadřuje, že při zvyšování hodnoty T2, klesá výsledná hodnota Q12. Co je na těchto výsledcích na první pohled zarážející, je nestejný citlivostní index u e1 a e2. V rovnici stojí na stejné úrovni, mají téměř stejnou střední hodnotu a standardní odchylku a přesto je zde relativně velký rozdíl. Další zvláštností bylo červené (tmavé) pozadí za všemi hodnotami u SRC, které vyjadřuje, že hodnoty neprošly tzv. Čebyševovým testem významnosti. Zelené (světlejší) pozadí znamená, že hodnoty testem prošly a bílé pozadí znamená, že test nebylo možné provést. Podle tohoto testu bychom vlastně neměli žádnou z hodnot u SRC indexu považovat za relevantní (statisticky významnou). U jiných hodnotících indexů byla situace obdobná a tak bylo třeba hledat jiné řešení. Nakonec se ukázalo, že výsledek se značně zlepší, pokud provedeme úpravu nastavení Simlabu ve smyslu navýšení počtu cyklů (výpočtů modelu). V tomto případě bylo zvoleno 65000, aby byl výsledný rozdíl co nejvýraznější. Pokud se nyní podíváme na výsledky (tab. 3), zjišťujeme, že SRC u e2 se srovnal s e1, S1 a f12 a i Čebyševův test dopadl mnohem lépe. T2 sice stále nevyhovuje, ale pro nás to jen znamená, že tuto proměnnou nebudeme ze statistického hlediska považovat pro hodnotu Q12 ani její nejistotu za významnou. Na základě tohoto příkladu můžeme shrnout, že je třeba:

K dokreslení možností, které Simlab nabízí, je zajímavé se ještě v obou případech podívat na analýzu nejistoty, tedy v jakých mezích se pohybují výstupní hodnoty. Z obr. 3 a 4 je stejně jako u citlivostní analýzy na první pohled zřejmý vliv počtu cyklů na celkové rozložení výsledných hodnot. Výhodou zejména u složitějších systémů je, že výsledky analýzy mohou být také zobrazeny ve formě statistických hodnot charakterizujících rozdělení závislé veličiny jako jsou střední hodnota, rozptyl, standardní odchylka a další. Pro obr. 3 a 4 bude uvedena střední hodnota se standardní odchylkou pod každým z obrázků v popisu. Výsledky i přes rozdílnost obou obrázků vykazují překvapivě velmi dobrou shodu. Můžeme konstatovat, že tabelární výstupy analýzy nejistoty se na rozdíl od grafických s narůstajícím počtem cyklů téměř nemění.

1. 2. 3.

vždy velmi uvážlivě vybírat, jaký citlivostní index zvolíme pro hodnocení; dávat pozor na závislost výsledků na počtu cyklů; zvážit, jak velká nepřesnost může vzniknout, pokud výsledek neprojde Čebyševovým testem.

6 5 4 3 2 1 0 0.10

0.20

0.30

0.40

0.50

Obr. 3 – Zobrazení frekvenční distribuční funkce závislé veličiny při 100 cyklech, osa x hodnota přeneseného výkonu [W], osa y četnost výskytu Q12 = (0,26 ± 0,06) W 3000 2500 2000 1500 1000 500 0 0.10

0.20

0.30

0.40

0.50

Obr. 4 – Zobrazení frekvenční distribuční funkce závislé veličiny při 65000 cyklech, osa x hodnota přeneseného výkonu [W], osa y četnost výskytu Q12 = (0,25 ± 0,06) W

CFD APLIKACE Až dosud jsme se zabývali pouze jednoduchým využitím citlivostní analýzy v příkladu, který byl snadno popsatelný. V praxi se však s takovými příklady téměř nesetkáme a vždy se jedná o model daleko složitější, málokdy popsatelný jednoduchými rovnicemi. A právě na takové případy lze využít externí model. V našem případě budeme hovořit o modelu vytvořeném pro simulační prostředí CFD (Gambit – Fluent). Jedná se o stejný případ jako je na obr. 1 s tím rozdílem, že přichází v úvahu více parametrů (tab. 4). Tab. 4 – Hodnoty parametrů a jejich nejistot pro případ citlivostní analýzy programem Simlab s externím modelem pro CFD Parametry

Popis

Hodnota

Odchylka*

a [m]

Šířka sálající plošky

0,069

± 0,004

b [m]

Délka sálající plošky

0,20

± 0,01

h [m]

Hloubka reflektoru

0,080

± 0,004

a [°]

Úhel podélného reflektoru

45

±2

b [°]

Úhel příčného reflektoru

45

±2

e1 [-]

Emisivita destičky

0,90

± 0,05

e2 [-]

Emisivita reflektoru

0,40

± 0,02

eo [-]

Emisivita stěn a stropu

0,85

± 0,04

ep [-]

Emisivita podlahy

0,85

± 0,04

T1 [K]

Povrchová teplota destiček

1173

± 59

T2 [K]

Povrchová teplota podlahy

283

± 14

*hodnoty standardních odchylek

Cílem celé analýzy je stanovení citlivostních indexů jednotlivých parametrů při přestupu tepla sáláním v uzavřeném prostoru z povrchu keramické destičky světlého plynového zářiče na plošku na podlaze. Jelikož se nacházíme zatím pouze na začátku této analýzy, byla zvolena pro první přiblížení mnohá zjednodušení, aby bylo dosaženo rychlé konvergence výsledků (při dlouhém výpočtu by nebylo možné provést takový počet cyklů, jaký by byl potřeba). Zejména se jednalo o: 1. 2. 3.

zjednodušení modelu (geometrické i fyzikální hledisko); snížení počtu parametrů (ne všechny ovlivňující hodnoty byly uvažovány jako parametry); odhad nejistot u všech parametrů na stejnou relativní hodnotu.

Obr. 5 – Schéma externího modelu citlivostní analýzy pro program Simlab a simulační prostředí CFD Externí model (obr. 5) se skládá ze tří základních prvků, programů Simlab 2.2, Gambit 2.3 a Fluent 6.2. Začátek postupu je stejný jako v případě prvního příkladu. Parametry se zadají do Simlabu (1). Ten podle nastavení připraví matici MxN a zapíše ji do externího souboru (2); v tomto případě byla využita možnost propojení s Excelem. Simlab vytvoří v zadaném souboru list se jménem Inputs a čeká do té doby, než bude soubor uložen a Excel uzavřen. Předpokládá přitom, že se ve stejném souboru v listu nazvaném Outputs objeví sloupec začínající jménem závislé proměnné a seznamem výsledných hodnot odpovídajících počtu cyklů N. Hodnoty z listu Inputs jsou mezitím vloženy do vstupního souboru pro Gambit, tzv. žurnálu a současně skriptu generujícím žurnál pro Fluent. Je spuštěn proces generace sítě (mesh) pro všech N kombinací vstupních hodnot (3) a zároveň vytvořen žurnál pro Fluent. Fluent je spuštěn, postupně si načítá jednotlivé soubory s hotovou sítí (4) a zapisuje po konvergenci každého případu hledané hodnoty dopadajícího radiačního toku do dalšího souboru (5). Výsledky jsou opět načteny do původního Excel souboru do listu Outputs a Excel se zavírá. Tímto krokem se opět dostáváme zpět k Simlabu, který si načte potřebné hodnoty (6) a celý výsledek vyhodnotí (7), jako tomu bylo v jednoduchém případě. Pokud se nyní opět podíváme na schéma celého externího modelu (obr. 5), můžeme si vytvořit představu, jaký může být rozdíl mezi jednoduchým interním modelem a modelem externím. Aby bylo možné zahrnout všechny parametry do úvahy, bylo třeba na různých úrovních modelu vyřešit celou řadu problémů. První problém se objevil s generováním velkého množství souborů se sítí v Gambitu s rozdílnými rozměry sálající destičky a rozměry reflektorů. Tento problém byl vyřešen vytvořením žurnálu, který postupně v každém cyklu načte příslušné rozměry a podle nich vytvoří síť. Pouhou změnou počtu cyklů pak lze generovat těchto souborů libovolný počet. Další problém se objevil při snaze zahrnout do analýzy jiné proměnné, které nebylo možné vložit do Gambitu

(emisivity, teploty). Snaha vytvořit žurnál do Fluentu podobným způsobem jako do Gambitu selhala a tak bylo třeba nalézt jinou cestu. Nakonec se podařilo vytvořit skript, který Fluent žurnál vytvoří sám. Podobně jako v případě Gambitu pozdější změna počtu cyklů je již velice snadnou záležitostí.

ních parametrů. Nemusí se dokonce jednat pouze o analýzu jedné závislé proměnné, ale takových proměnných může být v rámci jednoho výpočtu hned několik.

VÝSLEDKY Z výsledků prvního přiblížení, jak je vidět v tabulce 6, si již můžeme vytvořit představu, které veličiny jsou pro systém více důležité a které méně. Například váha některých veličin jako je emisivita reflektoru e1, emisivita okolních stěn a stropu eo, ale i úhel nastavení příčného reflektoru b se zdá být zanedbatelný. Naopak podle očekávání největší vliv má teplota povrchu 1, rozměry sálající destičky a a b, úhel podélného reflektoru a a překvapivě i teplota povrchu 2, jejíž vliv jsme v prvním příkladu mohli zanedbat. Zůstává ještě několik otázek, které je třeba vyřešit. Například jak je možné, že SRC u T2 je v tomto případě najednou kladné, když by mělo mít stejně jako v prvním příkladu záporný vliv? Nebo z jakého důvodu jsou různá znaménka u úhlů nastavení a a b? Je samozřejmě možné, že zvýšením počtu cyklů se ještě výsledek nějakým směrem posune, ale jak jsme viděli u prvního příkladu, posun nemůžeme očekávat v řádu desetin, ale spíše setin. Hlavně předpokládáme zvýšení počtu hodnot, které projdou Čebyševovým testem, protože podle těchto výsledků žádná z SRC tímto testem neprošla. Tab. 5 – Výsledky citlivostní analýzy prvního přiblížení pro 110 cyklů s programem Simlab a externím modelem v CFD

Kromě výsledků z citlivostní analýzy je také zajímavé se podívat na výsledné rozložení dopadajícího sálavého toku jednoho ze 110 případů (obr. 6). Místo, ve kterém snímáme hodnotící veličinu, je zde na první pohled velice dobře patrné a lze si i velice snadno vytvořit představu o rozložení sálavého toku v oblasti pod zářičem. Absolutní hodnoty nejsou důležité, neboť se nejedná o reálný případ, ale pouze modelový. Výhoda Fluentu jako externího modelu spočívá v širokých možnostech výběru vstupních a výstup-

Obr. 6 – Schéma externího modelu citlivostní analýzy pro program Simlab a simulační prostředí CFD

VERIFIKACE MODELU V průběhu řešení byla citlivostní analýza použita zároveň pro verifikaci CFD modelu. Bylo třeba určit optimální nastavení tzv. parametrů relaxace (under-relaxation factors URF), jejichž hodnota značně ovlivňuje rychlost konvergence. Jako proměnné byly tedy zvoleny tyto URF a za výstupní hodnotu číslo iterace při které dojde ke konvergenci. Jako měřítko konvergence byla brána v úvahu pouze konvergenční kritéria, nikoli další ukazatele jako např. výsledná tepelná bilance nebo průběh sledované veličiny. Důvodem bylo, že je ve Fluentu možné automaticky ukončit výpočet, pokud tento konvergenční ukazatel dosáhne konkrétní hodnoty. U ostatních ukazatelů by automatické ukončení výpočtu po dosažení určité hodnoty bylo jen velmi těžko proveditelné. Výsledkem bylo pořadí citlivosti jednotlivých faktorů. Jednalo se však pouze o pomocnou informaci, neboť právě u URF platí, že jejich citlivost není lineární. Na tento problém je třeba si dát pozor u každého analyzovaného systému a v žádném případě závěry automaticky nezobecňovat. Vždy záleží na tom, jaká hodnota parametru bude pro analýzu zvolena. Nicméně alespoň z hlediska získání představy o chování celého systému byla tato informace velmi cenná. Podobné analýzy by se dalo využít například pro zjištění vlivu různých okrajových podmínek na konvergenci řešení nebo na výslednou tepelnou bilanci. Aplikací je celá řada a určitě v mnoha případech vhodnějších než velmi rozšířená metoda pokus – omyl.

ZÁVĚR Hlavním výstupem je zdokumentování postupu, který umožňuje aplikovat citlivostní analýzu v programu Simlab s využitím simulačního prostředí Gambit – Fluent. Z uvedených případů vyplývá, že

uplatnění citlivostí analýzy s výše zmíněným simulačním prostředím je limitované velikostí modelu. Vzhledem k velkému množství cyklů výpočtu, které jsou nutné, je často problematické dosáhnout konvergence v rozumně krátkou dobu. Citlivostní analýza z toho důvodu není prakticky použitelná pro komplikované modely. Naproti tomu má celou řadu výhod a ani její použití se simulačním prostředím CFD by nemělo být v inženýrské praxi opomíjeno. V další fázi se předpokládá, jak již bylo konstatováno výše, rozšíření nejdříve počtu cyklů a následně zahrnutí rovnic popisujících proudění do výpočtu Fluentu. Očekává se, že bude poté možné výsledky lépe zobecnit a stanovit absolutní citlivostní stupnici pro ovlivňující parametry přestupu tepla od světlého plynového zářiče. Výsledky analýzy poté budou uplatněny pro optimalizaci světlého plynového zářiče a jeho geometrie sálání.

PODĚKOVÁNÍ Tento článek byl podpořen z výzkumného záměru MSM 6840770011 Technika životního prostředí.

LITERATURA [1] Ferson S., Tucker T.: Sensitivity in risk analyses with uncertain numbers. Sandia National Laboratories. SAND 2006-2801. Setauket, New York. 2006. [2] Saltelli A. Tarantola S., Campolongo, F. and Ratto, M.: Sensitivity Analysis in Practice. A Guide to Assessing Scientific Models, John Wiley & Sons publishers. 2004. [3] Simlab 2.2 Software and Reference manual. Dostupné na URL adrese [4] Fluent 6.2 Software and Manual. Dostupné na URL adrese [5] Gambit 2.3 Software and Manual. Dostupné na URL adrese

PŘEHLED OZNAČENÍ a b h S1 a b

e1 e2 eo ep f12 s

šířka sálající plošky [m] délka sálající plošky [m] hloubka reflektoru [m] plocha sálající plošky [m2] úhel nastavení podélného reflektoru [°] úhel nastavení příčného reflektoru [°] emisivita destičky [-] emisivita plechu reflektoru [-] emisivita okolních stěn a stropu [-] emisivita podlahy [-] poměr osálání mezi povrchy 1 a 2 [-] Stefan – Bolzmanova konstanta [W·m–2·K–4]

t1 t2 T1 T2

povrchová teplota sálajícího povrchu [°C] povrchová teplota přijímajícího povrchu [°C] absolutní povrchová teplota sálajícího povrchu [K] absolutní povrchová teplota přijímacího povrchu [K]