1 / 35

NAP 7

NUMERICAL ANALYSIS of PROCESSES. NAP 7. F inite difference method ( FD schemes and stability analysis). Hyperbolic equations - back to waterhammer Parabolic equation - temperature and concentration field (dr ying bean s ).

shaunam
Download Presentation

NAP 7

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. NUMERICAL ANALYSIS of PROCESSES NAP7 Finite difference method (FD schemes and stability analysis). Hyperbolic equations - back to waterhammer Parabolic equation - temperature and concentration field (drying beans) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

  2. Ti-1 Ti Ti+1 x PDE finite differences NAP7 All the types PDE (hyperbolic, parabolic and elliptic equations) can be solved by Finite Difference methods. Investigated region is covered with regular or irregular mesh of nodes, each node (tn,xi,yj,zk) is associated with calculated value Tnijk and in each node zero residual of PDE is required. First and second derivatives in PDE in node (tn,xi,yj,zk) are substituted by finite differences calculated from Tnijk in neigbouring nodes. Thus an algebraic equation is obtained in each node of mesh. For 1D equidistant nodes Asym. Diff. 1se order Centrál dif. 2ndorder error <K|x2|

  3. Ti-1 Ti Ti+1 x PDE finite differences NAP7 Finite difference formula (even for non-equidistant nodes) and order of accuracy can be derived from Taylor expansion of function T(t,x,y,z) in the vicinity of central node Tnijk or by using a polynomial interpolation. For example. Quadratic polynomial defined by 3 nodes x Contral formula of first and second derivatives follows from derivatives of polynomial at x=0 In other words: As soon as the T(x) is an arbitrary polynomial of the degree that is identical with the order of accuracy, finite differences are exact (identical with derivatives).

  4. Ti-1 Ti Ti+1 x PDE finite differences NAP7 Alternative derivation by using Taylor expansion. x Order of accuracy x2 Fourth derivative of any cubic polynomial is zero therefore this formula is exact not only for quadratic but also cubic polynomials

  5. PDE finite differences NAP7 Evolution problems

  6. PDE má exaktní řešení tn+1 t tn x tn-1 xk-1 xk xk+1 PDE finite differences NAP7 Stability analysis of explicit schemes (implicit schemes are always stable). Example: Wave equation in 1D describing transport of the function T(t=0,x) by velocity c Leap Frog method (multistep) Central approx.1st derivative (2Qnd order of accuracy O(t2)+O(x2)). Stable for t<x/c. Asymmetric 1st derivative (1st order of accuracyi O(t)+O(x2)). Always unstable!!!. Lax Friedrichs:Asymmetric substitution of 1st derivative (1st order of accuracy O(t)+O(x2)) Conditional stability for t<x/c Lax trick with modification of term Tkn in unstable Euler method

  7. PDE finite differences NAP7 Order of accuracy can be increased at single step methods in a similar way like Runge Kutta methods, by evaluating values of derivativbes at intermediate time steps. Example is explicite Lax-Wendroff method of the second order in time and space O(t2)+O(x2). The method is conditionally stable for t<x/c. First step – Lax Friedrichs with time step t/2 tn+1 t tn+1/2 Second step – central differencing tn x xk-1 xk xk+1

  8. MATLABLax (convective transport of a jump) NAP7 Lax Friedrichs method c=1; initial condition is the Heaviside unit step function 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 Discontinuity is smeared by numerical diffusion. CFL=1 Only at CFL=1 the solution is exact (the solution corresponds to MOC) CFL=1.1

  9. MATLABLax Wendroff (convective transport) NAP7 Lax Wendroff method c=1; initial condition is the Heaviside unit step function 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 Slightly oscillating but stable solution CFL=1 CFL=1.1 Stability is lost (see the order of magnitude)

  10. t t/2 k-1 k k+1 MATLABLax Wendroff (water hammer) NAP7 Lax Wendroff scheme can be used for solution of previous problem (water hammer) Bernoulli Continuity equation Formal transcript to matrixform Basic steps of explicit Lax Wendroff

  11. t t/2 p=0 v(t) xN=L x1=0 MATLABLax Wendroff (water hammer) NAP7 Boundary conditions: left (inlet) velocity v(t), pressure p(t)at outlet to je dáno parametry čerpadla v(t) Remark: the case with attached fixed tube at end means, that values of pressure and velocity must be solved from the system o equations

  12. MATLABLax Wendroff (water hammer) NAP7 Ramp function of velocity at inlet and constant pressure at outlet 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 MOC MATLABLax Wendroff (water hammer) NAP7 p(x,t) v(x,t) Comparison of 3 methods (Lax Wendroff-2nd orde,r Lax Friedrichs-1st order and MOC method of characteristics). 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 Remark: almost the same results can be obtained by the second order Lax Wendroff method 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 time course of pressure at inlet Remark: almost the same results can be obtained by the second order Lax Wendroff method 0.1 0.18 Results by MOC (method of characterisics)

  16. PDE stability analysis (von Neuman) NAP7 Have you testged the previous program at your computer? Are you still living? Your stability is stil undisturbed?

  17. PDE stability analysis (von Neuman) NAP7 Stability analysis is based on an approximation of PDE and differential equations by Fourier expansion, i.e. by terms of the type where km=m/L is a wave number (discrete frequency within the range of Nyquist frequency). Any initial condition T(0,x) can be substituted by Fourier expansion and to observe growth/decay of individual terms (|exp(at)|>1 growth, instability) or (|exp(at)|<1attenuation, stability). Why is the Euler‘s formula unstable? Courant number CFL criterion (Courant-Fridrichs-Levi) Gain factor G – absolute value is greater than 1 for arbitrary wave-number

  18. PDE stability analysis (von Neuman) NAP7 Hints

  19. PDE stability analysis (von Neuman) NAP7 Why is the Lax‘s formula conditionally stable? Substitute into Lax dif. formula It was 1 at Euler method Courant number  Gain factor G is less than 1 for any wave number if the courant number is less than 1

  20. PDE stability analysis (von Neuman) NAP7 The result of stability analysis the gain factor G = exp (at), which determines how the Fourier component corresponding to a specific wave number changesduring one time step. G is a complex quantity (coefficient is in general a complex number) and its absolute value (or G2) is an indicator of stability. The imaginary part of the gain factor G expresses the phase shift of the Fourier components (all components properly would be in this case transported at the same speed c). The phase shift and the imaginary part of the gain factor is a measure of dispersion errors numerical methods (manifested by unfair oscillations of the solution). The absolute value of G is a measure of damping (numerical diffusion). Value> 1 is absolutely inadmissible (exponential increase of errors), the value of <1 describes incorrect smoothing. The ideal would be 1.

  21. tn+1 t tn x tn-1 xk-1 xk xk+1 PDE finite differences NAP7 Analysis of PDE with second derivatives. Example: Diffusional equation in 1D for evolution of temperature profile controlled by conductionkondukcí Implicit scheme – solution of system of equations in each time step is necessary Implicit formula Laasonen (1st order O(t)+O(x2)). Always stable (arbitrary time step) Explicit scheme: Asym. substitution of 1st derivative (first order O(t)+O(x2)). Condit.stable for t<x2/a. Richardson: Central substitution of 1st derivative (2nd order O(t2)+O(x2)). Always unstable!! Du Fort Frankel: Modification of Richardsonmethod is always stable O(t2)+O(x2)+O((t/x)2).

  22. PDE stability analysis (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 stability analysis (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

More Related