NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP8 Metoda sítí s aplikací na parabolické a eliptické rovnice Metodicky de facto totéž, co v minulé přednášce jen jiné příklady: Neustálený rychlostní profil při oscilačním toku v trubce (parabolická rovnice pro cylindrický s.s.) Teplotní pole (parabolická rovnice) a elektrické pole (eliptická rovnice) při přímém ohmickém ohřevu. Identifikace parametrů modelu popisovaného metodou konečných diferencí použitím algoritmu Nelder Mead (fminsearch) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. Oscilační tok v trubce NAP8 Při řešení neustáleného toku stlačitelné tekutiny v trubce nebyl ověřován radiální rychlostní profil. Tiše jsme předpokládali, že rychlostní profil odpovídá ustálenému toku (a že se třeba v laminárním režimu vyvine parabolický profil). Při rychle se měnícím axiálním gradientu tlaku to tak není. Mírou odchylky je frekvence tlakových pulzací a Womersleyho číslo Mordechai Ardon

3. Oscilační tok v trubce NAP8 Prototypem oscilační toku je například tok krve v cévním systému. Pulzace tlaku vyvolává srdce (systola-diastola, opakující se s frekvencí cca 1Hz) a v cévách s větším průměrem (aorta) dochází k tomu, že krev u stěny teče opačným směrem než někde uprostřed. Je tomu tak proto, že největší rychlost a tedy i kinetická energie je největší u osy trubky, takže setrvačnost kapaliny překoná i opačně působící gradient tlaku při jeho náhlé změně. Naproti tomu, u stěny, kde je rychlost menší, nemusí gradient tlaku překonávat tak velkou hybnost kapaliny, a směr toku bude víceméně poslušně sledovat směr gradientu tlaku. Když zanedbáme stlačitelnost kapaliny a tuhou trubku, popisuje celý problém Navierova Stokesova rovnice pro složku rychlosti kapaliny ve směru osy trubky, zapsaná v cylindrickém souřadném systému V nekonečně dlouhé trubici musí být příčná složka rychlosti nulová a z NS rovnice v radiálním směru plyne, že gradient tlaku dp/dr=0, jinými slovy, tlak p nezávisí na poloměru. Axiální gradient tlaku pak může být pouze funkcí času t.

4. Oscilační tok v trubce NAP8 Pokud je gradient tlaku konstantní, vyvine se po jisté době parabolický rychlostní profil Pokud je gradient tlaku harmonicky oscilující funkcí času, budou i rychlosti proudění harmonickou funkcí času, se stejnou frekvencí , jen fázově posunuté. Radiální rychlostní profil je pro tento speciální případ vyjádřen Besselovou funkcí J0 s komplexním argumentem kde r je bezrozměrný poloměr,  je úhlová frekvence sinusových oscilací gradientu tlaku a Wo je Womersleyho číslo Když je Wo<<1 můžeme předpokládat, že rychlostní profil je parabolický a docela přesně sleduje časové (pomalé) změny gradientu tlaku. Pak ani není nutné použít numerické řešení a třeba i změny průtoku lze počítat z běžné formy Bernoulliho rovnice. Toto analytické řešení uvádím spíš jen pro zajímavost (v praxi se využívá při modelování toku krve v artériích), v dalším textu se zaměříme na numerické řešení, které je použitelné pro zcela libovolný časový průběh axiálního gradientu tlaku.

5. r <0,1> r 1 2 i-1 i i+1 N Oscilační tok v trubce NAP8 Navier Stokesovu rovnici upravme zavedením bezrozměrného poloměru a času, Diferenční přepis této rovnice na ekvidistantní síti uzlových bodů (pro jeden časový krok ) A to je soustava lineárních algebraických rovnic s tridiagonální maticí soustavy (vektor bi je diagonála, di je vektor pravé strany) Na ose symetrie Na stěně Zkuste si to odvodit: na ose je problém s dělením poloměrem r1 – musíte použít limitní přechod (třeba l’Hopitalovo pravidlo)

