NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP7 Metoda sítí a konečné diference (schémata a stabilitní analýza). Hyperbolické rovnice – znovu se vrátíme k hydraulickému rázu Parabolické rovnice - teplotní a koncentrační pole (sušení kávového zrna) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. Ti-1 Ti Ti+1 x PDE konečné diference NAP7 Všechny zmíněné typy PDE (hyperbolické, parabolické i eliptické rovnice) lze řešit metodou konečných diferencí. Vyšetřovaná oblast se pokryje pravidelnou či nepravidelnou sítí uzlů, každému uzlu (tn,xi,yj,zk) se přiřadí počítaná hodnota Tnijk a v každém uzlu se požaduje splnění PDE (nulové reziduum). První a druhé derivace PDE v uzlu (tn,xi,yj,zk) se nahradí aproximacemi, konečnými diferencemi, které jsou počítané z hodnot Tnijk v uzlech sousedních. Výsledkem je jedna algebraická rovnice pro každý uzel sítě. Např. v 1D pro ekvidistantní uzly Asymetrická diference 1 řádu přesnosti Centrální diference 2 řádu přesnosti chyba <K|x2| kde K je nějaká konstanta

3. Ti-1 Ti Ti+1 x PDE konečné diference NAP7 Diferenční formule (i pro neekvidistantní uzly) a řád přesnosti lze odvodit z Taylorova rozvoje funkce T(t,x,y,z) v okolí uzlu Tnijk nebo tak, že sousedními uzly proložíme interpolační polynom a vyčíslíme jeho derivace v uzlu. Např. hodnotami v trojici uzlů proložíme maximálně kvadratický polynom x a jeho první a druhá derivace pro x=0 dává předchozí centrální formule. Neboli: Pokud je T(x) libovolný polynom stupně, který je stejný jako řád přesnosti, dává konečná diference naprosto přesný výsledek.

4. Ti-1 Ti Ti+1 x PDE konečné diference NAP7 Alternativou odvození diferenční formule je Taylorův rozvoj. Uvažujme opět trojici uzlů a uzlové hodnoty vyjádříme Taylorovým rozvojem kolem centrálního uzlu i x Řád přesnosti x2 Čtvrtá derivace každého kubického polynomu je nulová, tudíž diferenční formule bude přesná nejen pro všechny kvadratické, ale i kubické polynomy

5. PDE konečné diference NAP7 Dále se zaměříme jen na evoluční problémy (PDE s derivacemi dle času), tj. na hyperbolické a parabolické rovnice. Tehdy je totiž nutné se zabývat i problémem stability řešení.

6. PDE má exaktní řešení tn+1 t tn x tn-1 xk-1 xk xk+1 PDE konečné diference NAP7 Nahrazení derivací PDE správnými konečnými diferencemi zaručí řád přesnosti (a garantuje velikost chyby jednoho časového kroku), ale vůbec ještě nezaručuje použitelnost navržené metody, která může být jen podmínečně stabilní nebo divergující. Vždy stabilní jsou pouze implicitní metody. Příklad: Vlnová rovnice v 1D popisující de facto jen transport funkce T(t=0,x) unášivou rychlostí c Metoda žabích skoků (leap frog) je ovšem vícekroková (a to je nevýhoda, je třeba použít nějakou jednokrokovou startovací metodu) Centrální náhrada prvních derivací (druhý řád přesnosti O(t2)+O(x2)). Podmínečně stabilní t<x/c. Asymetrická náhrada první derivace (první řád přesnosti O(t)+O(x2)). Vždy nestabilní!!!. Lax Friedrichs:Asymetrická náhrada první derivace (první řád přesnosti O(t)+O(x2)) umožňuje alespoň podmínečnou stabilitu t<x/c Laxův trik s modifikací členu Tkn nestabilní Eulerovy metody

7. PDE konečné diference NAP7 Zvýšení řádu přesnosti lze dosáhnout i u jednokrokových metod (počítajících hodnoty na hladině tn+1 jen z hodnot na hladině tn) podobně jako u metod typu Runge Kutta tím, že se vyčíslí hodnoty derivací v „mezičasech“. Příkladem je explicitní metoda Lax-Wendroff, která je druhého řádu přesnosti v čase i prostoru O(t2)+O(x2) a je podmínečně stabilní t<x/c. První krok – Lax Friedrichs s časovým krokem t/2 tn+1 t tn+1/2 Druhý krok – centrální náhrady derivací tn x xk-1 xk xk+1

