NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP4 Nelineární modely systémů popisované obyčejnými diferenciálními rovnicemi Počáteční problém (Runge Kutta, Adams) Dopravní zpoždění Podivné atraktory Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. Modely ODE počáteční problém NAP4 Systémy, které jsou popisované soustavou obyčejných diferenciálních rovnic, a jejichž řešení je plně popsáno výchozím stavem, jsou například integrální modely hmotnostních a entalpických bilancí elementárních jednotek („lumped parameter“, nebo „compartment“ modely). Cílem je zjistit vývoj stavových veličin (koncentrací, teplot) v čase, jde tedy o evoluční problémy. Těmi elementárními jednotkami mohou být konkrétní aparáty (reaktory, absorbéry, kolony, zásobníky, filtry,…) v „spread sheetech“ nebo takové zóny v aparátech, jejichž stav lze charakterizovat střední teplotou, a koncentrací složek (ideální mísič, zóna pístového toku). Příkladem je popis míchané nádoby soustavou ideálních mísičů pod a nad turbinovým míchadlem, nebo nahrazení kotlového výměníku tepla soustavou elementárních výměníků mezi přepážkami. Každé elementární jednotce odpovídají obyčejné diferenciální rovnice popisující bilanci hmoty (pro každou složku jedna diferenciální rovnice popisující časovou změnu koncentrace) a jedna bilanční rovnice entalpie (nebo celkové energie), která popisuje vývoj teploty kompartmentu objemový průtok ze sousedních kompartmentů objemový zdroj tepla entalpické toky ze sousedních kompartmentů

3. Modely ODE počáteční problém NAP4 Systém nemusí být (i když bývá) popisovánobyčejnými diferenciálními rovnicemi jen s prvními derivacemi. Např. v rovnicích kmitání budou i druhé derivace času (zrychlení) Jenomže každou takovou rovnici lze nahradit soustavou rovnic jen s prvními derivacemi, např. Počáteční problém je charakterizován tím, že každé počítané veličině odpovídá jedna diferenciální rovnice prvního řádu a každé diferenciální rovnici právě jedna počáteční podmínka (hodnota počítané veličiny) odpovídající výchozímu stavu, tedy pro stejnou hodnotu nezávisle proměnné, např. času. aneb, kolik rovnic tolik počátečních podmínek

4. Modely ODE počáteční problém NAP4 model zpětného promíchávání (disperze, difuze) Příklad modelu průtočného systému: (1+r+b)Q 6 5 (1+r)Q 4 2 1 Q,x(t) bQ Ci, Ti sQ rQ 3 7 model stagnačních zón model recirkulace s dopravním zpožděním Diferenciální rovnice koncentrací Parametry určující rozdělení průtoků (s-stagnantní zona, r-recirkulace, b-zpětné promíchávání) musí být zadané, r-třeba z měření průtokoměrem, b-ze vztahů pro disperzi. Stejně tak objemy jednotlivých kompartmentů. Nejčastěji se ale tyto parametry zjišťují optimalizací shody predikovaných charakteristik (třeba RTD-Residence Time Distribution) s experimentem. A dle těchto výsledků se hodnotí velikost recirkulačních nebo stagnačních zón, zkratové toky, velikost disperze, aktivní objem aparátu apod.

5. Modely ODE Fourierova transformace NAP4 Model lze řešit i metodou Fourierovy transformace (podobně jako v předchozí přednášce). Nestačí ovšem pouze aplikace konvoluce, protože model zahrnuje i kompartmenty se zpětným promícháváním, stagnantní zóny… Lze však vyjít z uvedených diferenciálních rovnic, a na každou z nich použít Fourierovu transformaci. Pro Fourierovu transformaci derivace platí (za předpokladu že koncentrace v čase nekonečno jsou nulové) Diferenciální rovnice tím přejdou na rovnice algebraické, kde místo času t bude frekvence f. Tato soustava algebraických rovnic (lineárních!) se už dá poměrně snadno řešit, a výsledkem budou Fourierovy transformace koncentrací. Průběhy koncentrací v čase se pak získají zpětnou FFT. Řešení je přesné a rychlé (při použití FFT), ale je omezené jen na lineární modely. Postup je výhodný tehdy, když jsou koeficienty soustavy (r,b,s) konstantní. Jinak by totiž bylo nutné počítat Fourierovu transformaci součinu dvou funkcí (a ta se nerovná součinu transformací). V obecném případě a s nelineárními členy je proto lepší zůstat v časové rovině a hledat řešení numerickou integrací ODE.