6. u>0 u<0 Oscilační tok v trubce NAP8 Tato tridiagonální soustava algebraických rovnic se dá velmi efektivně řešit faktorizací, algoritmizovanou např. v MATLABovské funkci tridag q()-průtok function u = tridag(a,b,c,r,n) bet=b(1); u(1)=r(1)/bet; for j=2:n gam(j)=c(j-1)/bet; bet=b(j)-a(j)*gam(j); u(j)=(r(j)-a(j)*u(j-1))/bet; end for j=n-1:-1:1 u(j)=u(j)-gam(j+1)*u(j+1); end Úplný výpis MATLABovského programu bude uveden až za okamžik. Zde je jen ukázka typického časového průběhu průtoku, který jak patrno není ve fázi se zadaným harmonickým průběhem gradientu tlaku. Druhý obrázek (radiální rychlostní profil) dokumentuje existenci zpětného toku v ose trubky – pro harmonicky oscilující gradient tlaku by to měla být Besselova funkce J0.

7. Bernoulliho rovnice NAP8 Z praktického hlediska má asi největší význam posouzení toho, do jaké míry ovlivní změna rychlostního profilu Bernoulliho rovnici. Otázka: Mohu použít následující tvar Bernoulliho rovnice pro výpočet časově proměnné střední rychlosti (nebo průtoku)? Rovnici si odvoďte z rovnováhy sil působících na váleček kapaliny o poloměru R a tloušťce dx. Jistá nesrovnalost takto napsané Bernoulliho rovnice je v tom, že zrychlení se odvozuje z předpokladu konstantního radiálního profilu rychlosti, zatímco třecí ztráty z představy parabolického profilu rychlosti. Tatáž Bernoulliho rovnice, ale zapsaná s bezrozměrným časem a gradientem tlaku ověřte, že to je analytické řešení (odvodí se metodou variace konstant z řešení homogenní rovnice) Odchylka predikovaných rychlostí nebo průtoků závisí na tom, jak velké jsou oscilace tlaku; srovnání s přesným (numerickým) řešením ukazuje, že amplituda průtoků je dle Bernoulliho rovnice, založené na tlakových ztrátách odpovídajících parabolickému rychlostního profilu, zhruba o 10% vyšší (parabolický profil má menší tlakovou ztrátu). Lepší přesnosti se dá dosáhnout aproximací rychlostního profilu polynomem 4 stupně, viz následující text...

8. Přibližné řešení NS rovnice NAP8 Vyjděme z Navier Stokesovy rovnice ve tvaru použitém pro numerické řešení (s bezrozměrným časem a radiální souřadnicí) a hledejme přibližné řešení ve tvaru, který zajistí automatické splnění okrajových podmínek v ose a na stěně (1) Tento rychlostní profil by měl v každém časovém kroku splnit počáteční podmínky (pro čas =0), vyjádřené např. dvojicí hodnot c1(0), c2(0), a také Navier Stokesovu rovnici a to pro libovolné hodnoty r a pro libovolnou zadanou funkci P()

9. Přibližné řešení NS rovnice NAP8 Porovnáním koeficientů u mocnin poloměru r (r0, r2, r4) získáme soustavu obyčejných diferenciálních rovnic Teď ovšem máme problém, protože jsme získali soustavu tří diferenciálních rovnic pro pouze dvě neznámé c1,c2. A ty tři rovnice jsou lineárně nezávislé, takže přesné řešení (c1(),c2()) asi existovat nebude, můžeme jen hledat přijatelný kompromis. • Lze zkusit například tyto varianty: • c2=0 (ignorujeme druhou a třetí rovnici), řešení • ignorujeme třetí rovnici a řešíme soustavu 2 rovnic metodou vlastních čísel • eliminujeme c2 z prvních dvou rovnic předpokládajíce přibližnou platnost třetí rovnice (přibližně konstantní c2 v řešeném časovém kroku) • V následujícím textu rozebereme variantu c), protože je nejjednodušší: Eliminace c2() z první rovnice a jejím řešením je pro konstantní gradient tlaku P Funkce c2() plyne z algebraické rovnice vzniklé sečtením předchozích 3 rovnic. Limitní případ pro  má řešení c1()=-P/4 a c2()=0 což správně odpovídá plně vyvinutému parabolickému rychlostnímu profilu.

10. Vážená residua–přibližné řešení NAP8 Na těchto šedivých stránkách je uvedeno alternativní řešení vypočtu koeficientů c1(t), c2(t) metodou vážených reziduí. Řešení metodou vážených reziduí je možná lepší (konzistentnější) než předchozí heuristické řešení.