8. MATLABLax (konvektivní přenos skoku) NAP7 Lax Friedrichs metoda c=1; počáteční podmínka je skoková funkce(tato skoková funkce by se správně měla jen posouvat vpravo rychlostí c) CFL=0.1 dt=0.011;c=1;l=1;n=101;dx=l/(n-1); cfl=c*dt/dx for i=1:n x(i)=(i-1)*dx; if x(i)<0.1 t0(i)=0; else t0(i)=1; end end tmax=1.;itm=tmax/dt; for it=1:itm for i=2:n-1 t(i)=(t0(i-1)+t0(i+1))/2-cfl*(t0(i+1)-t0(i-1))/2; end t(1)=t0(1);t(n)=t0(n); tres(it,:)=t(:); t0=t; end figure(1) hold off for ig=1:itm plot(x,tres(ig,:)) hold on end řešení je výrazně přetlumené, jako byste dali závodnímu koni prášek na spaní. Skoky jsou uměle rozmazány. To je důsledek prvního řádu přesnosti. CFL=1 jedině když nastavíte správný rytmus (časový krok) bude řešení přesné. Přesně se totiž respektuje oblast vlivu PDE a integruje se de facto dle jejícharakteristiky CFL=1.1

9. MATLABLax Wendroff (konvektivní přenos) NAP7 Lax Wendroffova metoda c=1; počáteční podmínka je skoková funkce CFL=0.1 dt=0.011;c=1;l=1;n=101;dx=l/(n-1); cfl=c*dt/dx for i=1:n x(i)=(i-1)*dx; if x(i)<0.1 t0(i)=0; else t0(i)=1; end end tmax=1.;itm=tmax/dt; for it=1:itm for i=2:n ta(i)=(t0(i-1)+t0(i))/2-cfl*(t0(i)-t0(i-1))/2; end ta(1)=t0(1); for i=2:n-1 t(i)=t0(i)-cfl*(ta(i+1)-ta(i)); end t(1)=t0(1);t(n)=t0(n); tres(it,:)=t(:); t0=t; end řešení se chvěje jako závodní kůň na startu, ale poměrně věrně modeluje realitu CFL=1 CFL=1.1 tady se řešení úplně zbláznilo, podívejte se na řád počítaných hodnot (1012).

10. t t/2 k-1 k k+1 MATLABLax Wendroff (vodní ráz) NAP7 Schema Lax Wendroff lze použít i pro dříve řešený problém vodního rázu Bernoulli Rovnice kontinuity a-rychlost zvuku Formální přepis do maticového tvaru Dva základní kroky explicitní metody Lax Wendroff

11. t t/2 p=0 v(t) xN=L x1=0 MATLABLax Wendroff (vodní ráz) NAP7 Problém bude (jako obvykle) s okrajovými podmínkami. Vlevo je zadaná rychlost v(t), vpravo tlak p(t) (a nedej bože Bernoulliho rovnice tuhé trubky). to je dáno parametry čerpadla v(t) Pozn.: Případ s napojenou tuhou trubicí znamená to, že hodnoty tlaku a rychlosti by se musí počítat ze soustavy rovnic