6. Modely ODE Laplaceova transformace NAP4 U Laplaceovy transformace je to podobné (stejný vztah pro konvoluci a korelaci), analogický pro derivace Zpětná Laplaceova transformace se obvykle hledá v tzv. slovnících, např.

7. ODEpočáteční problém numerické řešení c c ck ck t t tk+1 tk+1 tk tk NAP4 Jednokrokové metody aproximují řešení v čase tk+1=tk+ na základě odhadu střední derivace v intervalu <tk,tk+1> Eulerovy metody snadno naprogramujete i v Excelu první řád přesnosti znamená, že chyba numerického řešení je přímo úměrná kroku  když napíšete Taylorův rozvoj pro jeden časový krok, vyjde chyba O(2), tedy druhý řád přesnosti. Chyby se ale akumulují, takže po 1/ krocích bude výsledná chyba O(). Eulerova metoda explicitní řád přesnosti  (první) stabilní jen pro malý integrační krok  Eulerova metoda implicitní řád přesnosti  (první) stabilní, ale je nutné iterovat (protože ck+1 neznáme)

8. ODEpočáteční problém numerické řešení c ck t tk+1 tk NAP4 Jednokrokové metody aproximují řešení v čase tk+1=tk+ na základě odhadu střední derivace v intervalu <tk,tk+1> Metody typu Runge Kutta počítají střední hodnotu derivace jako vážený průměr derivací f(t,c). Nejčastěji používanáje varianta se čtyřmi vyhodnocovanými body (krajní body tk, tk+1 a dva body uprostřed) Tato varianta je implementována v MATLABu jako funkce ode45 (aspoň si to myslím, z manuálu to jasné není) Váhové koeficienty (1/6,1/3,1/3,1/6) a přírůstky c jsou navrženy tak, aby řád přesnosti byl co nejvyšší (tato varianta je čtvrtého řádu přesnosti). Obecně platí, že řád přesnosti metod RK je stejný jako počet vyhodnocovaných bodů.

9. ODEpočáteční problém numerické řešení c ck ck-1 t tk+1 tk-2 tk-1 tk NAP4 Vícekrokové metody aproximují řešení v čase tk+1 na základě několika předchozích hodnot, např. v časech tk, tk-1 (dvou), ev. tk-2 (tříkroková metoda). Zpravidla se jedná o kombinací prediktoru (explicitní) a korektoru (implicitní formule s vyšší přesností a stabilnější) Stačí vyčíslit jen jedinou hodnotu f Obě formule jsou druhého řádu přesnosti, ale korektor (implicitní formule) má v absolutní hodnotě chybu zhruba 9-krát menší (implicitní formule mívají chybu menší než explicitní, i když je řád chyby stejný). Nevýhodou je to, že pro první dva kroky se musí použít nějaká jednokroková metoda, třeba RK. Počet vyčíslovaných derivací je ale menší.