11. Vážená residua–přibližné řešení NAP8 Dosazením získáme residuum NS rovnice Funkce c1 c2 musí splňovat podmínku nulového váženého rezidua alespoň pro váhové funkce w1(r)=1 a w2(r)=r, tj. rovnice V maticovém zápisu A to už je standardní formulace počátečního problému pro 2 rovnice

12. Vážená residua–přibližné řešení NAP8 Z praktického hlediska je možná lepší řešit přímo průtoky Po dosazení původních rozměrových veličin a doplnění zrychlení, které působí proti směru osy x, získáme soustavu dvou obyčejných diferenciálních rovnic pro průtoky Pro konstantní gradient tlaku lze nalézt řešení evoluce vývoje průtoků v analytickém tvaru (řešením vlastních vektorů a vlastních čísel) jako lineární kombinaci exponenciálních funkcí.

13. Vážená residua–přibližné řešení NAP8 Speciální případ kdy gradient tlaku je konstantní a řešení je ustálené Je zřejmé, že řešení se redukuje na plně vyvinutý parabolický profil, stěnová složka mizí, a řešení splňuje Hagen Poisseuille zákon pro laminární tok v trubce s kruhovým průřezem.

14. Bernoulliho rovnice NAP8 Z rychlostního profilu získáme integrací průtok nebo střední objemovou rychlost Tento výsledek je možné srovnat s analogickým řešením Bernoulliho rovnice kde exponent m=8 odpovídá konstantnímu, zatímco m=4 parabolickému radiálnímu rychlostnímu profilu v nestacionárním členu (du/dt). Porovnání výsledků přesného numerického řešení, Bernoulliho rovnice a aproximace radiálního rychlostního profilu polynomem 4-tého stupně je naprogramováno jako MATLABovský program…

15. rychlostní profily pro stejnou střední objemovou rychlost=1 m/s bezrozměrný poloměr Bernoulliho rovnice NAP8 Na počátku každého časového kroku (v čase =0) je obecně třeba zadat kompletní radiální rychlostní profil jako počáteční podmínku. Za předpokladu, že rychlostní profil je parabolický, stačí zadat jen průtok (nebo střední objemovou rychlost), v případě s rychlostním profilem 4 stupně musí být podmínky dvě, např. c1(0) a c2(0), nebo střední objemová rychlost a relativní velikost stěnového průtoku uw, z nichž lze koeficienty c1,c2 snadno spočítat

16. Oscilační tok v trubce NAP8 Zadaný průběh gradientu tlaku , počítají se časové průběhy průtoků qfd (Finite Dif.), qb (Bernoulli), qode (Ordinary Dif.Eq.), včetně matic rychlostních profilů for i=1:nt %finite differences for ir=1:nr-1 d(ir)=-p(i)+u(ir)/dt; end u(1:nr)=tridag(a,b,c,d,nr); q=0; for ir=1:nr q=q+u(ir)*r(ir); end q=2*3.141*dr*q; qfd(i)=q; ufd(i,1:nr)=u(1:nr); %bernoulli umean0=(umean0+p(i)/8)*exp(-8*dt)-p(i)/8; qb(i)=umean0*3.141; for ir=1:nr ub(i,ir)=2*umean0*(1-r(ir)^2); end %ODE c1(i)=(c10+p(i)/4)*exp(-16*dt/3)-p(i)/4; c2(i)=(-p(i)-4*c1(i))/12; qode(i)=3.141/2*(c1(i)+c2(i)/3); c10=c1(i); for ir=1:nr uode(i,ir)=c1(i)*(1-r(ir)^2)+c2(i)*(r(ir)^2-r(ir)^4); end end nt=1501; nr=21; a=zeros(nr);b=zeros(nr);c=zeros(nr);d=zeros(nr);u=zeros(nr); ufd=zeros(nt,nr);uode=zeros(nt,nr);ub=zeros(nt,nr); qfd=zeros(nt);qb=zeros(nt);qode(nt)=0;t=zeros(nt); tmax=3;dt=tmax/(nt-1);t=linspace(0,tmax,nt); % defince casoveho prubehu gradientu tlaku p(1:nt)=-1-3*sin(10*t); r=linspace(0,1,nr);dr=1/(nr-1);dr2=dr^2; %pocatecni rychlostni profil stabilizovany for i=1:nr u(i)=-p(1)/4*(1-r(i)^2); end umean0=-p(1)/8; c10=-p(1)/4; for i=2:nr-1 a(i)=-(r(i)+r(i-1))/(2*r(i)*dr2); c(i)=-(r(i+1)+r(i))/(2*r(i)*dr2); b(i)=1/dt-a(i)-c(i); end b(1)=1/dt+4/dr2;c(1)=-4/dr2;d(1)=-p(1)+u(1)/dt; a(nr)=0;b(nr)=1;d(nr)=0;u(nr)=0;