12. MATLABLax Wendroff (vodní ráz) NAP7 Model se zadaným lineárním nárůstem rychlosti na vstupu a konstantním tlakem na výstupu function [vrel,dvdt]=pump(t) if t<.03 vrel=0.35;dvdt=0; elseif t<0.06 vrel=0.35*(1+(t-0.03)/0.03);dvdt=35/3; else vrel=0.7;dvdt=0; end cfl=.5;mju=0.001;l=.18;d=0.015;rho=1000;a=8;pe=0; v0=pump(0); ez=frict(v0,d,rho,mju)*rho;dp=l*ez; n=101;dx=l/(n-1); dt=dx/a*cfl;tmax=.1;itmax=tmax/dt; aa=[0 1/rho;a^2*rho 0] for i=1:n w(1:2,i)=[v0;pe+dp-(i-1)*dx*ez]; b(1:2,i)=[-ez/rho;0]; end for it=1:itmax t=it*dt; for i=1:n-1 wn(1:2,i)=0.5*(w(1:2,i+1)+w(1:2,i))-dt/(2*dx)*aa*(w(1:2,i+1)-w(1:2,i))+dt/4*(b(1:2,i)+b(1:2,i+1)); end for i=2:n-1 b(1,i)=-frict(w(1,i),d,rho,mju); wc(1:2,i)=w(1:2,i)-dt/dx*aa*(wn(1:2,i)-wn(1:2,i-1))+dt*b(1:2,i); end [v1,dvdt]=pump(t); wc(1,1)=v1;wc(2,1)=wc(2,2)+rho*dx*(dvdt+frict(v1,d,rho,mju)); wc(1,n)=wc(1,n-1);wc(2,n)=pe; b(1,n)=-frict(wc(1,n),d,rho,mju); b(1,1)=-frict(v1,d,rho,mju); vres(it,1:n)=wc(1,1:n); pres(it,1:n)=wc(2,1:n); w=wc; end x=linspace(0,l,n); time=linspace(0,tmax,itmax); figure(1) contourf(x,time,pres,30) figure(2) contourf(x,time,vres,30) function ez=frict(v,d,rho,mju) re=abs(v)*d*rho/mju; if re<1 f=64; elseif re<2100 f=64/re; else f=0.316/re^0.25; end ez=f*v*abs(v)/(2*d);

13. Lax Wendroff CFL=0.9 v [m/s] Lax Friedrichs CFL=0.9 0.7 0.35 0.06 t [s] 0.03 Metoda charakteristik MATLABLax Wendroff (vodní ráz) NAP7 p(x,t) v(x,t) Porovnání tří různých metod (varianty metody konečných diferencí Lax Wendroff-druhý řád přesnosti a Lax Friedrichs-první řád přesnosti s metodou charakteristik dle předchozí přednášky). 0.18

14. MATLABFridrichs (+fixed tube) v [m/s] 0.7 0.35 0.06 t [s] 0.03 NAP7 p(x,t) v(x,t) Lax Friedrichs Pozn.: prakticky stejné výsledky získáme použitím metody Lax Wendroff CFL=0.9 N=101 CFL=0.9 N=201 0.1 0.18 CFL=0.99 N=201

15. MATLABFridrichs (+fixed tube) v [m/s] 0.7 0.35 0.06 t [s] 0.03 NAP7 CFL=0.9 N=201 Lax Friedrichs časový průběh tlaků na vstupu Pozn.: prakticky stejné výsledky získáme použitím metody Lax Wendroff 0.1 0.18 takto to vyšlo metodou charakteristik

16. PDE stabilitní analýza (von Neuman) NAP7 Vyzkoušeli jste si předchozí program na svém počítači? Ještě žijete? Zatím nic nenarušilo vaši stabilitu?

17. PDE stabilitní analýza (von Neuman) NAP7 Stabilitní analýza je založena na aproximací PDE i diferenční rovnice Fourierovým rozvojem, tj. členy typu kde km=m/L je vlnové číslo (diskrétní frekvence pohybující se v rozsahu Nyquistových frekvencí). Libovolnou počáteční podmínku T(0,x) lze nahradit Fourierovým rozvojem a zkoumat to, jak se po dosazení do diferenční rovnice vyvíjí jednotlivé členy rozvoje, zda rostou (|exp(at)|>1 nestabilita) nebo klesají (|exp(at)|<1stabilita). Proč je Eulerova formule nestabilní? Dosaďme do diferenční formule Courantovo číslo, tzv. CFL kriterium (Courant-Fridrichs-Levi) Faktor zesílení G – jak vidno jeho absolutní hodnota bude větší než 1 dokonce pro libovolné vlnové číslo

18. PDE stabilitní analýza (von Neuman) NAP7 Pokud nerozumíte předchozím vztahům, vzpomeňte na definici sin, cos A také na to, jak se spočítá kvadrát modulu komplexního čísla

19. PDE stabilitní analýza (von Neuman) NAP7 A proč je Laxova formule podmínečně stabilní? Dosaďme do Laxovy diferenční formule U Eulerovy metody to byla jednička Courantovo číslo  Faktor zesílení G je menší než 1 pro každé vlnové číslo (každou frekvenci), pokud je ovšem Courantovo číslo menší než 1