10. ODEpočáteční problém numerické řešení NAP4 Dynamické nastavení velikosti integračního kroku  Nejjednodušší řešení spočívá v tom, že se každý krok provede dvakrát, jednou s krokem  a pak s krokem /2 (dvakrát po sobě). Když je rozdíl vyšší než zadaná tolerance, musí se krok  zmenšit. Není to zase tak neefektivní, protože když je znám řád přesnosti metody dá se z těch dvou výsledků (s plným a polovičním krokem) finální hodnota ještě zpřesnit (tzv.Aitkenova extrapolace). MATLAB používá dynamickou volbu časového kroku automaticky. Stiff problémy Když se řeší soustavy rovnic může se stát, že integrační krok je diktován rovnicí s extrémně krátkou časovou konstantou, a stabilita řešení pak vyžaduje extrémně krátký integrační krok (výpočet se neúnosně prodlužuje). Říkáme, že systém je stiff. Některé metody integrace ODE vliv těchto nepříjemných rovnic potlačí (příkladem je Gearova metoda). Dopravní zpoždění V ukázkovém příkladu bylo do recyklu zařazeno dopravní zpoždění, tj. posun signálu o určitý čas. Implementace řešiče si musí pamatovat celou historii počítaného řešení (a nejen poslední krok) a musí být definovány i hodnoty řešení pro záporné časy (tj. kompletní historie).

11. ODEpočáteční problém numerické řešení NAP4 Stabilitní analýza Uvažujme diferenciální rovnici popisující třeba rychlost chemické reakce 1. řádu stabilitní podmínka omezení velikosti časového kroku  (koeficient zesílení menší než 1) Eulerova explicitní integrační formule prvního řádu Eulerova implicitní integrační formule prvního řádu koeficient zesílení bude menší než 1 pro libovolné 

12. ODEpříklady řešení v MATLABu NAP4 ode45, ode23Runge-Kutta (jednokrokové metody s různým řádem přesnosti) ode113Adams (vícekroková metoda) ode15sStiff (Gearova vícekroková metoda s variabilním řádem přesnosti – vhodná pro stiff problémy) dde23 dopravní zpoždění

13. RTD série mísičů 1/4 x(t) y(t) (t) c1 c2 c3 c4 NAP4 Příklad, který byl řešen v předchozích přednáškách použitím FFT. Série 4 ideálně míchaných nádob (objem nádoby V, průtok Q). Cílem řešení je nalézt časové průběhy koncentrací c1,…,c4 pro libovolné počáteční podmínky a pro libovolný průběh koncentrace na vstupu x(t). Hmotnostní bilance tm=2 odpovídá poměru V/Q=0.5 (to je střední doba prodlení jedné nádoby) Jako vstupní koncentrace x(t) bude použita impulzní odezva soustavy dvou mísičů (se stejným objemem nádoby) a nulové počáteční koncentrace ve všech nádobách. Důvodem takového zadání je to, že existuje analytické řešení, které můžeme porovnat s výsledkem numerického řešení soustavy ODE. Vzruchová funkce x(t) bude impulsní odezva E(t) pro N=2 a tm=1, výstupní koncentrace y(t)=c4(t) bude impulzní odezva pro N=6 a tm=3

14. RTD série mísičů 2/4 NAP4 MATLAB nabízí pro řešení soustavy ODE řadu metod (ode45, ode113, ode15s, …), které si automaticky nastavují integrační krok tak, aby bylo dosaženo požadované přesnosti. Způsob zadání je stejný, nejprve se definuje funkce, jejímž výstupem je vektor prvních derivací pro danou hodnotu integrační proměnné (času) a vektoru řešení. V našem případě je třeba definovat funkci jejímž výsledkem je sloupcový vektor (4x1) pravých stran 4 diferenciálních rovnic. Napíšeme funkci, která bude předpokládat, že vzruchová funkce x(t) je impulzní odezvou serie nx mísičů, se střední dobou tmx. Parametr tm je střední doba modelovaného systému 4 mísičů. function dy = serie4(t,y,tm,tmx,nx) dy = zeros(4,1); xx=efun(t,tmx,nx); tm1=tm/4; dy(1)=(xx-y(1))/tm1; dy(2)=(y(1)-y(2))/tm1; dy(3)=(y(2)-y(3))/tm1; dy(4)=(y(3)-y(4))/tm1; Bylo by asi elegantnější, kdyby vzruchová funkce x(t) byla přímo dalším formálním parametrem funkce serie4. S tím mám ale v Matlabu trochu problém, třeba mi poradíte. Funkci efun bylo možné definovat jako subfunkci přímo v m-souboru serie4.m, ale možná je lepší ji definovat jako zvláštní m-soubor efun.m function y=efun(t,tm,n) y=n^n*t^(n-1)/(factorial(n-1)*tm^n)*exp(-t*n/tm);