17. Oscilační tok v trubce NAP8 Výsledky jsou znázorněny v grafech figure(1) plot(t(1:nt),qfd(1:nt),'r',t(1:nt),qb(1:nt),'b',t(1:nt),qode(1:nt),'g') figure(2) hold off for it=nt-100:50:nt plot(r(1:nr),uode(it,1:nr),'g',r(1:nr),ufd(it,1:nr),'r',r(1:nr),ub(it,1:nr),'b') hold on end q-diference q-Bernoulli q-aproximace

18. Přímý Ohmický ohřev NAP8 Vnitřní (objemové) zdroje tepla vyvolané průchodem elektrického proudu nebo absorpcí mikrovlnného záření Ardon

19. Napájení 50 Hz, konstantní napětí Rovinné nerezové elektrody plátek masa Teplotní sondy REFLEX R232 Nortech TS-RX4 MULTIPLEXER AGILENT 20CH 34901A TECTRA - LMG 95 Měření výkonu GPIB Ohmický ohřev masa NAP8 Na obrázku je schema aparatury pro experimentální přímý ohmický ohřev plátku masa (sekaná), který je sevřen mezi nerezovými elektrodami. Napájení 24V, frekvence 50Hz. Měří se elektrický výkon, teploty masa (optická vlákna) a teploty elektrod (termočlánky). Cílem je stanovení elektrické vodivosti masa a její závislost na teplotě. Kdyby byl ohřev naprosto rovnoměrný (teploty v každém místě stejné) stačil by jednoduchý integrální model a vyhodnocení elektrické vodivosti přímo z výkonu a teploty (viz dále). Nerovnoměrný ohřev je ale nutné modelovat řešením FK rovnice s proměnným objemovým zdrojem tepla. Je to podobný případ jako sušení zrna (podobné rovnice a metody řešení). Při elektrickém ohřevu ale musíme kromě parabolické rovnice vedení tepla (FK rovnice s akumulačním a zdrojovým členem) současně řešit eliptickou rovnici rozložení elektrického potenciálu.

20. E Ohmický ohřev masa integrální model NAP8 Pokud by byla teplota masa všude stejná a intenzita elektrického pole E konstantní, bylo by možné napsat jednoduchou integrální bilanci Zanedbání tepelné kapacity elektrod měrná elektrická vodivost rostoucí lineárně s teplotou S-povrch, V-objem vzorku jejímž řešením je exponenciální nárůst teploty Když je R>0 bude rychlost ohřevu klesat, pro R<0 bude exponenciálně růst (akcelerace ohřevu je ovlivněna vysokou teplotní citlivostí elektrické vodivosti , tj. vysokou hodnotou koeficientu 1). Následující numerické řešení je vlastně jen zpřesnění tohoto integrálního modelu a jeho analytického řešení.

21. Ohmický ohřevmasa numerický model NAP8 Schiele

22. he kontaktní vrstva W × W T0 z m-maso W hk 2H Metoda sítí – ohmický ohřev NAP8 Geometrie: ohřívaný vzorek masa má tloušťku 2H a příčné rozměry W x W Experimenty ukazují, že přímý ohmický ohřev je ovlivněn i elektrickým a termickým odporem kontaktní vrstvy, a že elektrické vlastnosti masa i této kontaktní vrstvy závisí nejenom na teplotě, ale i na přítlačné síle elektrod. Pro jednoduchost budeme vliv kontaktních odporů i efekt působení tlaku zanedbávat. • Předpoklady: • Kontaktní odpor je zanedbatelný, takže teplota masa je stejná jako teplota elektrody v místě kontaktu. Stejný je i elektrický potenciál (zanedbáme elektrickou dvouvrstvu, polarizační efekty, elektrolýzu). • Maso je homogenní, a termofyzikální vlastnosti jsou konstantní. Teplotně závislá je pouze měrná elektrická vodivost. • Povrch masa je chlazen přirozenou konvekcí vzduchem (konstantní teplota vzduchu Ta) a kontaktem s elektrodou (elektroda je rovněž chlazena přirozenou konvekci vzduchem). Elektrická vodivost elektrod (ocel) je tak velká, že potenciál U je v celé elektrodě konstantní a žádný elektrický výkon se v ní negeneruje.