20. PDE stabilitní analýza (von Neuman) NAP7 Výsledkem stabilitní analýzy je stanovení faktoru zesílení G=exp(at), který určuje jak se změní komponenta Fourierova rozvoje odpovídající určitému vlnovému číslu během jednoho časového kroku. G je komplexní veličina (koeficient a je obecně komplexní číslo) a její absolutní hodnota (nebo G2) je indikátorem stability. Imaginární část faktoru zesílení G pak vyjadřuje fázový posun Fourierovy komponenty (správně by všechny komponenty měly být v našem případě transportovány stejnou rychlostí c). Fázový posuv a imaginární složka faktoru zesílení jsou měřítkem disperzní chyby numerické metody (projevující se nekorektními oscilacemi řešení). Absolutní hodnota G je mírou tlumení (numerické difuze). Hodnota >1 je naprosto nepřípustná (exponenciální nárůst chyby), hodnota <1 popisuje nekorektní vyhlazování skutečných změn. Ideální by byla jednička.

21. tn+1 t tn x tn-1 xk-1 xk xk+1 PDE konečné diference NAP7 Podobná situace a podobná analýza je i u PDE s druhými derivacemi. Příklad: Difuzní rovnice v 1D popisující evoluci teplotního profilu vyvolanou kondukcí Implicitní schéma – v každém časovém kroku je třeba řešit soustavu rovnic Implicitní formule Laasonen (první řád přesnosti O(t)+O(x2)). Vždy stabilní (libovolný časový krok) Explicitní schémata: Asymetrická náhrada první derivace (první řád přesnosti O(t)+O(x2)). Podminečně stabilní t<x2/a. Richardson: Centrální náhrada prvních derivací (druhý řád přesnosti O(t2)+O(x2)). Vždy nestabilní!! Du Fort Frankel: Modifikace Richardsonovymetody je vždy stabilní O(t2)+O(x2)+O((t/x)2). Du Fort Frankelův trik s modifikací Tkn nestabilní Richardsonovy metody

22. PDE stabilitní analýza (von Neuman) NAP7 Proč je standardní explicitní schema podmínečně stabilní? Dosaďme do diferenční formule Stabilitní omezení časového kroku (když zmenšíte x dvakrát, musíte zkrátit časový krok 4 krát) Proč je Laasonenova implicitní formule nepodmínečně stabilní? Naprosto stejným postupem získáme Bude vždy kladné pro libovolné vlnové číslo Pozor na znaménko

23. PDE stabilitní analýza (von Neuman) NAP7 Úplně stejně se dá dokázat nepodmíněná stabilita schematu Crank Nicholson záporné číslo v intervalu (-2,0)

24. Sušení kávového zrna NAP7 Sušení, ohřev, chlazení, to jsou typické problémy, popisované diferenciálními rovnicemi druhého řádu. Pro jednoduché geometrie je optimálním nástrojem řešení právě metoda konečných diferencí (zpravidla implicitní varianta u nelineárních problémů) nebo metoda konečných prvků.

25. Sušení kávového zrna NAP7 Katrin Burmester, Rudolf Eggers: Heat and mass transfer during the coffee drying process, ChISA 2008

26. R Sušení kávového zrna – popis procesu NAP7 • Při sušení kulové částice probíhají souběžně dva procesy: • Ohřev částice (povrchem, konvekcí, teplota T(t,r) směrem k povrchu roste) • Přenos vlhkosti (difúze póry k povrchu, měrná vlhkost X(t,r) s poloměrem klesá) • O tom, co se v sušené částici děje, existují různé představy: • Voda difunduje k povrchu a teprve tam dojde k jejímu odpaření (často používaný model, který však předpokládá, že póry částice se během sušení smršťují a tím vymačkávají vodu k povrchu. Alternativně by bylo možné akceptovat tento model tehdy, když by póry byly částečně zaplněny vzduchem přicházejícím z povrchu částice a vyplňujícím prázdné místo odkud vydifundovala voda) • K odpařování dochází již uvnitř částice, difunduje vodní pára a povrchová voda Tato představa je realističtější, ale vede ke komplikacím matematických modelů. Měrná vlhkost X je vlastně součtem vlhkosti vody Xw a páry Xv, zvlášť se musí řešit transportní rovnice pro páru a vodu a navíc je třeba nějak stanovit rychlost vnitřního odpařování. Možný kompromis může vyjít z předpokladu, že rychlost změny dX/dt je přímo rychlost odpařování, přesněji Chápejte to tak, že veškerá pára odchází pryč a na vlhkosti vzorku se nepodílí kde s je hmotnostní koncentrace sušiny (kg sušiny/m3 částice), X=Xw+Xv, a m je rychlost odpařování (kg vody/m3 částice za sekundu). • Všeobecně však platí souhlas o tom, že na povrchu částice je rovnováha a měrnou vlhkost (na povrchu) lze stanovit ze sorpční izotermy.

