Mendelova univerzita v Brně Lesnická a dřevařská fakulta Ústav nauky o dřevě
Model vázaného pohybu vlhkostního a teplotního pole ve dřevě
Disertační práce
2012
Miroslav Trcala
2
Prohlášení Prohlašuji, že jsem disertační práci na téma: Model vázaného pohybu vlhkostního a teplotního pole ve dřevě zpracoval sám a uvedl jsem všechny použité prameny. Souhlasím, aby moje diplomová práce byla zveřejněna v souladu s § 47b Zákona č. 111/1998 Sb., o vysokých školách a uložena v knihovně Mendelovy univerzity v Brně, zpřístupněna ke studijním účelům ve shodě s Vyhláškou rektora MU v Brně o archivaci elektronické podoby závěrečných prací.
Autor kvalifikační práce se dále zavazuje, že před sepsáním licenční smlouvy o využití autorských práv díla s jinou osobou (subjektem) si vyžádá písemné stanovisko univerzity o tom, že předmětná licenční smlouva není v rozporu s oprávněnými zájmy univerzity a zavazuje se uhradit případný příspěvek na úhradu nákladů spojených se vznikem díla dle řádné kalkulace.
V Brně, dne: ........................................
Miroslav Trcala
3
Poděkování Děkuji mému školiteli doc. Dr. Ing. Petru Horáčkovi za rady a připomínky k této práci. Dále děkuji doc. Ing. Petru Koňasovi, Ph.D. za odborné diskuze, které pomohly nasměrovat mou práci. Velké poděkování patří prof. Ing. Janu Čermákovi, CSc. za dlouhodobou spolupráci a pomoc s publikacemi. V neposlední řadě chci touto cestou poděkovat lidem, kteří mi byli a stále jsou oporou, je to rodina, přítelkyně a všichni mí blízcí.
4
1 ÚVOD............................................................................................................................ 6 1.1 DETERMINISTICKÉ MODELOVÁNÍ FYZIKÁLNÍCH POLÍ VE DŘEVĚ ............................... 7 1.2 STOCHASTICKÉ MODELOVÁNÍ FYZIKÁLNÍCH POLÍ VE DŘEVĚ .................................... 7 1.3 MĚŘENÍ TRANSPIRAČNÍHO PROUDU A TEPELNÉ VODIVOSTI DŘEVA V KMENU STROMU ......................................................................................................................... 8 2 CÍL PRÁCE................................................................................................................ 11 3 LITERÁRNÍ PŘEHLED .......................................................................................... 12 3.1 DETERMINISTICKÉ MODELOVÁNÍ FYZIKÁLNÍCH POLÍ VE DŘEVĚ ............................. 12 3.2 STOCHASTICKÉ MODELOVÁNÍ FYZIKÁLNÍCH POLÍ VE DŘEVĚ .................................. 15 3.3 MĚŘENÍ TRANSPIRAČNÍHO PROUDU A TEPELNÉ VODIVOSTI DŘEVA V KMENU STROMU ....................................................................................................................... 17 4 METODIKA............................................................................................................... 19 4.1 MODELOVÁNÍ FYZIKÁLNÍCH POLÍ VE DŘEVĚ BĚHEM SUŠENÍ NEBO V DŘEVOSTAVBÁCH ..................................................................................................... 20 4.1.1 Deterministické modelování fyzikálních polí ve dřevě .................................. 20 4.1.1.1 Modelování vázaného teplotního a vlhkostního pole ve dřevě............... 20 4.1.1.2 Simulace teplotního a vlhkostního pole ve dřevě ................................... 23 4.1.1.3 Modelování vázaného teplotního, vlhkostního a deformačního pole ve dřevě.................................................................................................................... 26 4.1.2 Stochastické modelování fyzikálních polí ve dřevě........................................ 29 4.1.2.1 Stochastické modelování difuzních problémů při sušení dřeva.............. 29 4.1.2.2 Stochastické modelování teplotních problémů v konstrukci dřevostavby ............................................................................................................................ 31 4.2 MODELOVÁNÍ FYZIKÁLNÍCH POLÍ VE DŘEVĚ PŘI MĚŘENÍ TRANSPIRAČNÍHO PROUDU A TEPELNÉ VODIVOSTI V BĚLI KMENE STROMU ............................................................ 33 4.2.1 Měření s lineárním ohřevem .......................................................................... 33 4.2.1.1 Metoda tepelné bilance kolem lineárního ohřevu................................... 33 4.2.1.2 Výpočet transpiračního toku a tepelné vodivosti současně .................... 35 4.2.2 Měření s objemovým ohřevem ....................................................................... 36 4.2.2.1 Současný stav metody tepelné bilance (THB metody) .......................... 36 4.2.2.2 Vyjádření fiktivního toku v závislosti na tepelné vodivosti a okamžité hodnotě transpiračního toku................................................................................ 37 5 ZÁVĚR ....................................................................................................................... 40 SEZNAM POUŽITÉ LITERATURY......................................................................... 43 PUBLIKACE AUTORA DISERTAČNÍ PRÁCE ..................................................... 48
5
1 Úvod V inženýrských problémech se často řeší úlohy spojené s tepelným a vlhkostním chováním dřeva. Setkáváme se s nimi například při sušení řeziva nebo u dřevostaveb, ale také při měření transpiračního proudu nebo tepelné vodivosti dřeva stromů. Dřevo jako materiál se složitou vnitřní strukturou se často v těchto situacích chová nepředvídatelně a k analýze tohoto chování je výhodné použít kvalitních fyzikálních a matematických přístupů, které umožní daný problém vyřešit. Jde především o přístupy, které pro popis těchto fyzikálních jevů využívají soustav parciálních diferenciálních rovnic s příslušnými okrajovými a počátečními podmínkami. Tyto soustavy jsou pak řešeny pomocí numerických metod (nejčastěji metodou konečných prvků). Snahou vědecko-výzkumných pracovníků je sestavit matematické modely těchto složitých fyzikálních jevů spojených s teplotně-vlhkostním chováním dřeva v různých úlohách dřevařské a lesnické praxe. Matematické modelování se stalo důležitým nástrojem při simulacích, analýzách, předvídání a optimalizaci různých procesů, jevů a chování systémů. Systém je v této práci chápán jako určitá abstrakce reálného světa, kterou si lidé vytvářejí během jeho poznání. Použití matematického modelu systému umožňuje zjistit informace o chování systému, i když ze skutečného systému je to nemožné nebo obtížné, urychluje, usnadňuje a racionalizuje proces poznání objektivní reality, umožňuje variantní řešení a tím i nalezení optimálního řešení a v neposlední řadě identifikuje a kvantifikuje nejistoty v hledaném řešení vzhledem k nejistotám ve vstupních údajích modelu. Matematické modely lze pro účely této práce rozdělit do dvou skupin: deterministické a stochastické. Model je vždy jen zjednodušením reálného systému, a obsahuje tak prvky, kterým se říká neurčitosti nebo nejistoty. S analýzou nejistot je úzce spojena tzv. analýza citlivosti. Zatímco při analýze nejistot charakterizujeme veškeré nejistoty na vstupu a chceme kvantifikovat celkovou neurčitost na výstupu modelu, při analýze citlivosti zkoumáme, jak konkrétně ovlivňují výstup jednotlivé parametry modelu.
6
1.1 Deterministické modelování fyzikálních polí ve dřevě Náplní této práce je popsat teplotně-vlhkostní stavy dřeva při jeho sušení nebo zabudování do stavby. Tento popis lze zajistit sestavením adekvátního matematického modelu. Aby mohl být odvozený matematický model využíván osobami bez hlubší znalosti matematiky a programování, jsou k němu vyvíjeny (programovány) aplikace. Tyto aplikace mají příjemné grafické uživatelské rozhraní, kde uživatel jednoduše zadá jen základní materiálové charakteristiky (rozměry, fyzikální vlastnosti jako hustota v absolutně suchém stavu apod.), parametry okolního prostředí (vzduchu) a aplikace sama spočítá a v přehledné formě zobrazí údaje o vlhkosti a teplotě uvnitř zvolené konstrukce dřevostavby nebo sušeného řeziva. Sušení řeziva je dynamický proces, při kterém se podmínky působící na dřevo mění s časem a dřevo vysychá nerovnoměrně v celém svém objemu. Tyto vnější podmínky regulujeme dle sušících režimů. Pro srovnání různých sušících režimů je důležité znát nejen energetickou a časovou nákladnost, ale i kvalitu sušení (vysušit na požadovanou vlhkost s minimálními odchylkami, udržovat vznikající vnitřní napětí v dovolených mezích). Tyto požadavky jsou často protichůdné a důležitou činností je pak optimalizace sušících režimů, jejíž výsledek (rozhodnutí o tom, který sušící režim je v konkrétním případě vhodné použít) ovlivňuje ekonomiku sušení a rozhoduje tak o hospodářské úspěšnosti či neúspěšnosti daného sušícího podniku. Pro získání představy o tom, který sušící režim je či není agresivní pro daný materiál, lze v současné době využít výše zmíněných počítačových simulací. Toto lze provést sestavením vhodného modelu, který co nejvěrněji popíše stavové proměnné, kterými jsou vlhkost a deformace sušeného materiálu. Z tohoto důvodu jsou k teplotně-vlhkostním rovnicím přidány rovnice popisující deformace vznikající v důsledku změn vlhkosti a v menší míře (často zanedbatelné) v důsledku změn teploty. 1.2 Stochastické modelování fyzikálních polí ve dřevě Cílem simulace je zjistit, jak se bude systém chovat pro zadaná vstupní data. Například jak se bude měnit vlhkostní pole (nebo deformačně-napěťové pole) ve dřevě během sušení pro zadaný sušící režim. Pokud tuto simulaci provedeme opakovaně pro různé sušící režimy, tak lze podle stanoveného kritéria (nejen fyzikálního, ale i 7
ekonomického apod.) rozhodnout o tom, který sušící režim je nejvhodnější pro sušení daného druhu řeziva. Ovšem toto řezivo je pro účely simulace bráno jako fyzikální kontinuum, které je jednoznačně definováno pomocí svých fyzikálních vlastností jako je hustota, měrné teplo, tepelná vodivost, vlhkostní vodivost apod. Tyto vlastnosti jsou stejně jako sušící režim vstupními daty, na kterých závisí výsledky simulace a tedy i přesnost, se kterou aproximujeme reálný systém (reálné vlhkostní pole). Tyto vlastnosti však mnohdy nejsou u dřeva známy s dostatečnou přesností a hodnota vlastnosti reálného systému je dost odlišná od tabelovaných průměrných hodnot. Je to dáno složitou strukturou dřeva, nepřesností měření materiálových vlastností a také lidskou neschopností přesně popsat všechny fyzikální zákonitosti. Tyto nejistoty ve vstupních parametrech modelu lze statistickými nástroji kvantifikovat a počítat s nimi při simulaci. Nejčastěji se tento problém řeší metodou Monte Carlo, která používá generátor pseudonáhodných čísel. Princip metody je jednoduchý. Pomocí generátoru náhodných čísel opakovaně generujeme hodnoty náhodných vstupních parametrů modelu a pro každý takto vygenerovaný vstup provedeme numerický výpočet, který je založen na standardním deterministickém modelu. Výsledky těchto simulací se ukládají a po provedení dostatečného počtu simulací se výsledky statisticky vyhodnotí, konkrétně se zjistí pravděpodobnostní rozložení výsledků (histogram) a vypočítá se průměr, směrodatná odchylka apod. V této práci je věnována pozornost speciální třídě metod, které jsou založené na sestavení tzv. spektrálního stochastického modelu, který je pomocí vhodných výpočetních programů řešen a výsledky jsou srovnány se statistickými výstupy metody Monte Carlo. Ukazuje se, že použití těchto spektrálních metod je u některých úloh výhodnější, protože dají stejný výsledek za mnohem kratší dobu výpočtu a tím se výrazně ušetří náklady na testování a vývoj nových konstrukcí nebo technologií.
1.3 Měření transpiračního proudu a tepelné vodivosti dřeva v kmenu stromu Vodní
provoz
rostlin
je
všeobecně
spojen
s největšími
toky
energie
v ekosystémech. Voda je nejčastější přírodní limitující faktor růstu a nezbytným základem umožňujícím existenci dalších fyziologických procesů. Lesy jsou největším kontrolním mechanismem koloběhu vody na kontinentech a klimatu. Největší objem 8
struktur stromů a porostů představuje jejich vodo-vodivý systém, zahrnující jak dřevo, tak lýko. Změny velikosti tohoto vodivého systému představují růst. Současná technika nám umožňuje měřit mobilní soupravou přístrojů nezávislou na stacionárních objektech kvantitativní parametry vodního provozu, struktury a růstu na úrovni celých stromů a porostů a jejich kombinací i na úrovni povodí či větších lesních celků. Metod pro měření transpiračního proudu je několik (viz. kapitola 3 Literární přehled). Tato práce se zabývá metodami termodynamickými, které jsou zde rozděleny do dvou skupin. První skupinu tvoří metody, které pro měření transpiračního proudu využívají lineárního ohřevu měřeného místa (pomocí vyhřívané jehle) a druhou skupinu tvoří metody využívající objemového ohřevu v místě měření (například pomocí elektrod). Cíl obou metod je stejný - ze zákonů termodynamiky z naměřených teplotních diferencí odvodit transpirační proud. V prvním případě (lineární ohřev) je příkladem metoda deformace teplotního pole (HFD - heat field deformation) a ve druhém případě metoda tepelné bilance části kmene (THB - trunk heat balance). Metoda HFD je založena na měření teplotního pole kolem lineárního ohřevu (kolem vyhřívané jehly) vloženého do kmene v radiálním směru. Tato metoda je založena na poměru teplotních diferencí v axiálním a tangenciálním směru a umožňuje měřit v řadě bodů napříč bělí a tím získat radiální profily transpiračního proudu. Tyto diference jsou měřeny dvěma páry termočlánků vloženými (podobně jako ohřev) v injekční jehle z nerezové oceli o průměru 1.5 mm. Reverzní toky spojené s vodní redistribucí mohou být touto metodou také měřeny stejně jako noční resaturační toky, které poskytují cenné informace o vodním stavu rostlin. U velkých stromů může být HFD metoda použita, jen pokud aplikujeme multibodová teplotní čidla, která umožňují měření radiálního profilu transpiračního proudu a lze pak provést přibližnou integraci proudu přes celou tloušťku běle. Jednobodová teplotní čidla tedy mohou být použita pouze v případě malých kmenů, kořenů nebo větví. Nevýhodou této metody je (jako u všech metod s lineárním ohřevem) nejistota při měření velkých transpiračních proudů. Metoda THB je založena na výpočtu tepelné bilance rozměrově definované, mírně zahřívané části kmene. Přiváděné množství tepelné energie je odváděno kondukcí skrze rostlinná pletiva a také vodním (transpiračním) proudem, který těmito pletivy (xylemem) protéká. Většina tepla je odváděna proudící vodou a kondukcí (zahříváním pletiv v okolí měřiště) je ztracena menší část tohoto tepla (asi 10-20%). Tyto tepelné 9
ztráty zahrnuté do aplikovaných rovnic jsou částečně technicky eliminovány (tepelnou izolací měřiště např. polyuretanovou pěnou a jeho odstíněním proti přímé sluneční radiaci reflexním štítem), avšak přesto k nim dochází. Mění se spolu s tepelným polem v závislosti na intenzitě vodního proudu. Naproti tomu vliv změn obsahu vody v xylemu je zanedbatelný. Tepelné ztráty jsou zřetelně viditelné na záznamech vodního proudu jako určitá hodnota “fiktivního toku”, který je zaznamenáván i když skutečný proud se rovná nule. Jestliže vypočítáváme skutečnou hodnotu vodního toku, je nutné tento fiktivní tok ze zaznamenávaných dat odečíst. Jediným problémem, který je řešen v této práci, je jak tento fiktivní tok vypočítat. Je zde také řešen problém měření nulového a záporného toku pomocí nového zapojení se symetrickým uspořádáním termočlánků nad a pod vyhřívanou oblastí (elektrodami).
10
2 Cíl práce Cílem práce je v podobě matematických modelů a numerických simulací aplikovat současné teoretické poznatky z oblasti studie fyzikálních polí ve dřevě, zejména pole teplotního a vlhkostního, a tyto teoretické poznatky podle potřeb v daných oblastech dále rozvinout. Práce je zaměřena na studium fyzikálních polí ve dřevě v oblasti sušení řeziva nebo stavební fyziky u dřevostaveb. Dále se práce zabývá měřením transpiračního proudu a tepelné vodivosti dřeva v kmenu stromu. Přehledně jsou dílčí cíle práce shrnuty do následujících bodů: • sestavení a řešení multifyzikálního modelu popisujícího teplotní a vlhkostní pole ve dřevě v mezích vody vázané s obecně-ortotropním charakterem materiálu • návrh a vývoj aplikací pro simulaci vlhkostního a teplotního pole v řezivu během jeho teplovzdušného sušení nebo ve vybraných konstrukcích dřevostavby tak, aby tyto aplikace mohli využívat i uživatelé bez hlubších znalostí matematiky a programování • rozšíření modelu o mechanickou analýzu (deformace a napětí), která je nezbytná při optimalizaci parametrů sušícího prostředí • stochastické modelování vlhkostního pole ve dřevě během sušení • stochastické modelování teplotního pole v rámové konstrukci dřevostavby • měření transpiračního proudu a tepelné vodivosti dřeva lineárním ohřevem v bělové části kmene stromu • měření transpiračního proudu a tepelné vodivosti dřeva objemovým ohřevem v bělové části kmene stromu
11
3 Literární přehled Pohyb vody ve dřevě lze rozdělit na objemový tok (popisován nejčastěji Darcyho zákonem) a molekulární nebo tzv. difuzní tok (popisován nejčastěji Fickovým zákonem). Teplo se šíří vedením (Fourierův zákon), prouděním nebo sáláním. Všechny zmiňované způsoby šíření tepla a vody ve dřevě mohou nastat ve všech možných kombinacích, tedy i všechny najednou. V této práci je u tzv. „mrtvého“ dřeva pro účely sušení řeziva nebo stavební fyziky pozornost zaměřena na popis fyzikálních polí ve dřevě s vlhkostí pod mezí hygroskopicity (tedy v oblasti vody vázané, kdy dochází pouze k difuzi vody ve dřevě). U tzv. „živého“ dřeva stromu je pro účely měření transpiračního proudu a tepelné vodivosti v běli kmene stromu brána do úvahy i konvekce (proudění tepla), ale pouze z makroskopického (fenomenologického) hlediska. 3.1 Deterministické modelování fyzikálních polí ve dřevě Vývoj matematických modelů popisující teplotní a vlhkostní pole (pod mezí hygroskopicity) ve dřevě byl založen na makroskopickém popisu těchto dvou fyzikálních polí (Avramidis et al. 1994). Tento makroskopický popis byl odvozen ze zmiňovaného Fourierova zákona pro teplotu a Fickova zákona pro vlhkost (Babiak 1995) a také ze zákonů termodynamiky (Stanish et al. 1985, Plumb et al. 1984, Beard et al. 1983, Thomas et al. 1980, Adesanya et al. 1988, Bramhall 1979a, 1979b, Liu 1990, Kayihan 1986, Avramidis a Siau 1987, Avramidis et al. 1992, Barbour a Johnson 1989). Termodynamický
popis
obou
fyzikálních
polí
je
vhodný
při
sestavování
multifyzikálních (vázaných) matematických modelů. Jako první pro modelování pohybu vody ve dřevě byly použity difuzní modely (Sherwood 1929) a jsou stále běžně uplatňovány (Simpson 1993, Bramhall 1995). U listnatých dřev, lze k vůli jejich nepropustnosti pro tekutiny v příčném směru využít difuzních modelů pro šíření vlhkostního pole u dřev s vlhkostí pod i nad mezí hygroskopicity (Horáček 2004). U jehličnatých dřev tento předpoklad zaujmout nemůžeme (nelze zanedbat objemový tok tekutin v makrokapilárním systému dřeva), tedy difuzní modely se mohou při sušení použít u těchto dřev pouze s počáteční vlhkostí pod mezí hygroskopicity. Difuzní modely vycházející z gradientu vlhkosti (koncentrace 12
vody ve dřevě), jako konkrétní hybné síly celého difuzního děje, jsou nejrozšířenější (Droin et al. 1988, Droin-Josserand et al. 1988, 1989a, 1989b, Vergnaud 1991). Další hybnou silou difuze může být parciální tlak vodních par ve dřevě (Nelson 1986a, 1986b, 1986c, Siau 1983b, 1984a, 1984b, Cloutier et al. 1992). Pro teoretický popis multifyzikálních jevů se termodynamické modely zdají být nejužitečnější. Luikov (1966) a také Whitaker (1977) objevili způsob, jak popsat sdružený tepelný a vlhkostní pohyb při sušení dřeva, který je založen na ireverzibilních termodynamických procesech. Potíže při matematickém popisu vázaného šíření vlhkosti a tepla jsou vyvolány počtem fyzikálních mechanismů vysvětlující tento jev, vzájemnými závislosti mezi nimi a mnoha proměnnými, které s tím souvisejí. V případě šíření vlhkostního a teplotního pole bylo provedeno několik experimentů (Avramidis et al. 1994), ale doposud nebyly všechny navržené modely dostatečně verifikovány. Mnoho prací se zabývalo modelováním vlhkostních a tepelných toků v materiálech jako jsou polymery, dřevo nebo zemědělské produkty (Salin 1991, Kamke and Vanek 1994). Modely vlhkostních toků v závislosti na nestacionárních neisotermálních podmínkách nebyly experimentálně verifikovány, ačkoliv teorie se vyskytla už v práci (Siau 1983a). Vlhkostní a tepelný tok by měl být uvažován jako vázaný proces zohledňující: gradientem teploty indukovaný vlhkostní tok tzv. Soretův efekt (Siau 1984a, Avramidis et al. 1994) a tepelnou spotřebu vycházející z difuze vlhkosti tzv. Dufourův efekt (Siau 1992). Zmíněný termodynamický přístup tedy umožňuje v rovnici pro vlhkostního pole popsat fyzikální jev zvaný termodifuze. Další multifyzikální analýzy teplotního a vlhkostního pole jsou popsány v pracích (Voigt et al. 1940, Siau a Babiak 1983, Siau et al. 1986, Horáček 2004). Rychlost a kvalita teplovzdušného sušení jsou především ovlivňovány teplotou, vlhkostí a rychlostí proudění vzduchu, což jsou primární faktory, které rozhodují o průběhu sušení dřeva. Průběh sušení je však ovlivňován i vlastnostmi dřeva, které jsou většinou nezanedbatelně závislé na teplotě a vlhkosti. Za předpokladu konstantní hodnoty materiálových koeficientů mají rovnice lineární charakter a jejich řešení se podstatně zjednoduší. Ovšem je to mnohdy na úkor přesnosti řešení a dosažené výsledky se nemusí shodovat se skutečnými stavovými hodnotami. U dřeva je potřeba tyto materiálové koeficienty vyjádřit v závislosti na hodnotách stavových veličin (vlhkosti, teplotě). K popisu koeficientů difuze vody ve dřevě je možno použít vztahy z 13
Siau (1995) a Skaar (1988), pro koeficienty tepelné vodivosti MacLean (1941), pro měrné teplo Skaar (1988) a pro hustotu Kollmann (1951). Hodnoty teplotních a vlhkostních koeficientů pro numerické simulace sušícího procesu byly diskutovány v Soderstrom a Salin (1993), Avramidis et al. (1994), Siau (1995), Pang (1996), Plumb et al. (1985). Samotný materiál je nejčastěji považován za homogenní kontinuum, protože reálná stavba dřeva analytické i numerické postupy znesnadňuje. U tohoto homogenního kontinua se musí zohlednit anizotropní charakter dřeva. Materiálové koeficienty, které specifikují konkrétní materiál, musí být určeny pro každý anatomický směr dřeva. Tyto anatomické (hlavní) koeficienty se dají použít pouze u speciálně ortotropních těles. U těles s obecně ortotropním charakterem, jakými většinou reálná tělesa jsou, je nutné pomocí transformace další koeficienty dopočítat (Trcala 2012). Sušení dřeva je stále ve velké míře řízeno empirickými zkušenostmi. Z tohoto důvodu se mnozí výzkumníci už několik let snaží vyvinout model, kterým by kompletně popsali sušící proces na teoretické i praktické úrovni. Tyto modely se využívají k popisu sušícího procesu a nalezení způsobu, kterým by šla zvýšit kvalita sušení při současném snížení energetické spotřeby a času. Vývoj takových modelů je založen na studii tří fundamentálních jevů (přenos tepla, šíření vlhkosti a mechanické deformace). Řada teoretických modelů byla navržena a analyzována (např. Chen a Pei 1989, Irudayaraj et al. 1990), ale modely se nedají kompletně aplikovat tak, aby simulovaly reálný sušící proces. Důvodem je nízká komplexnost předkládaných modelů. Tato komplexnost je dána
širokým
rozsahem
vázaných
polí,
specifikací
okrajových
podmínek,
experimentálním vyhodnocováním materiálových charakteristik atd. (Koňas 2008). Dřevo při změnách vlhkosti mění své rozměry a také tvar. Těmto změnám se říká vlhkostní deformace a lze je popsat pomocí experimentálně zjišťovaných koeficientů vlhkostní deformace. V důsledku těchto deformací vznikají ve dřevě napětí, které lze pomocí konstitutivních vztahů z těchto deformací vypočítat. Pro malé deformace a předpoklad lineární závislosti napětí na těchto deformacích lze využít obecně známého Hookova zákona. V případě nelineární závislosti napětí na deformaci (plastické chování, viskoplastické chování apod.) je situace mnohem komplikovanější, ale numerická matematika již má nástroje, které tento problém dokáží vyřešit (například pomocí Newton-Raphsonovy metody). Výpočty značně komplikují i geometrické 14
nelinearity, které se vyskytují při velkých deformacích bez ohledu na to, zda je vztah mezi napětími a deformacemi lineární nebo nelineární. I tento problém je po matematické stránce vyřešen a existují softwary, které mají tyto matematické algoritmy implementovány. 3.2 Stochastické modelování fyzikálních polí ve dřevě Přesnost matematických modelů je dána, jak už bylo řečeno, zahrnutím všech fyzikálních jevů (multifyzikální úlohy), ale také přesností zadávaných materiálových vlastností. Hodnoty materiálových vlastností je obtížné přesně určit a v některých případech je to přímo vyloučeno. Bohužel materiálové vlastnosti dřeva patří k těmto případům, kdy nelze přesně předvídat hodnotu fyzikální nebo mechanické vlastnosti dřeva. Většinou je k dispozici jen průměrná hodnota materiálové vlastnosti a zřídka na internetu nebo v literatuře najdeme další důležité charakteristiky jako je směrodatná odchylka, typ pravděpodobnostního rozložení dané materiálové vlastnosti, přestože lze tyto údaje získat ze stejných dat, ze kterých byla vypočítána průměrná hodnota. Tyto informace
jsou
velmi
důležité,
protože
umožňují
kvantifikovat
nejistotu
v materiálových vlastnostech a tím umožňují tuto nejistotu v matematických modelech zohlednit a jako výsledek z těchto modelů získat nejenom informaci o „průměrném chování reálného systému“, ale navíc informaci o pravděpodobnostním chování systému pro všechny možné stavy výsledků simulací (předpovědí), které mohou nastat. Tímto přístupem je pak celá analýza systému zpřesněna a lze tak dospět k důkladnějšímu zhodnocení výstupů tohoto systému. Kvantifikace nejistot a jejich šíření skrze fyzikální systémy hraje důležitou roli při zlepšování vypovídací schopnosti (predikci) těchto systémů (Nouy 2009). Standardním přístupem pro řešení těchto problémů byla a stále je metoda Monte Carlo (MC). Tato metoda vzorkuje vstupní náhodné veličiny, pro každý vzorek vstupu počítá výstup systému a tyto výstupy pak statisticky vyhodnotí (Helton a Davis 2003). Tato metoda je stále široce využívána vzhledem k její robustnosti a jednoduché implementaci. MC metoda má však také nevýhody, hlavní z nich je pomalá konvergence. Pokud vyhodnocení výstupu pro každý vstupní vzorek je náročné (numerické řešení parciálních diferenciálních rovnic pro každý vstupní vzorek), pak vyhodnocení pro tisíce vzorků může být dosti časově nákladné. V těchto případech je vhodné použít 15
jinou metodu, která by srovnatelného výsledku dosáhla za mnohem kratší čas a byla by tedy efektivnější (Xiu a Karniadakis 2002). Bylo vyvinuto a použito několik alternativních metod, které dokáží vypočítat nejistoty ve stochastických modelech (Hosder et al. 2006, Debusschere et al. 2004, Reagan et al. 2003, Isukapalli 1999, Debusschere et al. 2004, Mathelin et al. 2005). Jednou z těchto metod je spektrální metoda, která dokáže zohlednit nejistoty na vstupu a efektivně vypočítat nejistoty ve výstupech daného matematického modelu tvořeného parciálními diferenciálními rovnicemi (Constantine 2009). V posledních dvou desetiletí je zaznamenán rostoucí zájem o tyto spektrální metody, kterými se řeší tzv. stochastické diferenciální rovnice (obyčejné i parciální). V této práci je zaměřena pozornost na spektrální metody, které využívají tzv. „polynomial chaos expansion“ hledaných náhodných výstupních veličin. Hlavní výhodou této metody jsou především nízké výpočetní náklady (hlavně časové, ale i paměťové). Ghanem a Spanos (1990) a Ghanem (1999) aplikovali tuto metodu ve svých pracích. Mathelin et al. (2004) tuto metodu použil pro stochastický popis turbulentních proudění. Xiu and Karniadakis (2003) analyzovali tok kolem válce a rozšířili metodu od původní Wienerovi formulace (1938), aby mohli zahrnout různé bázové funkce, které tvoří daný pravděpodobnostní prostor. Walters (2003) aplikoval „polynomial chaos expansion“ metodu na dvoudimenzionální problém vedení tepla, ve kterém se nejistoty vyskytovaly v geometrii. Stochastické Eulerovy rovnice byly řešeny v pracích (Perez a Walters 2005). Polynomial chaos expansion metoda je založena na rozkladu náhodných veličin vyskytující
se
v matematickém
modelu
do
bázových
funkcí
daného
pravděpodobnostního prostoru. Koeficienty tohoto spektrálního rozkladu jsou neznámé funkce, které budou nalezeny pomocí Galerkinovy projekce. Po provedení známé Galerkinovy projekce postupně na každou z těchto bázových funkcí dostaneme ze stochastické parciální diferenciální rovnice soustavu vázaných deterministických parciálních diferenciálních rovnic, kterou lze řešit standardními numerickými metodami (například metodou konečných prvků).
16
3.3 Měření transpiračního proudu a tepelné vodivosti dřeva v kmenu stromu K měření transpiračního proudu jsou většinou používány metody založené na různých variantách termodynamického principu. Mnohem méně často jsou používány metody založené na jiných principech, např. magnetohydrodynamickém, elektrickém, nukleární magnetické rezonanci, využívající radioaktivních izotopů apod. Mezi nejčastěji ve světě používané termodynamické metody patří např. metoda tepelného pulzu (HPV), kdy je měřena postupná rychlost proudu dle pohybu tepelné vlny v krátkodobě zahřáté části kmene (Huber 1932, Huber a Schmidt 1936, Marshall 1958, Swanson 1971, Green a Clothier 1988) a metoda dissipace tepla (HD), vycházející z úměry mezi teplotou zahřívaného čidla a hustotou proudu (Granier 1985). Mezi metody, které poskytují kvantitativní údaje aniž by vyžadovaly jakoukoli kalibraci patří metoda tepelné bilance kmene (Daum 1967, Sakuratani 1981) a metoda deformace teplotního pole v kmeni, HFD (Naděždina a Čermák 1998 b, Naděždina et al. 1998 a). Metoda deformace teplotního pole (heat field deformation - HFD) je založena na měření teplotního pole kolem lineárního ohřevu (kolem vyhřívané jehly) vloženému do kmene v radiálním směru (Nadezhdina et al. 1998a, 2002a, 2006; Čermák et al. 2004). Přední pohled na HFD zapojení s teplotním polem vypadá jako symetrická elipsa vzhledem k různým tepelným vodivostem kmene v příčném a podélném směru za podmínky nulového toku. Kmen je brán jako fyzikální kontinuum, ale při stanovení jeho tepelných vlastností je zohledňována skutečnost, že je běl kmene tvořena dřevní substancí, vodou a vzduchem. Vzorec HFD metody pro výpočet transpiračního proudu je odvozen ze závislosti tvaru teplotních isočar na aktuální hodnotě transpiračního proudu (na různé deformaci teplotního pole v důsledku různého transpiračního proudu). Tedy tvar elipsy, který má isočára teplotního pole za předpokladu nulového toku se v případě rostoucího toku více a více deformuje. Vzorec vyhodnocuje transpirační proud ze dvou měřených teplotních diferencí (symetrické vertikální dTsym a horizontální asymetrické dTasym) a vychází z jejich poměru. Tyto diference jsou měřeny dvěma páry termočlánků vložené (podobně jako ohřev) do injekční jehly z nerezové oceli o průměru 1.5 mm. Reverzní toky spojené s vodní redistribucí (Daum 1967; Burgess et al. 1998) mohou být touto metodou také měřeny (Nadezhdina et al. 2009, 2010) stejně jako noční resaturační toky (Kunia 1955; Daum 1967; Burgess et al. 2000a; Brooks et al. 2002), které poskytují cenné informace o vodním stavu rostlin. U velkých stromů může být 17
HFD metoda použita, jen pokud aplikujeme multibodová teplotní čidla, které umožňují měření radiálního profilu transpiračního proudu (Čermák a Nadezhdina 1998; Nadezhdina a Čermák 1998a, 1998b; Nadezhdina et al. 2009, 2010). Nevýhodou této metody je (jako u všech metod s lineárním ohřevem) nejistota při měření velkých transpiračních proudů. Multibodová teplotní čidla umožňují vytvoření 3D obrazu transpiračního proudu v závislosti na umístění v běli (Čermák et al. 2004). Jednobodová teplotní čidla mohou být použita u metod založených na lineárním ohřevu pouze v případě malých kmenů, kořenů nebo větví (nad 10 mm v průměru). Metoda tepelné bilance kmene s přímým elektrickým ohřevem pletiv a vitřním měřením teploty (metoda THB) byla na lesnické fakultě v Brně původně vyvinuta pro velké stromy (Čermák a Deml 1974; Čermák et al., 1973, 1976, 1982; Kučera 1977; Kučera et al., 1977). Vychází z principu, že určité místo kmene, které je známým výkonem ohříváno je současně transpiračním proudem ochlazováno. Vybraná sekce kmene stromu je zahřívána zevnitř jouleovým teplem, které se uvolňuje při průchodu od země odděleného střídavého (aby nedocházelo k elektrolýze) elektrického proudu xylemem. Teplo je uvolňováno uvnitř kmene poměrně rovnoměrně a nedochází ke ztrátě energie ohřevem skrze krycí pletiva (především borku). Elektronické obvody jsou schopny udržovat konstantní přiváděný výkon P (který je přímo úměrný intenzitě proudu) nebo teplotní rozdíl dT mezi zahřívanou a nezahřívanou částí měřiště (která je nepřímo úměrná intenzitě proudu), zatímco změny druhé veličiny jsou zaznamenávány. (jiné způsoby regulace výkonu, (např. jeho změny s denní dobou) by mohly vést k významným chybám v důsledku tepelné inerce pletiv. THB metoda je robustní, poskytuje spolehlivá data po dlouhou dobu měření u stromů s průměrem kmene přes 15 cm pro různé dřeviny a podmínky prostředí a má dobré dynamické vlastnosti (Čermák a Kučera 1991), ikdyž v malé míře může výsledek ovlivnit teplotní setrvačnost. Tato metoda přibližně integruje radiální profil transpiračního proudu a byla používána jako standard pro testování jiných metod (Offenthaler and Hietz 1998; Schubert 1999).
18
4 Metodika Celou práci lze shrnout do následujících bodů:
1. Modelování fyzikálních polí ve dřevě během sušení nebo v dřevostavbách
a) deterministické modelování •
sestavení matematického modelu popisující vázané šíření vlhkosti a teploty ve dřevě v mezích vody vázané s obecně-ortotropním charakterem materiálu (Trcala a Koňas, 2012)
•
návrh a vývoj aplikací vycházející z výše odvozeného modelu a popisující vlhkostní a teplotní pole v řezivu během jeho teplovzdušného sušení nebo ve dřevě zabudovaném ve stavbě tak, aby je mohl použít uživatel bez hlubších znalostí matematiky a programování
•
rozšíření modelu o mechanickou analýzu (deformace a napětí), která je nezbytná pro optimalizaci parametrů sušícího prostředí (Trcala a Koňas, submitted 2012)
b) stochastické modelování •
stochastické modelování difuzních problémů při sušení dřeva (Trcala, submitted 2012a)
•
stochastické
modelování
teplotních
problémů
v rámové
konstrukci
dřevostavby (Trcala, submitted 2012b)
2. Modelování fyzikálních polí ve dřevě při měření transpiračního proudu a tepelné vodivosti v běli kmene stromu •
měření s lineárním ohřevem (Trcala a Čermák, submitted 2011)
•
měření s objemovým ohřevem (Trcala a Čermák, submitted 2012)
19
4.1 Modelování fyzikálních polí ve dřevě během sušení nebo v dřevostavbách 4.1.1 Deterministické modelování fyzikálních polí ve dřevě 4.1.1.1 Modelování vázaného teplotního a vlhkostního pole ve dřevě Jde o úlohy modelování a simulace teplotního a vlhkostního pole ve dřevě při jeho sušení nebo při jeho zabudování do obvodového pláště dřevostavby. Důraz je kladen na provázanost těchto dvou fyzikálních polí prostřednictvím Soretova efektu, Duforova efektu a závislosti materiálových vlastností na teplotě a vlhkosti. Je zohledněna anizotropie dřeva a jsou odvozeny transformační vztahy pro obecně-ortotropní případ modelu. Tímto způsobem je odvozena soustava parciálních diferenciálních rovnic s příslušnými okrajovými a počátečními podmínkami pro popis vázaného šíření vlhkostního a teplotního pole ve dřevě. Pro reálný popis distribuce zmíněných fyzikálních polí není možné použít Fickova a Fourierova zákona izolovaně. Difuzní přenos tepla a hmoty je možné, za předpokladu zavedení teorie kontinua a zohlednění anisotropie vlastností dřeva, popsat soustavou parciálních diferenciálních rovnic. V rovnicích vystupují jednak koeficienty z Fickova a Fourierova zákona, také ale koeficienty související například s termodifuzí (Horáček 2004). Tyto koeficienty byly odvozené z teorie termodynamické rovnováhy chemické (koncentrace vody-vlhkosti) a fyzikální (teplotní). Vycházelo se z abstraktní soustavy Luikovových rovnic se zahrnutím termodifuze a tepelných ztrát daných potřebou aktivační energie pro pohyb vody vázané. Odvození obou koeficientů není náplní této práce a je podrobněji popsáno v práci Horáčka (2004). Výsledný tvar koeficientu termodifuze vody ve dřevě (tzv. Soretův efekt) (Avramidis et al. 1994, Horáček 2004) lze zapsat v následujícím tvaru:
s=
ϕ ∂EMC Eb RT ∂ϕ T ,
(4.1)
kde EMC [-] je rovnovážná vlhkost dřeva (equilibrium moisture content), která je funkcí
teploty
T
[K]
a
relativní
vlhkosti
vzduchu
φ
[-],
EMC (ϕ , T ) = 1 / B ln( A / ln(1 / ϕ )) , ∂EMC / ∂ϕ je tedy směrnice sorpční izotermy,
Eb = 38500 − 29000 M
[Jmol-1]
je 20
aktivační
energie
vody
vázané, A = 7 .731706 − 0 .014348 T ,
B = 0 .008746 + 0 .000567 T
jsou
de Boer-
Zwickerovy koeficienty a R [Jmol-1K-1] je univerzální plynová konstanta. Obecně, tyto izotermy jsou specifické pro daný druh dřeva, protože diferenciální teplo sorpce se pro různé druhy dřev liší (Babiak 1990).
Příslušná soustava rovnic potom vypadá následovně: ∂M − ∇ ⋅ ( D∇ M + sD∇T ) = 0 ∂t
ρC
(4.2)
E ∂M ∂T − b − ∇ ⋅ λ∇T = 0 ∂t 1.8C ∂t
kde D je matice difuzních koeficientů, λ je matice koeficientů tepelné vodivosti, s je Soretův efekt [-/K], sD je matice difuzních koeficientů D vynásobená Soretovým efektem s, M je vlhkost [-], T je teplota [K], C je měrná tepelná kapacita [Jkg-1K-1], ρ je hustota dřeva [kgm-3]. V modelu jsou zahrnuty okrajové podmínky 3. řádu. Tyto podmínky zachycují tok vlhkosti nebo tepla na povrchu tělesa. Zachycují tok ve směru kolmém na tento povrch a v podstatě popisují vypařování vody z povrchu tělesa nebo ohřev povrchu tělesa. Hustota difuzního toku je počítána z rozdílu mezi povrchovou vlhkostí a rovnovážnou vlhkostí dřeva odpovídající okamžitým okolním parametrům vzduchu. Matematicky jde o lineární vztah mezi hustotou vlhkostního toku a zmíněným rozdílem povrchové a rovnovážné vlhkosti. Vztah vypadá následovně: −n ⋅ ( D∇M + sD∇T ) = α M ( M ∂Ω − EMC ) ,
(4.3)
kde α M je koeficient přestupu vlhkosti [ms-1], M ∂Ω je vlhkost na povrchu tělesa [-], EMC je rovnovážná vlhkost dřeva [-]. Nejdůležitějším a také nejtěžším úkolem je přesné stanovení koeficientu přestupu vlhkosti α M (Siau 1992, Avramidis et al. 1994, Horáček 2004). Ten závisí na spoustě parametrů (vlhkosti, teplotě, tvaru a rozměrech tělesa, apod.).
21
Hustota tepelného toku je počítána analogicky (vychází z Newtonova zákona ochlazování) a má následující tvar: −n ⋅ λ∇T = αT (T∂Ω − Tair ) ,
(4.4)
kde αT je koeficient přestupu tepla [Wm-2K-1], T∂Ω je teplota na povrchu tělesa [-], Tair je teplota okolního vzduchu [K]. I zde je nejobtížnější stanovit koeficient přestupu. Ten závisí na spoustě jiných parametrech a všechny tyto závislosti jsou v této práci zohledněny a zahrnuty do modelu. Počáteční podmínky jsou zadány jako funkce polohy v tělese. V této práci jsou uváženy počáteční podmínky udávající rovnoměrně rozloženou vlhkost i teplotu po celém objemu tělesa.
M=M0
(4.5)
T=T0
Mohli by být uváženy i nekonstantní podmínky. Potom by šlo o počáteční podmínky zadané jako funkční závislosti.
M = f M ( x, y , z )
(4.6)
T = fT ( x, y, z )
Běžně je teplota dřeva před sušením v celém objemu rozložena rovnoměrně, ale vlhkost nikoliv.
22
4.1.1.2 Simulace teplotního a vlhkostního pole ve dřevě Pro výše uvedený model jsou vyvíjeny dvě aplikace, které uživatel může použít pro simulaci teplotního a vlhkostního pole ve dřevě při jeho sušení nebo ve dřevě zabudovaném ve stavbě. Obě aplikace jsou vytvářeny tak, aby je mohl použít uživatel bez hlubších znalostí matematiky a programování. Použití aplikací je velmi intuitivní a uživatel musí znát jen základní vstupní údaje, jako jsou rozměry, okolní teploty a vlhkosti vzduchu, počáteční teploty a vlhkosti a v případě stavební fyziky jen základní fyzikální vlastnosti (hustota, měrné teplo, tepelná vodivost, difuzní vodivost, apod.) materiálů dané konstrukce. Zde je základní jednoduché metodické schéma při vytváření aplikací:
Fig. 4.1 Schéma postupu při vytváření aplikací pro simulaci teplotního a vlhkostního pole ve dřevě při jeho sušení nebo ve dřevě zabudovaném do stavební konstrukce.
23
Simulace teplotního a vlhkostního pole ve dřevě během sušení Aplikace je zaměřena především na analýzu procesu teplovzdušného sušení dřeva. Jde o aplikaci, kterou může uživatel využít bez větších znalostí fyziky, matematiky a programování k tomu, aby simuloval rozložení a průběh vlhkosti v čase uvnitř dřeva, které je vystaveno nějakému sušícímu režimu nebo obecně jakýmkoli parametrům okolního vzduchu. Aplikace je vytvořena v programu MATLAB a COMSOL Multiphysics. V programu COMSOL Multiphysics je zajištěno řešení matematického modelu popisující daný fyzikální problém sušení, v programu MATLAB je vytvořeno grafické uživatelské rozhraní (GUI) aplikace, zajištěna parametrizace úloh a obohacuje aplikaci o nezbytné programátorské prvky.
Simulace teplotního a vlhkostního pole v konstrukci dřevostavby Je
vyvíjena
aplikace
pro
simulaci
teplotního
a
vlhkostního
pole
ve
vybraných konstrukcích dřevostavby. Je zde snaha o sestavení a implementaci modelu do aplikace, který bude popisovat vázané šíření tepla a vlhkosti v konstrukci dřevostavby obecně zapsaném v následujícím tvaru
C pp
∂p ∂T + C pT − ∇ ⋅ (K pp ∇p + K pT ∇T ) = 0 ∂t ∂t ,
CTp
∂p ∂T + CTT − ∇ ⋅ (K Tp ∇p + K TT ∇T ) = 0 ∂t ∂t
(4.7)
kde T je teplota, p je parciální tlak vodních par. Pro dřevo je problém vázaného šíření tepla a vlhkosti řešen rovnicemi (4.2) a je otázkou, jak tento model použít v rámci celé konstrukce, kde se vyskytují i nedřevěné materiály s velkými vzduchovými póry, ve kterých se šíření vlhkosti popisuje pomocí gradientu parciálního tlaku vodních par a sleduje se zároveň i tlak nasycených vodních par. Pokud hodnota parciálního tlaku vodních par přesáhne v některých místech konstrukce hodnotu tlaku nasycených vodních par, tak v tomto místě dojde ke kondenzaci vodních par. Kondenzace vodních par je nežádoucím jevem ve stavebních konstrukcích a je třeba, aby byla při návrzích konstrukcí vzhledem k daným 24
exteriérovým a interiérovým podmínkám posouzena a vyloučena pomocí co nejdůvěryhodnějších simulací. Simulace mohou být založeny na matematickém modelu, který principiálně vychází z následující teorie. Hustota difuzního toku může být vyjádřena v závislosti na gradientu koncentrace a také na gradientu parciálního tlaku vodních par. V případě gradientu koncentrace se ve vztahu vyskytuje tzv. difuzní koeficient Dx a v případě gradientu parciálního tlaku vodních par tzv. součinitel difuzní vodivosti δx. U dřeva je nám známa funkční závislost difuzního koeficientu Dx na teplotě a vlhkosti a bylo vhodné toto zohlednit při simulaci teplot a vlhkostí uvnitř konstrukce obsahující dřevo jako konstrukční materiál. U ostatních stavebních materiálů jsou však především známy hodnoty součinitelů difuzní vodivosti δx, protože kvůli možnosti kondenzace vodních par se v konstrukci sleduje parciální tlak vodních par p a jako hybné síly difuze se využívá gradientu této veličiny. Cílem této části práce tedy bylo vyjádřit u dřeva součinitel difuzní vodivosti ze známého difuzního koeficientu. Pro hustotu difuzního toku jx platí následující vztahy, ze kterých lze vyjádřit součinitel difuzní vodivosti dřeva δx ze známého difuzního koeficientu Dx: ∂M ∂M ∂p ∂M ∂p = δx , jx = −δ x , ρ r Dx , δ x = ρ r Dx ∂p ∂x ∂x ∂x ∂x
(4.8)
p ∂ ∂M ∂M ∂ϕ ∂M psat ρ ∂M = ρ r Dx = ρ r Dx = Dx r δ x = ρ r Dx ∂p ∂ϕ ∂p ∂ϕ ∂p psat ∂ϕ
(4.9)
jx = − ρ r Dx
a následně hledaný vztah pro hustotu difuzního toku
j x = − Dx
ρ r ∂M ∂p ρ ∂M j= − r D ∇p . , psat ∂ϕ ∂x psat ∂ϕ
(4.10)
Základní diferenciální rovnice pro popis vlhkostního pole ve dřevě vypadá následovně
∂c −∇⋅ j = 0 ∂t
(4.11)
kde 25
p ∂ ∂p ∂c ∂M ∂M ∂ϕ ∂M psat ∂M 1 ∂p (4.12) = ρr = ρr = ρr = ρr psat − sat p 2 ∂t ∂t ∂ϕ ∂t ∂ϕ ∂t ∂ϕ psat ∂t ∂t Redukovaná hustota se vykrátí a dostaneme následující tvar nestacionární rovnice
1 ∂M ∂p ∂M 1 ∂p psat − sat p − ∇ ⋅ D∇p = 0 , 2 ∂ϕ psat ∂t ∂t psat ∂ϕ
(4.13)
po úpravě a dopsání rovnice popisující teplotní pole dostaneme následující soustavu dvou rovnic ∂M 1 ∂p ∂M 1 ∂psat − 2 ∂ϕ psat ∂t ∂ϕ psat ∂t
ρC
1 ∂M D ∇p = 0 p − ∇ ⋅ p ∂ ϕ sat
(4.14)
∂T − ∇ ⋅ λ∇T = 0 . ∂t
Problém lze také vyřešit tak, že u dřeva se bude řešit M a jinak se bude řešit p a na hranicích se tato dvě fyzikální pole provážou následujícími okrajovými podmínkami:
j = −n ⋅ δ∇p = −n ⋅ D∇c = −n ⋅ D∇ ( ρr M )
.
(4.15)
4.1.1.3 Modelování vázaného teplotního, vlhkostního a deformačního pole ve dřevě V tomto bodě jsou ke stávajícím teplotně-vlhkostním rovnicím přidány rovnice popisující deformace a napětí ve dřevě, které vznikají při změnách vlhkosti a teploty. Vliv gradientu teploty na deformace dřeva je většinou zanedbatelný. Simulace je provedena na konkrétním experimentálním vzorku v programu Comsol Multiphysics s podporou programu MATLAB. Pro popis deformací a napětí je třeba definovat složky posunutí ve směru os x, y, z vektorem ( u , v, w ) . Deformace jsou potom popsány symetrickým tensorem obsahující tři normálové složky ( ε x , ε y , ε z ) a v případě symetrie tři smykové složky ( ε xy , ε yz , ε xz ) .
26
Napětí jsou analogicky k deformacím popsány symetrickým tensorem o třech normálových (σ x , σ y , σ z ) a třech smykových složkách (τ xy ,τ yz ,τ xz ) . Vztah mezi napětími a deformacemi je v této práci popisován tzv. Hookovým zákonem, který předpokládá malé deformace a lineární vztah mezi napětími a deformacemi. Zaveďme následující značení
σ x σ y σ σ = z , τ xy τ yz τ xz
∂u ∂x ∂v ∂y εx ∂ w ε y ∂z εz ε = = 1 ∂u ∂v , ε xy + ε yz 2 ∂y ∂x ∂ v ∂ w 1 ε xz + 2 ∂z ∂y 1 ∂w ∂u + 2 ∂x ∂z
ε Mx ε Tx CMET CTET ε My ε Ty CMER CTER ε ε Tz CMEL CTEL ε M = Mz = ( M − M 0 ) = (T − T0 ) , εT = ε Mxy ε Txy 0 0 ε Myz ε Tyz 0 0 0 0 ε Mxz ε Txz
(4.16)
(4.17)
kde εM je vektor vlhkostních deformací, ε T je vektor teplotních deformací, CMET, CMER,
CMEL jsou koeficienty vlhkostní deformace a CTET , CTER , CTEL jsou
koeficienty teplotní deformace. Potom podle Hookova zákona platí následující vztah σ = C(ε − ε M − ε T ) ,
(4.18)
kde 27
E x (1 − µ yz µ zy ) ∆ E y ( µ xy + µ xz µ zy ) ∆ C = ( Cij ) = Ez ( µ xz + µ xy µ yz ) i , j =1,..,6 ∆ 0 0 0
E x ( µ yx + µ zx µ yz )
Ex ( µ zx + µ yx µ zy )
∆ E y (1 − µ zx µ xz )
∆ E y ( µ zy + µ zx µ xy )
∆ Ez ( µ yz + µ xz µ yx )
Ez (1 − µ xy µ yx )
0
0
∆
0
0
∆ 0
∆ 0
0
0
2Gxy
0
0
0
0
2G yz
0
0
0
0
∆ = 1 − µ xy µ yx − µ yz µ zy − µ zx µ xz − 2 µ xy µ yz µ zx
0 0 0 0 0 2Gxz
(4.19)
je Hookova matice s mechanickými vlastnostmi material jako jsou moduly pružnosti Ex, Ey, Ez a Poissonovy čísla µ kl, kde indexy k,l=x, y, z. Výsledná soustava tří parciálních diferenciálních rovnic pro neznámé tři složky posunutí ( u , v, w ) vypadá následovně: ∂ 2u ∂σ x ∂τ xy ∂τ xz − − = Fx ρ 2 − ∂t ∂x ∂y ∂z ∂ 2 v ∂τ xy ∂σ y ∂τ yz − − = Fy ρ 2− ∂t ∂x ∂y ∂z , 2 ∂τ ∂ w ∂τ ∂σ ρ 2 − xz − yz − z = Fz ∂t ∂x ∂y ∂z
(4.20)
kde F=(Fx, Fy, Fz) značí vnější silové působení na těleso.
28
4.1.2 Stochastické modelování fyzikálních polí ve dřevě 4.1.2.1 Stochastické modelování difuzních problémů při sušení dřeva Parciální diferenciální rovnice popisující nestacionární vlhkostní pole ve dřevě vypadá následovně: ∂M − ∇ ⋅ D∇ M = 0 , ∂t
(4.21)
kde D je matice difuzních koeficientů [m2s-1] and M je vlhkost [-]. Okrajové podmínky jsou dány přestupem vlhkosti z povrchu materiálu do okolního vzduchu. Pro numerické řešení byla použita následující okrajová podmínka: −n ⋅ D∇M = α M ( M ∂Ω − EMC ) ,
(4.22)
kde α M koeficient přestupu vlhkosti [ms-1], M ∂Ω je vlhkost povrchu sušeného materiálu [-], EMC je rovnovážná vlhkost dřeva [-].
Použitím polynomial chaos expansion je vlhkostní pole vyjádřeno následujícím způsobem (je rozloženo do bázových funkcí Φ i (ξ ) pravděpodobnostního prostoru) Q
M ( x, t , ξ (ω )) = ∑ M i ( x, t )Φ i (ξ ) .
(4.23)
i=0
Matice difuzních koeficientů jako náhodná veličina: D → D + σξ ,
(4.24)
kde ξ je náhodná proměnná s rovnoměrným rozdělením pravděpodobnosti ξ ∈ U ( −1,1).
Dosazením výrazů (4.23, 4.24) můžeme rovnice (4.21, 4.22) přepsat do tvaru
29
Q
∂ ∑ M i (x, t )Φ i (ξ ) i =0
∂t
Q
− ∇ ⋅ ( D + σξ ) ∑ ∇M i (x, t )Φ i (ξ ) = 0
(4.25)
i =0
Q
Q
i =0
i =0
−n ⋅ ( D + σξ ) ∑ ∇M i (x, t )Φ i (ξ ) = α M (∑ M i ( x, t )Φ i (ξ ) − EMC ). Galerkinovou projekcí rovnic (4.25, 4.26) do bázových funkcí
(4.26)
{Φ
j
(ξ )} pro každé
j = {0,1,..., Q} se obdrží soustava (Q + 1) vázaných parciálních diferenciálních rovnic:
Q
∑ i =0
∂M i ( x, t ) Φ i (ξ ), Φ j (ξ ) − ∇ ⋅ D Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇M i ( x, t ) = 0 ∂t
(
)
(4.27) Q Q −∑ n ⋅ D Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇M i ( x, t ) = α M ∑ M i (x, t ) Φ i (ξ ), Φ j (ξ ) − EMC 1, Φ j (ξ ) i =0 i =0
(4.28)
Tato soustava je numericky pomocí metody konečných prvků řešena a výsledkem jsou koeficienty polynomial chaos expansion, ze kterých lze podle následujících vztahů přibližně spočítat průměr Q
Q
i=0
i =0
E[ M ( x, t , ξ (ω ))] = E[∑ M i ( x, t )Φ i (ξ )] = ∑ M i ( x, t )E[Φ i (ξ )] = M 0 ( x, t ) ,
(4.29)
rozptyl
Q
Q
Var[ M (x, t , ξ (ω ))] = E[(∑ M i (x, t )Φ i (ξ ) − M 0 (x, t )) 2 ] = E[(∑ M i (x, t )Φ i (ξ )) 2 ] = i =0
i =1
Q
Q
i =1
i =1
= ∑ M i 2 (x, t )E[Φ i (ξ )2 ] = ∑ M i 2 (x, t ) Φ i (ξ ), Φ i (ξ ) (4.30) a směrodatnou odchylku výsledného vlhkostního pole
Std [ M (x, t , ξ (ω ))] =
Q
∑M i =1
2
i
(x, t ) Φ i (ξ ), Φ i (ξ ) .
30
(4.31)
.
4.1.2.2 Stochastické dřevostavby
modelování
teplotních
problémů
v
konstrukci
Parciální diferenciální rovnice popisující stacionární teplotní pole vypadá následovně:
−∇ ⋅ λ∇T = 0
(4.32)
kde λ je matice koeficientů tepelné vodivosti dřeva [Wm-1K-1] and T je teplota [K]. Okrajové podmínky jsou dány přestupem tepla mezi povrchem stěny a okolní vzduchem. Vztah mezi tepelným tokem z (do) stěny a rozdílem mezi povrchovou teplotou a teplotou okolního vzduchu je lineární: −n ⋅ λ∇T = α T (T∂Ω − Tair ) ,
(4.33)
kde αT je koeficient přestupu tepla [Wm-2K-1], T∂Ω je povrchová teplota [K], Tair je teplota okolního vzduchu [K]. Použitím polynomial chaos expansion je teplotní pole vyjádřeno následujícím způsobem (je rozloženo do bázových funkcí Φ i (ξ ) pravděpodobnostního prostoru) Q
T (x, ξ (ω )) = ∑ Ti (x)Φ i (ξ ) .
(4.34)
i=0
Matice koeficientů tepelné vodivosti je brána jako náhodná veličina λ → λ + σξ ,
(4.35)
kde ξ je náhodná veličina rovnoměrného rozdělení pravděpodobnosti ξ ∈ U ( −1,1). Dosazením výrazů (4.34, 4.35) můžeme rovnice (4.32, 4.33) přepsat do tvaru
Q −∇ ⋅ ( λ + σξ ) ∑ ∇Ti (x)Φ i (ξ ) = 0 i =0
(4.36)
Q Q −n ⋅ ( λ + σξ ) ∑ ∇Ti (x)Φ i (ξ ) = α T ∑ Ti (x)Φ i (ξ ) − Text . i =0 i =0
(4.37)
31
Galerkinovou projekcí rovnic (4.36, 4.37) do bázových funkcí
{Φ
j
(ξ )} pro každé
j = {0,1,..., Q} se obdrží soustava (Q + 1) vázaných parciálních diferenciálních rovnic:
(
Q
)
−∑ ∇ ⋅ λ Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇Ti (x) = 0 i =0
(4.38)
Q Q −∑ n ⋅ λ Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇Ti ( x) = α T ∑ Ti (x) Φ i (ξ ), Φ j (ξ ) − Text 1, Φ j (ξ ) i =0 i =0
(4.39) a pro části konstrukce s deterministickými koeficienty tepelné vodivosti
Q
(
)
−∑ ∇ ⋅ λ Φ i (ξ ), Φ j (ξ ) ∇Ti (x) = 0 i =0
(4.40)
Q −∑ n ⋅ λ Φ i (ξ ), Φ j (ξ ) ∇Ti (x) = αT ∑ Ti (x) Φ i (ξ ), Φ j (ξ ) − Text 1, Φ j (ξ ) i =0 i=0 Q
(4.41)
Tato soustava je numericky pomocí metody konečných prvků řešena a výsledkem jsou koeficienty polynomial chaos expansion, ze kterých lze podle následujících vztahů přibližně spočítat průměr
Q
Q
i=0
i =0
E[T (x, ξ (ω ))] = E[∑ Ti (x)Φ i (ξ )] = ∑ Ti (x)E[Φ i (ξ )] = T0 (x) ,
(4.42)
rozptyl Q
Q
Var[T (x, ξ (ω ))] = E[(∑ Ti (x)Φ i (ξ ) − T0 (x)) 2 ] = E[(∑ Ti (x)Φ i (ξ )) 2 ] = i=0
i =1
Q
Q
i =1
i =1
= ∑ Ti 2 (x)E[Φ i (ξ ) 2 ] = ∑ Ti 2 (x) Φ i (ξ ), Φ i (ξ )
(4.43)
a směrodatnou odchylku výsledného teplotního pole
Std [T (x, ξ (ω ))] =
Q
∑T i =1
i
2
(x) Φ i (ξ ), Φ i (ξ ) . 32
(4.44)
4.2 Modelování fyzikálních polí ve dřevě při měření transpiračního proudu a tepelné vodivosti v běli kmene stromu Modelování a simulace jevů při měření transpiračního proudu a tepelné vodivosti termodynamickými metodami jsou založeny na řešení známé parciální diferenciální rovnice popisující teplotní pole pro kondukčně-konvekční šíření tepla v anisotropním kontinuu (Tatarinov et al., 2005):
ρc
∂T = qconduction + qconvection + P , ∂t
(4.45)
kde
qconduction = ∇⋅ λ∇T , qconvection = cwQw
∂T ∂y
(4.46)
a potom dostáváme výsledný tvar rovnice
ρc
∂T ∂T = ∇ ⋅ λ∇T + cwQw +P, ∂t ∂y
(4.47)
kde T je teplota (K), Qw je hustota transpiračního proudu (kg.m-2.s-1), λ je matice koeficientů tepelné vodivosti čerstvého bělového dřeva (W.m-1.K-1), ρ je hustota
čerstvého bělového dřeva (kg.m-3), c je měrné teplo čerstvého bělového dřeva (J.kg-1.K1
), cw je měrné teplo vody (J.kg-1.K-1), P je zdroj tepla (W.m-3). Tuto rovnici lze použít pro simulaci lineárního nebo objemového ohřevu části běle
s libovolně definovanou hodnotou hustoty transpiračního proudu Qw. Tyto simulace jsou provedeny v softwarech COMSOL Multiphysics a MATLAB a jsou založeny na
řešení výše zminěné rovnice pomocí metody konečných prvků (finite element method – FEM).
4.2.1 Měření s lineárním ohřevem 4.2.1.1 Metoda tepelné bilance kolem lineárního ohřevu
Byla vyvinuta nová tzv. „Metoda tepelné bilance kolem lineárního ohřevu“, anglicky „Linear Heat Balance method – LHB method“, která při měření sap flow 33
lineárním ohřevem zohledňuje vliv koeficientu tepelné vodivosti na výsledek měření. Tato metoda je založena na parciální diferenciální rovnici, která popisuje teplotní pole v případě kondukce a konvekce tepla v běli kmene stromu. Za předpokladu ustáleného teplotního pole (přibližně lze tento předpoklad pro reálný systém termodynamického měření transpiračního proudu zaujmout) se rovnice zjednoduší na tvar ∂ 2T ∂ 2T ∂T 0 = λx 2 + λ y 2 + cwQw +P, ∂x ∂y ∂y
(4.48)
ve kterém je analyticky řešitelná, po substituci λ y = 2λx a jednoduchých úpravách dostaneme hledaný vzorec
Pl T ( x, y ) = e π 8λx
Qw cw y 4 λx
Qw cw (2 x 2 + y 2 ) ), K0 ( 4λ x
(4.49)
kde T(x,y) je teplota v ustáleném stavu (K) v bodě [x, y] v souřadnicové soustavě, která má počátek v místě lineárního ohřevu, K0(-) je modifikovaná Besselova funkce druhého
řádu. Odvozený vzorec lze použít následujícím způsobem
Tupper −
Tright + Tleft 2
Qw4cwλ0.01 Q c 0.01 Q c 0.01 2 = ) − K0 ( w w ) , e x K 0 ( w w 4λx 4λx π 8λx Pl
(4.50)
kde je již dosazena měřená teplotní diference mezi horní teplotní sondou Tupper a průměrem bočních teplotních diferencí Tright, Tleft
ve vzdálenostech 0.01 m od
lineárního ohřevu. Nyní je připraven vzorec, ze kterého lze, za předpokladu známé tepelné vodivosti, vypočítat sap flow a naopak.
34
4.2.1.2 Výpočet transpiračního toku a tepelné vodivosti současně
Z výše uvedeného máme k dispozici jednu rovnici o dvou neznámých (transpirační tok, tepelná vodivost), což vede k myšlence, zda nelze najít další rovnici, která by také dávala do vztahu tyto dvě neznámé. Tím bychom dostali soustavu dvou rovnic o dvou neznámých, která v případě nezávislosti a vzájemné konzistenci rovnic by se dala využít pro výpočet obou neznámých současně. Obecně lze tuto myšlenku matematicky zapsat následujícím způsobem: F1 (Qw , λx ) = 0 F2 (Qw , λx ) = 0
,
(4.51)
kde F1(-) a F2(-) jsou obecně nelineární funkce obsahující hledané dvě neznámé. První rovnice může být dána výše odvozeným vzorcem pro LHB metodu
F1 (Qw , λx ) = Tupper −
Tright + Tleft 2
Qw c4wλ0.01 Q c 0.01 Q c 0.01 2 − ) − K0 ( w w ) . e x K 0 ( w w 4λ x 4λ x π 8λx
Pl
(4.52) Druhá rovnice by mohla být dána modifikovanou verzí již známého vzorce pro výpočet sap flow na principu deformace teplotního pole kolem lineárního ohřevu tzv. HFD metoda. Modifikovaná verze proto, protože nezahrnuje koeficient tepelné vodivosti. Modifikace vzorce z HFD metody byla provedena pomocí numerických simulací, tak, aby výsledná rovnice obsahovala koeficient tepelné vodivosti a zároveň byla konzistentní s rovnicí z LHB metody.
F2 (Qw , λx ) = Qw − f (Qw , λx )
(dTasym − dTsym )0 + dTsym − dTasym dTasym
= 0.
(4.53)
Lze použít vzorce z HFD metody ve standardní nemodifikované verzi pro výpočet sap flow, který se dosadí do LHB rovnice, ze které se dopočítá druhá neznámá, tedy tepelná vodivost. Tento přístup je však riskantní z pohledu nekonzistence standardního
35
HFD vzorce s LHB vzorcem a problém tak nemusí mít řešení nebo řešení nemusí odpovídat realitě.
4.2.2 Měření s objemovým ohřevem 4.2.2.1 Současný stav metody tepelné bilance (THB metody) Rovnice tepelné bilance vyhřívané části kmene je založena na objemovém ohřevu pomocí elektrod (přímý vnitřní elektrický ohřev objemu mezi elektrodami). Část takto generované tepelné energie (10-20%) je odvedena do okolního prostředí kondukcí a větší část je pohlcena a transportována proudem vody (transpiračním proudem). Pro tento případ byla odvozena následující diferenciální rovnice (Kučera et al., 1977)
P = QdTcw + dT λ + ρ cdT ′ ,
(4.54)
kde P je množství tepla vstupující do měřené části kmene (W), Q je transpirační proud (kg.s-1), ρ c = mwcw + mx cx , cw a cx je měrné teplo vody (4186.8 J.kg-1.K-1) a měrné teplo suchého dřeva (J.kg-1.K-1), mw a mx je hmotnost vody a je hmotnost suchého dřeva v měřené části kmene (kg), λ je množství tepla odvedené do okolního prostředí kondukcí (W.K-1), dT je teplotní diference (K) mezi vyhřívanou a nevyhřívanou částí a t je čas (s). Analytické řešení výše uvedené rovnice dává vzorec pro teplotní diferenci dT ve tvaru dT (t ) =
−t ( Qcw + λ ) P 1 − exp . Qcw + λ ρc
(4.55)
Pro již ustálený stav se rovnice (4.54) zjednoduší na tvar P = QdTcw + dT λ ,
(4.56)
ze kterého se jednoduše vyjádří transpirační proud Q
Q=
P λ − . dTcw cw
(4.57) 36
Tepelné ztráta (λ) není přesně specifikována, ale je zahrnuta v tzv. fiktivním toku (fictitious flow - Qfic), který je přibližně odhadován za podmínky nulového toku. S hodnotou tepelných ztrát koreluje obsah vody ve dřevě, takže znalostí jedné z této veličin se lze pak dopočítat ke druhé (za předpokladu znalosti hustoty dřeva v absolutně suchém stavu). Když počítáme proud Q, je nezbytné odečíst zmiňovaný fiktivní tok Qfic, který byl doposud zjišťován (1) uvážením minimálních hodnot Qrec pro každý jednotlivý den, toto zapříčiňuje podhodnocení nízkých toků a resaturační toky zcela eliminuje, (2) je brána minimální hodnota Qrec za více dnů (pevně stanovená časová perioda), tento přístup způsobuje většinou nadhodnocení denních minim toku. Rovnici (4.57) lze přepsat ve smyslu nyní zmíněného
Q = Qrec − Q fic ,
(4.58)
kde
P dTcw ,
Qrec =
Q fic =
λ
(4.59)
cw
kde Q fic je zmiňovaný fiktivní tok s neznámou tepelnou ztrátou λ z vyhřívaného objemu do okolí. Tento problém je v následujícím bodě řešen, jsou v něm odvozeny dvě nové varianty vzorců pro výpočet transpiračního proudu založené na přesnějším výpočtu fiktivního toku. Díky novému uspořádání termočlánků v okolí elektrod (tzv. horizontálnímu a symetrickému) lze měřit i nulový a záporný tok.
4.2.2.2 Vyjádření fiktivního toku v závislosti na tepelné vodivosti a okamžité hodnotě transpiračního toku
Fyzikální formulace Q fic odvozená ze standardních koeficientů tepelné vodivosti je následující:
Q fic =
λ cw
=
λ x S x + λ y S y + λz S z dycw
,
(4.60)
37
kde dy je nespecifikovaná teoretická vzdálenost, na které se realizuje teplotní diference dT a λx , λ y , λz resp. S x , S y , S z jsou koeficienty tepelné vodivosti (W.m-1.K-1) resp. plochy ohraničující vyhřívaný objem ve směrech x, y, z, které také nejsou přesně specifikovány (jejich hodnota je závislá na aktuální hodnotě transpiračního proudu). Vzorec je obtížné aplikovat v této podobě, protože obsahuje mnoho neznámých parametrů, ale zamyšlením se nad těmito parametry a nad fyzikálním principem vzorce (4.60) nás inspiruje k zavedení fiktivního toku jako funkce tepelné vodivosti a transpiračního proudu. Navíc pokud předpokládáme λy = 2λx a tepelné ztráty v radiálním směru jsou malé (zanedbatelné), pak lze psát
Q fic = λx
S x + 2S y dycw
.
(4.61)
Vzhledem k již zmíněnému může být fiktivní tok brán jako funkce tepelné vodivosti λx (W.m-1.K-1) a transpiračního proudu Q (kg.s-1) Q fic = f (λx , Q ) .
(4.62)
Hustota fiktivního toku Qw.fic pro hustotu transpiračního proudu Qw (kg.m-2.s-1) může být přibližně linearizována, jak ukazují numerické simulace
Qw. fic = a + bλx + cQw .
(4.63)
Potom fiktivní tok Qfic je roven
Q fic = F1 + F2Q ,
(4.64)
kde F1 = A ( a + bλx ) = (Q fic )0 a F2 = Ac je míra závislosti Q fic on Qrec , kde A přepočtová plocha měřiště, kterou při praktickém měření dále znát nemusíme (byla jen třeba v simulacích pro přepočet hustoty toku na celkový tok přes plochu vyhřívaného 38
objemu). Výsledky ze simulací ukázaly, že hodnota F2 je přibližně 0.5. Tato hodnota má vliv především na vyšší toky a určuje amplitudy odpoledních toků. Hodnota F1 je odvozena z teplotní diference dT v případě nulového toku, značena jako (dT)0 a odvozena extrapolací lineárního trendu teplotních diferencí na vztahovém grafu mezi dT (na ose y) a dTa-dTb (na ose x). Po zjištění (dT)0 pro každý jednotlivý den (od západu po východ slunce) se pro výpočet F1 použije následující vztah
F1 = (Q fic )0 = ( Qrec )0 =
P ( dT )0 cw .
(4.65)
Po substituci Q = Qrec − Q fic do (15) dostáváme rovnici Q fic = F1 + F2 ( Qrec − Q fic ) ,
(4.66)
ze které lze vyjádřit Qfic
Q fic =
F1 F + 2 Qrec . 1 + F2 1 + F2
(4.67)
Pokud nezohledníme závislost Q na Qfic (F2=0) tak dostaneme Variantu 1 Q = Qrec − F1 = Qrec − (Q fic ) 0 ,
(4.68)
a pokud tuto závislost zohledníme, tak dostaneme Variantu 2
F F Q = Qrec 1 − 2 − 1 . 1 + F2 1 + F2
(4.69)
Varianta 2 je ze všech způsobů měření transpiračního proudu nejlepší, protože fiktivní tok stanovuje nejpřesněji.
39
5 Závěr V práci bylo postupováno podle metodických bodů a stanovené dílčí cíle byly splněny. Těchto bodů bylo celkem 7 a jsou níže podrobně shrnuty. První tři se zabývaly deterministickým modelováním teplotního, vlhkostního a deformačního pole ve dřevě. Další dva body se zabývaly stochastickým modelováním teplotního a vlhkostního pole. Těchto pěti bodů lze využít při sušení řeziva nebo ve stavební fyzice u dřevostaveb. Poslední dva body se zabývaly problémem měření transpiračního proudu a tepelné vodivosti v běli kmene stromu. Oba body se skládaly z experimentální a teoretické části. V obou šlo o to, jak co nejpřesněji změřit transpirační proud a tepelnou vodivost běle pomocí zákonů termodynamiky. V práci byl sestaven matematický model popisující vázané šíření vlhkosti a teploty ve dřevě v mezích vody vázané s obecně-ortotropním charakterem materiálu. Byl zohledněn Soretův, Duforův efekt a závislost materiálových vlastností na teplotě a vlhkosti. Pro popis pohybu vody bylo využito poznatku, že se voda ve dřevě s vlhkostí pod mezí hygroskopicity šíří tzv. difuzí (molekulárním tokem). Aplikovány byly dva fyzikální zákony. První se nazývá Fickův (udává lineární vztah mezi gradientem koncentrace vody a hustotou vlhkostního toku ve dřevě prostřednictvím koeficientu difuze) a druhý se nazývá Fourierův (udává lineární vztah mezi gradientem teploty a hustotou tepelného toku ve dřevě prostřednictvím koeficientu tepelné vodivosti). V dalším výzkumu by bylo dobré model rozšířit pro dřeva s vlhkostí nad mezí hygroskopicity. Při šíření vlhkosti ve dřevě nad mezí hygroskopicity dochází k objemovému toku vody v makrokapilárách dřeva a tento tok není vhodné kvantifikovat na základě gradientu vlhkosti nebo koncentrace vody ve dřevě, ale je lepší využít gradientu tlaku, tedy Darcyho zákona, který udává lineární vztah mezi gradientem tlaku a hustotou vlhkostního toku prostřednictvím koeficientu propustnosti. Byly a stále ještě jsou vyvíjeny (programovány) dvě aplikace vycházející z výše odvozeného multifyzikálního modelu popisující vlhkostní a teplotní pole. První aplikace je ve fázi, kdy ji může uživatel použít pro simulaci vlhkostního pole v řezivu během jeho teplovzdušného sušení pro různé sušící režimy a výsledky těchto simulací porovnat. Tato aplikace je tedy zatím funkční jen pro testování různých sušících režimů a mohla by být obohacena o další prvky. Autor v současné době vyvíjí druhou aplikaci, 40
která bude podobného charakteru jako ta první, ale bude zaměřena na stavební fyziku, tedy na simulaci teplotního a vlhkostního pole ve vybraných konstrukcích dřevostavby za různých exteriérových a interiérových podmínek. Model vázaného teplotního a vlhkostního pole ve dřevě byl rozšířen o mechanickou analýzu, která je důležitá při optimalizaci parametrů sušícího prostředí. Tato mechanická analýza spočívá v popisu deformací a s tím souvisejících napětí. Deformace vyvolané během sušení jsou spojené se změnou vlhkosti, deformace vyvolané změnou teploty jsou zanedbatelné. Pro popis mechanického chování dřeva během sušení byly využity základní rovnice rovnováhy v nichž se vyskytují derivace napětí podle prostorových proměnných a časové derivace posunutí. Pro vyjádření složek napětí v závislosti na deformaci tělesa bylo využito Hookova zákona, který platí pouze v lineární oblasti závislosti napětí na deformacích, tedy jen do určitých mezních hodnot napětí, při jejichž překročení je třeba počítat s nelineárním vztahem mezi těmito dvěma veličinami. Jde o tzv. materiálové nelinearity (plasticita, viskoplasticita, apod.). Při velkých deformacích v matematickém modelu vznikají další nelinearity, které nazýváme geometrické nelinearity. Oba typy nelinearit (materiálové, geometrické) si zaslouží samostatnou studii a jsou pro autora této práce velmi zajímavým tématem, kterému se chce v dalších letech věnovat. Byl sestaven spektrální stochastický model jako účinný nástroj k analýze nejistot v nestacionárních difuzních úlohách sušení dřeva. Stejná úloha byla řešena metodou Monte Carlo a výsledky byly v dobré shodě s výsledky spektrální metody. Pro další analýzu nejistot v difuzních problémech sušení dřeva je tedy vhodné použít testovanou spektrální metodu, protože je efektivnější (ušetří spoustu výpočetního času a také paměti počítače). Podobně byl sestaven spektrální stochastický model pro řešení teplotního pole v konstrukci obvodového pláště dřevostavby. Opět byly výsledky srovnány s výsledky metody Monte Carlo a opět bylo dosaženo shody. Tedy i v tomto případě je výhodné pro kvantifikaci nejistot na výstupu použít spektrálního přístupu. Byla provedena teoretická analýza, numerické simulace a experimentální měření teplotního pole v okolí lineárního ohřevu (elektrickým odporem vyhřívaná jehla) pro účely měření transpiračního proudu a tepelné vodivosti v bělové části kmene stromu. Podařilo se vyvinout a aplikovat nové zapojení včetně nového vzorce pro 41
vyhodnocování transpiračního proudu v závislosti na tepelné vodivosti. Dále byl teoreticky analyzován a numerickými simulacemi verifikován nový přístup pro současné stanovení hustoty transpiračního proudu a tepelné vodivosti. Podobná práce byla provedena pro měření transpiračního proudu a tepelné vodivosti objemovým ohřevem (elektrodami). Zde byla zlepšena stávající THB metoda. Zlepšení bylo založeno na zpřesnění fiktivního toku, konkrétně šlo o zohlednění závislosti tohoto toku na okamžité hodnotě transpiračního proudu. Tato závislost byla odhalena pomocí numerických simulací a byla následně aplikována na naměřená reálná data. Měření bylo provedeno novým zapojením, které dokáže měřit i nulový a záporný tok.
42
Seznam použité literatury Adesanya BA, Nanda AK, Beard JN. 1988. Drying rates during high temperature drying of yellow poplar. Drying Technol. 6(1): 95-112 Avramidis S, Englezos R, Papathanasiou T. 1992. Dynamic nonisothermal transport in hygroscopic porous media: moisture diffusion in wood. AIChE Journal 38(8): 12791287 Avramidis S, Hatzikiriakos SG, Siau JF. 1994. An irreversible thermodynamics model for unsteadystate nonisothermal moisture diffusion in wood. Wood Science and Technology 28: 349-358 Avramidis S, Siau JF. 1987. An investigation of the external and internal resistance to moisture diffusion in wood. Wood Sci. Technol. 21: 249-256 Babiak M. 1990. Wood water system. TU Zvolen. 63 pp. Babiak, M. 1995. Is Fick’s law valid for the adsorption of water by wood? Wood Science and Technology 29: 227-229 Barbour JR, Johnson JA. 1989. A microstructurally based model for moisture movement in wood below the fiber saturation point. In: Kayihan E, Johnson JA, Smith RW. (ed.) Wood Drying Symposium. IUFRO. 247-254, Seattle, WA Beard JN, Rosen HN, Adesanya BA. 1983. Temperature distributions and heat transfer during the drying of lumber. Drying Technology 1(1): 117-140 Bramhall G. 1979 a. Mathematical model for lumber drying. I. Principles involved. Wood Sci. 12(1): 14-21 Bramhall G. 1979 b. Mathematical model for lumber drying. II. The model. Wood Sci. 12(1): 22-31 Bramhall G. 1995. Diffusion and the drying of wood. Wood Science and Technology 29: 209-215 Brooks JR, Meinzer FC, Coulombe R, Gregg J. 2002. Hydraulic redistribution of soil water during summer drought in two contrasting Pacific Northwest coniferous forests. Tree Physiol 22:1107–1117 Burgess SSO, Adams MA, Turner NC, Ong CK. 1998. The redistribution of soil water by tree root systems. Oecologia 115:306–311 Burgess SSO, Adams MA, Bleby TM. 2000a. Measurement of sap flow in roots of woody plants: a commentary. Tree Physiol 20:909–913 Chen P, Pei DCT. 1989. A mathematical model of drying processes. Int. J. Heat Mass Transfer 32(2): 297-310 Cloutier A, Fortin Y, Dhatt G. 1992. A wood drying finite element model based on the water potential concept. Drying Technology 10(5): 1151-1181 Constantine PG. 2009.Spectral methods for parameterized matrix equations. A dissertation Čermák J, Deml M, Penka M. 1973. A new method of sap flow rate determination in trees. Biol. Plant. (Praha) 15(3): 171-178. Čermák J, Deml M. 1974. Method of water transport measurements in woody species, especially in adult trees (in Czech). Patent (Certification of authorship) CSFR, No.155622 (P.V.5997-1972) Čermák J, Kučera J, Penka M. 1976. Improvement of the method of sap flow rate determination in adult trees based on heat balance with direct electric heating of xylem. Biol Plant (Praha) 18: 105–110. 43
Čermák J, Úlehla J, Kučera J, Penka M. 1982. Sap flow rate and transpiration dynamics in the full-grown oak (Quercus robur L.) in floodplain forest exposed to seasonal floods as related to potential evapotranspiration and tree dimensions. Biol. Plant. (Praha) 24(6): 446-460. Čermák J, Nadezhdina N. 1998. Sapwood as the scaling parameter - defining according to xylem water content or radial pattern of sap flow? Ann.Sci.For.55: 509521. Čermák J, Kučera J. 1991. Extremely fast changes of xylem water flow rate in mature trees, caused by atmospheric, soil and mechanical factors. 181-190pp. In: Proc.CEC International Workshop "Methodologies to assess the impacts of climatic changes on vegetation: Analysis of water transport in plants and cavitation of xylem transport in plants and cavitation of xylem conduits". Raschi A. Borghetti M. (eds.), May 2931.1991. Firenze, Italy. Čermák J, Kučera J, Nadezhdina N. 2004. Sap flow measurements with two thermodynamic methods, flow integration within trees and scaling up from sample trees to entire forest stands. Trees, Structure and Function 18: 529-546. Debusschere BJ, Najm HN, Pebay PP, Knio OM, Ghanem RG, Maitre OPL. 2004. Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes. SIAM Journal on Scientific Computing 26(2):698–719. Daum CR. 1967. A method for determining water transport in trees. Ecology 48: 425– 431. Droin A, Taverdet IL, Vergnaud JM. 1988. Modeling the kinetics of moisture adsorption by wood. Wood Sci. Technol. 22: 11-20 Droin-losserand A, Taverdet JL, Vergnaud JM. 1988. Modelling the absorption and desorption of moisture of wood in an atmosphere of constant and programmed relative humidity. Wood Sci. Technol. 22: 299-310 Droin-losserand A, Taverdet JL, Vergnaud JM. 1989 a. Modelling the process of moisture absorption in three dimensions by wood samples of various shapes: cubic, parallelepipedic. Wood Sci. Technol. 23: 259-271 Droin-losserand A, Taverdet JL, Vergnaud JM. 1989 b. Modelling of moisture absorption within a section of parallelepipedic sample of wood by considering longitudinal and transversal diffusion. Holzforschung. 43(5): 297-302 Ghanem RG, Spanos PD. 1990. Polynomial Chaos in Stochastic Finite Elements. Journal of Applied Mechanics 57: 197–202. Ghanem RG. 1999. Stochastic Finite Elements with Multiple Random Non-Gaussian Properties. Journal of Engineering Mechanics 26–40. Ghanem RG. 1999. Ingredients for a General Purpose Stochastic Finite Element Formulation. Computer Methods in Applied Mechanics and Engineering 168: 19–34 Granier A. 1985. Une nouvelle methode pour la mesure dy flux de seve brute dans le trons des arbres. Annales Scientifique Forestiere 22: 193-200. Green SR, Clothier BE. 1988. Water use of kiwifruit vines and apple trees be the heatpulse technique. Journal of Experimental Botany 39: 115-123. Helton JC, Davis FJ. 2003. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliability Engineering & System Safety 81(1):23–69. Horáček, P. 2004. Model vázaného šíření vlhkostního a teplotního pole při sušení dřeva, MZLU Brno, 126 s, ISBN 80-86386-59-7 44
Hosder S, Walters RW, Perez R. 2006. A Non-Intrusive Polynomial Chaos Method For Uncertainty Propagation in CFD Simulations. 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada Huber B. 1932. Beobachtung und Messung pflanzlicher Saftstrome. Berliner Deutsche Botanisches Geselschaft 50: 89-109. Huber B, Schmidt E. 1936. Weitere thermoelektrische Untersuchungen uber den Transpirationsstrom der Baume. Tharandter Forstliche Jahrsblad 87: 369-412. Irudayaraj J, Haghighi K, Stroshine RL. 1990. Nonlinera finite element analysis of coupled heat and mass transfer problems with application to timber drying. Drying technology 8 (4): 731-749. Isukapalli SS. 1999. Uncertainty Analysis of Transport-Transformation Models. PhD Dissertation Kayihan F. 1986. Adaptive optimal scheduling and control of batch lumber kilns using a population balance model. In: Mujnmdar, A. S. (ed.): Drying '86: 391-400 Kamke PA, Vanek M. 1994. Comparison of wood drying models. Proc. 4th IUFRO Wood Drying Conference, Rotorua, New Zeland, 1-21 Kollmann F, 1951. Technologie des Holzes und der Holzwerkstoffe. 1. Band, Springer Vrelag Berlin, 1050 pp. Koňas P. 2008. General model of wood in typical coupled tasks Part I. Phenomenological approach. Acta Universitatis agriculturae et silviculturae Mendelianae Brunensis, vol. LVI: 95-102 Kučera J. 1977. A system for water flux measurements in plants (in Czech). Patent (Certificate of authorship) CSFR No. 185039 (P.V. 2651-1976). Bureau for Inventions and Discoveries, Prague. Kučera J, Čermák J, Penka M. 1977. Improved thermal method of continual recording the transpiration flow rate dynamics. Biol. Plant. (Praha) 19(6): 413-420. Kunia Y. 1955. Studies on the sap streaming in plants by the thermoelectric method. Sci Rep Tokyo Univ. Biol 21:153–178 Liu IY. 1990. Lumber drying in a medium with variable potentials. In: General papers: Phase change and convective heat transfer: Proceedings of AIAA/ASME thermophysics and heat transfer conference, Vafai, K., et. al. (eds.): 149-156, Seattle, WA. Luikov AV. 1966. Heat and mass transfer in capillary-porous bodies. Pergamon Press, New York, 523 pp. MacLean JD. 1941. Thermal conductivity of wood. Heat. Pipe Air Cond. 13: 380-391. Marshall DC. 1958. Measurements of sap flow in conifers by heat transport. Plant Physiology 33: 385-396. Mathelin L, Hussaini MY, Zang TA, Bataille,F. 2004. Uncertainty Propagation for Turbulent. Compressible Nozzle Flow Using Stochastic Methods. AIAA Journal 42(8):1669–1676. Mathelin L, Hussaini MY, Zang TA. 2005. Stochastic Approaches to Uncertainty Quantification in CFD Simulations. Numerical Algorithms 38(1): 209–236. Nadezhdina N, Čermák J, Ceulemans R. 2002a. Radial patterns of sap flow in woody stems of dominant and understory species: scaling errors associated with positioning of sensors. Tree Physiology 22: 907-918. Nadezhdina N, Čermák J, Nadezhdin V. 1998a. Heat field deformation method for sap flow measurements. Proc. 4th. International Workshop on Measuring Sap Flow in Intact Plants. Židlochovice, Czech Republic, Oct.3-5, 1998. 72-92. IUFRO Publications. Publishing house of Mendel Univ.Brno. 45
Nadezhdina N, Čermák J. 1998b. "The technique and instrumentation for estimation the sap flow rate in plants". Patent No.286438 (PV-1587-98). Nadezhdina N, Čermák J, Neruda J, Prax A, Ulrich R, Nadezhdin V, Gašpárek J, Pokorný E. 2006. Roots under the load of heavy machinery in spruce trees. European J.For.Res. 125: 111-128. Nadezhdina N, Steppe K, De Pauw DJW, Bequet R, Čermák J, Ceulemans R. 2009. Stem-mediated hydraulic redistribution in large roots on opposing sides of a Douglas-fir tree following localized irrigation. New Phytologist 184:932–943 Nadezhdina N, David TS, David JS, Ferreira MI, Dohnal M, Tesar M, Gartner K, Leitgeb E, Nadezhdin V, Čermák J, Jimenez MS, Morales M. 2010. Trees never rest: the multiple faces of hydraulic redistribution. Ecohydrology 3: 431-444. Nelson RM. Jr. 1986 a. Diffusion of bound water in wood. Part 1: The driving force. Wood Sci. Technol. 20: 125-135 Nelson RM. Jr. 1986 b. Diffusion of bound water in wood. Part 2: A Model for isothermal diffusion. Wood Sci. Technol. 20: 235-251 Nelson RM. Jr. 1986 c. Diffusion of bound water in wood. Part 3: A model for nonisothermal diffusion. Wood Sci. Technol. 20: 309-328 Nouy A. 2009. Recent Developments in Spectral Stochastic Methods for the Numerical Solution of Stochastic Partial Differential Equations. Archives of Computational Methods in Engineering 16: 251-285. Offenthaler I, Hietz P. 1998. A comparison of different methods to measure sap flow in spruce. In: Čermák J, Nadezhdina N (eds) Measuring sap flow in intact plants. Proceedings of 4th International Workshop, Židlochovice, Czech Republic, IUFRO Publ. Mendel University, Brno, Czech Republic, 55–64. Pang, S. 1996. Moisture content gradient in softwood board during drying: simulation from a 2-D model and measurement. Wood Science and Technology 30: 165-178. Perez R, Walters R. 2005. An Implicit Compact Polynomial Chaos Formulation for the Euler Equations AIAA-Paper 2005-1406, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, CD-ROM. Plumb OA, Brown CA, Olmstead BA. 1984. Experimental measurements of heat and mass transfer during convective drying of southern pine. Wood Science Technology 18: 187-204 Plumb OA, Spolek GA, Olmstead BA. 1985. Heat and mass transfer in wood during drying. Intern. J. Heat Mass Transfer 28(9): 1669-1678 Reagan M, Najm HN, Ghanem RG, Knio OM. 2003. Uncertainty Quantification in Reacting Flow Simulations through Non-Intrusive Spectral Projection. Combustion and Flame 132: 545–555. Sakuratani T. 1981. A heat balance method for measuring water flux in the stem of intact plants. J Agric Meteorol (Jpn) 37: 9–17. Salin JG. 1991. Modeling of wood drying: a bibliography. Drying Technology 9(3): 775-793 Sherwood TK. 1929. The drying of solids, Ind. Eng. Chem., 21: 12-16. Schubert H. 1999. Qualitative and quantitative Untersuchung verschiedener Methoden der Xylemflussmessung an Baumen. Diploma Thesis, Forstwissentschaftlichen Fakultat der Ludwig- Maxmilians-Universitats Munchen, Lehrstuhl fur Forstbotanik, Freising Siau JF. 1983 a. A proposed theory for nonisothermal unsteady-state transport of moisture in wood. Wood Science and Technology 17: 75-77. 46
Siau JF. 1983 b. Chemical potential as a driving force for nonisothermal moisture movement in wood. Wood Sci. Technol. 17: l01-105 Siau, JF. 1984 a. Chemical potential and nonisothermal diffusion. Letter to the editor. Wood Fiber Sci. 16(4): 628-629 Siau JF. 1984 b. Transport processes in wood. Springer - Verlag, Berlin, Heidelberg, New York. 245 pp. Siau JF, Babiak M. 1983. Experiments on nonisothermal moisture movement in wood. Wood Fiber Science 15(1): 40-46. Siau JF, Bao F, Avramidis S. 1986. Experiments in nonisothermal diffusion of moisture in wood. Wood Fiber Science 18(1): 84-89 Siau, JF. 1992. Nonisothermal diffusion model based on irreversible thermodynamics. Wood Science and Technology 26(5): 325-328. Siau JF. 1995. Wood. Influence of moisture on physical properties. Virginia Polytechnic Institute and State University, NY, 227 pp. Simpson W. 1993. Specific gravity, moisture content, and density relationships for wood. U. S. Department of Agriclture Gen Skaar Ch. 1988. Wood-Water relations. Springer- Verlag Berlin, New York, Tokio. 263 pp. Söderström O, Salin JG. 1993. On determination of surface emission factors in wood drying. Holzforschung 47: 391-397 Stanish MA, Schajer GS, Kayihan E. 1985. Mathematical modeling of wood drying from heat and mass transfer fundamentals. Drying 1985: 360-367 Swanson RH. 1971. Velocity distribution patterns in ascending xylem sap during transpiration. In: Symposium on flow - its measurement and control in science and industry. Canadian Forestry Service Paper No.4/2/171: 11p. Tatarinov FA, Kučera J, Cienciala E. 2005. The analysis of physical background of tree sap flow measurements based on thermal methods. Measurements Science and Technology 16: 1157-1169. Thomas, HR, Lewis RW, Morgan K. 1980. An application of the finite element method to the drying of timber. Wood and Fiber. 11(4): 237-243 Vergnaud IM. 1991. Liquid transport processes in polymeric materials: modelling and industrial applications. Prentice Hall, Englewood Cliffs, N.J., 362 pp. Voight H, Krisher O, Schauss H. 1940. Die Feuchtigkeitsbewegung bei der Verdunstungstrocknung von Holz. Holz Roh- Werkst. 3: 305-321. Walters R. 2003. Towards stochastic fluid mechanics via Polynomial Choas-invited AIAA-Paper 2003-0413, 41st AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, CD-ROM. Whitaker S. 1977. Simultaneous heat, mass and momentum transfer in porous media: A theory of drying. In: Hartnett J.P., Irvine T.F. Jr. (eds.): Advances in heat transfer, Academic Press, Vol. 13: 119-203 Wiener N. 1938. The Homogeneous Chaos. American Journal of Mathematics 60(4): 897–936. Xiu D, Karniadakis G E. 2003. Modeling Uncertainty in Flow Simulations via Generalized Polynomial Chaos. Journal of Computational Physics 187(1): 137–167. Xiu D, Karniadakis G. 2002. The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal of Scientific Computing 24: 619–644.
47
Publikace autora disertační práce Trcala M, Koňas P. 2012. Transformation relations and matrix implementation of multiphysics model for temperature and moisture fields in wood. Wood Research 57(1). Trcala M, Koňas P. (submitted 2012) Numerical solution of coupled temperature, moisture and strain fields in wood. Wood Research. Trcala M. (submitted 2012a). Stochastic modelling of uncertainties in diffusion problems of wood drying. European Journal of Wood and Wood Products. Trcala M. (submitted 2012b) Spectral stochastic finite element approach to effective uncertainty quantification in heat conduction problems. Applied Numerical Mathematics. Trcala M, Čermák J. (submitted 2011) Sap flow and heat conductivity measurement and their interactions. Acta Horticulturae. Trcala M, Čermák J. (submitted 2012). Improvement of the trunk heat balance method including measurement of zero and reverse sap flows. Agricultural and Forest Meteorology.
48
TRANSFORMATION RELATIONS AND MATRIX IMPLEMENTATION OF MULTIPHYSICS MODEL FOR TEMPERATURE AND MOISTURE FIELDS IN WOOD Miroslav Trcala, Petr Koňas Mendel University in Brno, Brno, Czech Republic
ABSTRACT This paper deals with the numerical simulation of 3D heat and moisture transfer in wood. The appropriate mathematical model is based on a set of coupled temperature and mass transfer equations. The main contributions of this paper are: first, the transformation relations between matrixes of material coefficients (diffusion and thermal conductivity coefficients) in anatomical and geometrical directions in three-dimensional case; second, the matrix implementation of a numerical solution for the three-dimensional problem of nonisothermal moisture transfer in the anisotropic structure of wood in convective drying with transient phenomena (Soret effect, Duffour effect). The simulation is based on the unsteady-state nonisothermal nonlinear diffusion of moisture and heat with respect to the orthotropic nature of wood. It deals with a transport model which makes use of the gradient of moisture and temperature for motive power and respects the dependence of material coefficients on temperature and moisture. KEY WORDS: temperature and moisture, multiphysics modelling, matrix transformations
INTRODUCTION Wood drying is still controlled by experience to a great degree. For this reason, many scholars have been trying to develop a drying model to describe wood drying processes completely, both in the theory and practice for years. There is an attempt to use the model to guide, design and optimize the wood drying process and to find out the way to increase the drying quality, decrease the energy consumption and decrease the drying time. The development of the wood drying model is based on the study of three fundamental phenomena (heat transfer, moisture transfer and strain-stress analysis). The theory of transport phenomena in porous materials was summarised by Luikov (1975, 1980) and Whitaker (1977). Luikov (1975) developed a set of coupled partial differential equations to describe the heat and mass transport in capillary porous media by assuming that the transfer of moisture is analogical to heat transfer and that capillary transport is proportional to moisture and temperature gradient. This approach was used by Thomas et al. (1980) and Irudayaraj (1990). But the solution, even numerical, is complicated, involving complex eigenvalues (Liu, Cheng 1991). Malan and Lewis 49
(2003) and Lewis et al. (1996) demonstrated the efficacy of the finite element method in solving highly non-linear drying systems. A number of theoretical models have been suggested in wood drying (e.g. Chen and Pei 1989) but they are not completely applicative to simulate real drying processes as seen in industry. It is caused by a low level of generality of the proposed models that is brought about by a wide range of coupled physical fields included in the model, the specification of boundary conditions, the evaluation of experimental characteristics, etc. Koňas (2008, part I) assembled a coefficient form of a coupled stress-strain task with moisture/temperature dependency of wood with orthotropic, viscoelastic material properties suitable for current numerical solvers. The weak solution of this task and the subgrid upscaling homogenization method for the large scale hierarchical structure which is typical for wood was obtained in Koňas and Přemyslovská (2008, part II). The modified Ritz-Galerkin method was derived for a simple solution of this problem. The difficulty with the mathematical formulation of the wood drying process lies in the number of combined transfer mechanisms, the interdependencies among these mechanisms, and the different variables controlling them. Many external factors can also affect the wood drying speed, of which temperature, humidity (equilibrium moisture content) and airflow speed are the primary factors, and these three main factors determine the change of the wood moisture content. Therefore, the establishment of a phenomenological model based on the relation between the macroscopical characteristics of wood and drying conditions is another important focus of research in modelling. The use of Fick’s law for the description of the moisture content has been questioned by (Babiak 1995). The process of wood drying can be interpreted as a simultaneous heat and moisture transfer with local thermodynamic equilibrium at each point within the timber (Horáček 2004). Drying of wood is in its nature an unsteadystate non-isothermal diffusion of heat and moisture, where temperature gradients may counteract with the moisture gradient. For a theoretical description of these phenomena, thermodynamic models seem to be the most suitable. Although unsteady-state nonisothermal experiments have been conducted (e.g. Avramidis et al. 1994), we are convinced that practical solutions of mathematical models that couple moisture and heat transfer with inhomogeneous material, combined with an experimental verification, are still lacking (Irudayaraj et al. 1990). A considerable volume of research has been carried out regarding the modelling of the moisture and heat transfer in materials like polymers, wood, or agricultural products (Salin 1991, Kamke and Vanek 1994). The modelling of moisture fluxes under unsteady-state non-isothermal conditions has been noticeably absent in literature (Avramidis et al. 1994), although the theory behind was proposed twenty years ago (Siau 1983). The heat and moisture transfer should be considered as coupled processes: a thermally induced mass transfer, the Soret effect (Siau 1984, Avramidis et al. 1994, Horáček 2004) and the heat flux resulting from the moisture diffusion, the Duffour effect (Siau 1992), should be all taken into account. Luikov (1966), and in much details Whitaker (1977) developed a unique approach that describes the simultaneous heat and moisture transfer in drying processes, based on irreversible thermodynamic processes. The determination of the moisture diffusion coefficient is based on models derived by Siau (1995) and Skaar (1988), the thermal conductivity by MacLean (1941), the specific heat of wood by Skaar (1988), and the density of wood by Kollmann (1951). In contrast to the usually applied assumption by Crank (1975), the initial values of 50
temperature and moisture do not have to be uniformly distributed within the specimen. The numerical value used for the heat and mass transfer coefficient in the simulation of wood drying has been discussed in Soderstrom and Salin (1996a, 1996b), Avramidis et al. (1994), Siau (1995) and Pang (1996) together with the correction coefficient introduced by Plumb et al. (1985) and Siau (1995).
Theory The coupled and highly non-linear nature of the transport equations that govern the drying process highlights the applicability of numerical simulation in this field (Perré, Turner 1999). The mathematical model used in the present work is formed by a system of two partial differential equations and corresponding third-order boundary conditions. It is a purely mathematic treatment of wood as a continuous and homogenous medium, without explaining those physical mechanisms associated with the moisture transport that takes place simultaneously inside wood. The moisture content and the temperature gradient are set as driving forces, since all the other possible factors related to the moisture content are applicable only in the hygroscopic region. The proposed equations are inspired by Siau (1992) and Avramidis (1994), who with these equations found the best agreement between the experimental fluxes and thermodynamic model. The threedimensional transfer of heat and moisture is generally described as follows:
∂w − ∇(D∇w + sD∇T ) = 0 ∂t E ∂w ∂T Cρ + b − ∇λ∇T = 0 ∂t 1.8C ∂t
(1) (2)
where D is the matrix of moisture diffusion coefficients of wood [m2s-1], λ the matrix of thermal conductivity coefficients of wood [Wm-1K-1], s is the Soret effect [-/K], w is the moisture content [-], T is the temperature [K], C is the specific heat of heat [Jkg-1K-1], ρ is the wood density [kgm-3]. The relationship s=
ϕ ∂EMC Eb RT ∂ϕ T
(3)
refers to the Soret effect (thermodiffusion) based on the slope of the sorption isotherms ∂EMC and the activation energy for water diffusion Eb, where φ is the relative air ∂ϕ humidity and R is the universal gas constant. In general, these isotherms are specific to wood specimens, since differential heat of sorption for different wood differs noticeably (Babiak 1990). Boundary conditions are obtained from the heat and moisture transfer between wood surfaces and external air. The following boundary conditions were applied for the numerical solution: −n(D∇w + sD∇T ) = α w ( w∂Ω − EMC )
(4) 51
where α w is the mass transfer coefficient [ms-1], w∂Ω is the surface moisture [-], EMC is the equilibrium moisture content of wood [-] −nλ∇T = α T (T∂Ω − Tair )
(5)
where α T is the heat transfer coefficient [Wm-2K-1], T∂Ω is the surface temperature [K], Tair is the air temperature [K]. The initial conditions can be defined as functions of the position in the body of wood. In this work the initial conditions are taken as constants: w=fw(x,y,z)=w0 (6) T=fT(x,y,z)=T0
MATERIAL AND METHODS During physical analysis and mathematical modelling it is necessary to consider the anisotropic wood structure. Often the values of the diffusion and heat conductivity coefficients in geometrical directions are not known. These coefficients are established (measured, analytically determined) only in anatomical directions and it is necessary to transform them into geometrical directions before solving the system of partial differential equations (PDEs) describing the simultaneous, nonisothermal, nonlinear, nonstationary moisture and heat transfer in a generally orthotropic wood body. The transformation relations and the matrix form of the mentioned system of PDEs are described below. The numerical simulations are performed for pinewood 0.045 m thick, 0.2 m wide and 5 m long with an initial moisture of 0.25 and an initial temperature of 20 °C. The air temperature is kept at a constant value of 80 °C and the relative humidity (φ) and the corresponding equilibrium moisture content (EMC) depend on time as featured in Tab.1. The air flow is constant - 2 ms-1. Tab. 1: Time drying schedule
t (hours) Tair (°C)
φ (-)
EMC (-)
4.5 9 17 29 45.5 54.5
0.77 0.67 0.59 0.46 0.34 0.3
0.108 0.087 0.073 0.056 0.042 0.039
80 80 80 80 80 80
52
Transformation relations and implementation into COMSOL Multiphysics Let L: V→V be a linear operator on a vector space V. Let Lα be the matrix of L relative to a basis α for V. Let Lβ be the matrix of L relative to another basis β for V. Let Pβα be the transition matrix from the basis α to β. Then
Lβ = Pβα L α Pβα-1 .
(7)
Thus, if the matrix of diffusion coefficients is known relative to a basis α, it can be expressed relative to any other basis β. Then the matrix Lα has a diagonal form:
DR Lα = 0 0
0 DT 0
0 0 DL
(8)
For generally orthotropic material as wood is this matrix Lα must be transformed to a matrix of diffusion and heat-conductive coefficients Lβ in general axes of the model at the deviation of anatomic directions from these axes
Dxx Lβ = Dyx D zx
Dxy Dyy Dzy
Dxz DR Dyz = Pβα 0 0 Dzz
0 DT 0
0 0 Pαβ . DL
(9)
The transition matrixes are assembled from the following relations: Pαβ = Pβα -1 ,
Pβα = β -1α , where α, β are matrixes with columns formed by vectors of bases α, β. In this work β is assumed to be the standard basis (thus it is formed by unit vectors and then β is the identity matrix) and consequently Pβα = α . The previous information is important, because it enables us to gain transition matrixes Pβα and Pαβ easily. Consequently, by using relation (9) the matrix of diffusion (or thermal conductivity coefficients) can be obtained and be used at Fick´s (or Fourier´s) laws for the mathematical description of moisture (or heat) transfer inside generally orthotropic material. This theory is valid for an arbitrary basis α and therefore the above mentioned relations are applicable for an arbitrary fiber lean. The simulations of wood drying were provided for wood with a fiber lean as illustrated in Fig.1. The angle u denotes the transverse fiber lean and v denotes the longitudinal fiber lean.
53
Fig. 1: Deviation of anatomic directions (R, T, L) from geometric axes (x, y, z) of the model (u=35°, v=10°)
For the fiber lean as illustrated in Fig.1 the transition matrix Pβα is following:
0 cos(u ) -sin(u ) Pβα = sin(u ) cos(u ) sin(v) 0 0 cos(v)
(10)
and the searched matrix Lβ : Dxx Lβ = Dyx D zx
Dxy Dyy Dzy
Dxz cos(u ) -sin(u ) 0 DR Dyz = sin(u ) cos(u ) sin(v) 0 0 cos(v) 0 Dzz 0
0 DT 0
0 cos(u ) -sin(u ) 0 0 sin(u ) cos(u ) sin(v) 0 cos(v) DL 0
−1
DR cos 2 (u ) + DT sin 2 (u ) sin(u ) cos(u )( DR − DT ) sin(u ) cos(u )( DR − DT ) tan(v) 2 2 2 2 Lβ = sin(u ) cos(u )( DR − DT ) DT cos (u ) + DR sin (u ) −( DT cos (u ) + DR sin (u ) − DL ) tan(v) 0 0 DL
(11) Regarding the notations in the previous theory section, the matrix Lβ is written as D ( Lβ = D) so that it is obvious that this is the matrix of diffusion coefficients. The system of PDEs (1, 2) can be written into the FEM environment of COMSOL Multiphysics as follows: ea
∂ 2u ∂u + da + ∇(−c∇u − αu + γ ) + au + β∇u = f 2 ∂t ∂t
where:
54
(12)
1 w D sD u = , c = , da = Eb T 0 λ 1.8C Dxx D = Dyx D zx
Dxy Dyy Dzy
sDxx Dxz Dyz , sD = sDyx sD Dzz zx
0 C ρ sDxy sDyy sDzy
(13) λxx sDxz sDyz , λ = λ yx λ sDzz zx
λxy λ yy λzy
λxz λ yz λzz
and others input parameters are zero matrixes 0. Thus: D sD ∇ w c∇u = 0 λ ∇T
(14),
1 ∂u da = E ∂t b 1.8C
∂w 0 ∂t C ρ ∂T ∂t
D sD ∇w α w ( w∂Ω − EMC ) −n . = 0 λ ∇T αT (T∂Ω − Tair )
(15)
(16)
RESULTS AND DISCUSSION Mathematical relations for the transformation of the diffusion and thermal conductivity matrixes for generally orthotropic material were developed. The above mentioned diffusion model of the coupled moisture and heat transfer was solved numerically in program COMSOL Multiphysics (the solver is based on the finite element method) with the moisture initial condition below the fiber saturation point (w0=0.25) and temperature initial condition (T0=293.15 K). During numerical calculations the dependencies of material coefficients on moisture and temperature were taken into account. The simulations were performed for wood species with an absolutely dry wood density ρ0=0.4 gcm-3 (corresponding to pinewood).
a)
b)
Fig. 2: Influence of temperature T and moisture w on a) diffusion coefficient Dxy and b) thermal conductivity coefficient lxy (=λxy) with fiber lean as in figure 1
55
Fig. 2 a) depicts the dependence of the diffusion coefficient Dxy on moisture and temperature. The coefficient together with the other eight diffusion coefficients are the results of the above mentioned transformation of diffusion coefficients from anatomical to geometric directions. The pictured dependence is conditioned by the dependence of diffusion coefficients in anatomical directions on moisture and temperature according to Siau (1995). Fig. 2 b) shows the dependence of the thermal conductivity coefficient lxy on moisture and temperature. This coefficient is, as the diffusion coefficient, a result of the above mentioned transformation. The pictured dependence is conditioned by the dependence of thermal conductivity coefficients in anatomical directions on moisture and temperature according to MacLean (1941) and Siau (1995). As for wood density, its dependence on moisture is considered according to Kollmann (1951) and the analytical formula for calculating of the specific heat was taken from Skaar (1972, 1988). The numerical values used for the heat and mass transfer coefficient in the simulation of wood drying were discussed in Soderstrom and Salin (1996), Avramidis et al. (1994), Siau (1995) and Pang (1996) together with the correction coefficient introduced by Plumb et al. (1985) and Siau (1995). This work assumes the dependence of the mass transfer coefficient on moisture and temperature according to Avramidis et al. (1994).
Fig. 3: Influence of wood moisture w and relative air humidity phi (=φ) on the Soret effect coefficient ST(=s)
The heat and moisture transfer should be considered as coupled processes: the thermally induced mass transfer, the Soret effect (Siau 1984, Avramidis et al. 1994) and the heat flux resulting from the moisture diffusion, the Duffour effect (Siau 1992), should all be taken into account. The Soret effect is displayed in Fig. 3 and we can see that both the temperature and moisture influence it. The Soret effect decreases with the increase in moisture and increases with the increase in temperature. There is a noticeable high dependence of the effect on the relative air humidity within the range of the relative air humidity from 0.6 to 1. The dependence of the Soret effect on the moisture of wood is smaller and linear. Luikov (1966), and in much details Whitaker 56
(1977), developed a unique approach that describes the simultaneous heat and moisture transfer in drying processes, based on irreversible thermodynamic processes. Numerical solution of the model
a)
b)
c)
Fig. 4: a) The affect of fiber deviation on the progress of average moisture content; b) the Soret effect on the progress of moisture content in the medium part of timber; c) the Soret effect on the progress of moisture content in the surface of timber
Figure a) shows that influence of the fiber lean (Fig. 1) on the drying speed is almost insignificant. When comparing the simulated progess of moisture content at the beginning of drying, Fig. b) shows that the Soret effect causes an increase in moisture content in the middle of the timber; Fig. c) shows that it causes a faster drying of timber surface, which corresponds to the theory according to Siau (1984), Avramidis et al. (1994)
a)
b)
Fig. 5: Simulated moisture content profiles in the middle of the thickness (y=0.0225 m) and along the timber (along the x axis) in a) the surface plane (z=0) and b) the medium plane (z=2.5)
Fig. 5 displays the influence of the Soret effect on the moisture distribution within the body at the beginning of drying. We can see that from the very beginning of the drying process up to the time between 20000s and 50000s (about 10 hours) the effect of the temperature gradient on the moisture diffusion decreases and even vanishes.
57
a)
b)
Fig. 6: Moisture distribution of cross section z=2.5 (in the middle of the length) a) for selected times in the range 0 to 100000 seconds ( 0 to 27.8 hours); b) in time 100000 seconds (27.8 hours)
a)
b)
Fig. 7: Differentiations of moisture with respect to spatial variables in the middle of the wood body (0.1,0.0225,2.5) a) and at surface point (0.1,0.0225,0); b) with the Soret effect and without the Soret effect
The graphs (Fig. 6) of moisture distribution during the drying process show that the wooden board is actually subjected to relatively high differences between the surface moisture and the inner moisture after vanishes the Soret effect. The numerical results for the moisture gradient curves are displayed in Fig.7. The largest gradients are at surface points at the beginning of the drying process. Thus the largest stresses (depending on the moisture gradients) are assumed at these points and times. The predicted values by the assembled theoretical model are compared with norm data taken under actual drying conditions to show the difference between the predictive model and time drying schedules values.
58
0,3
average MC (-)
0,25 0,2 model
0,15
norm
0,1 0,05 0 4,5
9
17
29
45,5
54,5
time (hours)
Fig. 8: The differences between simulated values of average moisture content (w) and the values stipulated by the norm (wn)
Fig. 8 shows the prediction of the average moisture content simulated by the model together with norm values according to the drying schedule depicted in Tab. 1. The maximum difference between the values does not exceed 2.5% of moisture content. At the beginning of the drying process the model predicts a higher drying speed than the norm features. Next, the trend reverts and the model drying speed is smaller and after about 17 hour of drying process the curves intersect. At the end of drying the model predicts average MC of 10.5 % and the norm features desired 8 % of average MC.
CONCLUSION A mathematical model describing the coupled heat and moisture transfer in a wood body during the drying process was assembled and solved. The transformation relations between matrixes of material coefficients (diffusion and thermal conductivity coefficients) in anatomical and geometrical directions in three-dimensional case were developed. There were respected the dependencies of material coefficients on temperature and moisture. The numerical solution of the model was performed because an analytical solution to the mentioned coupled unsteady-state nonisothermal diffusion phenomena is difficult to find. The matrix form of the mentioned mathematical model was developed and implemented into FEM solver.
ACKNOWLEDGEMENT The Institutional research plan MSM6215648902 – Forest and Wood: the support of functionally integrated forest management and use of wood as a renewable raw material (2005 – 2010, Ministry of Education, Youth and Sport, Czech Republic) supported this project. And the project was also supported by the Internal Grant Agency (IGA) of the Faculty of Forestry and Wood Technology Mendel University in Brno project no. 59/2010.
59
REFERENCES 1. Avramidis, S., Hatzikiriakos, S.G., Siau, J.F. 1994: An irreversible thermodynamics model for unsteadystate nonisothermal moisture diffusion in wood. Wood Science and Technology 28, Pp. 349-358 2. Babiak, M., 1990: Wood water system. TU Zvolen. 63 pp. 3. Babiak, M., 1995: Is Fick’s law valid for the adsorption of water by wood? Wood Science and Technology 29: 227-229 4. Chen, P., Pei, D.C.T., 1989: A mathematical model of drying processes. Int. J. Heat Mass Transfer 32(2), Pp. 297-310 5. Crank, J., 1975: Mathematics of diffusion. Cranendon Press, New York - Oxford. 6. Horáček, P., 2004: Modelling of coupled moisture and heat transfer during wood drying, MZLU Brno, 126 s, ISBN 80-86386-59-7 7. Irudayaraj, J., Haghighi, K., Stroshine, R.L. 1990: Nonlinear finite element analysis of coupled heat and mass transfer problems with application to timber drying. Drying technology 8 (4), Pp. 731-749 8. Kamke, P.A., Vanek, M., 1994: Comparison of wood drying models. Proc. 4th IUFRO Wood Drying Conference, Rotorua, New Zeland, Pp. 1-21 9. Kollmann, F., 1951: Technologie des Holzes und der Holzwerkstoffe. 1. Band, Springer Vrelag Berlin, 1050 pp. 10. Koňas, P., 2008: General model of wood in typical coupled tasks Part I. Phenomenological approach, Acta Universitatis agriculturae et silviculturae Mendelianae Brunensis, vol. LVI, Pp. 95-102 11. Koňas, P., Přemyslovská, E., 2008: General model of wood in typical coupled tasks Part II. Weak solution, Acta Universitatis agriculturae et silviculturae Mendelianae Brunensis, vol. LVI, Pp. 103-108 12. Lewis, R.W., Morgan, K., Thomas H.R., 1996: The Finite Element Method in Heat Transfer Analysis, Wiley, 1996 13. Liu, J.Y., Cheng S., 1991: Solution of Luikov equations of heat and mass transfer in capillary porous bodies. International Journal of Heat and Mass Transfer 34 (7), 1747–1754. 14. Luikov, A.V., 1966: Heat and mass transfer in capillary-porous bodies. Pergamon Press, New York, 523 p. 15. Luikov, A.V., 1975: Systems of differential equations of heat and mass transfer in capillary-porous bodies. International Journal of heat and Mass Transfer 8, 1-14 16. Luikov, A.V., 1980: Heat and Mass Transfer, Mir, Moscow, 623 p. 17. MacLean J.D., 1941: Thermal conductivity of wood. Heat. Pipe Air Cond. 13, Pp. 380-391 18. Malan, A.G., Lewis, R.W., 2003: Modelling coupled heat and mass transfer in drying non-hygroscopic capillary particulate materials, Communications in Numerical Methods in Engineering 19, 669–677 19. Pang, S., 1996: Moisture content gradient in softwood board during drying: simulation from a 2-D model and measurement. Wood Science and Technology 30, Pp. 165-178 20. Perré, P., Turner I.W., 1999: A 3-D version of TransPore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media, International Journal of Heat and Mass Transfer 42, 4501±4521 60
21. Plumb, O.A., Spolek, G.A., Olmstead, B.A., 1985: Heat and mass transfer in wood during drying. Intern. J. Heat Mass Transfer 28(9), Pp. 1669-1678 22. Söderström, O., Salin, J-G., 1993: On determination of surface emission factors in wood drying. Holzforschung 47, Pp. 391-397 23. Salin, J-G., 1991: Modeling of wood drying: a bibliography. Drying Technology 9(3), Pp. 775-793 24. Salin, J-G., 1996a: Mass Transfer from Wooden Surfaces and International Moisture Non-Equilibrium. Drying technology 14 (10), Pp. 2213-2224 25. Salin, J-G., 1996b: Prediction of heat and mass transfer coefficients for individual boards and board surfaces. A review. Proceedings 5th IUFRO International Wood Drying Conference, Quebec City, Canada, Pp. 49-58 26. Siau, J.F., 1983: A proposed theory for nonisothermal unsteady-state transport of moisture in wood. Wood Science and Technology 17, Pp. 75-77 27. Siau, J.F., 1984: Transport processes in wood. Springer- Verlag, Berlin, Heidelberg, New York. 245 pp. 28. Siau, J.F., 1992: Nonisothermal diffusion model based on irreversible thermodynamics. Wood Science and Technology 26(5), Pp. 325-328 29. Siau, J.F., 1995: Wood. Influence of moisture on physical properties. Virginia Polytechnic Institute and State University, NY, 227 pp. 30. Skaar, Ch., 1988: Wood-Water relations. Springer- Verlag Berlin, New York, Tokio. 263 pp. 31. Thomas, H.R., Morgan K., Lewis, R.W., 1980: A fully nonlinear analysis of heat and mass transfer problems in porous media. International Journal of Numerical Methods in Engineering 15, 1381-1393 32. Whitaker, S., 1977: Simultaneous heat, mass and momentum transfer in porous media: A theory of drying. In: Hartnett J.P., Irvine T.F. Jr. (eds.): Advances in heat transfer, Academic Press, Vol. 13, Pp. 119-203
61
Numerical solution of coupled temperature, moisture and strain fields in wood Miroslav Trcala, Petr Koňas Department of Wood Science, Mendel University in Brno, Czech Republic
Abstract This paper deals with the multiphysics modelling and numerical simulations of unsteady-state nonlinear diffusion of heat and moisture and corresponding moisture deformations in convective drying of wood. The model respects the dependence of material coefficients on temperature and moisture, the Soret effect, the Duffour effect and the anisotropic nature of wood. The model was implemented into FEM solver and numerical simulations were performed. Keywords: multiphysics modelling; temperature; moisture; moisture deformations; wood drying 1. Introduction Thermo-hygro-mechanical problems exist in many engineering areas including drying applications, soil and building applications etc. In this paper, the problem of wood drying is solved. The drying of wood is one of the most important production processes in the technology of wood processing. The complexity of the drying of lumber is attributable to the variety of physico-mechanical properties of different species, anisotropy of wood structure, interrelation between the phenomena of heat and moisture transfer, strains and stresses occurring in drying, and large differences in the dimensions of assortments. There is an attempt to use the model to guide, design and optimize the wood drying process and to find out the way to increase the drying quality, decrease the energy consumption and decrease the drying time. The development of the wood drying model is based on the study of three fundamental phenomena (heat transfer, moisture transfer and stress-strain analysis). The theory of transport phenomena in porous materials was summarised by Luikov (1975, 1980) and Whitaker (1977). Luikov (1975) developed a set of coupled partial differential equations to describe the heat and mass transport in capillary porous media by assuming that the transfer of moisture is analogical to heat transfer and that capillary transport is proportional to moisture and temperature gradient. The use of Fick’s law for the description of the moisture content has been questioned by (Babiak 1995). The process of wood drying can be interpreted as a simultaneous heat and moisture transfer with local thermodynamic equilibrium at each point within the timber. Process of drying of wood can be simplified into an unsteady-state non-isothermal diffusion of heat and moisture, where temperature gradients may counteract with the moisture gradient. For a theoretical description of these phenomena, thermodynamic models seem to be the most suitable. Although unsteady-state non-isothermal experiments have been conducted (e.g. Avramidis et al. 1994), we are convinced that practical solutions of 62
mathematical models that couple moisture and heat transfer with inhomogeneous material, combined with an experimental verification, are still lacking (Irudayaraj et al. 1990). A considerable volume of research has been carried out regarding the modelling of the moisture and heat transfer in materials like polymers, wood, or agricultural products (Salin 1991, Kamke and Vanek 1994). The modelling of moisture fluxes under unsteady-state non-isothermal conditions has been noticeably absent in literature (Avramidis et al. 1994), although the theory behind was proposed twenty years ago (Siau 1983). The heat and moisture transfer should be considered as coupled processes: a thermally induced mass transfer, the Soret effect (Siau 1984, Avramidis et al. 1994) and the heat flux resulting from the moisture diffusion, the Duffour effect (Siau 1992), should be all taken into account. Luikov (1966), and in much details Whitaker (1977) developed a unique approach that describes the simultaneous heat and moisture transfer in drying processes, based on irreversible thermodynamic processes. The determination of the moisture diffusion coefficient is based on models derived by Siau (1995) and Skaar (1988), the thermal conductivity by MacLean (1941), the specific heat of wood by Skaar (1988), and the density of wood by Kollmann (1951). In contrast to the usually applied assumption by Crank (1975), the initial values of temperature and moisture do not have to be uniformly distributed within the specimen. The numerical value used for the heat and mass transfer coefficient in the simulation of wood drying has been discussed in Soderstrom and Salin (1996), Avramidis et al. (1994), Siau (1995) and Pang (1996) together with the correction coefficient introduced by Plumb et al. (1985) and Siau (1995). Because of wood shringing during drying process, strains and stresses develop and it can lead to unusuable products. Realistic modelling of thermo-hygro-mechanical behaviour of wood undergoing moisture changes is essential for optimization of drying processes. Accurate prediction of such behaviour requires not only mathematical models capable of describing multi-physics phenomena of heat and moisture, hygroexpansion, creep as well as thermo-mechanical and mechano-sorptive effects as they develop in time under changing external conditions, but also, clearly defined, meaningful, and reliable material data. The mechanical behaviour of wood is complex and combines different behaviours that include elasticity, plasticity, viscoelasticity and mechano-sorption. Understanding the drying behaviour using experimental approach alone, poses practical difficulties due to the extremely large number of variables that must be considered in addition to the large time involved and the cost. Drying is fundamentally a problem of simultaneous heat and mass transfer under transient conditions resulting in a system of coupled non-linear partial differential equations. The simulation of wood drying has been the subject of numerous works, including physical and mechanical formulations and numerical solutions (Thibeault et al. 2010). A number of theoretical models have been suggested in wood drying (e.g. Chen and Pei 1989) but they are not completely applicative to simulate real drying processes as seen in industry. It is caused by a low level of generality of the proposed models that is brought about by a wide range of coupled physical fields included in the model, the specification of boundary conditions, the evaluation of experimental characteristics, etc. (Koňas 2008).
63
2. Theory The coupled and highly nonlinear nature of transport equations that govern the drying process highlights the applicability of numerical simulation in this field (Perré, Turner 1999). The mathematical model used in the present work is formed by a system of five partial differential equations for temperature, moisture and displacements. It is a purely mathematic treatment of wood as a continuous and homogenous medium, without explaining those physical mechanisms associated with the moisture transport that takes place simultaneously inside wood. The proposed equations are inspired by Siau (1992) and Avramidis (1994), who with these equations found the best agreement between the experimental fluxes and thermodynamic model. The three-dimensional transfer of heat and moisture is generally described as follows: ∂M − ∇ ⋅ (D∇M + sD∇T ) = 0 ∂t E ∂M ∂T ρC − b − ∇ ⋅ λ∇T = 0 ∂t 1.8C ∂t
(1) (2)
where D is the matrix of moisture diffusion coefficients of wood [m2s-1], λ the matrix of thermal conductivity coefficients of wood [Wm-1K-1], s is the Soret effect [-/K], M is the moisture content [-], T is the temperature [K], C is the specific heat of heat [Jkg-1K-1], ρ is the wood density [kgm-3]. The relationship s=
ϕ ∂EMC Eb RT ∂ϕ T
(3)
refers to the Soret effect (thermodiffusion) (Siau 1984, Avramidis et al. 1994) based on the slope of the sorption isotherms ∂EMC / ∂ϕ and the activation energy for water diffusion Eb = 38500 − 29000 M , where EMC is the equilibrium moisture content and it is a function of both relative humidity and temperature of surrounding air EMC (ϕ , T ) = 1 / B ln( A / ln(1 / ϕ )) , where φ is the relative air humidity, A = 7.731706 − 0.014348T , B = 0.008746 + 0.000567T and R is the universal gas constant. In general, these isotherms are specific to wood specimens, since differential heat of sorption for different wood differs noticeably (Babiak 1990). Boundary conditions are obtained from the heat and moisture transfer between wood surfaces and external air. The following boundary conditions were applied for the numerical solution: −n ⋅ (D∇M + sD∇T ) = α M ( M ∂Ω − EMC )
(4)
where α M is the mass transfer coefficient [ms-1], M ∂Ω is the surface moisture [-], EMC is the equilibrium moisture content of wood [-] 64
−n ⋅ λ∇T = αT (T∂Ω − Tair )
(5)
where αT is the heat transfer coefficient [Wm-2K-1], T∂Ω is the surface temperature [K], Tair is the air temperature [K]. The initial conditions can be defined as functions of the position in the body of wood. In this work the initial conditions are taken as constants: M=fM(x,y,z)=M0 (6) T=fT(x,y,z)=T0 It is possible to completely describe the strain conditions at a point with the displacement components ( u, v, w ) in 3D and their derivatives. The strain in a material
is described by the symmetric tensor consisting of three normal strains ( ε x , ε y , ε z ) and
six, or if symmetry is used, three shear strains ( ε xy , ε yz , ε xz ) . The stress in a material is
described by the symmetric tensor consisting of three normal stresses (σ x , σ y , σ z ) and three shear stresses (τ xy ,τ yz ,τ xz ) .
The stress-strain relationship for the small-displacement assumption and for linear conditions are given as follows:
σ = C(ε - ε M ) ,
(7)
where
σ x σ y σ σ = z , τ xy τ yz τ xz
∂u ∂x ∂v ∂y ε Mx εx CMET ∂w ε My εy CMER ∂z ε ε CMEL ε = z = 1 ∂u ∂v , ε M = Mz = ( M − M 0 ) ε xy + ε Mxy 0 ε yz 2 ∂y ∂x ε Myz 0 1 ∂ v ∂ w ε ε 0 xz Mxz + 2 ∂ z ∂ y 1 ∂w ∂u + 2 ∂x ∂z
65
(8)
where ε M is hygroscopic strain and CMET, CMER, CMEL are coefficients of moisture expansion. Ex (1 − µ yz µ zy ) ∆ E y ( µ xy + µ xz µ zy ) ∆ C = ( Cij ) = Ez ( µ xz + µ xy µ yz ) i , j =1,..,6 ∆ 0 0 0
Ex ( µ yx + µ zx µ yz )
Ex ( µ zx + µ yx µ zy )
∆ E y (1 − µ zx µ xz )
∆ E y ( µ zy + µ zx µ xy )
∆ Ez ( µ yz + µ xz µ yx )
Ez (1 − µ xy µ yx )
0
0
∆
0
0
∆ 0
∆ 0
0
0
2Gxy
0
0
0
0
2Gyz
0
0
0
0
0 0 0 0 0 2Gxz (9)
where ∆ = 1 − µ xy µ yx − µ yz µ zy − µ zx µ xz − 2 µ xy µ yz µ zx
The equilibrium equations expressed in the stresses for 3D are ∂ 2u ∂σ x ∂τ xy ∂τ xz − − − = Fx ∂t 2 ∂x ∂y ∂z ∂σ ∂τ ∂ 2 v ∂τ ρ 2 − xy − y − yz = Fy ∂t ∂x ∂y ∂z 2 ∂τ ∂ w ∂τ ∂σ ρ 2 − xz − yz − z = Fz ∂t ∂x ∂y ∂z
ρ
(10)
where F denotes the volume forces (body forces). 3. Methods The above mentioned model of the coupled temperature, moisture and strain fields was solved numerically in program COMSOL Multiphysics (the solver is based on the finite element method) with the moisture initial condition below the fiber saturation point (M0=0.25) and temperature initial condition (T0=293.15 K). During numerical calculations the dependencies of material coefficients on moisture and temperature were taken into account as well as the Soret and the Duffour effects. The simulations were performed for wood species with an absolutely dry wood density ρ0=0.4 gcm-3 (corresponding to pinewood). The material is 0.045 m thick, 0.2 m wide and 2 m long. The air temperature is kept at a constant value of 80 °C and the relative humidity (φ) and the corresponding equilibrium moisture content (EMC) depend on time as featured in Tab.1. The air flow is constant - 2 ms-1. 66
Tab. 1: Time drying schedule
t (hours) Tair (°C)
φ (-)
EMC (-)
4.5 9 17 29 45.5 54.5
0.77 0.67 0.59 0.46 0.34 0.3
0.108 0.087 0.073 0.056 0.042 0.039
80 80 80 80 80 80
4. Results and discussion 4.1.
Material properties
a)
b)
Fig.1: Influence of temperature T and moisture M on diffusion coefficient for a) transverse direction (tangential) DT and b) longitudinal direction DL
a)
b)
Fig.2: Influence of temperature T and moisture M on thermal conductivity coefficient for a) transverse direction (tangential) λT and b) longitudinal direction λL
Fig. 1 depicts the dependence of the diffusion coefficients DT and DL on moisture and temperature according to Siau (1995). Fig. 2 shows the dependence of the thermal conductivity coefficients λT and λL on moisture and temperature according to MacLean (1941) and Siau (1995). 67
a)
b)
Fig.3: a) Influence of moisture M on density ρ; b) Influence of temperature T and moisture M on specific heat C
As for wood density, its dependence on moisture is considered according to Kollmann (1951) and the analytical formula for calculating of the specific heat was taken from Skaar (1988). The numerical values used for the heat and mass transfer coefficient in the simulation of wood drying were discussed in Soderstrom and Salin (1996), Avramidis et al. (1994), Siau (1995) and Pang (1996) together with the correction coefficient introduced by Plumb et al. (1985) and Siau (1995). This work assumes the dependence of the mass transfer coefficient on moisture and temperature according to Avramidis et al. (1994). The heat and moisture transfer should be considered as coupled processes: the thermally induced mass transfer, the Soret effect (Siau 1984, Avramidis et al. 1994) and the heat flux resulting from the moisture diffusion, the Duffour effect (Siau 1992), should all be taken into account. The Soret effect decreases with the increase in moisture and increases with the increase in temperature. 4.2.
Numerical solutions
Fig. 4: Moisture field in the wood body at 54.5 hours of drying time
68
Fig. 5: Deformed shape and moisture field in the middle of the length for selected times
The numerical results for the moisture and displacements changes are displayed in Fig.5. The difference between the average deformations in the tangential and radial directions is significant.
a)
b)
Fig. 6: Simulated moisture content profiles in the middle of the thickness (y=0.0225 m) and along the timber (along the x axis) in a) the surface plane (z=0) and b) the medium plane (z=1)
Fig. 6 displays the influence of the Soret effect on the moisture distribution within the body at the beginning of drying. We can see that from the very beginning of the drying process up to the time about 10 hours the effect of the temperature gradient on the moisture diffusion decreases and even vanishes. At the beginning of drying the Soret effect causes an increase in moisture content in the middle of the timber and a faster drying of timber surface.
69
5. Conclusion A mathematical model of coupled temperature, moisture and strain fields was assembled and solved by FEM solver. A theory section of the paper introduced the basic equations describing drying process of wood. A methods section presented the time drying schedule and the dependencies of material coefficients on temperature and moisture. A results section presented the numerical solution of the coupled unsteadystate nonlinear model. This model can be used to predict the distribution of temperature, moisture and strain fields in the hygroscopic moisture content and below 100°C temperature ranges for single boards. Acknowledgements The project was supported by the Internal Grant Agency (IGA) of the Faculty of Forestry and Wood Technology Mendel University in Brno project no. 7/2011. References [1] Avramidis, S., Hatzikiriakos, S.G., Siau, J.F. 1994: An irreversible thermodynamics model for unsteadystate nonisothermal moisture diffusion in wood. Wood Science and Technology 28, 349-358 [2] Babiak, M., 1990: Wood water system. TU Zvolen, 63 pp. [3] Babiak, M., 1995: Is Fick’s law valid for the adsorption of water by wood? Wood Science and Technology 29, 227-229 [4] Chen, P., Pei, D.C.T., 1989: A mathematical model of drying processes. Int. J. Heat Mass Transfer 32(2), 297-310 [5] Crank, J., 1975: Mathematics of diffusion. Cranendon Press. New York – Oxford, 414 pp. [6] Irudayaraj, J., Haghighi, K., Stroshine, R.L. 1990: Nonlinear finite element analysis of coupled heat and mass transfer problems with application to timber drying. Drying technology 8 (4), 731-749 [7] Kamke, P.A., Vanek, M., 1994: Comparison of wood drying models. Proc. 4th IUFRO Wood Drying Conference, Rotorua, New Zeland, 1-21 [8] Kollmann, F., 1951: Technologie des Holzes und der Holzwerkstoffe. 1. Band, Springer Verlag Berlin, 2233 pp. [9] Koňas, P., 2008: General model of wood in typical coupled tasks Part I. Phenomenological approach, Acta Universitatis agriculturae et silviculturae Mendelianae Brunensis, vol. LVI, 95-102 [10] Lewis, R.W., Morgan, K., Thomas H.R., 1996: The Finite Element Method in Heat Transfer Analysis, Wiley [11] Liu, J.Y., Cheng S., 1991: Solution of Luikov equations of heat and mass transfer in capillary porous bodies. International Journal of Heat and Mass Transfer 34 (7), 1747–1754 [12] Luikov, A.V., 1966: Heat and mass transfer in capillary-porous bodies. Pergamon Press. New York, 523 pp. [13] Luikov, A.V., 1975: Systems of differential equations of heat and mass transfer in capillary-porous bodies. International Journal of heat and Mass Transfer 8, 1-14 [14] Luikov, A.V., 1980: Heat and Mass Transfer, Mir Publishers. Moscow, 623 pp. 70
[15] MacLean J.D., 1941: Thermal conductivity of wood. Heat. Pipe Air Cond. 13, 380391 [16] Malan, A.G., Lewis, R.W., 2003: Modelling coupled heat and mass transfer in drying non-hygroscopic capillary particulate materials, Communications in Numerical Methods in Engineering 19, 669–677 [17] Pang, S., 1996: Moisture content gradient in softwood board during drying: simulation from a 2-D model and measurement. Wood Science and Technology 30, 165-178 [18] Perré, P., Turner I.W., 1999: A 3-D version of TransPore: a comprehensive heat and mass transfer computational model for simulating the drying of porous media, International Journal of Heat and Mass Transfer 42, 4501-4521 [19] Plumb, O.A., Spolek, G.A., Olmstead, B.A., 1985: Heat and mass transfer in wood during drying. Intern. J. Heat Mass Transfer 28(9), 1669-1678 [20] Söderström, O., Salin, J-G., 1993: On determination of surface emission factors in wood drying. Holzforschung 47, 391-397 [21] Salin, J-G., 1991: Modeling of wood drying: a bibliography. Drying Technology 9(3), 775-793 [22] Salin, J-G., 1996a: Mass Transfer from Wooden Surfaces and International Moisture Non-Equilibrium. Drying technology 14 (10), 2213-2224 [23] Salin, J-G., 1996b: Prediction of heat and mass transfer coefficients for individual boards and board surfaces. A review. Proceedings 5th IUFRO International Wood Drying Conference, Quebec City, Canada, 49-58 [24] Siau, J.F., 1983: A proposed theory for nonisothermal unsteady-state transport of moisture in wood. Wood Science and Technology 17, 75-77 [25] Siau, J.F., 1984: Transport processes in wood. Springer- Verlag. Berlin, Heidelberg, New York, 245 pp. [26] Siau, J.F., 1992: Nonisothermal diffusion model based on irreversible thermodynamics. Wood Science and Technology 26(5), 325-328 [27] Siau, J.F., 1995: Wood: Influence of moisture on physical properties. Virginia Polytechnic Institute and State University, 227 pp. [28] Skaar, Ch., 1988: Wood-Water relations. Springer- Verlag. Berlin, New York, Tokio, 283 pp. [29] Thibeault F., Marceau, Younsi R., Kocaefe D., 2010: Numerical and experimental validation of thermo-hygro-mechanical behaviour of wood during drying proces. International Communications in Heat and Mass Transfer 37, 756–760 [30] Thomas, H.R., Morgan K., Lewis, R.W., 1980: A fully nonlinear analysis of heat and mass transfer problems in porous media. International Journal of Numerical Methods in Engineering 15, 1381-1393 [31] Whitaker, S., 1977: Simultaneous heat, mass and momentum transfer in porous media: A theory of drying. In: Hartnett J.P., Irvine T.F. Jr. (eds.): Advances in heat transfer, Academic Press, Vol. 13, 119-203
71
Stochastic modelling of uncertainties in diffusion problems of wood drying Miroslav Trcala Department of Wood Science, Mendel University in Brno, Czech Republic Abstract This paper deals with the stochastic numerical analysis of moisture transfer in wood with the random diffusion coefficient of wood. The spectral solution of this problem is based on discretization the resulting random field (moisture) in the stochastic dimension by the orthogonal polynomials (generalized polynomial chaos algorithm). A Galerkin projection is applied in the stochastic dimension to obtain the deterministic set of partial differential equations that is solved by finite element method. The main purpose of this paper is to demonstrate that the stochastic spectral method based on polynomial chaos expansion can be more efficient in modelling uncertainties associated with moisture transfer in wood than Monte Carlo method. This spectral approach has a big advantage over the Monte Carlo method (statistical approach) in terms of computer time. Numerical example of diffusion of moisture in convective drying of wood is given and there is shown that the results (mean and the standard deviation) obtained with the stochastic spectral method are in good agreement with the results of the Monte Carlo simulations. Keywords: spectral stochastic finite element method; polynomial chaos; random diffusion coefficient; Monte Carlo method 1. Introduction Diffusion problems exist in many engineering areas including drying applications, soil and building applications etc. In this paper, the stochastic diffusion problem of wood drying is solved. The drying of wood is one of the most important production processes in the technology of wood processing. The complexity of the drying of lumber is attributable to the variety of physico-mechanical properties of different species, anisotropy of wood structure, interrelation between the phenomena of heat and moisture transfer, strains and stresses occurring in drying, and large differences in the dimensions of assortments. There is an attempt to use the model to guide, design and optimize the wood drying process and to find out the way to increase the drying quality, decrease the energy consumption and decrease the drying time. The development of the wood drying model is based on the study of three fundamental phenomena (heat transfer, moisture transfer and stress-strain analysis). The theory of transport phenomena in porous materials was summarised by Luikov (1975, 1980) and Whitaker (1977). Luikov (1975) developed a set of coupled partial differential equations to describe the heat and mass transport in capillary porous media by assuming that the transfer of moisture is analogical to heat transfer and that capillary 72
transport is proportional to moisture and temperature gradient. The use of Fick’s law for the description of the moisture content has been questioned by (Babiak 1995). Numerical simulations in diffusion problems rely on obtaining accurate physical properties of materials. In many situations, it is difficult to obtain the accurate values of material properties. Most of building envelopes are porous media and their structure is complex and material properties are difficult to determine. Thus, uncertainty quantification and propagation in physical systems appear as a critical path for the improvement of the prediction of their response (Nouy 2009). In the last two decades, a growing interest has been devoted to a new family of methods, called spectral stochastic methods, for the propagation of uncertainties through physical models governed by stochastic partial differential equations. In this paper we focus on uncertainty propagation with the polynomial chaos (PC) method, which is based on the spectral representation of the uncertainty. The biggest advantage of the spectral PC method is given mainly by low computational costs. Several methods have been used to model and propagate uncertainty in stochastic computational simulations (Hosder et al. 2006). Several researchers have studied and implemented the PC approach for a wide range of problems. Ghanem and Spanos (1990) and Ghanem (1999) applied the PC method to several problems of interest to the structures community. Mathelin et al. (2004) studied uncertainty propagation for a turbulent, compressible nozzle flow by this technique. Xiu and Karniadakis (2003) analyzed the flow past a circular cylinder and incompressible channel flow by the PC method and extended the method beyond the original formulation of Wiener (1938) to include a variety of basis functions. In 2003, Walters applied the PC method to a twodimensional steady-state heat conduction problem for representing geometric uncertainty. Following this effort, an implicit compact PC formulation was implemented for the stochastic Euler equations (Perez and Walters 2005). The PC method for the propagation of uncertainty in computational simulations involves the substitution of uncertain variables and parameters in the governing equations with the polynomial expansions. In general, an intrusive approach will calculate the unknown polynomial coefficients by projecting the resulting equations onto basis functions (orthogonal polynomials) for different modes. As its name suggests, the intrusive approach requires the modification of the deterministic code and this may be inconvenient for many complex computational problems since the modification of the existing code can be difficult, expensive, and time consuming. To alleviate this inconvenience, non-intrusive methods have been investigated by many researchers. Most of the non-intrusive PC approaches in the literature are based on sampling (Debusschere et al. 2004; Reagan et al. 2003; Isukapalli 1999) or quadrature methods (Debusschere et al. 2004; Mathelin et al. 2005) to determine the projected polynomial coefficients. Before their introduction, the standard was to employ Monte Carlo (MC) methods to sample the inputs, compute the output for each sample, and aggregate statistics of the output (Helton and Davis 2003) In fact, this is still the most widely used method in practice due to its robustness and ease of implementation. However, the MC methods suffer from a dreadfully slow convergence rate proportional to the inverse of the square root of the number of samples. And if each sample evaluation is expensive, such as the solution of a partial differential equation, then obtaining hundreds of thousands of samples may be entirely infeasible. Thus, the initial applications of spectral methods showed orders-of-magnitude reduction in the work needed to estimate statistics with comparable accuracy (Xiu and Karniadakis 73
2002). Such results spurred interest in applying spectral methods to differential equations with stochastic (i.e. parameterized) inputs (Constantine 2009). This paper presents a spectral polynomial chaos method for the propagation of input uncertainty in diffusion coefficient of wood. 1.1 Generalized polynomial chaos Generalized polynomial chaos (gPC) is employed to represent stochastic processes. Stochastic process can be seen as a process involving some randomness. They can be represented by a stochastic mathematical model, often expressed in terms of stochastic differential equations. Stochastic mathematical models are based on a probability space (Ω, F , P) where Ω denotes the space of elementary events, F ⊂ 2Ω its σ-algebra of events, and P its probability measure. In addition considering some physical domain D ⊂ R d × T (d=1, 2, 3), which can be a combination of spatial and temporal dimensions, a stochastic process can be seen as a scalar or vector valued function u(x, t , ω ) : D × Ω → Rb where x is an element of the physical space, t denotes the time and ω is a point in the sample space. Furthermore, because of the infinite-dimensional nature of the probability space, we discretize this space by characterizing it by a finite number of random variables
{ξ (ω )} j
N j =1
. This can be seen as assigning a finite number of coordinates {ξ j }
N j =1
to the
probability space reducing it to a finite dimensional space Λ ⊂ R N . Consequently, the stochastic process u becomes a mapping u(x, t , ξ ) : D × Λ → Rb . It is important to note that in this work, we assume that the occurring stochastic processes are already characterized by a known set of random variables. Wiener was the first to represent stochastic processes by orthogonal polynomial expansions (Wiener 1938). To accomplish this, he used Hermite polynomials in terms of Gaussian random variables to represent Gaussian processes, which is referred to as homogeneous chaos. In this way, the stochastic process is represented in the form: ∞
u(x, t , ξ (ω )) = ∑ ui (x, t ) H i (ξ (ω ))
(1)
i=0
in which H i are Hermite polynomials and ξ is a vector of Gaussian random variables with zero mean and unit variance. It is a spectral expansion in the random dimensions employing deterministic coefficients. According to the Cameron-Martin theorem (Cameron and Martin 1947), for a fixed value of x and t, this expansion converges to any L2 (Ω) functional in the L2 (Ω) sense. This implies that the application of polynomial chaos is restricted to those stochastic processes yielding
∫ω
∈Ω
2
u(x, t , ω ) dP(ω ) < ∞.
(2)
As a result, polynomial chaos is restricted to second-order stochastic processes, i.e. processes with finite second-order moments. These are processes with finite variance, 74
and this applies to most physical processes. Although Wiener’s original polynomial chaos expansion converges to any second-order stochastic process, it is most suitable to represent Gaussian processes, due to the random variable’s Gaussian nature, yielding a fast convergence. In order to deal with a broader range of stochastic processes, the Wiener-Hermite chaos has been generalized to the generalized polynomial chaos (Xiu and Karniadakis 2002), also referred to as Wiener-Askey polynomial chaos. Analogously, gPC is a means of representing second-order stochastic processes through the expansion: ∞
u(x, t , ξ (ω )) = ∑ ui (x, t )Φ i (ξ (ω )).
(3)
i=0
Here the random trial base {Φ i (ξ(ω ))} exists out of orthogonal polynomials from the Askey-scheme, of which the Hermite polynomials are a subset, in terms of a random vector ξ = {ξi (ω )} j =1 . The combination of random vector and polynomials is carefully N
selected based on the distribution of the random input. It seems that for certain random variables, their probability density function (PDF) uniquely corresponds to one of the weighting functions ω (ξ ) in the orthogonality relation of the different orthogonal polynomials of the Askey-scheme. Choosing the corresponding combination leads to a proper gPC expansion. Since each of the polynomials of the Askey-scheme forms a complete basis in the Hilbert space determined by their corresponding support, it is expected, according to Xiu and Karniadakis (2002), that each type of Wiener Askey expansion converges to any L2 (Ω) functional in the L2 (Ω) sense in the corresponding Hilbert functional space as a generalized result of Cameron-Martin theorem. This generalized polynomial chaos has later been further generalized for arbitrary random inputs (Wan and Karniadakis 2006; Wittveen and Bijl 2006). In this way, gPC is not restricted anymore to random inputs of standard types. In order to expand the solution in a polynomial chaos expansion in terms of this arbitrary random variables, a proper random trial base {Φ j (ξ )} should be created, and this according to the gPC rules. This means that a set of orthogonal polynomials is created in such a way that they are orthogonal with respect to a weighting function which equals the probability density function of the random input. Different algorithms can be employed for this orthogonalization, among others the Stieltjes procedure, the Lanczos algorithm (Wan and Karniadakis 2006) or the Gram-Schmidt orthogonalization (Wittveen and Bijl 2006). Taking this further generalization into account, one can truly speak of generalized polynomial chaos. The polynomials of the gPC’s random trial base satisfy following orthogonality relation Φ i Φ j = Φ i2 δ ij where δ ij is the Kronecker delta and
⋅ denotes the ensemble average. This ensemble average is the inner product in the Hilbert space determined by the measure of the random variables f (ξ ) g (ξ ) = ∫ f (ξ ) g (ξ )dP(ω ) = ∫ f (ξ ) g (ξ ) w(ξ )dξ with w(ξ ) the corresponding ω∈Ω
weighting function and with integration in the last integral taken over the support of ξ . With the introduction of the gPC, optimal performance can now be achieved for a broader range of stochastic processes. In order to accomplish this, it is important to choose the appropriate type of chaos depending on the stochastic process’ nature. 75
Although every type of chaos converges to any second-order process, choosing the appropriate polynomial chaos expansion leads to an optimal solution. This has been presented in (Xiu and Karniadakis 2002), where it is computationally shown that exponential convergence is achieved when applying the optimal gPC through a Galerkin projection in random space on a stochastic ordinary differential equation (Vos 2006). 1.2 Galerkin projection The stochastic problem can be formulated as: find the stochastic process u ( x, t , ξ (ω )) , which solves the stochastic differential equation:
L ( x, t , ξ; u ) = f ( x, t , ξ ) ,
(4)
where f ( x, t , ξ ) is a source term and L is a differential operator involving differentiation in space and time. In the spectral expansion of the solution u ( x, t , ξ (ω )) Q
u(x, t , ξ (ω )) = ∑ u i (x, t )Φ i (ξ (ω ))
(5)
i=0
the unknown deterministic coefficients ui (x, t ) , which represent the different modes of the solution, need to be determined. To accomplish this, start with substituting the expansion for the solution (5) into the governing equation (4) Q L x, t , ξ ; ∑ u i Φ i = f ( x, t , ξ ) i=0
(6)
Next, successively multiplying this equation by the different orthogonal polynomials of the finite expansion, and taking the statistical average Q L x, t , ξ; u ∑ ui Φ i , Φ j = f ( x, t , ξ ) , Φ j i=0
j=0, 1, …, Q
(7)
results in a set of (Q + 1) coupled equations for the different random modes ui (x, t ) . The total number of expansion terms (Q+1) is determined by the dimension N of the random vector ξ and the highest degree P of the orthogonal polynomials {Φ i } Q +1 =
( N + P )!.
(8)
N !P!
As a result, the governing set of equations for the expansion coefficients ui (x, t ) is completely deterministic. This system should further be solved in space and time by finite element method or other numerical solvers.
76
2. Governing equations The steady-state of temperature field is generally described as follows:
∂M − ∇ ⋅ D∇M = 0 ∂t
(9)
where D is the matrix of moisture diffusion coefficients of wood [m2s-1] and M is the moisture content [-]. Boundary conditions are obtained from the heat transfer between envelope surfaces and external air. The following boundary conditions were applied for the numerical solution:
−n ⋅ D∇M = α M ( M ∂Ω − EMC )
(10)
where α M is the mass transfer coefficient [ms-1], M ∂Ω is the surface moisture [-], EMC is the equilibrium moisture content of wood [-]. By using the generalized polynomial chaos expansion, the temperature can be expand as Q
M (x, t , ξ (ω )) = ∑ M i (x, t )Φ i (ξ ) .
(11)
i=0
The random heat conductivity has the form D → D + σξ ,
(12) where ξ is an uniform random variable ξ
U ( −1,1).
After substituting the terms (11, 12) we can rewrite the above equations (9, 10) as Q
∂ ∑ M i (x, t )Φ i (ξ ) i =0
∂t
Q
− ∇ ⋅ ( D + σξ ) ∑ ∇M i (x, t )Φ i (ξ ) = 0
(13)
i=0
Q
Q
i =0
i =0
−n ⋅ ( D + σξ ) ∑ ∇M i (x, t )Φ i (ξ ) = α M (∑ M i (x, t )Φ i (ξ ) − EMC ).
77
(14)
By Galerkin projection of the above equation onto polynomial basis {Φ j (ξ )} for each
j = {0,1,..., Q} we obtain this set of (Q + 1) coupled partial differential equations Q
∑ i =0
∂M i (x, t ) Φ i (ξ ), Φ j (ξ ) − ∇ ⋅ D Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇M i (x, t ) = 0 ∂t (15)
(
)
Q Q −∑ n ⋅ D Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇M i (x, t ) = α M ∑ M i (x, t ) Φ i (ξ ), Φ j (ξ ) − EMC 1, Φ j (ξ ) i =0 i =0
(16) Once we have the coefficients of the expansion, we can compute approximate statistics of the solution with the formulas for mean value Q
Q
i=0
i =0
E[ M (x, t , ξ (ω ))] = E[∑ M i (x, t )Φ i (ξ )] = ∑ M i (x, t )E[Φ i (ξ )] = M 0 (x, t ) ,
(17)
for variance Q
Q
Var[ M (x, t , ξ (ω ))] = E[(∑ M i (x, t )Φ i (ξ ) − M 0 (x, t )) 2 ] = E[(∑ M i (x, t )Φ i (ξ )) 2 ] = i =0
Q
i =1
Q
(18)
= ∑ M i (x, t )E[Φ i (ξ ) ] = ∑ M i (x, t ) Φ i (ξ ), Φ i (ξ ) 2
2
2
i =1
i =1
and finally for standard deviation
Std [ M (x, t , ξ (ω ))] =
Q
∑M i =1
2 i
(x, t ) Φ i (ξ ), Φ i (ξ ) .
(19)
3. Example Numerical simulations are performed for example of moisture diffusion problem in convective drying of wood. The material is 0.05 m thick, 0.2 m wide. Geometry and mesh of the material are shown in figure 1. Boundary conditions are obtained from the moisture transfer between wood surfaces and the ambient air corresponding to the equilibrium moisture content of wood EMC=0.05 [-]. The moisture transfer coefficient is α M = 10 −7 [ms-1].
78
.
Figure 1: Geometry and triangular mesh (Lagrange-Quadratic elements, 976 elements, Number of degrees of freedom 10145)
Now replace some deterministic model inputs with stochastic quantities that better represent my understanding of the problem. One example of this may be an equation parameter (diffusion coefficient) that admits some measurement error, which we can model with a stochastic representation. Consider the uniform random diffusion coefficients, tangential DT defined on the interval DT ∈ 7 ⋅10−10 ,13 ⋅10−10 and radial
DR defined on the interval DR ∈ 12 ⋅10−10 ,18 ⋅10 −10 . For the matrix of moisture diffusion coefficients of wood
D D= T 0
0 is the following relationship D = E [ D] + σξ DR
where
10 ⋅10−10 0 σ 0 3 ⋅10−10 E [ D] = , σ = = 0 15 ⋅10−10 0 σ 0
0 3 ⋅10−10
and ξ is an uniform random variable ξ ∈ U (−1,1) with the corresponding weighting 1 function w(ξ ) = . Thus the Legendre polynomials must be used in this work. 2 Convergent results of spectral solution were obtained by fourth-order polynomial chaos expansion with one-dimensional random vector ξ = ξ . The total number of (1 + 4 )! = 5 . expansion terms is Q + 1 = 1!4! For orthogonal Legendre polynomials
79
Φ 0 (ξ ) = 1 Φ1 (ξ ) = ξ 3 1 Φ 2 (ξ ) = ξ 2 − 2 2 5 3 3 Φ 3 (ξ ) = ξ − ξ 2 2 35 4 30 2 3 Φ 4 (ξ ) = ξ − ξ + 8 8 8 the inner products in the Hilbert space are following: 1 1 Φ i (ξ ), Φ j (ξ ) = ∫ Φ i (ξ )Φ j (ξ ) w(ξ )d ξ = ∫ Φ i (ξ )Φ j (ξ )d ξ 2 −1
Φ 0 (ξ ), Φ 0 (ξ ) = 1
Φ1 (ξ ), Φ1 (ξ ) = 0.3333 Φ 2 (ξ ), Φ 2 (ξ ) = 0.2 Φ 3 (ξ ), Φ 3 (ξ ) = 0.1429 Φ 4 (ξ ), Φ 4 (ξ ) = 0.1111
ξΦ 0 (ξ ), Φ1 (ξ ) = 0.6667 ξΦ1 (ξ ), Φ 2 (ξ ) = 0.2667 ξΦ 2 (ξ ), Φ 3 (ξ ) = 0.1714 ξΦ 3 (ξ ), Φ 4 (ξ ) = 0.127
Many programming languages have the ability to generate pseudo-random numbers which are effectively distributed according to the standard uniform distribution U(0,1). If u is a value sampled from the standard uniform distribution, then the value a+(b−a)u follows the uniform distribution parameterized by a and b, denoted as U(a,b). This is used for Monte Carlo solution of our problem. The main purpose of this example is to perform spectral solution of the mentioned diffusion problem and to illustrate that it is more efficient approach than the Monte Carlo method in terms of computational time. 4. Results and discussions
The above derived system of PDEs was solved numerically by FEM solver.
(a)
(b)
Figure 2: Spectral solution of stochastic problem for (a) mean and (b) standard deviation in the middle of the wood body
80
The graph in Figure 2 (a) shows time curve of mean and the graph in Figure 2 (b) shows time curve of standard deviation in the middle of the wood body [0.1,0.025]. The influence of input uncertainty (random diffusion coefficient) on the progress of standard deviation is interesting. We can see that between 300000s and 400000s the influence is the most significant.
81
(a)
(b)
Figure 3: Spectral solution of stochastic problem for (a) mean and (b) standard deviation at times 2, 4, 6,8 ⋅105 seconds
1
2
3
4 Figure 4: Monte Carlo solution of stochastic problem (200 samples) for standard deviation
at times 2, 4, 6,8 ⋅105 seconds
The graphs show the mean and the standard deviation distribution during the drying process. These graphs were obtained by spectral method (polynomial chaos expansion of moisture, Galerkin projection in the stochastic dimension and solving the system of coupled PDEs by FEM solver). Figure 3 (a) shows standard moisture distribution (mean field) in wood for selected times. The standard deviation field is shown in Figure 3 (b) and it is more interesting than Figure 3 (a) because it informs us about the level of the output uncertainties (uncertainty in moisture field). It is seen that the values of the standard deviation is small and thus the output uncertainties in this case (for the input uncertainties-random diffusion coefficients) is also small. The variation of the standard 82
deviation across the domain changes during drying time. At the beginning of the drying process the standard deviation has the two local extremes that are symmetric about the vertical center line and then this variation decreases and even vanishes. Convergent results of spectral solution were obtained by fourth-order polynomial chaos expansion with one-dimensional random vector ξ = ξ . Monte Carlo simulations were conducted to validate the results by spectral method. For case of the unsteady diffusion problem a 200 realizations were needed to obtain converged Monte Carlo results that agree well with those of chaos expansion as is shown in Figure 4. Figure 4 shows Monte Carlo results only for standard deviation field at the selected times because for mean field are the results the same as is shown in Figure 3 (a).
5. Conclusion
A stochastic spectral approach to model uncertainty in wood during drying has been developed. The generalized polynomial chaos and Galerkin projection were applied to the solution of these problems. Monte Carlo simulations were also conducted to verify the spectral method and there was shown that the results are in good agreement with the results of the Monte Carlo simulations. This spectral approach has advantages over the Monte Carlo approach mainly in terms of computational costs (mainly time). This method can be used for a wider range of the problems and deserves further research. Acknowledgement
This work was supported by the Internal Grant Agency (IGA) of the Faculty of Forestry and Wood Technology, Mendel University in Brno, project no. 7/2011. References
Babiak M (1995) Is Fick’s law valid for the adsorption of water by wood? Wood Science and Technology 29, 227-229 Cameron R, Martin W (1947) The orthogonal development of nonlinear functionals in series of Fourier-Hermite functionals. Ann. Math.;48(2):385-392. Chen P, Pei DCT (1989) A mathematical model of drying processes. Int. J. Heat Mass Transfer;32(2):297-310. Constantine PG (2009) Spectral methods for parameterized matrix equations. A dissertation. Crank J (1975) Mathematics of diffusion. Cranendon Press. Debusschere BJ, Najm HN, Pebay PP, Knio OM, Ghanem RG, Maitre OPL (2004) Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes. SIAM Journal on Scientific Computing; 26(2):698–719. Ghanem R, Spanos PD (1990) Polynomial Chaos in Stochastic Finite Elements. Journal of Applied Mechanics;57:197–202. Ghanem R (1999) Stochastic Finite Elements with Multiple Random Non-Gaussian Properties. Journal of Engineering Mechanics:26–40.
83
Ghanem R (1999) Ingredients for a General Purpose Stochastic Finite Element Formulation. Computational Methods in Applied Mechanical Engineering;168:19–34 Helton JC, Davis FJ (2003) Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliability Engineering & System Safety;81(1):23–69. Hosder S, Walters RW, Perez R (2006) A Non-Intrusive Polynomial Chaos Method For Uncertainty Propagation in CFD Simulations. 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada 2006. Isukapalli SS (1999) Uncertainty Analysis of Transport-Transformation Models. PhD Dissertation. Lewis RW, Morgan K, Thomas HR (1996) The Finite Element Method in Heat Transfer Analysis. Luikov AV (1975) Systems of differential equations of heat and mass transfer in capillary-porous bodies. International Journal of heat and Mass Transfer 8, 1-14 Luikov AV (1980) Heat and Mass Transfer, Mir Publishers. Moscow, 623 pp. Mathelin L, Hussaini MY, Zang TA, Bataille,F (2004) Uncertainty Propagation for Turbulent. Compressible Nozzle Flow Using Stochastic Methods. AIAA Journal;42(8):1669–1676. Mathelin L, Hussaini MY, Zang TA (2005) Stochastic Approaches to Uncertainty Quantification in CFD Simulations. Numerical Algorithms;38(1):209–236. Nouy A (2009) Recent Developments in Spectral Stochastic Methods for the Numerical Solution of Stochastic Partial Differential Equations. Arch. Comput. Methods Eng.;16:251-285. Perez R, Walters R (2005) An Implicit Compact Polynomial Chaos Formulation for the Euler Equations; AIAA-Paper 2005-1406, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, CD-ROM. Reagan M, Najm HN, Ghanem RG, Knio OM (2003) Uncertainty Quantification in Reacting Flow Simulations through Non-Intrusive Spectral Projection. Combustion and Flame;132:545–555. Vos P (2006) Time-Dependent Polynomial Chaos. Master of Science Thesis. Walters R (2003) Towards stochastic fluid mechanics via Polynomial Choas-invited AIAA-Paper 2003-0413, 41st AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, CD-ROM. Wan X, Karniadakis GE (2006) Beyond Wiener-Askey expansions: Handling arbitrary PDFs, J. Sci. Comp.;27(1-3):455–464. Whitaker S (1977) Simultaneous heat, mass and momentum transfer in porous media: A theory of drying. In: Hartnett J.P., Irvine T.F. Jr. (eds.): Advances in heat transfer, Academic Press, Vol. 13, 119-203 Wiener N (1938) The Homogeneous Chaos. American Journal of Mathematics;60(4):897–936. Wittveen JAS, Bijl H (2006) Modeling arbitrary uncertainties using Gram-Schmidt polynomial chaos, 44th AIAA Aerospace Sciences Meeting and Exhibit. Xiu D, Karniadakis GE (2003) Modeling Uncertainty in Flow Simulations via Generalized Polynomial Chaos. Journal of Computational Physics;187(1):137– 167. Xiu D, Karniadakis GE (2002) The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal of Scientific Computing;24:619–644. 84
Spectral stochastic finite element approach to effective uncertainty quantification in heat conduction problems Miroslav Trcala Department of Wood Science, Mendel University, Brno, Czech Republic Abstract
This paper deals with the stochastic analysis of heat transfer in building envelope with the random heat conductivity of building materials. The spectral solution of this problem is based on discretization the resulting temperature field in the stochastic dimension by the orthogonal polynomials. A Galerkin projection is applied in the stochastic dimension to obtain the deterministic set of partial differential equations that is solved by FEM. The main purpose of this paper is to demonstrate that the stochastic spectral method based on polynomial chaos expansion has a big advantage over the Monte Carlo method in terms of computational time. Keywords: spectral stochastic finite element method; polynomial chaos; random heat conductivity; Monte Carlo method
1. Introduction
Numerical simulations in building physics rely on obtaining accurate physical properties of building materials. In many situations, it is difficult to obtain the accurate values of material properties. Most of building envelopes are porous media and their structure is complex and material properties are difficult to determine. Thus, uncertainty quantification and propagation in physical systems appear as a critical path for the improvement of the prediction of their response [20]. In the last two decades, a growing interest has been devoted to a new family of methods, called spectral stochastic methods, for the propagation of uncertainties through physical models governed by stochastic partial differential equations. In this paper we focus on uncertainty propagation with the polynomial chaos (PC) method, which is based on the spectral representation of the uncertainty. The biggest advantage of the spectral PC method is given mainly by low computational costs. Heat conductivity is one of the major thermodynamic coefficients in the study of energy transport through complex systems due to its increasing importance in agronomy, geothermal and geo environmental engineering [23]. Uncertainties in computer models arise from a number of sources, and can be defined in different ways depending on the scale and application of the model [1]. Stochastic analysis in relation to building energy modelling has mainly been restricted to the scale of an individual building [11, 15, 16, 17, 24, 29]. Several methods have been used to model and propagate uncertainty in stochastic computational simulations [12]. Several researchers have studied and implemented the PC approach for a wide range of problems. Ghanem and Spanos [7] and Ghanem [9] 85
applied the PC method to several problems of interest to the structures community. Mathelin et al. [18] studied uncertainty propagation for a turbulent, compressible nozzle flow by this technique. Xiu and Karniadakis [31] analyzed the flow past a circular cylinder and incompressible channel flow by the PC method and extended the method beyond the original formulation of Wiener [28] to include a variety of basis functions. In 2003, Walters applied the PC method to a two-dimensional steady-state heat conduction problem for representing geometric uncertainty. Following this effort, an implicit compact PC formulation was implemented for the stochastic Euler equations [21]. The PC method for the propagation of uncertainty in computational simulations involves the substitution of uncertain variables and parameters in the governing equations with the polynomial expansions. In general, an intrusive approach will calculate the unknown polynomial coefficients by projecting the resulting equations onto basis functions (orthogonal polynomials) for different modes. As its name suggests, the intrusive approach requires the modification of the deterministic code and this may be inconvenient for many complex computational problems since the modification of the existing code can be difficult, expensive, and time consuming. To alleviate this inconvenience, non-intrusive methods have been investigated by many researchers. Most of the non-intrusive PC approaches in the literature are based on sampling [6, 22, 13] or quadrature methods [6, 19] to determine the projected polynomial coefficients. Before their introduction, the standard was to employ Monte Carlo (MC) methods to sample the inputs, compute the output for each sample, and aggregate statistics of the output [10]. In fact, this is still the most widely used method in practice due to its robustness and ease of implementation. However, the MC methods suffer from a dreadfully slow convergence rate proportional to the inverse of the square root of the number of samples. And if each sample evaluation is expensive, such as the solution of a partial differential equation, then obtaining hundreds of thousands of samples may be entirely infeasible. Thus, the initial applications of spectral methods showed orders-of-magnitude reduction in the work needed to estimate statistics with comparable accuracy [32]. Such results spurred interest in applying spectral methods to differential equations with stochastic (i.e. parameterized) inputs [4]. This paper presents a spectral polynomial chaos method for the propagation of input uncertainty in heat conductivity of building materials. 1.2 Generalized polynomial chaos
Generalized polynomial chaos (gPC) is employed to represent stochastic processes. Stochastic process can be seen as a process involving some randomness. They can be represented by a stochastic mathematical model, often expressed in terms of stochastic differential equations. Stochastic mathematical models are based on a probability space (Ω, F , P) where Ω denotes the space of elementary events, F ⊂ 2Ω its σ-algebra of events, and P its probability measure. In addition considering some physical domain D ⊂ d × T (d=1, 2, 3), which can be a combination of spatial and temporal dimensions, a stochastic process can be seen as a scalar or vector valued function u(x, t , ω ) : D × Ω → b where x is an element of the physical space, t denotes the time and ω is a point in the sample space. Furthermore, because of the infinite-dimensional nature of the probability space, we discretize this space by characterizing it by a finite number of random variables 86
{ξ (ω )} j
. This can be seen as assigning a finite number of coordinates {ξ j }
N
, N∈ j =1
N j =1
to the probability space reducing it to a finite dimensional space Λ ⊂ N . Consequently, the stochastic process u becomes a mapping u(x, t , ξ ) : D × Λ → b . It is important to note that in this work, we assume that the occurring stochastic processes are already characterized by a known set of random variables. Wiener was the first to represent stochastic processes by orthogonal polynomial expansions [28]. To accomplish this, he used Hermite polynomials in terms of Gaussian random variables to represent Gaussian processes, which is referred to as homogeneous chaos. In this way, the stochastic process is represented in the form: ∞
u(x, t , ξ (ω )) = ∑ ui (x, t ) H i (ξ (ω ))
(1)
i=0
in which H i are Hermite polynomials and ξ is a vector of Gaussian random variables with zero mean and unit variance. It is a spectral expansion in the random dimensions employing deterministic coefficients. According to the Cameron-Martin theorem [2], for a fixed value of x and t, this expansion converges to any L2 (Ω) functional in the L2 (Ω) sense. This implies that the application of polynomial chaos is restricted to those stochastic processes yielding
∫ω
∈Ω
2
u(x, t , ω ) dP(ω ) < ∞.
(2)
As a result, polynomial chaos is restricted to second-order stochastic processes, i.e. processes with finite second-order moments. These are processes with finite variance, and this applies to most physical processes. Although Wiener’s original polynomial chaos expansion converges to any second-order stochastic process, it is most suitable to represent Gaussian processes, due to the random variable’s Gaussian nature, yielding a fast convergence. In order to deal with a broader range of stochastic processes, the Wiener-Hermite chaos has been generalized to the generalized polynomial chaos [32], also referred to as Wiener-Askey polynomial chaos. Analogously, gPC is a means of representing second-order stochastic processes through the expansion: ∞
u(x, t , ξ (ω )) = ∑ ui (x, t )Φ i (ξ (ω )).
(3)
i=0
Here the random trial base {Φ i (ξ(ω ))} exists out of orthogonal polynomials from the Askey-scheme, of which the Hermite polynomials are a subset, in terms of a random vector ξ = {ξi (ω )} j =1 . The combination of random vector and polynomials is carefully N
selected based on the distribution of the random input. It seems that for certain random variables, their probability density function (PDF) uniquely corresponds to one of the weighting functions ω (ξ ) in the orthogonality relation of the different orthogonal polynomials of the Askey-scheme. Choosing the corresponding combination leads to a 87
proper gPC expansion. Since each of the polynomials of the Askey-scheme forms a complete basis in the Hilbert space determined by their corresponding support, it is expected, according to Xiu and Karniadakis [32], that each type of Wiener Askey expansion converges to any L2 (Ω) functional in the L2 (Ω) sense in the corresponding Hilbert functional space as a generalized result of Cameron-Martin theorem. This generalized polynomial chaos has later been further generalized for arbitrary random inputs [27, 30]. In this way, gPC is not restricted anymore to random inputs of standard types. In order to expand the solution in a polynomial chaos expansion in terms of this arbitrary random variables, a proper random trial base {Φ j (ξ )} should be created, and this according to the gPC rules. This means that a set of orthogonal polynomials is created in such a way that they are orthogonal with respect to a weighting function which equals the probability density function of the random input. Different algorithms can be employed for this orthogonalization, among others the Stieltjes procedure, the Lanczos algorithm [27] or the Gram-Schmidt orthogonalization [30]. Taking this further generalization into account, one can truly speak of generalized polynomial chaos. The polynomials of the gPC’s random trial base satisfy following orthogonality relation Φ i Φ j = Φ i2 δ ij where δ ij is the Kronecker delta and ⋅ denotes the ensemble average. This ensemble average is the inner product in the Hilbert space determined by the measure of the random variables f (ξ) g (ξ) = ∫ f (ξ ) g (ξ )dP(ω ) = ∫ f (ξ ) g (ξ)ω (ξ)dξ with ω (ξ ) the ω∈Ω
corresponding weighting function and with integration in the last integral taken over the support of ξ . With the introduction of the gPC, optimal performance can now be achieved for a broader range of stochastic processes. In order to accomplish this, it is important to choose the appropriate type of chaos depending on the stochastic process’ nature. Although every type of chaos converges to any second-order process, choosing the appropriate polynomial chaos expansion leads to an optimal solution. This has been presented in [32], where it is computationally shown that exponential convergence is achieved when applying the optimal gPC through a Galerkin projection in random space on a stochastic ordinary differential equation [25]. 1.2 Galerkin projection
The stochastic problem can be formulated as: find the stochastic process u ( x, t , ξ (ω )) , which solves the stochastic differential equation:
L ( x, t , ξ; u ) = f ( x, t , ξ ) ,
(4)
where f ( x, t , ξ ) is a source term and L is a differential operator involving differentiation in space and time. In the spectral expansion of the solution u ( x, t , ξ (ω )) Q
u(x, t , ξ (ω )) = ∑ u i (x, t )Φ i (ξ (ω ))
(5)
i=0
88
the unknown deterministic coefficients ui (x, t ) , which represent the different modes of the solution, need to be determined. To accomplish this, start with substituting the expansion for the solution (5) into the governing equation (4) Q L x, t , ξ ; ∑ u i Φ i = f ( x, t , ξ ) i=0
(6)
Next, successively multiplying this equation by the different orthogonal polynomials of the finite expansion, and taking the statistical average Q L x, t , ξ; u ∑ ui Φ i , Φ j = f ( x, t , ξ ) , Φ j i=0
j=0, 1, …, Q
(7)
results in a set of (Q + 1) coupled equations for the different random modes ui (x, t ) . The total number of expansion terms (Q+1) is determined by the dimension N of the random vector ξ and the highest degree P of the orthogonal polynomials {Φ i } Q +1 =
( N + P )!.
(8)
N !P!
As a result, the governing set of equations for the expansion coefficients ui (x, t ) is completely deterministic. This system should further be solved in space and time by finite element method or other numerical solvers.
2. Governing equations
The steady-state of temperature field is generally described as follows:
∇λ∇T = 0
(9)
where λ is the matrix of thermal conductivity coefficients of wood [Wm-1K-1] and T is the temperature [K]. Boundary conditions are obtained from the heat transfer between envelope surfaces and external air. The following boundary conditions were applied for the numerical solution: −nλ∇T = αT (T∂Ω − Tair )
(10)
where αT is the heat transfer coefficient [Wm-2K-1], T∂Ω is the surface temperature [K], Tair is the air temperature [K]. By using the generalized polynomial chaos expansion, the temperature can be expand as 89
Q
T (x, ξ (ω )) = ∑ Ti (x)Φ i (ξ ) .
(11)
i=0
The random heat conductivity has the form λ → λ + σξ ,
(12)
where ξ is an uniform random variable ξ U ( −1,1). After substituting the terms (11, 12) we can rewrite the above equations (9, 10) as Q ∇ ( λ + σξ ) ∑ ∇Ti (x)Φ i (ξ ) = 0 i =0
(13)
Q Q n ( λ + σξ ) ∑ ∇Ti (x)Φ i (ξ ) = αT ∑ Ti (x)Φ i (ξ ) − Text . i =0 i=0
(14)
By Galerkin projection of the above equation onto polynomial basis {Φ j (ξ )} for each
j = {0,1,..., Q} we obtain this set of (Q + 1) coupled partial differential equations
∑ ∇ ( λ Q
i =0
)
Φ i (ξ ), Φ j (ξ ) + σ ξΦ i (ξ ), Φ j (ξ ) ∇Ti ( x) = 0
(15)
Q n λ Φ ( ξ ), Φ ( ξ ) + σ ξ Φ ( ξ ), Φ ( ξ ) ∇ T ( x ) = α ∑ i j i j i T ∑ Ti ( x ) Φ i (ξ ), Φ j (ξ ) − Text 1, Φ j (ξ ) i =0 i =0 (16) and for subdomain with deterministic heat conductivity Q
∑ ∇ ( λ Q
i =0
)
Φ i (ξ ), Φ j (ξ ) ∇Ti (x) = 0
(17)
Q n λ Φ i (ξ ), Φ j (ξ ) ∇Ti (x) = αT ∑ Ti (x) Φ i (ξ ), Φ j (ξ ) − Text 1, Φ j (ξ ) ∑ i =0 i =0 Q
. (18)
Once we have the coefficients of the expansion, we can compute approximate statistics of the solution with the formulas for mean value Q
Q
i=0
i =0
E[T (x, ξ (ω ))] = E[∑ Ti (x)Φ i (ξ )] = ∑ Ti (x)E[Φ i (ξ )] = T0 (x) ,
for variance 90
(19)
Q
Q
Var[T (x, ξ (ω ))] = E[(∑ Ti (x)Φ i (ξ ) − T0 (x)) 2 ] = E[(∑ Ti (x)Φ i (ξ )) 2 ] = i=0
Q
i =1
(20)
Q
= ∑ Ti (x)E[Φ i (ξ ) ] = ∑ Ti (x) Φ i (ξ ), Φ i (ξ ) 2
2
2
i =1
i =1
and finally for standard deviation
Std [T (x, ξ (ω ))] =
Q
∑T i =1
i
2
(x) Φ i (ξ ), Φ i (ξ ) .
(21)
3. Example
Numerical simulations are performed for example of steady-state heat transfer through building envelope construction. Geometry and materials of the building envelope construction are shown in figure 1 and table 1. Boundary conditions are obtained from the heat transfer between envelope surfaces and external air. The heat transfer coefficient is 20 Wm-2K-1, exterior temperature is -15 °C and interior temperature 20 °C.
λ
layers
Dim
x Dim
(W.m1.K1) (cm)
(cm)
R1-plasterboard R2-chipboard R3-mineral wool R4-spruce timber R5-chipboard R6-polystyrene
0.22 0.11 0.042 0.15 0.11 0.04
30 30 23 7 30 30
Figure 1: Geometry
1.3 1.3 14 14 1.3 10
y
Table 1: Materials
Now replace some deterministic model inputs with stochastic quantities that better represent my understanding of the problem. One example of this may be an equation parameter (heat conductivity coefficient) that admits some measurement error, which we can model with a stochastic representation. Consider the three uniform random variables in our example. First is heat conductivity of mineral wool defined on the interval [0.032,0.052], second is heat conductivity of spruce timber defined on the interval [0.1,0.2] and third is heat conductivity of polystyrene defined on the interval [0.03,0.05].
91
Convergent results of spectral solution were obtained by third-order polynomial chaos expansion with three-dimensional random vector ξ . The total number of ( 3 + 3)! = 20. expansion terms is Q + 1 = 3!3! Many programming languages have the ability to generate pseudo-random numbers which are effectively distributed according to the standard uniform distribution U(0,1). If u is a value sampled from the standard uniform distribution, then the value a+(b−a)u follows the uniform distribution parameterized by a and b, denoted as U(a,b). This is used for Monte Carlo solution of our problem. The main purpose of this example is to perform spectral solution of the mentioned heat conduction problem and to illustrate that it is more effective approach than the Monte Carlo method. 4. Results and discussions
The above derived system of PDEs was solved numerically by FEM solver.
Figure 2: Triangular mesh, Lagrange-Quadratic elements, Number of degrees of freedom 162980
The figures below compare the mean and standard deviation of the solution using stochastic spectral method (polynomial chaos expansion, Galerkin projection) and straightforward Monte Carlo method.
92
(a)
(b)
Figure 3: Spectral solution of stochastic problem for (a) mean and (b) standard deviation
(a)
(b)
Figure 4: Monte Carlo solution of stochastic problem (2000 samples) for (a) mean and (b) standard deviation
A 2000 realizations were needed to obtain converged Monte Carlo results that agree well with those of chaos expansion as is shown in Figure 4. Convergent results of spectral solution were obtained by third-order polynomial chaos expansion with threedimensional random vector ξ .
93
(a)
(b)
Figure 5: Contours of difference between spectral and Monte Carlo solution for (a) mean and (b) standard deviation
A more precise comparison between spectral solution and Monte Carlo solution is shown in Figure 5. It is seen that the differences of the mean as well as standard deviation across the domain are small and the maximum values are in the area of spruce timber. The mean and the standard deviation obtained with the stochastic spectral method are in good agreement with the mean and the standard deviation of the Monte Carlo simulations. 5. Conclusion
A stochastic spectral approach to model uncertainty in building heat conduction problems has been developed. The generalized polynomial chaos and Galerkin projection were applied to the solution of these problems. Monte Carlo simulations were also conducted to verify the spectral method and there was shown that the results are in good agreement with the results of the Monte Carlo simulations. This spectral approach has advantages over the Monte Carlo approach mainly in terms of computational time. This method can be used for a wider range of the building physics problems and deserves further research. Acknowledgement
This work was supported by the Internal Grant Agency (IGA) of the Faculty of Forestry and Wood Technology, Mendel University in Brno, project no. 7/2011.
94
References
[1] Booth AT, Choudhary R, Spiegelhalter DJ. Handling uncertainty in housing stock models. Building and Environment 2011;48:35-47. [2] Cameron R, Martin W. The orthogonal development of nonlinear functionals in series of Fourier-Hermite functionals. Ann. Math. 1947;48(2):385-392. [3] Chen P, Pei DCT. A mathematical model of drying processes. Int. J. Heat Mass Transfer 1989;32(2):297-310. [4] Constantine PG. Spectral methods for parameterized matrix equations. A dissertation 2009. [5] Crank J. Mathematics of diffusion. Cranendon Press 1975. [6] Debusschere BJ, Najm HN, Pebay PP, Knio OM, Ghanem RG, Maitre OPL. Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes. SIAM Journal on Scientific Computing 2004; 26(2):698– 719. [7] Ghanem R, Spanos PD. Polynomial Chaos in Stochastic Finite Elements. Journal of Applied Mechanics 1990;57:197–202. [8] Ghanem R. Stochastic Finite Elements with Multiple Random Non-Gaussian Properties. Journal of Engineering Mechanics 1999:26–40. [9] Ghanem RG. Ingredients for a General Purpose Stochastic Finite Element Formulation. Computational Methods in Applied Mechanical Engineering 1999;168:19–34 [10] Helton JC, Davis FJ. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliability Engineering & System Safety 2003;81(1):23–69. [11] Hopfe C. Uncertainty and sensitivity analysis in building performance simulation for decision support and design optimization, PhD thesis, Technische Universiteit Eindhoven 2009. [12] Hosder S, Walters RW, Perez R. A Non-Intrusive Polynomial Chaos Method For Uncertainty Propagation in CFD Simulations. 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada 2006. [13] Isukapalli SS. Uncertainty Analysis of Transport-Transformation Models. PhD Dissertation 1999. [14] Lewis RW, Morgan K, Thomas HR. The Finite Element Method in Heat Transfer Analysis 1996. [15] Lomas KJ, Eppel H. Sensitivity analysis techniques for building thermal simulation programs. Energ Build 1992;19:21-44. [16] Macdonald I. Quantifying the effects of uncertainty in building simulation. PhD thesis. University of Strathclyde 2002. [17] Macdonald I, Clarke J. Applying uncertainty consideration to energy conservation equations. Energ Build 2007;39:1019-26. [18] Mathelin L, Hussaini MY, Zang TA, Bataille,F. Uncertainty Propagation for Turbulent. Compressible Nozzle Flow Using Stochastic Methods. AIAA Journal 2004;42(8):1669–1676. [19] Mathelin L, Hussaini MY, Zang TA. Stochastic Approaches to Uncertainty Quantification in CFD Simulations. Numerical Algorithms 2005;38(1):209–236.
95
[20] Nouy A. Recent Developments in Spectral Stochastic Methods for the Numerical Solution of Stochastic Partial Differential Equations. Arch Comput Methods Eng 2009;16:251-285. [21] Perez R, Walters R. An Implicit Compact Polynomial Chaos Formulation for the Euler Equations 2005; AIAA-Paper 2005-1406, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, CD-ROM. [22] Reagan M, Najm HN, Ghanem RG, Knio OM. Uncertainty Quantification in Reacting Flow Simulations through Non-Intrusive Spectral Projection. Combustion and Flame 2003;132:545–555. [23] Singh R, Bhoopal RS, Kumar S. Prediction of effective thermal conductivity of moist porous materials using artificial neural network approach. Building and Environment 2011;46:2603-2608. [24] Struck C, Hensen J. On supporting design decisions in conceptual design addressing specification uncertainties using performance simulation;10th International IBPSA Conference, Tshingua University Beijing/China. [25] Vos P. Time-Dependent Polynomial Chaos. Master of Science Thesis. 2006. [26] Walters R. Towards stochastic fluid mechanics via Polynomial Choas-invited 2003 AIAA-Paper 2003-0413, 41st AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, CD-ROM. [27] Wan X, Karniadakis GE. Beyond Wiener-Askey expansions: Handling arbitrary PDFs, J. Sci. Comp. 2006;27(1-3):455-464. [28] Wiener N. The Homogeneous Chaos. American Journal of Mathematics 1938;60(4):897–936. [29] de Wit S, Augenbroe G. Analysis of uncertainty in building design evaluations and its implications. Energ Buid 2002;34:951-8. [30] Wittveen JAS, Bijl H. Modeling arbitrary uncertainties using Gram-Schmidt polynomial chaos 2006 44th AIAA Aerospace Sciences Meeting and Exhibit. [31] Xiu D, Karniadakis G E. Modeling Uncertainty in Flow Simulations via Generalized Polynomial Chaos. Journal of Computational Physics 2003;187(1):137–167. [32] Xiu D, Karniadakis G. The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal of Scientific Computing 2002;24:619–644.
96
Sap Flow and Heat Conductivity Measurement and Their Interaction Trcala M., Čermák J. Faculty of Forestry and Wood Technology Mendel University Brno, Czech Republic Keywords: sap flow, heat conductivity, linear heat balance, heat field deformation, trunk heat balance Abstract There are many different methods for sap flow measurement. Quantitative output of such methods is a precondition of their application. However it is not easy to quantify, if not all needed input parameters are known. Heat conductivity is one of them. Introducing a constant value into corresponding equations has usually solved this problem, although this means introducing an unknown error. We have been working with two approaches: trunk heat balance method, THB, giving total flow rate per tree and the heat field deformation method, HFD, giving total, but also particular flows in different sapwood depths. There was a tendency to unify both these useful approaches, which led to development of intermediate (stationary at present), linear heat balance method, LHB (2010). The applied main equation (LHB formula) is derived from fundamental partial differential equation describing convective-conductive heat flow. A new approach was developed most recently, capable to calculate both parameters (sap flow and heat conductivity) at the same time, which is being applied and gradually improved now. This approach is based on two equations with the above-mentioned parameters. One includes heat balance around a heater; the other one includes deformation the heat field around a heater. The LHB method is illustrated on the example of a series of diurnal courses of sap flow calculated on the basis of different principles. INTRODUCTION This paper introduces new approach to sap flow measurement underlining influence of the heat conductivity of fresh sapwood on these measurements. There was an attempt to unify the trunk heat balance method and heat field deformation method. The aim of this work was to find a way how to calculate sap flow by means of heat balance around linear heater inserted into fresh sapwood with respect to the heat conductivity of fresh sapwood. This new approach, the linear heat balance method, LHB, combines both above mentioned methods. MATERIAL AND METHODS Experimental site and species The experimental site (stand 933/7A), was situated at Forest district Sobesice, University Training Forest Enterprise in Krtiny. Altitude 360 m, slight southward slope. Mean annual temperature 7.8 oC, mean annual precipitation 583 mm (360 mm over growing season). Well-illuminated sites with soil subjected to desiccation, illimerized soil type with an admixture of loess to a depth of 30-40 cm, medium supply of nutrients,
97
pH (H2O) 7.1 at surface, lower down 5.0 Fraction of physical clay falls in a downward direction to the subsoil from 19 to 9 %. Carpini querceta forest-type group in weathered granodiorite with a cover of loess in the anthropically strongly influenced suburban area around Brno. The herb layer does not exhibit clearly dominant species, prominent species being Poa nemoralis, Poa pratensis ssp.angustifolia Brachypodium sylveticum and Festuca heterophylla, with a large fraction of thermophilous and oligotrophic species (Vašíček et al. 1984). Mature Douglas-fir tree (Pseudotsuga menziessii (Mirb.) Franco) was selected as the sample tree for experiments (diameter at brest height 411 mm, bark thickness 11 mm). Common sap flow measurement methods Our study was based on evaluation of the trunk heat balance THB (Čermák et al., 1973, 2004, Kučera et al., 1977, Tatarinov et al. 2005) and the heat field deformation HFD (Nadezdhina et al., 1998, 2006, 2010), methods. Sensors were inserted at the same height at the tree stem. Linear heat balance method (LHB), theory The new LHB approach was developed (stationary at present) to analyse the sap flow with respect to the heat conductivity. The appropriate formula was derived from general partial differential equation (PDE) for heat conduction and convection
∂T ∂T = ∇λ∇T + cwQw + P, ∂t ∂y where T is temperature (K), Qw is sap flow (kg.m-2.s-1), λ is matrix of heat conductivity coefficients of fresh sapwood (W.m-1.K-1), ρ is density of fresh sapwood (kg.m-3), cx is specific heat of fresh sapwood (J.kg-1.K-1), cw is specific heat of water (J.kg-1.K-1), P is heat power (W.m-3). The mentioned stationary approach is more simple than nonstationary because the partial derivative of temperature with respect to time is equal to zero and there is no reason to know density and specific heat of sapwood here. Then the PDE has the following form
ρ cx
∂ 2T ∂ 2T ∂T 0 = λx 2 + λ y 2 + cwQw +P ∂x ∂y ∂y and it is possible to solve this equation, after substituting λ y = 2λx and modifications we gets Qw cw y 4 λx
Qw cw (2 x 2 + y 2 ) ), 4λ x π 8λx where T(x,y) is temperature (K) at point [x, y], cw is specific heat of water (J.kg-1.K-1), K0 is the modified Bessel function of the second kind. This formula was used by following way (LHB formula): T ( x, y ) =
Tupper −
Pl
e
Tright + Tleft 2
K0 (
Qw4cwλ0.01 Q c 0.01 Q c 0.01 2 = ) − K0 ( w w ) e x K 0 ( w w 4λ x 4λ x π 8λx Pl
98
Linear heat balance method (LHB), experiments The LHB sensor was partially similar to HFD ones, but differed in several points. It consisted of the linear heater 45 mm of operating length below bark surface (a resistance wire), inserted in the stainless hypodermic needle 1.5 mm in outer diameter (inner 1.0 mm) situated in the sapwood and four thermocouples copper-constantan, with measuring points 5 mm from needle tip in similar needles and situated at the distance of 10 in axial direction (up and down the heater) and 10 mm in tangential one (right and left side of the heater). Insulated reference ends of all thermocouples were situated at the same point (brass tube 5 mm in diameter, 15 mm long, surrounded by a ceramic block 12 mm thick, insulated by a layer of bubbled polyethylene) – therefore when temperature of the reference point is known, thermocouples measured absolute temperature. The sensor was installed at breast height from North stem side, so that the depth of sensing ends was 15 mm beneath the cambium in sapwood (needle tips 20 mm). The heater tip was inserted 15mm deeper. Total resistance of the heater wire was 60.5 Ω (204.8 Ω m-1) and applied direct current voltage of 4.24 (low heating) or 5.83 (high). this gives the power (U2/R) = 0.297 or 0.652 W, respectively. Simultaneous calculation of sap flow and heat conductivity An idea of simultaneous calculation of sap flow and heat conductivity from measured temperature differences around linear heater is introduced. This idea is inspired by mathematical theory saying that two unknown parameters (here sap flow and heat conductivity) is possible to obtain by solving of system of two independent equations including these two unknowns. Generally it looks like this:
F1 (Qw , λx ) = 0 F2 (Qw , λx ) = 0 In the above mentioned work, this system was designed Tupper −
Tright + Tleft 2
Qw − f (Qw , λx )
=
Pl
π 8λx
e
Qw cw 0.01 4 λx
Q c 2 × 0 + 0.012 ) Pl K0 ( w w e )− 4λ x π 8λx
Qw cw 0 4 λx
K0 (
Qw cw (2 × 0.012 + 02 ) ) 4λx
(dTas − dTsym )0 + dTsym − dTas
= 0. dTas First equation is from above mentioned LHB approach and second is modified well known HFD formula. HFD formula had to be modified because it doesn’t include heat conductivity. RESULTS AND DISCUSSION Numerical simulations According to generally know PDE and commonly used simulations based on this PDE, the above LHB method can be verified by these simulations. The influence of the sap flow on temperature field around linear heater is shown on these figures. These figures were obtained by simulations in COMSOL Multiphysics with MATLAB. These simulations are based on solving above PDE by finite element method (FEM). At the first figure you can see temperature field for sap flow equal to zero, at second 1.8, third for 5.4 and at the last fourth figure for 21.6 g.cm-2h-1. 99
Linear heat balance method (LHB) The derived LHB formula was verified by the simulations. The result of this step is at the graph (Fig. 2). The graph shows that the results from LHB formula are consistent across whole testing range of sap flux densities for higher heat conductivity but for low heat conductivity the formula is applicable only from 0 to 22 g.cm-2h-1. On the graph (Fig. 3) there is an application of LHB on real data. The graph shows the time curve of sap flow for different values of heat conductivity. The LHB method delimitates the range of the heat conductivities from 0.4 to 0.5 W.m-1.K-1 and thus most probably mean value of the heat conductivity for our real data is 0.45 W.m-1.K-1. Simultaneous calculation of sap flow and heat conductivity The system of two equations with two unknowns (sap flow and heat conductivity) was designed. First equation is from LHB approach and second is modified HFD formula by including heat conductivity (this is the topic of the following paper). Thermal diffusivity in classical HFD formula lead us to introduce a function of heat conductivity and sap flow f (Qw , λx ) . This function was found by simulations (Fig. 4). The modified HFD formula and subsequently the whole system were verified by simulations (Fig. 5). As you can see on the graph, the system is applicable only in the range of 7.2 and 21.6 g.cm-2h-1 at present (Fig. 5). CONCLUSIONS The paper presented two new approaches to sap flow measurement based on interaction sap flow and heat conductivity of fresh sapwood. First, a new LHB method was developed for sap flow measurement underlining the importance of heat conductivity. This was based on theoretic analysis, simulations and application on real measured data. The advantage of this approach is exact formula derived from physical law of eat conservation and reflecting the influence of the heat conductivity during sap flow measurement. The disadvantage is a difference between the real and theoretical radius of the linear heater. Second, an idea of simultaneous calculation of sap flow and heat conductivity was introduced. This is based only on theoretic analysis and simulations. The advantage is possibility of estimating sap flow together with heat conductivity. The disadvantage is that it is not easy to find a system of equations applicable to most sap flow values. ACKNOWLEDGEMENTS The study was supported by the project Research Intention of Czech MSM 6215648902. Literature Cited Čermák J. Deml M. Penka M. 1973. A new method of sap flow rate determination in trees. Biol. Plant. (Praha) 15(3): 171-178. Čermák J., Kučera J. and Nadezhdina N. 2004. Sap flow measurements with two thermodynamic methods, flow integration within trees and scaling up from sample trees to entire forest stands. Trees, Structure and Function 18:529-546 Kučera J. Čermák J. Penka M. 1977. Improved thermal method of continual recording the transpiration flow rate dynamics. Biol. Plant. (Praha) 19(6): 413-420. 100
Nadezhdina N. and Čermák J. 1998. "The technique and instrumentation for estimation the sap flow rate in plants". Patent No.286438 (PV-1587-98). Nadezhdina N. Čermák J. Nadezhdin V. 1998. Heat field deformation method for sap flow measurements. Proc. 4th. International Workshop on Measuring Sap Flow in Intact Plants. Židlochovice, Czech Republic, Oct.3-5, 1998. 72-92 pp. IUFRO Publications. Publishing house of Mendel Univ.Brno. Nadezhdina N, Čermák J, Neruda J, Prax A, Ulrich R, Nadezhdin V, Gašpárek J, Pokorný E. 2006: Roots under the load of heavy machinery in spruce trees. European J.For.Res. 125: 111-128. Nadezhdina N., David T.S., David J.S., Ferreira M.I., Dohnal M., Tesar M., Gartner K., Leitgeb E., Nadezhdin V., Čermák J., Jimenez M.S. and Morales M. 2010. Trees never rest: the multiple faces of hydraulic redistribution. Ecohydrology 3: 431-444. (DOI: 10.1002/eco.148) Tatarinov F.A., Kucera J., Cienciala E. 2005. The analysis of physical background of tree sap flow measurements based on thermal methods. Measurements Science and Technology 16: 1157-1169. Vašíček (ed.) 1984: Ecophysiological and ecomorphological studies of individual trees in the spruce ecosystems of the Drahanská vrchovina uplands (Czechoslovakia). Folia Universitatis Agriculturae, facultas Silviculturae, Brno, Czechoslovakia.
Figures
Zero
Low
Medium
High
Fig. 1. Numerical simulations of temperature field for different values of sap flow Qw in COMSOL Multiphysics with MATLAB (FE analysis)
101
Fig. 2. Verification of LHB formula
Fig. 3. Sap flow measurement for different values of heat conductivity
Fig. 4. Function in the modified HFD formula 102
Fig. 5. Verification of the modified HFD formula (left) and verification of system (right)
103
Improvement of the trunk heat balance method including measurement of zero and reverse sap flows Trcala M.1, Čermák J.2 1
Department of Wood Science, 2Department of Forest Botany, Dendrology and Geobiocenology, Mendel University in Brno, Brno, Czech Republic
[email protected];
[email protected]
Abstract The trunk heat balance (THB) method for tree sap flow measurement has been successfully applied for many years in different field conditions and tree species. When using the classical equations, calculation of flow was done considering heat losses expressed as so-called “fictitious flow” (Qfic), which is visible on recorded data. There were two weak points of this approach: (1) Qfic was assessed only approximately from the minimum values of sap flow and (2) dependence of the heat loses from the measuring point, expressed on the actually recorded flow data (Qrec) was not taken into account till now. The aim of this work was (1) to improve the THB method by more accurate evaluating of fictitious flow similarly as it is done when working with the heat field deformation (HFD) method (2) In addition by considering of Qfic dependence on sap flow. This way the sap flow sensor with modified geometry (from vertical to horizontal) will also allow measurement of zero and reverse sap flows.
Keywords: sap flow; trunk heat balance method; fictitious flow; heat field deformation method; heat conductivity
104
Introduction Sap flow was found as a very important indicator of tree behavior, therefore a series of measurement methods for its measurements was developed. These methods are mostly based on thermodynamic principles, either pulsing systems or heat balance systems (Daum, 1967; Swanson, 1994; Čermák et al., 1973, 1982, 2004; Sakuratani, 1981, Campbell, 1991; Smith & Allen, 1996; Čermák & Nadezhdina, 1998; Wullschleger et al., 1998; Kostner et al., 1998; Burgess et al., 2001). Heat balance methods giving quantitative data and not requiring any calibration or measurements of additional variables (e.g., xylem water content) appeared very suitable for this purpose. Variant of the trunk heat balance (THB) method with internal heating and sensing of temperature, designated especially for large trees was developed in early seventies (Čermák & Deml, 1972/74, Čermák et al., 1973, 1976, 1982, 2004; Kučera, 1977, Kučera et al., 1977), being theoretically analyzed later (Tatarinov et al., 2005). A certain section of a large tree trunk is heated from inside by the alternating electric current separated from the ground and passing through the tissues. Heat (Joule's heat) is released more uniformly within the bulk xylem tissue between electrodes situated at distances d=20 mm. Heat does not leave much through thick (insulating) bark neither enters much the heartwood, providing that it is dry enough (although e.g. wet rotten heartwood causes problems). Measured stem segments were 8 cm wide (when using older system of 5 electrodes) or 4 cm (when using 3 electrodes). Power (P) or temperature difference (dT) between the heated and reference points can be held constant. Application of constant power means higher and variable diurnal temperature differences (usually 2 to 5 K). Application of constant temperature difference instead of constant power (Kučera et al., 1977) improved sensor sensitivity and especially its dynamic properties and lead to decreasing of temperature differences at the measuring point (down about 1 to 2 °C). The less power consumption also means substantial spare of energy for heating, which is especially significant for long-term measurements in large trees. Proper installment of sensors can be done only on the basis of good knowledge of the anatomy and physiology of the conductive systems in trees. The temperature differences between heated and non-heated part of the measuring point were measured with a battery of thermocouples arranged between the electrodes, below them, at their sides or both in a way compensating the influence of naturally 105
changing temperature from the outside as well as inside stems this way (Čermák & Kučera, 1981). Thermocouples inserted within stainless steel hypodermic needles with prolonged metal ends measured temperature at short (several mm) abscissa. They were located in radial direction within the sapwood approximately through places where the minimum and maximum sap flow takes place in order to represent its mean value (Čermák et al., 1982). Improving temperature integration over sapwood depth was done in newer sensors by placing thermocouples in stainless steel needles 1 mm in diameter into direct contact with the metal of electrodes, when situated in the center of plate electrodes. The sap flow rate in the measuring point, Qm was expressed per 1 cm wide stem section (dividing by the distance of parallel outer electrodes) over the entire sapwood depth. In order to calculate the sap flow for the whole tree, Qt, the Qm was multiplied with a dimensional constant, particularly by the length of xylem circumference of the tree at the same height Lm. The measuring device was made by the Ecological Measuring Systems (EMS), Brno, Czech Republic. The THB method was tested on both shoots and large trees and found quantitative when comparing the results obtained with other methods such as gravimetric (Penka et al., 1979), gas-exchange (Schulze et al. 1985), volumetric and other techniques (Čermák & Kučera, 1990; Cienciala et al., 1992; Kostner et al., 1998). It was also applied as a standard for testing other methods (Offenthaler & Hietz, 1998; Schubert 1999; Lundblad et al., 2001). The method was found sensitive and suitable for any tree species with sapwood located in the outermost part of xylem. There is almost no limit of tree size, only too deep sapwood can cause problems. The sensors do not significantly interfere with growth and therefore the method can be used for measurement over whole growing seasons. The measured points represent significant parts of tree stems, almost an order larger than any only needle heater systems. Thermocouples occurring in direct contact with the metal of electrodes better integrate the radial pattern of sap flow across sapwood. Keeping low temperature differences and constant temperature gradients decreases impact of temperature on stem tissues and improves the dynamic properties of sensors. The dT const variant of the THB method may be applied to detect quantitatively rapid changes in sap flow up to about 30%.min-1 (Čermák & Kučera, 1990). The drawback of the method is that the internal heating of tissues may be 106
applied with difficulties on some species (e.g. tropical) with sapwood located deep inside of stems and in palms. The electrodes used in measuring points cause minor injuries to trees and should not be used after once loosing contact with tissues (e.g., due to freezing). The varying temperature or heat input requires sophisticated electronic control of each individual gauge. Application of needles in direct contact with the metal of electrodes increases substantially the value of Qfic. The important problem is, that usually applied vertical sensor geometry does not allow exact estimates of heat loses from the measured section caused by wood heat conductivity and also does not allow measurement of zero and reverse flows, which can be done by other methods (Daum, 1967; Nadezhdina et al., 1998; Sakuratani et al., 1999). Therefore the present study was focused at overcoming of the above-mentioned drawbacks and improvement of sap flow calculations through better involvement the heat loses.
Theoretical analysis Present state of the arts According to generally known partial differential equation describing temperature field for conduction-convection heat transfer in anisotropic medium, the simulations of volume heating for THB method can be performed. These simulations were obtained by COMSOL Multiphysics with MATLAB and are based on solving partial differential equation by the finite element method. The partial differential equation is following (Tatarinov et al., 2005):
ρc
∂T = qconduction + qconvection + P , ∂t
(1)
where
qconduction = ∇ ⋅ λ∇T , qconvection = cwQw
∂T ∂y
(2)
and thus
ρc
∂T ∂T = ∇ ⋅ λ∇T + cwQw +P, ∂t ∂y
(3)
where T is temperature (K), Qw is sap flow density (kg.m-2s-1), λ is matrix of heat conductivity coefficients of fresh sapwood (W.m-1.K-1), ρ is density of fresh sapwood (kg.m-3), C is specific heat of fresh sapwood (J.kg-1.K-1), cw is specific heat of water 107
(J.kg-1.K-1), P is heat power (W.m-3). For steady state the derivative of a temperature with respect to time is equal to zero and the basic equation for the calculation of the sap flux density is following
0 = ∇ ⋅ λ∇T + cwQw
∂T + P. ∂y
(4)
The heat balance equation is based on a volume (three-dimensional) heating of stem segment. Basically, the input energy has to be split between the conductive heat losses and the warming of water passing through, according to the equation (5). Most heat goes up with streaming water, but a certain part of it (about 10 to 20%) is lost by heating of stem tissues surrounding the measuring point, but nevertheless the heat loss takes place there. The fundamental differential equation of heat balance for a given measuring point, if considering a homogenous heat field in the plant part being measured can be expressed in the form (Kučera et al., 1977)
P = QdTcw + dT λ + ρ cdT ′
(5)
where P is the heat input into the measured part of xylem (W), Q is sap flow rate per measuring point (kg.s-1), ρ c = mwcw + mx cx , cw and cx is the specific heat of water (4186.8 J.kg-1.K-1) and of dry xylem (J.kg-1.K-1), mw and mx is the mass of water and of dry xylem in the measured plant part (kg), λ is a coefficient of heat transfer from measured plant part to its environment (W.K-1), dT (K) is temperature difference between heated and nonheated measuring point and t is time (s). Solution of the above equation gives the temperature difference dT in the form dT (t ) =
−t ( Qcw + λ ) P 1 − exp . ρc Qcw + λ
(6)
For steady state case described by the equation P = QdTcw + dT λ (7) we can derive the basic relation for the calculation of the sap flow rate through the measured part of xylem Q=
P λ − . dTcw cw
(8)
108
The heat loss magnitude (λ) is not clearly specified but is included in the socalled “fictitious flow” Qfic, which is recorded even when the actual flow is zero. The influence of possible changes in water content of tissues, dMw (a relation of dMw to the total mass of tree tissue) is of little significance. When calculating the actual sap flow rate Q, it is necessary to subtract Qfic. This can be approached e.g. (1) considering minimum Qrec values for each separate day (this causes underestimating error because saturating flows are eliminated), (2) periodically occurring (more days should be taken into account) early morning values after prolonged rain from the recorded flow data Qrec (this causes mostly overestimating errors). According to equation (8) we can write
Q = Qrec − Q fic ,
(9)
where
Qrec =
P λ , Q fic = dTcw cw
(10)
where dT = dTa + dTb and Q fic is for now unspecified quantity corresponding with unknown heat losses (λ) from heated volume by conduction.
New contribution The physical formulation of Q fic from standard heat conductivity coefficients follows: Q fic =
λ cw
=
λ x S x + λ y S y + λz S z dycw
,
(11)
where dy is unspecified theoretic distance on which is realized temperature difference dT and λx , λy , λz resp. S x , S y , S z are standard heat conductivity coefficients resp. cross-section of heated volume in axis x, y, z that are not also exactly known. The formula is difficult to apply because there are several unknown parameters, but it inspires us to introduce the idea of expressing the fictitious flow as a function of heat conductivity and sap flow. If we assume that λy = 2λx and that heat loss in z direction are small (zero) so we can write Q fic = λx
Sx + 2S y dycw
.
(12)
109
Therefore, the fictitious flow can be expressed as function of heat conductivity λx and sap flow rate Q (kg.s-1). Q fic = f (λx , Q ) .
(13)
The Qw.fic for sap flow density Qw (kg.m-2.s-1) can be approximately linearized (Fig. 3)
Qw. fic = a + bλx + cQw .
(14)
Thus the fictitious flow Qfic is equal to
Qfic = F1 + F2Q ,
(15)
where F1 = A ( a + bλx ) = (Q fic )0 and F2 = Ac is measure of dependence of Q fic on Qrec . Results obtained by simulations showed this value about 0.5 when compared to a classical variant of the THB method. It has an impact presumably on the afternoon flow values, reaching the daily maximum at this period of time. Value of F1 is derived from figures (Fig. 5), from which we can get the temperature difference dT, which is marked as (dT)0 in case of zero flow, described in the equation F1 = (Q fic )0 = ( Qrec )0 =
P ( dT )0 cw
(16)
this way we will get the required coefficient F1. After substitution Q = Qrec − Q fic into (15) we get this formula
Q fic = F1 + F2 ( Qrec − Q fic )
(17)
and when we express Qfic so we get
Q fic =
F1 F + 2 Qrec . 1 + F2 1 + F2
(18)
If we do not consider dependence of Qfic on Q (F2=0) we have in mind the Variant 1 Q = Qrec − F1 = Qrec − (Q fic )0
(19)
and if we do consider this dependence we get the best expression of sap flow rate called Variant 2 F Q = Qrec 1 − 2 1 + F2
F1 . − 1 + F2
(20)
110
Material and methods The impact of improved THB sensor using horizontal geometry was tested on the adult specimen of Douglas-fir (Pseudotsuga menziessii (Mirbel) Franco) with diameter at breast height, DBH = 41.0 cm, bark thickness = 2.0 cm, (smooth 1.2 cm) top height = 22 m, crown base height = 11 m, age about 60 years. The forest stand was situated at the altitude of 360 m, very slight southward slope. Mean annual temperature 7.8 oC, mean annual precipitation 583 mm (360 mm over growing season). Wellilluminated sites with soil subjected to desiccation, illimerized soil type with an admixture of loess to a depth of 30-40 cm, medium supply of nutrients, pH (H2O) 7.1 at surface, lower down 5.0. Fraction of physical clay falls in a downward direction to the subsoil from 19 to 9 %. The stand was characterized as Carpini querceta forest-type group in weathered granodiorite with a cover of loess in the anthropically strongly influenced suburban area around Brno. The herb layer does not exhibit clearly dominant species, prominent species being Poa nemoralis, Poa pratensis ssp.angustifolia Brachypodium sylvaticum and Festuca heterophylla, with a large fraction of thermophilous and oligotrophic herbaceous species. Dominance of species, which are semisciophytes to semiheliophytes, and a little fraction of heliophylous plants, together with the low incidence of sciophytes characterize a relatively high stand illumination. Dominance of plant species preferring low nitrogen content with a minority of species bound to medium-rich nitrogen soils. A fraction of about 10 % of species entirely nitrophyllous occurring mainly due to weeds penetration from the neighboring agrocenoses (Vašíček et al., 1984). Classical THB sensor (Čermák et al., 1973, 1982; Kučera et al., 1977) equipped with compensation series of four pairs of thermocouples (Čermák & Kučera, 1981) were installed at breast height North-facing side of a Douglas fir stem. Sap flow depth estimated with multi-point sensors (Čermák & Nadezhdina, 1998) of the heat field deformation method (Nadezhdina et al., 1998, 2006, 2010; Čermák et al., 2004) reached maximum about 60 mm, showing rather asymmetric flow density, therefore thermocouples were placed into the depth of 15 mm beneath cambium in order to better represent the flow. The heat loss (λ) considered in the applied equations was eliminated partially by the technical design of the measuring point (insulation by polyurethane foam and shielded against direct sun radiation and eventually application of 111
compensating pairs of thermocouples at a certain (about 10 cm) horizontal distance from the heated stem section. Several variants of sensor horizontal geometry were tested, of which one as shown on (Fig. 1) gave the best results. In this horizontal geometry of the trunk heat balance (THB) method the sensor is symmetrical in respect of central vertical axis and temperature references are at the level of electrodes center. Average values of left and right thermocouple positions situated above (dTa) and below (dTb) electrodes were applied for sap flow calculation. The axial symmetrical difference dTa-dTb characterizes the direction of the flow and dT = dTa+dTb is used to calculate sap
flow.
Results and discussion Numerical simulations The influence of the sap flow on temperature field around heated volume for THB method with internal (direct electric) heating (Fig. 2), for heat conductivity of wood in the traverse direction λx =0.3 W.m-1.K-1 and for sap flow density Qw=0.02 kg.m-2.s-1, side (radial) view (left panel) and frontal view (right panel). Situation for acropetal sap flow is shown and for basipetal flow the temperature field would be vice verza. Frontal view shows part of the stem circumference and corresponds to the width of measuring point and maximum temperature difference of heated volume is about 3.5 K. Side view characterizes the radial pattern of temperature field reaching of the edge of the sap low (max. 60 mm). The dependence of fictitious flow Qw.fic on the actual value of sap flow density Qw is noticeable but the dependence of fictitious flow Qw.fic on actual value of heat
conductivity in transverse direction λx is minor (Fig. 3). Influence of sap flow density Qw on fictitious flow Qw.fic is shown for different heat conductivities λx (upper panel)
and influence of heat conductivity λx on fictitious flow Qw.fic is shown for different sap flow densities Qw (lower panel).A certain point at the upper panel of the above figure correspond to a certain line in the lower panel of the figure. There is almost linear dependence (bit a more curvilinear at the upper panel).
112
Experimental Classical ways of sap flow measurement are based on vertical sensor arrangement with temperature reference of thermocouples situated about 10 cm below electrodes (Čermák, 1973, 1982, 2004) with eventual applications of compesating pairs of thermocouples (Čermák & Kučera, 1981). Sap flow is expressed e.g. per one meter of stem xylem circumference. Values of Q shown for Classic 1 (Fig. 4 - upper panel) have been calculated by subtracting of Qfic taken as Qmin each separate day (this causes underestimating error because saturating flows are eliminated). Values of Q shown in Classic 2 (Fig. 4 - lower panel) have been calculated by subtracting Qmin taken in predown sap flow records in relatively rare days with rain (supposing full water saturation of tissues) more days should be taken into account of course here (this causes mostly overestimating errors). These errors are important especially in the range of low flows. Figures applied to derive temperature differences corresponding to zero sap flow, for introducing into the equation (16) in order to calculate the value of fictitious flow under such conditions (Fig. 5) modified from the heat field deformation (HFD) method (Nadezhdina et al., 1998, 2006, 2010). In cases, when (dT)0 cannot be estimated this way, we can approximate it by taking values of maximum dT over the period from sunset to sunrise. This can happen when we cannot put a line through the data, which would well fit to them and would have an increasing trend toward y axis. New variants of sap flow measurement are based on horizontal sensor arrangement with temperature reference of thermocouples situated symmetrically on both sides of electrodes (Fig. 1). Sap flow is calculated by formula based on more accurate detemination of (Qfic)0 calculated from the temperature difference (dT)0 under zero sap flow from the relationship relating x = dTa-dTb to y = dT similarly like it is applied at the HFD method – see (Nadezhdina et al., 1998, 2006, 2010). This represents Variant 1 (Fig. 6 - upper panel) (Eq. 19). Variant 2 (Fig. 6 - lower panel) applies the same sensor arrangement and also calculates (Qfic)0. However it includes also the impact of Qrec on Qfic (Eq. 20). We can see from the relationship of sap flow rate calculated according to Variant 2 with another Variant 1 and Classical 1, 2 (Fig. 7) that the Classical 1 method underestimates a bit lower sap flow rates, while the Classical 2 method partially overestimate the lower sap flow rate. This is because all estimated daily minima are not 113
absolutely exact and therefore the systems cannot recognize true zero and also reverse flows (if any). Both classical methods underestimates very high flow rates. Both new variants can recognize zero and reverse flows thanks to symmetrical arrangements of thermocouples above and below the electrodes. Variant 2 is the best, because as the single one considers the impact of actual sap flow on the fictitive flow, which makes the results of flow measurements more accurate. That is why it was also taken as a reference (i.e., x axis). If considering Variant 2 as the closed aproximation of true sap flow it means that the Variant 1 makes a linearization and Variant 2 corrected mostly the amplitude which is visible on the higher values (see Fig. 6, 7).
Conclusions The paper presented the new approach focused on improving of sap flow measurement with the THB method based on volume heating of trunk segment of trees. The modified variant of the method based on a new arrangement of temperature probes was developed based on theoretical analysis, simulations and their application on real measured data. This involves the dependence of the so-called “fictitious flow” (Qfic) on the recorded flow data (Qrec). Comparing the results of the modified method with the results of the classical THB method indicates possibilities of measurements of zero and reverse sap flow.
Acknowledgements The study was supported by the project Research Intention of Czech MSM 6215648902 and by the Internal Grant Agency 2011 (IGA 2011) of the Faculty of Forestry and Wood Technology Mendel University in Brno.
114
References Burgess SO., Adams MA, Turner NC, Beverly CR, Ong CK, Ghan AAH, Bleby TM 2001. An improved heat pulse method to measure low and reverse rates of sap flow in woody plants. Tree Physiology 21: 589-598
Campbell GS. 1991. An overview of methods for measuring sap flow in plants. In: Collected summaries of papers at the 83rd annual meeting of the American Society of Agronomy, Division A-3. Agroclim and Agron Modeling, Denver, Colorado, 2–3. Čermák J, Deml M, Penka M. 1973. A new method of sap flow rate determination in
trees. Biol. Plant. (Praha) 15(3): 171-178. Čermák J, Deml M. 1974. Method of water transport measurements in woody species,
especially in adult trees (in Czech). Patent (Certification of authorship) CSFR, No.155622 (P.V.5997-1972) Čermák J, Kučera J, Penka M. 1976. Improvement of the method of sap flow rate
determination in adult trees based on heat balance with direct electric heating of xylem. Biol Plant (Praha) 18: 105–110. Čermák J, Kučera J. 1981. The compensation of natural temperature gradient in the
measuring point during the sap flow rate determination in trees. Biol. Plant. (Praha)
23(6): 469-471. Čermák J, Úlehla J, Kučera J, Penka M. 1982. Sap flow rate and transpiration
dynamics in the full-grown oak (Quercus robur L.) in floodplain forest exposed to seasonal floods as related to potential evapotranspiration and tree dimensions. Biol. Plant. (Praha) 24(6): 446-460. Čermák J, Kučera J. 1990a. Scaling up transpiration data between trees, stands and
watersheds. Silva Carelica 15: 101-120. Čermák J, Kučera J. 1990b. "Transpiration, its limiting factors and hydro logically
important biometric parameters of forest tree species" (in Czech). In: Proc."
Československý príspevek do MHP-UNESCO 1985-89", 10p, Modra u Bratislavy, CSFR, Sept.11-13.1989. Čermák J, Kučera J. 1991. Extremely fast changes of xylem water flow rate in mature
trees, caused by atmospheric, soil and mechanical factors. 181-190pp. In: Proc.CEC International Workshop "Methodologies to assess the impacts of climatic changes on vegetation: Analysis of water transport in plants and cavitation of xylem transport in 115
plants and cavitation of xylem conduits". Raschi A. Borghetti M. (eds.), May 2931.1991. Firenze, Italy. Čermák J, Nadezhdina N. 1998. Sapwood as the scaling parameter - defining
according to xylem water content or radial pattern of sap flow? Ann.Sci.For.55: 509521. Čermák J, Kučera J, Nadezhdina N. 2004. Sap flow measurements with two
thermodynamic methods, flow integration within trees and scaling up from sample trees to entire forest stands. Trees, Structure and Function 18: 529-546.
Cienciala E, Lindroth A, Čermák J, Hallgren J-E, Kučera J. 1992. Assessment of transpiration estimates for Picea abies trees during a growing season. Trees 6: 121127.
Daum CR. 1967. A method for determining water transport in trees. Ecology 48: 425– 431.
Kostner B, Granier A, Čermák J. 1998. Sap flow measurements in forest standsmethods and uncertainties. Ann.Sci.For. 55: 13-27.
Kučera J. 1977. A system for water flux measurements in plants (in Czech). Patent (Certificate of authorship) CSFR No. 185039 (P.V. 2651-1976). Bureau for Inventions and Discoveries, Prague.
Kučera J, Čermák J, Penka M. 1977. Improved thermal method of continual recording the transpiration flow rate dynamics. Biol. Plant. (Praha) 19(6): 413-420.
Lundblad M, Lagergren F, Lindroth A. 2001. Evaluation of heat balance and heat dissipation methods for sap flow measurements in pine and spruce. Ann.Sci.For. 58: 625–638.
Nadezhdina N, Čermák J, Nadezhdin V. 1998a. Heat field deformation method for sap flow measurements. Proc. 4th. International Workshop on Measuring Sap Flow in Intact Plants. Židlochovice, Czech Republic, Oct.3-5, 1998. 72-92. IUFRO Publications. Publishing house of Mendel Univ.Brno.
Nadezhdina N, Čermák J. 1998b. "The technique and instrumentation for estimation the sap flow rate in plants". Patent No.286438 (PV-1587-98).
Nadezhdina N, Čermák J, Neruda J, Prax A, Ulrich R, Nadezhdin V, Gašpárek J, Pokorný E. 2006. Roots under the load of heavy machinery in spruce trees. European J.For.Res. 125: 111-128.
116
Nadezhdina N, David TS, David JS, Ferreira MI, Dohnal M, Tesar M, Gartner K, Leitgeb E, Nadezhdin V, Čermák J, Jimenez MS, Morales M. 2010. Trees never rest: the multiple faces of hydraulic redistribution. Ecohydrology 3: 431-444.
Offenthaler I, Hietz P. 1998. A comparison of different methods to measure sap flow in spruce. In: Čermák J, Nadezhdina N (eds) Measuring sap flow in intact plants. Proceedings of 4th International Workshop, Židlochovice, Czech Republic, IUFRO Publ. Mendel University, Brno, Czech Republic, 55–64.
Penka M, Čermák J, Štepánek V, Palát M. 1979. Diurnal courses of transpiration rate and transpiration flow rate as determined by the gravimetric and thermometric methods in a full-grown oak tree (Quercus robur L). Acta Univ. Agric. Brno, Ser C,
48 (1-4): 3-30. Sakuratani T, Aoe T, Higuchi H. 1999. Reverse flow in roots of Sesbania rostrata measured using the constant power heat method. Plant Cell Environ 22: 1153–1160.
Sakuratani T. 1981. A heat balance method for measuring water flux in the stem of intact plants. J Agric Meteorol (Jpn) 37: 9–17.
Schubert H. 1999. Qualitative and quantitative Untersuchung verschiedener Methoden der Xylemflussmessung an Baumen. Diploma Thesis, Forstwissentschaftlichen Fakultat der Ludwig- Maxmilians-Universitats Munchen, Lehrstuhl fur Forstbotanik, Freising
Schulze ED, Čermák J, Matyssek R, Penka M, Zimermann R, Vašícek F, Gries W, Kučera J. 1985. Canopy transpiration and flow rate fluxes in the xylem of the trunk of Larix and Picea trees - a comparison of xylem flow, porometer and cuvette measurements. Oecologia (Berlin) 66: 475-483.
Smith DM, Allen SJ. 1996. Measurement of sap flow in plant stems. J Exp Bot 47: 1833–1844.
Swanson RH. 1994. Significant historical development in thermal methods for measuring sap flow in trees. Agric For Meteorol 72: 113–132.
Tatarinov FA, Kučera J, Cienciala E. 2005. The analysis of physical background of tree sap flow measurements based on thermal methods. Measurements
Science and
Technology 16: 1157-1169.
Trcala M, Čermák J. 2011. Sap flow and heat conductivity measurement and their interactions 8th International Workshop on Sap Flow Voltera, Italy, 8-12 May 2011. 117
Vašíček F. (ed.) 1984. Ecophysiological and ecomorphological studies of individual trees in the spruce ecosystems of the Drahanská vrchovina uplands (Czechoslovakia). Folia Universitatis Agriculturae, facultas Silviculturae, Brno, Czechoslovakia.
Wullschleger SD, Meinzer FC, Vertessy RA. 1998. A review of whole-plant water use studies in trees. Tree Physiology 18: 499–512.
118
Fig. 1: Horizontal geometry of the trunk heat balance (THB) sensor. Electrodes are marked by thick vertical lines (dotted for 5 electrode variant and full for 3 electrode ones), position of thermocouples is marked by large dots. Distances between electrodes and thermocouples are given in mm. The sensor is symmetrical in respect of central vertical axis (the reference is shown here on the left for 5 electrodes, on the right for 3 electrode system). Average values of left and right thermocouple positions situated above (dTa) and below (dTb) electrodes were applied for sap flow calculation.
119
Fig. 2: Example of simulated temperature field [K] for THB method with internal (direct electric) heating, for heat conductivity of wood in the traverse direction λx =0.3 W.m-1.K1 and for sap flux density Qw=0.02 kg.m-2.s-1, side (radial) view (left panel) and frontal view (right panel). Situation for acropetal sap flow is shown.
120
Fig. 3: Influence of sap flow density Qw on fictitious flow Qw.fic for different heat conductivities λx (upper panel). Influence of heat conductivity λx on fictitious flow Qw.fic for different sap flux densities Qw (lower panel).
121
Fig. 4: Sap flow dynamics measured by the classic ways expressed per one meter of stem xylem circumference. Values of Q shown in Classic 1 (upper panel) have been calculated by subtracting of Qfic taken as Qmin each sparate day. Values of Q shown in Classic 2 (upper panel) have been calculated by subtracting Qmin taken in pre-down sap flow records in (relatively rare) days with rain (supposing full water saturation of tissues).
122
Fig. 5: Estimation of (dT)0 (dT for zero sap flow) – applied values are over the period from sunset to sunrise of each separate day. Analogue of the approach applied in the heat field deformation method.
123
Fig. 6: Sap flow dynamics measured by the variants 1, 2 of the new sensor arrangement expressed per one meter of stem xylem circumference. Both variant consider (Qfic)0 taken from dTa-dTb to dT relationships (Fig. 5). However variant 2 consider in addition the impact of sap flow rate on Qfic.
124
Fig. 7: Comparison the results obtained by the new variant 1 and the classic THB method (y axis) with the new (optimal) variant 2 (x axis).
125