23. Metoda sítí – FK rovnice NAP8 Teplotní pole v mase popisuje FK rovnice se zdrojem tepla, který je úměrný kvadrátu intenzity elektrického pole (což je gradient potenciálu U).  je měrná elektrická vodivost masa, která závisí lineárně na teplotě Tuto rovnici je třeba doplnit počátečními podmínkami (konstantní teplota T0 v čase 0) a okrajovými podmínkami (přestup tepla  na povrchu vzorku a teplota elektrod na ploškách z=0 a z=2H). Rozložení teplot na povrchu elektrod je třeba stanovit současným řešením FK rovnice v elektrodách, ale tentokrát bez objemového zdroje tepla.

24. x z W y Metoda sítí – FK rovnice NAP8 3D řešení s 3D sítí uzlových bodů je časově náročné (cílem je identifikace parametrů 0 a 1 a optimalizační algoritmus bude tudíž vyžadovat opakované řešení této diferenciální rovnice). Je tudíž žádoucí formulaci zjednodušit převedením na 1D problém. Zbavit se derivací x,y v příčném směru lze integrací FK rovnice přes celý průřez (teplota T(t,x,y,z) je po integraci nahrazena střední teplotou v místě z T(t,z)). Zanedbáme teplotní závislost tepelné vodivosti  a složky intenzity elektrického pole v příčném směru (na povrchu, který je elektricky izolovaný, i v ose musí být U/x i U/y nulové). Tento člen vznikne integrací druhých derivací ve směrech x a y a použitím okrajových podmínek typu Poznámka: Techniku redukce dimenze problému integrací 3D rovnic přes méně důležité rozměry jste použili v předmětu Tepelné procesy nebo Přenos hybnosti, tepla a hmoty např. při stanovení teplotního profilu podél chladicího žebra nebo při převodu Prandtlových rovnic mezní vrstvy na Karmánovy integrální rovnice. W

25. Metoda sítí – FK rovnice NAP8 Objemový zdroj tepla je vyjádřen intenzitou elektrického pole (dominantní složkou U/z), která se podél souřadnice z mění (mění-li se teplota a tedy i elektrická vodivost). Alternativně lze tento člen vyjádřit pomocí hustoty elektrického proudu I [A/m2], která se dle souřadnice z nemění, závisí jen a pouze na čase V takto formulované FK rovnici se vyskytuje pouze teplota (a ne elektrický potenciál). Nelinearitu zdrojového členu lze odstranit linearizací (ověřte Taylorovým rozvojem vzhledem k teplotě T0): Hustotu elektrického proudu I (i když je konstantní) je třeba stanovit v každém časovém kroku řešením Laplaceovy rovnice elektrického potenciálu

26. Metoda sítí – program NAP8 Rámcové schéma numerického řešení (vývojový diagram) Hlavní program Podprogram MODEL Počáteční podmínky (teploty a potenciál ve všech uzlech) Čtení geometrie a měřených dat (časy, výkony, teploty) Cyklus časových kroků Použití algoritmu Nelder Mead pro minimalizaci součtu čtverců odchylek numerické predikce (podprogram MODEL) a měřených teplot a výkonů. Optimalizované parametry jsou 0 1 Cyklus iterací Sestavení a řešení soustavy rovnic pro teplotní pole Sestavení a řešení soustavy rovnic pro elektrický potenciál Výpočet hustoty proudu I Konec iterací Konec cyklu časových kroků

27. H h 1 2 3 W P E N-1 N Metoda sítí – diskretizaceFK NAP8 Jednoduchá geometrie je dobrý předpoklad pro použití metody sítí nebo metody kontrolních objemů. Vzhledem k symetrii stačí řešit teplotní a potenciálové pole v intervalu z:<0,H>, který pokryjeme sítí N ekvidistantních uzlů. Každému uzlu přiřadíme teplotu TP a potenciál UP. V metodě sítí požadujeme, aby v každém uzlu byla splněna příslušná diferenciální rovnice (počet podmínek je potom stejný jako počet neznámých). První a druhé derivace v diferenciálních rovnicích se nahradí konečnými diferencemi (numerickou aproximací derivací v bodě). Pro uzly P=2,3,…,N tudíž píšeme algebraické rovnice pro triádu sousedních teplot TW TP TE na nové hladině času