27. R Sušení kávového zrna – rovnice X NAP7 Ať již přijmeme hypotézu o povrchovém nebo vnitřním odpařování, lze transportní problém vody (páry) popsat difúzní rovnicí a okrajovými podmínkami X sorpční izotermy pro kávu kde Xeq je sorpční isoterma GAB tady by se hodil referát Lenky Dundáčkové   je relativní vlhkost vzduchu, , ,, jsou tři parametry modelu GAB, závislé na teplotě dle Arrheniova vztahu. Do efektivního difúzního součinitele je třeba zahrnout všechny nejistoty o mechanizmech transportu vody a páry, kontrakci pórů a především teplotní závislost difuzivity Aktivační energie Ea=29.5 kJ/mol, b=2.5 kg/kg. Dwa=2.26e-5 m2/s je součinitel binární difúze vody ve vzduchu. X je vlhkost materiálu (kg vody/kg sušiny). Platí pro zelená kávová zrna Robusta.

28. R Sušení kávového zrna – rovnice T NAP7 • Teplotní pole v částici bude nepochybně ovlivněno přijatými hypotézami: • Varianta povrchový odpar Teplo na odpar povrchové vody • Varianta vnitřní odpar Teplo na odpar vnitřní vody Teplo dodané ze vzduchu Součinitel přenosu tepla na povrchu kuličky lze stanovit z korelace Součinitel tepelné vodivosti závisí na vlhkosti materiálu X s tepelná vodivost sušiny (0.2 W/m/K káva), w vody

29. r1=0 r2=r rN=R Sušení kávového zrna – diskretizace X NAP7 Rozepsání derivací na pravé straně difúzní rovnice vlhkosti materiálu Použitím asymetrické náhrady časové derivace a centrálních náhrad prostorových derivací získáme soustavu rovnic pro uzlové vlhkosti Xk na nové časové hladině n k=2,3,…,N-1 Implicitní varianta Laasonen (O(r2)+O(t)) Okrajová podmínka na povrchu (k=N) je silná – přímo hodnota sorpční izotermy pro souběžně počítanou teplotu povrchu. Okrajová podmínka v ose vyplývá ze symetrie: buď jednoduše X1=X2, tím bychom ovšem metodu degradovali jen na první řád přesnosti, nebo z předchozí formule

30. Sušení kávového zrna – diskretizace T NAP7 Rovnice pro transport tepla je téměř stejná a stejný je i diferenční přepis k=2,3,…,N-1 Okrajová podmínka na povrchu (k=N) je dána konvektivním přenosem tepla Okrajová podmínka v ose vyplývá ze symetrie: buď jednoduše T1=T2, tím bychom ovšem metodu degradovali jen na první řád přesnosti, nebo

31. Sušení kávového zrna – konečné prvky NAP7 Alternativní metodou řešení je metoda konečných prvků (Galerkin) Aplikace integrace per partes Bázové a váhové funkce Ni(r) vedoucí na soustavu obyčejných diferenciálních rovnic Integrály matice hmot [[M]] i difuze [[K]] je možné například pro lineární bázové funkce odvodit v analytickém tvaru, ale někdy se spokojíme s přibližnou numerickou integrací Diagonalizovaná matice hmot (hmotnosti sférických slupek o tloušťce r) Tridiagonální matice difuze Když navíc nahradíme časovou derivaci vlhkosti X první diferencí, získáme podobnou soustavu rovnic jako v předchozím