15. RTD série mísičů 3/4 NAP4 Provedení integrace zajistí volání standardní funkce se třemi vstupními parametry Sol = ode45(funkce popisující pravou stranu, rozsah času, počáteční podmínky) v našem případě např. sol=ode45(@(t,y)serie4(t,y,2,1,1),[0 10],[0 0 0 0]) Vektor počátečních podmínek (nulové koncentrace) solver: 'ode45' extdata: [1x1 struct] x: [1x29 double] y: [4x29 double] stats: [1x1 struct] idata: [1x1 struct]. Integrace od t=0, do t=10. Zdá se, že Matlabu se nedá vnutit konstantní integrační krok. Výsledkem je struktura sol v níž je schován vektor časových kroků sol.x (v tomto případě řešení pro zadanou přesnost vyšlo na 29 kroků) a korespondující matice sol.y výsledných koncentrací (čtyři řádky odpovídají čtyřem koncentracím) Použití funkce jako skutečného parametru vypadá v Matlabu trochu divně. Nestačí napsat prostě jméno této funkce (serie4), protože ode45 požaduje, aby tento parametr byla funkce jen se dvěma parametry (t,y) a funkce serie4 má těch parametrů 5. Proto se použije tzv. anonymní funkce @(t,y)výraz.

16. RTD série mísičů 4/4 NAP4 Numerické řešení metodou Runge Kutta (ode45) můžeme porovnat s analytickým řešením E(t) pro N=6 a tm=3 Ukázka toho jak definovat jednopříkazovou funkci (bez použití funkce v samostatném m-souboru) e=@(t,tm,n) n^n.*t.^(n-1)./(factorial(n-1)*tm^n).*exp(-n.*t/tm) t=linspace(0,10.23,1024); y=e(t,3,6); Standardní funkce linspace umožní vygenerovat vektor 1024 ekvidistantních časů Sy=deval(sol,t) tyto body jsou výsledkem řešení funkce ode45 s automatickým nastavováním integračního kroku funkce deval interpoluje neekvidistantní výsledky řešení sol.y dle předepsané časové škály t

17. RTD dopravní zpoždění NAP4 Následující příklad dokumentuje to, že když se v systému objeví dopravní zpoždění (dlouhé potrubní úseky, recykly), nelze použít standardní integrační funkce (ode23,ode45,…), ale jinou „rodinu“ MATLABovských funkcí, které si umí uchovávat předchozí výsledky (s předepsaným vektorem různých zpoždění 1 2 …)

18. RTD dopravní zpoždění 1/2 x(t) (1+r)Q 1 2 Q rQ NAP4 Příklad: průtočný systém, tvořený dvěma mísiči a recyklem s dopravním zpožděním . Je to dost běžné: každé potrubí propojující aparáty lze modelovat dopravním zpožděním. Cílem řešení mají být časové průběhy koncentrací c1(t), c2(t), při známých hodnotách průtoků (Q), poměru recirkulace (r) a objemů nádob (V1,V2), a pro zadaný časový průběh koncentrace x(t) na vstupu. Systém popisují ODE hmotnostních bilancí: dopravní zpoždění

19. RTD dopravní zpoždění 2/2 NAP4 V MATLABu je třeba použít speciální řešiče ODE (dde23 pro konstantní zpoždění, resp ddesd pro variabilní zpoždění), např. sol=dde23(@derivace, [1, 2…],historie,[tmin,tmax]) function dcdt=dclag(t,c,z) q=1; v1=2; v2=1; r=10; if t<.2 x=1; else x=0; end dcdt=zeros(2,1); clag=z(:,1); dcdt(1)=q/v1*(x-c(1)); dcdt(2)=q/v2*(c(1)+r*clag(2)-(1+r)*c(2)); z je matice hodnot řešení v čase t-k. Index sloupce (1) je indexem vektoru posunutí nulová historie c1=c2=0 pro t<0 sol=dde23(@dclag,1,[0;0],[0,20]) plot(sol.x,sol.y) vektor posunutí 1=1 (zde má vektor jen jeden prvek)