28. Metoda sítí – diskretizaceFK NAP8 Sdružením členů u neznámých TW TP TE vyjádříme koeficienty soustavy rovnic s tridiagonální maticí bP aP cP dP Tuto rovnici píšeme i pro uzel na ose, tj. pro P=N, s tím, že TE=TW=TN-1, čímž aplikujeme i okrajovou podmínku symetrie v ose (T/z=0). Okrajovou podmínku na elektrodě (P=1) je ovšem nutné vyjádřit jinak, ne metodou sítí (konečných diferencí), ale kontrolních objemů. Metodu konečných diferencí lze použít v takových uzlech, kde je víceméně hladký průběh řešení, a kde nejsou výrazné nespojitosti koeficientů řešené diferenciální rovnice. Pokud to tak není, je vhodnější použít metodu kontrolních objemů, která vyšetřovanou oblast rozdělí na menší podoblasti (kontrolní objemy) z nichž každý obklopuje „svůj“ uzlový bod. Místo požadavku na splnění diferenciální rovnice v izolovaném uzlu se požaduje splnění diferenciální rovnice v „průměru“, ve smyslu integrální bilance kontrolního objemu.

29. 1 2 3 he Metoda sítí – diskretizace FK NAP8 Kontrolní objem, který přiřadíme uzlu číslo 1 zahrnuje nalevo elektrodu, v pravé části maso. Z pravé části vtéká do kontrolní objemu entalpický tok T/z, z levé části odtéká tok povrchem elektrody. Jistá komplikace je v tom, že povrch elektrody Ae je větší než je průřez ohřívaného masa A=W2 a kontrolní objem má i svou tepelnou kapacitu, v níž převažuje tepelná kapacita elektrody. Velmi přibližnou integrální bilanci tohoto kontrolního objemu můžeme vyjádřit takto tepelná kapacita elektrody vztažená k průřezu masa tepelné ztráty povrchem elektrody Nahrazením derivací diferencemi získáme analogon (integrální bilanci kontrolního objemu) d1 c1 b1

30. Metoda sítí – diskretizace FK NAP8 Okrajová podmínka na elektrodě bude asi největším zdrojem nepřesností numerické predikce, protože je založena na předpokladech, které zcela určitě nejsou splněny: například teplota elektrody není stejná, protože součinitel teplotní vodivosti není nekonečně velký (u nerezové oceli je docela malý cca 15 W/m/K). Asi by tedy nebylo správné uvažovat celkovou tepelnou kapacitu elektrody, ale jen její část (redukcí povrchu elektrody Ae). Nejistoty se týkají i součinitele přenosu tepla z povrchu elektrody do okolního vzduchu – jde o přirozenou konvekci. Nejistoty související s odhadem  se ovšem týkají i chlazení povrchu masa,  na horizontálním a na vertikálním povrchu určitě nebude stejné. V následujícím programu bude použita korelace, kterou navrhl Churchill

31. Metoda sítí – diskretizace FK NAP8 Diskretizace rovnice pro elektrický potenciál i okrajové podmínky jsou mnohem jednodušší WPE kde Na elektrodě předepíšeme přímo polovinu napájecího napětí (střední hodnotu u napětí střídavého) a v ose symetrie nulu

32. Program MATLAB NAP8 Feininger