32. Sušení kávového zrna – diskretizace a testy NAP7 Výsledné algebraické rovnice lze přepsat do tvaru s tridiagonální maticí soustavy Bývá užitečné ověřit zda je výsledná formule správná pro konstantní řešení X=1 Také je docela snadné upravit rovnice pro explicitní variantu řešení a testovat, kdy jsou všechny koeficienty hledaných parametrů X kladné To je docela zajímavá podmínka, která zaručuje splnění minimax principu, který praví, že řešení transportní rovnice bez zdrojových členů nesmí mít žádná lokální minima nebo maxima, ta mohou být jen na hranici oblasti řešení. Praktický důsledek je ten, že schéma, které má všechny koeficienty kladné, nebude nikdy oscilovat a bude stabilní. Z toho lze např. dovodit maximální velikost časového kroku, nutnou pro stabilitu Zcela stejné úpravy a předběžné testy lze provést i s rovnicí pro teploty

33. Sušení kávového zrna – diskretizace a testy NAP7 r=0.005;fi=0.05;ta=80;x0=0.3;t0=30;n=51;rho=1100;cp=2700;hlg=2.2e6; dr=r/(n-1);dr2=dr^2;dt=100;tmax=30000;ntim=tmax/dt; vel=5;alf=alpha(vel,r); for i=1:n to(i)=t0;xo(i)=x0;tn(i)=t0;xn(i)=x0;rk(i)=(i-1)*dr; end for it=1:ntim time(it)=(it-1)*dt/3600; for i=1:n dif(i)=def(to(i),xo(i)); lak(i)=lam(xo(i)); end for i=2:n-1 a(i)=-(dif(i)/dr2-dif(i)/(rk(i)*dr)-(dif(i+1)-dif(i-1))/(4*dr2)); c(i)=-(dif(i)/dr2+dif(i)/(rk(i)*dr)+(dif(i+1)-dif(i-1))/(4*dr2)); b(i)=1/dt+2*dif(i)/dr2; d(i)=xo(i)/dt; end a(1)=0;b(1)=1/dt+2*dif(1)/dr2;c(1)=-2*dif(1)/dr2;d(1)=xo(1)/dt; a(n)=0;b(n)=1;c(n)=0;d(n)=gab(fi,tn(n)); xn=tridag(a,b,c,d,n); for i=2:n-1 a(i)=-(lak(i)/dr2-lak(i)/(rk(i)*dr)-(lak(i+1)-lak(i-1))/(4*dr2)); c(i)=-(lak(i)/dr2+lak(i)/(rk(i)*dr)+(lak(i+1)-lak(i-1))/(4*dr2)); b(i)=rho*cp/dt+2*lak(i)/dr2; d(i)=rho/dt*(cp*to(i)+(xn(i)-xo(i))*hlg); end a(1)=0;b(1)=rho*cp/dt+2*lak(1)/dr2;c(1)=-2*lak(1)/dr2; d(1)=(rho*cp*to(1)+(xn(1)-xo(1))*hlg)/dt; a(n)=-lak(n)/dr; b(n)=lak(n)/dr+alf; c(n)=0; d(n)=alf*ta; tn=tridag(a,b,c,d,n); tres(:,it)=tn;xres(:,it)=xn; xo=xn;to=tn; end Laasonen Xn Laasonen Tn vnitřní odpar

34. Sušení kávového zrna – diskretizace a testy NAP7 Difúzní součinitel function d=def(t,x) r=8.314; d=2.26e-5*exp(-29500/(r*(t+273))-2.5*x); function xeq=gab(fi,temp) r=8.314;t=temp+273; kappa=1.3668*exp(-894/(r*t)); zeta=3.66e-17*exp(1.077e5/(r*t)); omega=0.0051*exp(-5725/(r*t)); xeq=100*zeta*kappa*omega*fi/((1-kappa*fi)*(1-kappa*fi+zeta*kappa*fi)); function l=lam(x) l=(x*0.6+0.2)/(1+x); function alf=alpha(vel,r) pr=0.71;ny=20e-6;lam=0.03;re=2*r*vel/ny; nulam=0.664*re^0.5*pr^(1/3); nut=0.037*re^0.8*pr/(1+2.443*(pr^(2/3)-1)/re^0.1); nu=2+(nulam^2+nut^2)^0.5; alf=nu*lam/(2*r); 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 GAB sorpční izoterma lambda Součinitel přenosu tepla TDMA

35. Sušení kávového zrna – diskretizace a testy NAP7 teplota povrchu povrchový odpar Porovnání variant s vnitřním a vnějším odparem vnitřní odpar