20. Trajektorie částic NAP4 Počáteční problém se týká i popisu pohybu částic, na které působí síly. Příklady: trajektorie kapiček v komoře rozprašovací sušárny, pohyb částic ve fluidním loži, hořících uhelných částic ve spalovací komoře, částic v míchadle, usazovací nádrži, apod. Pohybové rovnice jsou ODE typu zrychlení vztlakové, elektrické či odstředivé síly síla odporu při pohybu částice v tekutině s rychlostí v

21. Trajektorie částic NAP4 V posledních přednáškách tohoto kurzu budeme věnovat pozornost i moderním metodám typu DEM (Discrete Element Method), které provádí simulace dynamického chování partikulárních systémů jako jsou výsypky, sila, mlýny, mísiče, separátory, pevná a fluidní lože. Výpočty jsou založené na řešení obyčejných diferenciálních rovnic pohybu a na uvažování vzájemných interakcí jednotlivých částic. Verletova integrace pohybových rovnic

22. Teorie chaosu podivný atraktor NAP4 Trajektorie pohybu částic počítal i meteorolog E.Lorenz (1963) a snažil se namodelovat přirozenou konvekci v atmosféře, vrstvě vzduchu, která je zespodu ohřívána a shora chlazena. Drastickým zjednodušení získal soustavu tří obyčejných diferenciálních rovnic pro souřadnice x,y,z částice v atmosféře Rayleigh-Bénárdova nestabilita Model má pouze 3 parametry, Rayleighovo číslo Ra, Prandtlovo číslo Pr a parametr b je koeficient štíhlosti rotujících buněk. Tyto 3 rovnice v MATLABu vyřešíme snadno, viz dále…

23. Teorie chaosu podivný atraktor NAP4 trajektorie částice pro Ra=10, Pr=28, b=8/3 function dy = chaos(t,y) dy = zeros(3,1); dy(1) = 10*(y(2)-y(1)); dy(2) = -y(1)*y(3)+28*y(1)-y(2); dy(3) = y(1)*y(2)-8/3*y(3); sol=ode45(@chaos,[0 30],[.10001 .1 .1]) plot(sol.y(2,:),sol.y(3,:)) tt=linspace(0,30,3000); sy=deval(sol,tt); plot(sy(2,:),sy(3,:)) počáteční souřadnice x,y,z integrace do 30 s vykreslení grafu y-z z automaticky vypočtených integračních kroků

24. Teorie chaosu podivný atraktor NAP4 Numericky vypočtené trajektorie se při hodnotách Pr>1 chovají chaoticky (nekonvergují k bodu x=y=z=0, který by samozřejmě splňoval rovnice soustavy), ale vytváří ve fázovém prostotu (např. y-z) něco, co vypadá jako protisměrně rotující víry konečných rozměrů, a tyto trajektorie zdánlivě náhodně „přeskakují“ z levého do pravého křídla, nikde se neprotínají, jsou nesmírně citlivé na počáteční podmínky a na jakékoliv poruchy (tedy i aproximační chyby numerické integrace). Limitnímu stavu se říká „podivný atraktor“ parametrická trajektorïe částice z Takové chování (deterministický chaos) je podmíněno nelinearitou soustavy a překročením meze stability (např. Prantlova, Reynoldsova nebo Rayleighova čísla). Deterministický chaos je typický pro turbulentní proudění. y

25. Co je třeba si pamatovat NAP4 Prosím, zapamatujte si z této přednášky věnované metodám numerické integrace diferenciálních rovnic (pro počáteční problém) alespoň to, co je na následující folii

26. Co je třeba si pamatovat NAP4 Fourierova transformace derivací Co je to počáteční problém Co je to řád přesnosti O(n) Explicitní a implicitní Eulerova metoda Princip RK metoda a metod vícekrokových (prediktor, korektor) Princip stabilitní analýzy Princip optimalizace integračního kroku