33. Ohmický ohřev MATLAB NAP8 Funkce fmeat, která popisuje numerickou predikci teplot i výkonů function [vtime,vtemean,vtmmean,vpower] = fmeat ... (h,w,u,he,ae,rkm0,rkm1,rlamm,n,tmax,ntime,ta,tinit) % % h - tloustka vzorku [m] % w - prurez vzorku W x W [m] % u - napeti [V] (rozdil potencialu na elektrodach) % he - tloustka elektrody % ae - pomer plochy elektrody a prurezu vzorku % rkm0,rkm1 - koeficienty el.vodivosti masa kappa=rkm0+rkm1*T % rlamm - tepelna vodivost masa % N -celkovy pocet uzlu (uzel N na ose symetrie) % tmax -max.cas % ntime-pocet casovych kroku % Ta - teplota vzduchu (okoli) ve stupnich Celsia % Tinit(N) - pocatecni rozlozeni teplot % % Vystupy: % vtime(ntime) - casy simulace % vtemean(ntime),... - stredni teplot v elektrode, a v mase % vpower(ntime),... - dissipovane vykony [W] % % numericka simulace % prostorove kroky a fixni hodnoty dz=double(h/(n-1)); dt=double(tmax/(ntime-1)); rhom=980.; cpm=2700.; rlame=15.0; rhoe=7800.0; cpe=460.0; % pocatecni profil teplot a vodivosti for i=1:n tv0(i)=tinit(i); end % pocatecni profil rozlozeni napeti resenim tridiagonalni soustavy av(1:n)=0; bv(1:n)=1; cv(1:n)=0; dv(n)=0; dv(1)=u/2; for i=2:n-1 tm2=(tv0(i-1)+tv0(i))/2; tp2=(tv0(i+1)+tv0(i))/2; dv(i)=0; rm2=(rkm0+rkm1*tm2); rp2=(rkm0+rkm1*tp2); av(i)=-rm2; cv(i)=-rp2; bv(i)=rm2+rp2; end uv0=tridag(av,bv,cv,dv,n); % Prekopirovani pocatecniho rozlozeni teplot a napeti uv=uv0; tv=tv0; % hustota elektrickeho proudu [A/m^2] cur=(uv(3)-uv(2))/dz*(rkm0+rkm1*(tv(3)+tv(2))/2); av-vektor pod diagonálou, bv-vektor diagonálních prvků a cv vektor prvků nad diagonálou matice soustavy řešení soustavy n-rovnic s tridiagonální maticí soustavy

34. Ohmický ohřev MATLAB NAP8 … pokračování funkce fmeat % Nove rozlozeni potencialu av(1:n)=0; bv(1:n)=1; cv(1:n)=0; dv(1)=u/2; dv(n)=0; for i=2:n-1 tm2=(tv(i-1)+tv(i))/2; tp2=(tv(i+1)+tv(i))/2; dv(i)=0; rm2=(rkm0+rkm1*tm2); rp2=(rkm0+rkm1*tp2); av(i)=-rm2; cv(i)=-rp2; bv(i)=rm2+rp2; end uv=tridag(av,bv,cv,dv,n); cur=(uv(2)-uv(3))/dz*(rkm0+rkm1*(tv(3)+tv(2))/2); % konec cyklu iterace na hladine IT end % vyhodnoceni strednich teplot v casovem kroku itime tmm=0; % vyhodnoceni strednich vykonu (kontakt, maso / Watty) for i=1:n tmm=tmm+tv(i); end vtemean(it)=tv(1); vtmmean(it)=tmm/n; vpower(it)=cur*u*w^2; % konec cyklu casovych kroku - nova casova hladina tv0=tv; uv0=uv; end end elegantnější by bylo omezit počet iterací testováním rozdílu vypočtených proudových hustot (while…). Mně se zdá, že 2 iterace úplně stačí… % Cyklus casovych kroku for it=1:ntime vtime(it)=double(dt*(it-1)); % cyklus iteraci na jedne hladine casu for iterat=1:2 % Teplotni profil for i=2:n-1 % maso alpha=(0.02+1.3*(w^3*(tv(i)-ta))^0.25)/w; av(i)=-rlamm/dz^2; cv(i)=-rlamm/dz^2; bv(i)=rhom*cpm/dt+2*rlamm/dz^2+4/w*alpha+cur^2*rkm1/(rkm0+rkm1*tv(i))^2; dv(i)=rhom*cpm*tv0(i)/dt+4/w*alpha*ta+cur^2*(rkm0+2*rkm1*tv(i))/(rkm0+rkm1*tv(i))^2; end % Rovnice v ose alpha=(0.02+1.3*(w^3*(tv(n)-ta))^0.25)/w; av(n)=-2*rlamm/dz^2; cv(n)=0; bv(n)=rhom*cpm/dt+2*rlamm/dz^2+4/w*alpha+cur^2*rkm1/(rkm0+rkm1*tv(n))^2; dv(n)=rhom*cpm*tv0(n)/dt+4/w*alpha*ta+cur^2*(rkm0+2*rkm1*tv(n))/(rkm0+rkm1*tv(n))^2; % Rovnice elektrody - korekce na odpor (zda se nevyznamna) alphae=(0.02+1.3*(w^3*(tv(1)-ta))^0.25)/w; alpha=1/(1/alphae+he/rlame); cv(1)=-rlamm/dz; bv(1)=rhoe*cpe*he*ae/dt+rlamm/dz+alpha*ae; av(1)=0; dv(1)=tv0(1)*rhoe*cpe*he*ae/dt+alpha*ae*ta; % Nova iterace teplot na hladine it tv=tridag(av,bv,cv,dv,n); s tou rovnicí teplot elektrody je možné si pohrát a zkoušet různé varianty…

35. Ohmický ohřev MATLAB NAP8 … pro řešení soustavy rovnic s tridiagonální maticí soustavy byla použita pomocná procedura tridag function u = tridag(a,b,c,r,n) bet=b(1); u(1)=r(1)/bet; for j=2:n gam(j)=c(j-1)/bet; bet=b(j)-a(j)*gam(j); u(j)=(r(j)-a(j)*u(j-1))/bet; end for j=n-1:-1:1 u(j)=u(j)-gam(j+1)*u(j+1); end

36. Ohmický ohřev MATLAB NAP8 Funkce fdev počítající odchylku predikce modelu a experimentálních dat function sumsq =fdev(p,h,w,u,he,ae,rlamm,n,ta,wel,wmeat,wpower,ndata,timdata,teldata,tmdata,powerdata) % vypocet odchylek modelu - teplot i vykonu pro konkretni experiment % sumsq - soucet ctvercu odchylek % vstupy % p(1)=rkm0, p(2)=rkm1 % wel,wmeat,wpower - vahove koeficienty pro soucet ctvercu % ndata pocet casovych kroku v datech % timdata(1:ndata) casy (predpoklad ekvidistantni vzorkovani) % teldata(1:ndata) teploty elektrod % tmdata(1:ndata) teploty masa (nebo sondy uprostred) % powerdata(1:ndata) vykon ve watech celkovy tmax=timdata(ndata)-timdata(1); dtdata=tmax/(ndata-1); nmult=10; dt=dtdata/nmult; ntime=tmax/dt+1; % provedeni simulace pro pocatecni podminky teplot tinit(1)=teldata(1); for i=2:n tinit(i)=tmdata(1); end [vtime,vtemean,vtmmean,vpower] = fmeat(h,w,u,he,ae,p(1),p(2),rlamm,n,tmax,ntime,ta,tinit); sumsq=0; iv=0; for i=1:nmult:ntime iv=iv+1; sumsq=sumsq+wel*(vtemean(i)-teldata(iv))^2; sumsq=sumsq+wmeat*(vtmmean(i)-tmdata(iv))^2; sumsq=sumsq+wpower*(vpower(i)-powerdata(iv))^2; end end numerické řešení bude prováděno s 10x menším časovým krokem než je vzorkování dat. To je otázka přesnosti-třeba vyzkoušet.

37. Ohmický ohřev MATLAB NAP8 Skript čtení dat a volání optimalizační funkce fminsearch aplikované na fdev % Fixni data h=0.01; w=0.02; u=24; he=0.002; ae=2; rlamm=0.71; % parametry site, a casovy krok n=101; tmax=50; ntime=101; % teploty ta=30; t0=30; % experimentalni data t, te,tm,power rdata=load('fmeat.load'); % vahy wel=.1; wmeat=.1; wpower=1; [pnew,sums,exitflag]=fminsearch(@(x)fdev(x,h,w,u,he,ae,rlamm,n,ta,wel,wmeat,wpower,length(rdata),rdata(:,1),rdata(:,2),rdata(:,3),rdata(:,4)),[0.2,0.03]) Vzorek dat souboru fmeat.load (čas, teploty elektrod a masa, příkony…) 0.000e+000 3.029e+001 3.075e+001 1.762e+001 5.000e-001 3.073e+001 3.140e+001 1.785e+001 1.000e+000 3.058e+001 3.198e+001 1.805e+001 1.500e+000 3.078e+001 3.157e+001 1.805e+001 2.000e+000 3.045e+001 3.203e+001 1.825e+001 2.500e+000 3.046e+001 3.331e+001 1.820e+001 3.000e+000 3.118e+001 3.337e+001 1.876e+001 3.500e+000 3.119e+001 3.334e+001 1.871e+001 …. menší váha je přiřazena teplotám – věříme víc predikovaným výkonům použití anonymní funkce, viz. druhá přednáška NAP2 počáteční odhady pnew je vektor identifikovaných parametrů masa 0 1 indikátor úspěšnosti fminsearch (překročení max. počtu kroků…)