NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP6 Parciální diferenciální rovnice (PDE), klasifikace na hyperbolické, parabolické a eliptické Hyperbolické PDE (kmitání táhel, nosníků, ráz v potrubí) MOC-metoda charakteristik (stlačitelné proudění) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. PDE parciální diferenciální rovnice NAP6 L.Wagner

3. PDE parciální diferenciální rovnice NAP6 • Když je pro popis systému nutný větší počet nezávisle proměnných (třeba prostorové souřadnice a čas) musíme řešit místo obyčejných parciální diferenciální rovnice (PDE). Většinou jsou to diferenciální rovnice s maximálně druhými derivacemi závisle proměnné(výjimkou je např. biharmonická rovnice deformace membrán se čtvrtými derivacemi): • Hyperbolickérovnice (Typické pro popis kmitání a vln, třeba elastických vln při rázové deformaci, a pro nadzvukové proudění. Klíčovou roli hraje konečná rychlost šíření tlakové vlny.) • Parabolickérovnice(Typické pro evoluční problémy, třeba vývoj teplotního nebo koncentračního profilu v čase, ale třeba i vývoj profilů v mezní vrstvě v závislosti na vzdálenosti od vstupu. Můžeme to chápat i jako mezní případ hyperbolických rovnic s nekonečnou rychlostí šíření poruchy tlaku.) • Eliptickérovnice(Typické pro stacionární problémy popisu rozložení teplot, koncentrací, deformací. Často jde o problém popisovaný parabolickými rovnicemi, ale až po ustálení, tj. pro nekonečně dlouhý čas. I většina pružnostních problémů je popisována eliptickými rovnicemi. ) Příklad: vlnová rovnice kmitání pružné tyče, ráz v potrubí Příklad: Fourierova rovnice vedení tepla Příklad: Poissonova rovnice

4. PDE parciální diferenciální rovnice NAP6 Typ PDE je dán pouze koeficienty druhých derivací. Každou PDE druhého řádu můžeme napsat ve tvaru kde a,b,c,f mohou být libovolné funkce x,y i hledaného řešení , včetně jeho prvních derivací. Koeficienty a,b,curčují rovnice tzv. charakteristik, čar závislostí y(x), které splňují rovnici b2-4ac>0 hyperbolická rovnice (existuje dvojice reálných kořenů, tudíž dvojice charakteristik) b2-4ac=0 parabolická rovnice (existuje jen jeden kořen a jedna charakteristika) b2-4ac<0 eliptická rovnice (kořeny jsou imaginární a žádná reálná charakteristika neexistuje) hyperbolická PDE parabolická PDE eliptická PDE y y y x x x znázornění oblasti vlivu (modrá oblast pod charakteristikami) a okrajových podmínek (červeně) na hodnotu řešení (x,y)

5. PDE Hyperbolické PDE NAP6 Všude tam, kde se projeví setrvačné síly a stlačitelnost, popisují problém hyperbolické parciální diferenciální rovnice, charakterizované reálnými charakteristikami (jejichž vlastnosti se někdy využívají při řešení hyperbolických PDE metodou charakteristik). Tato metoda ale není jediná, zhusta se používá metoda konečných prvků, zvláště při analýze kmitání pružných konstrukcí (nosníků, skořepin,…). Lze využít výsledky získané v předchozí přednášce.

6. fx(x), ux(x) PDE kmitání táhla NAP6 Rovnici popisující kmitání táhla (elastické tyče) získáme doplněním rovnice rovnováhy o zrychlující sílu vidíte, že je to opravdu hyperbolická rovnice, znaménka druhých derivací jsou opačná Metoda vážených residuí založená na aproximaci Nj v roli bázové funkce Ni v roli váhové funkce Integrace per partes aplikovaná na první člen matice hmot[[M]] budicí síly [[b]] matice tuhosti [[K]]

7. PDE kmitání táhla NAP6 Uvažujme volné kmitání, bez budících sil což je soustava jen obyčejných diferenciálních rovnic druhého řádu. Ale s konstantními koeficienty, takže existuje analytické řešení jehož dosazením do předchozí rovnice dostaneme soustavu algebraických rovnic pro koeficienty amplitud kmitů Protože je soustava homogenní, bude mít netriviální (nenulové) řešení jen když je matice soustavy singulární, tj. když Čtvercové matice [[K]] i [[M]] jsou dané a mají N řádků. Determinant je de facto polynom N-tého stupně proměnné 2a algebraická rovnice má N-kořenů, N-vlastních frekvencí. Tyto frekvence  spolu s vektory amplitud lze nalézt řešením vlastního problému v MATLABU lze tento problém řešit funkcí [u,omega]=eig(m\k)

8. PDE kmitání táhla NAP6 Pro lineární bázové funkce popisu amplitud v obecném elementu snadno stanovíme matici tuhosti elementu (stejná jako dřív) račte si povšimnout, že součet všech prvků se rovná hmotnosti celého elementu stejně jako matici hmot Matice hmot se často nahrazuje nahrazuje diagonalizovanou maticí, která odpovídá rozdělení hmotnosti táhla do jeho uzlů Výhoda diagonalizované matice hmot spočívá v tom, že inverze [[M]] je triviální.

9. PDE kmitání táhla MATLAB NAP6 Příklad: Soustava ocelových táhel (test pro jedno táhlo) x=[0 1];con=[1 2];a=[1e-4]; nu=length(x);ne=length(a);n=nu; k=zeros(n,n); m=zeros(n,n); for e=1:ne [kl,ml,ig]=kloc(e,x,con,a); k(ig(1:2),ig(1:2))=k(ig(1:2),ig(1:2))+kl(1:2,1:2); m(ig(1:2),ig(1:2))=m(ig(1:2),ig(1:2))+ml(1:2,1:2); end [u,omega]=eig(m\k); function [kl,ml,ig]=kloc(ie,x,con,a) E=200e9;R=7000; ae=a(ie); i1=con(ie,1);i2=con(ie,2); ig=[i1 i2]; le=abs(x(i1)-x(i2)); kl=E*ae/le*[1 -1;-1 1]; ml=R*ae*le/6*[2 1;1 2]; V tomto případě můžeme snadno zkontrolovat výsledek (pro diagonalizovanou matici hmot vychází vlastní frekvence menší )

10. PDE kmitání nosníků NAP6 Možná důležitější je případ příčného kmitání nosníků, např. trubek výměníku tepla, podepřených v několika místech přepážkami a na koncích vetknutých do trubkovnice. Použijte výsledky předchozí přednášky, kde spojité příčné zatížení nosníku nahradíte setrvačnými silami matice hmot[[M]] matice tuhosti [[K]] Dál už je to stejné jako u táhla, použijete jenom dříve definované matice tuhosti a hmot odvozené pro nosníkový prvek. A nemusíte si dělat starosti s pravou stranou, protože místa podepření trubek i jejich vetknutí do trubkovnice jsou silné okrajové podmínky.

11. PDE Ráz v potrubí NAP6 Nestacionární tok stlačitelné tekutiny v potrubních sítích, fungování kardiovaskulárního oběhu, to jsou typické příklady hyperbolických rovnic. Zvláštní případ reprezentuje hydraulický ráz, ke kterému dojde při náhlém uzavření ventilu v potrubí… projeví se náhlým vzrůstem tlaku (rázovou vlnou), který se šíří rychlostí zvuku v tekutině a stěně potrubí. Víte jak tu rychlost spočítat? Víte, jaký maximální tlak může vzniknout? Barbara Wagner M. Rohani, M.H. AfsharSimulation of transient flow caused by pump failure: Point-Implicit Method ofCharacteristics Annals of Nuclear Energy, Volume 37, Issue 12, December 2010, Pages 1742-1750

12. PDE Ráz v potrubí NAP6 Formulace problému: Potrubím s proměnným průřezem A(t,x), protéká tekutina střední rychlostí v(t,x).Tekutina je stlačitelná a souvislost její hustoty a tlaku je charakterizována objemovým modulem stlačitelnosti K [Pa] Śouvisí s rychlostí zvuku v tekutině a Je to vlastně stejný problém jako dříve uvažované potrubní sítě, kdy cílem bylo stanovit průběhy tlaku a průtoky. Tentokrát však tlaky i rychlosti závisí i na čase, p(t,x), v(t,x). Tyto veličiny popisuje stejně jako ve stacionárním případě dvojice rovnic, rovnice kontinuity a bilance hybnosti (Bernoulliho rovnice). Ve zjednodušeném tvaru Rovnici kontinuity v této poněkud zvláštní podobě vysvětlím na následující folii takto se zohlední i pružnost potrubí (závislost průřezu na tlaku) Tuto rovnici poznáváte: Bernoulliho rovnice kde je zanedbán konvektivní člen (kinetická energie)

13. ,A,v dx PDE Ráz v potrubí NAP6 Rovnice kontinuity pro případ, kdy průřezová plocha trubkyA(t,x) se mění po délce i v čase (elastická, nebo elasticko-plastická trubka, která se „nafukuje“ působením vnitřního přetlaku) hmotnostní bilance elementu konstantní délky dx akumulace mění se hustota, průřez i střední rychlost . souvislost hustoty a tlaku přítok-odtok Rozderivováním levé i pravé strany a využitím definice objemového modulu stlačitelnosti získáme zcela obecnou rovnici kontinuity (to je totéž, jen zapsané v materiálových derivacích) Tento člen je v předchozí rovnici zanedbán (dělá se to často, viz např.Khamlichi, Wave motion 1995)

14. PDE Ráz v potrubí NAP6 To, že je výše uvedený systém rovnic hyperbolický, vyplyne z eliminace tlaku nebo rychlosti (stačí derivovat rovnici kontinuity dle času, Bernoulliho rovnici dle souřadnice x, a odečíst). Výsledné rovnice pro tlak p(t,x) nebo rychlost v(t,x) Rychlost zvuku jsou stejné a dokonce stejné jako hyperbolická rovnice kmitání táhla. Je tedy možné použít i stejnou metodu řešení (metodu vážených residuí). Objeví se ovšem malý technický problém: I když jsou rovnice stejné, mají různá řešení, protože jsou různé okrajové podmínky, např. časový průběh průtoku na jednom konci a časový průběh tlaku na druhém konci. V případě rázu v potrubí je to třeba skoková změna průtoku na jednom konci (zavření ventilu) a konstantní tlak na druhém konci (zásobník). Metodou vážených reziduí je tedy nutné řešit soustavu obou PDE.

15. PDE NAP6 Pro hyperbolické rovnice se ale často využívá jiná metoda, která využívá specifickou vlastnost hyperbolických rovnic a tou je existence dvou reálných charakteristik. Nazývá se metoda charakteristik (MOC). Doporučená literatura Wylie, Streeter: Fluid Transients. McGraw Hill, 1978 Zábranský

16. Metoda charakteristik MOC NAP6 Co vlastně chceme řešit? Proudění stlačitelné tekutiny v trubce nebo nestlačitelné kapaliny v elastické trubce (která se může „nafukovat“ při zvýšení tlaku) – výsledný efekt stlačitelnosti je úplně stejný. Ve hře jsou dva vlivy: setrvačnost tekutiny a tlak vyvolaný jejím stlačením. Cílem je stanovení střední rychlosti v (tudíž i průtoku q) a tlaku p, přičemž obě tyto veličiny (v,p) se mění v čase i po délce trubky. Budeme uvažovat různé okrajové podmínky na koncích trubky: buď zadávaný (a časově proměnný) průtok, nebo zadávaný (a časově proměnný) tlak.

17. Metoda charakteristik MOC t x() t= x NAP6 Vraťme se k výchozí soustavě rovnice kontinuity a rovnice Bernoulliho spoustu věcí jsme zanedbali, například unášivou rychlost, takže následující analýza je správná jen tehdy, když rychlost zvuku je mnohem vyšší než rychlost proudění. V řadě uváděných příkladů tento předpoklad splněn není. a tyto rovnice (první z nich vynásobme libovolnou nenulovou konstantou ) sečteme Zvolme libovolnou křivku v rovině t-x (parametrem křivky x() může být přímo čas t). Úplné diferenciály tlaku a rychlosti podél této křivky jsou Jestliže bude křivka x(t=) splňovat rovnice (současně) získáme finální rovnici vyjádřenou pomocí úplných diferenciálů podél křivky x(t=) (a tuto křivku nazveme charakteristikou diferenciální rovnice)

18. C x2=-at x1=at h/a B A h Metoda charakteristik NAP6 Předchozí rovnice má dvě řešení a jim odpovídají dvě charakteristické křivky a diferenciální rovnice, které je třeba na nich integrovat Třecí ztráty vyjádřené součinitelem třecích ztrát f Přibližná integrace těchto rovnic ve směru charakteristik vede na soustavu dvou algebraických rovnic pro neznámé pC, vC např. přibližnost integrace spočívá jen v zanedbání proměnnosti třecích ztrát podél charakteristiky

19. Metoda charakteristik NAP6 Řešení této soustavy vyjádříme v explicitním tvaru a můžeme ho ihned využít pro numerické řešení rychlostí a tlaků na časové hladině C ze známých hodnot v bodech A,B na staré časové hladině.

20. Metoda charakteristik NAP6 Alternativní (a možná názornější) formulace je analogií postupu, který jsme použili při řešení soustavy diferenciálních rovnic výměníků tepla, viz Whitham: Linear and nonlinear waves, Wiley 1974. Soustavu dvou diferenciálních rovnic pro tlak a rychlost zapíšeme v maticovém tvaru a použijeme transformaci vektoru W takovou, že se ve výsledných rovnicích objeví vždy jen jedna neznámá pokud je [[Q]] maticí vlastních vektorů matice [[A]] bude tato matice diagonální Pro náš konkrétní případ jsou vlastní vektory a vlastní hodnoty matice [[A]] funkcí Z1 získáme integrací rovnice podél charakteristikydx/dt=-a podél charakteristikydx/dt=a funkcí Z2 získáme integrací rovnice

21. MOC hydraul.ráz MATLAB NAP6 Následující příklad ukazuje řešení konkrétního problému stlačitelného proudění v trubce, kde na vstupu je zadaný tlak (např. konstantní p0) a na výstupu je zadaný časový průběh průtoku, např. rychle (téměř skokově) klesající průtok, což je situace odpovídající rychlému zavírání ventilu. Otázkou je především to, jaké zvýšení tlaku uzavírání ventilu způsobí (hydraulický ráz v potrubí).

22. MOC hydraul.ráz MATLAB C x2=-at x1=at h/a B A h NAP6 Vraťme se k předchozímu výsledku integrace dle dvojice charakteristik který popisuje rychlosti a tlaky ve vnitřních uzlech sítě. Do okrajových uzlů se lze dostat jen integrací dle jedné charakteristiky. Druhou potřebnou rovnicí je okrajová podmínka, např. předepsaný průběh tlaků na vstupu a předepsaná (třeba nulová) rychlost na výstupu pro případ hydraulického rázu. Zavírání ventilu lze v MATLABU popsat jako funkci času, např. function vrel=valve(t) if t<.1 vrel=1; else vrel=exp(-10*(t-.1)); end C v bodě C je tlak daný okrajovou podmínkou, dopočítá se jen rychlost integrací přes charakteristiku vlevo p=p0 v(t) B A

23. Průběhy tlaku (pro různé hustoty sítě), vždy je dobré takto testovat chybu aproximace n=301 n=101 čas (do 3 s) MOC hydraul.ráz MATLAB NAP6 Trubice L=1m, D=0.01 m, rychlost zvuku a=1 m/s, ustálená dopředná rychlost v=0.6325 m/s. l=1;d=0.01;rho=1000;f=0.1;a=1;p0=2e3; v0=(p0*2*d/(l*rho*f))^0.5 n=101;h=l/(n-1);v(1:n)=v0;p(1)=p0; for i=2:n p(i)=p(i-1)-f*rho*v0^2*h/(2*d); end dt=h/a;tmax=3;itmax=tmax/dt;fhr=f*h/(2*a*d); for it=1:itmax t=it*dt; for i=2:n-1 pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1); pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+fhr*(vb*abs(vb)-va*abs(va))); vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-fhr*(vb*abs(vb)+va*abs(va))); end pc(1)=p0; vb=v(2);pb=p(2); vc(1)=vb+(pc(1)-pb)/(a*rho)-fhr*vb*abs(vb); vc(n)=v0*valve(t); va=v(n-1);pa=p(n-1); pc(n)=pa-rho*a*(vc(n)-va)+f*h*rho/(2*d)*va*abs(va); vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n); p=pc;v=vc; end pmax=max(max(pres))/p0 x=linspace(0,1,n); time=linspace(0,tmax,itmax); contourf(x,time,pres,30) V čase t=1.1 už dorazila zpětná vlna na vstup V čase t=0.1 začíná zavírat výstupní ventil

24. váleček se po nárazu zastaví a zkrátí (a pak se zase odrazí, ale o to teď nejde) MOC hydraul.ráz MATLAB NAP6 Vliv rychlosti zavírání ventilu na maximální přetlak v uzávěru (není příliš výrazný) vrel=exp(-50*(t-.1)); velmi rychlé zavření v0 =0.6325 m/s a=1 m/s=1000 kg/m2 pmax =2180 Pa n=101 Přibližná analýza maximálního tlaku: kinetická energie „válečku“ tekutiny pohybujícího se rychlostí v se po nárazu promění v deformační energii vrel=exp(-10*(t-.1)); rychlé zavření cca 0.1 s v čas (do 6 s) L vrel=exp(-0.1*(t-.1)); pomalé zavírání cca 10s To odpovídá Žukovského rovnici, viz další folie rychlost zvuku Pro náš případ tedy pmax=632 Pa, což je jen cca 30% skutečné hodnoty. Ta analýza je jen orientační, nestlačuje se rovnoměrně celá trubice.

25. MOC hydraul.ráz MATLAB NAP6 Souvislost mezi rychlostí zvuku, rychlostí sloupce kapaliny a přetlakem p vyjadřuje několik století stará rovnice hydraulického rázu (Young 1808, v modernější podobě pak Žukovského rovnice 1898 pro ráz v elastické trubce) Rovnice je tak známá, že se už většinou ani neuvádí odvození. dv rychlost pístu arychlost čela rázové vlny Odvození vychází z hmotnostní a hybnostní bilance kontrolního objemu, který se pohybuje konstantní rychlostí rázové vlny a směrem vpravo kontrolní objem pohybující se rychlostí zvuku vpravo Předpoklad konstantní rychlosti až do místa rázové vlny bilance hmoty v tok hybnosti je hmotnostní průtok krát rychlost rychlost tekutiny vzhledem k rozhraní, které se pohybuje rychlostí zvuku a p, bilance hybnosti Předpoklad skokového poklesu hustoty i tlaku v zoně rázové vlny Po dosazení rovnice bilance hmotnosti plyne rovnice hydraulického rázu

26. MOC hydraul.ráz MATLAB NAP6 Tyto černé stránky můžete bez obav přeskočit. Jsou jen ukázkou toho, jak numerický model zrealističtit (výpočtem realistického součinitele třecích ztrát) a současně demonstrovat skutečnost, že pro numerické výpočty se MATLAB zase tak moc nehodí, je to totiž jen interpret a tudíž pomalý. Následující program byl beze změny přepsán do FORTRANU (místo MATLABovských cyklů for i=2:n-1, se napíše do i=2,n-1, místo rozhodovacích příkazů typu if re<2100se napíše if(re.lt.2100)then atd.). Výpočet následující úlohy zavírání ventilu trval na mém počítači 2s ve FORTRANu a 125 s v MATLABu (cca 60 krát pomalejší).

27. MOC hydraul.ráz MATLAB NAP6 Program uvažující závislost součinitele třecích ztrát na Re, a identifikující laminární/turbulentní režim toku. Viz funkce frict. Blasiův vztah pro třecí ztráty function f=frict(v,d,rho,mju) re=abs(v)*d*rho/mju; if re<2100 f=64/re; else f=0.316/re^0.25; end for it=1:itmax t=it*dt; for i=2:n-1 pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1); fa=frict(va,d,rho,mju); fb=frict(vb,d,rho,mju); ea=fa*h/(2*a*d)*va*abs(va); eb=fb*h/(2*a*d)*vb*abs(vb); pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+eb-ea); vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-eb-ea); end pc(1)=p0; vb=v(2);pb=p(2); fb=frict(vb,d,rho,mju); eb=fb*h/(2*a*d)*vb*abs(vb); vc(1)=vb+(pc(1)-pb)/(a*rho)-eb; vc(n)=v0*valve(t); va=v(n-1);pa=p(n-1); fa=frict(va,d,rho,mju); pc(n)=pa-rho*a*(vc(n)-va)+fa*h*rho/(2*d)*va*abs(va); vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n); p=pc;v=vc; end konstantní tlak na vstupu p0=1e3;l=1;d=0.01;rho=1000;mju=0.001;a=4; vlam=p0*d^2/(32*mju*l); re=vlam*d*rho/mju if re<2100 v0=vlam f=64/re; else v0=(p0*d^(5/4)/(0.158*mju^0.25*rho^0.75*l))^(4/7) re=v0*d*rho/mju; f=0.316/re^0.25 end n=401;h=l/(n-1);v(1:n)=v0;p(1)=p0; for i=2:n p(i)=p(i-1)-f*rho*v0^2*h/(2*d); end dt=h/a;tmax=3;itmax=tmax/dt; zadaná rychlost na výstupu

28. MOC hydraul.ráz MATLAB NAP6 Voda, a=4 m/s, L=1m, D=0.01 m, p0=1 kPa, v0=1.14 m/s, Re=31000 Vývoj rychlosti Rozložení tlaku průběhy tlaku na zavíracím ventilu To, že se pulzace zrychlily, je způsobeno vyšší zvolenou rychlostí zvuku (tužší elastická trubice) a=4 m/s. Průlet tlakové vlny pak odpovídá času 0.25 s.

29. MOC hydraul.ráz MATLAB NAP6 Téměř stejným způsobem lze řešit případ, kdy je na vstupu zadaný průtok (objemové čerpadlo) a na výstupu je konstantní tlak (např. nulový) function vrel=pump(t) vrel=0.5*sin(3*t); v(t) p=0 B A vc(1)=pump(it*dt);vb=v(2);pb=p(2); pc(1)=pb+rho*a*(vc(1)-vb)+f*h*rho/(2*d)*vb*abs(vb); va=v(n-1);pa=p(n-1); pc(n)=0; vc(n)=va-(pc(n)-pa)/(rho*a)-fhr*va*abs(va);

30. MOC elastická a tuhá trubice NAP6 Řešení, které je asi úplně špatně

31. MOC elastická a tuhá trubice NAP6 Uvažujme případ toku nestlačitelné tekutiny elastickou trubicí (např. cévou) na níž je napojena trubice (třeba se stejným průměrem), ale dokonale tuhá. Tato zdánlivě nevinná kombinace okrajových podmínek představuje pro numeriku problém. V tuhé trubce a pro nestlačitelnou kapalinu je rychlost zvuku nekonečná, charakteristiky jsou rovnoběžné s osou trubky a časový krok vychází nekonečně malý. Tuhá trubice. Rychlost zvuku je nekonečně velká. Elastická trubice. Rychlost zvuku a je dána její elasticitou. v(t) pe(t) L Le skoková změna impedance (odporu) vyvolá odraz vln, které interferují s vlnami, které postupují směrem toku vpravo

32. MOC elastická a tuhá trubice Le NAP6 Nekorektní (?) okrajová podmínka na konci (podmínka s první derivací) D C A e Vraťme se k původní Bernoulliho rovnici, která by měla platit i ve vyústění elastické trubice a odečtěme od ní Bernoulliho rovnici pro tuhou trubici

33. MOC elastická a tuhá trubice NAP6 Program v MATLABu le=0.001;de=0.1;pe=0;l=1;d=0.01;rho=1000;f=0.1;a=1;p0=0;v0=0; n=401;h=l/(n-1);v(1:n)=v0;p(1)=p0; for i=2:n p(i)=p(i-1)-f*rho*v0^2*h/(2*d); end dt=h/a;tmax=3;itmax=tmax/dt; fhr=f*h/(2*a*d); for it=1:itmax t=it*dt; for i=2:n-1 pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1); pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+fhr*(vb*abs(vb)-va*abs(va))); vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-fhr*(vb*abs(vb)+va*abs(va))); end vc(1)=pump(it*dt);vb=v(2);pb=p(2); pc(1)=pb+rho*a*(vc(1)-vb)+f*h*rho/(2*d)*vb*abs(vb); va=v(n-1);vb=v(n);pa=p(n-1); pc(n)=(pc(n-1)/h+pe/le-f*rho*vb*abs(vb)/(2*d*de)*(de-d))/(1/h+1/le); vc(n)=va+1/(rho*a)*(pa-pc(n)-f*h*rho/(2*d)*va*abs(va)); vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n); p=pc;v=vc; end pmax=max(max(pres)) x=linspace(0,1,n); time=linspace(0,tmax,itmax); contourf(x,time,pres,30)

34. MOC elastická a tuhá trubice function vrel=pump(t) vrel=0.5*sin(3*t); n=401,Le=2, pmax=803 n=401,Le=1, pmax=709 n=401,Le=0, pmax=669 n=401,Le=0.1, pmax=669 n=401,Le=0.5, pmax=669 NAP6 Tlakové profily v elastické trubici (a=1 m/s, L=1 m, D=De=0.01 m)

35. MOC elastická a tuhá trubice NAP6 Experimenty prováděné v naší laboratoři neoperovaly s harmonickým průběhem průtoku, ale jen s jeho dočasným zvýšením. I skutečná rychlost zvuku byla vyšší. A o to větší jsou problémy s numerickým řešením.

36. MOC elastická a tuhá trubice NAP6 Elastická trubice (AORTA, rychlost šíření pulzní vlny=rychlost zvuku cca 8 m/s) a na ni napojená tuhá trubice (plexisklo). Puls průtoku na vstupu generuje elektronický řízená SUPERPUMP (Varimex). Výstup do uklidňovací nádrže. Ukázka experimentů z naší laboratoře (pulzní vlna způsobuje deformační vlnu aorty, a ta je monitorována dvojicí vysokorychlostních kamer, DIC=Digital Image Correlation). Pro řešení byla použita metoda sítí Lax Wendroff a MUSCL (Monotone Upstream Scheme for Conservation Laws, programováno ve Fortranu). Problém je v tom, že metoda nefunguje, když je na elastickou trubici napojená dokonale tuhá trubice (v modelech bylo nutné modelovat tuto nástavnou trubici také jako elastickou, jen s vysokou tuhostí). Matematici říkají, že je to tím, že v této hyperbolické rovnici nelze použít okrajovou podmínku s derivací. Pokud mají pravdu, je následující řešení chybné…

37. MOC elastická a tuhá trubice NAP6 Lineární změna průtoku na vstupu (rampová funkce). Laminární i turbulentní režim. le=0.1;de=0.015;pe=0;mju=0.001;l=.18;d=0.015;rho=1000;a=8; v0=pump(0);ve=v0*(d/de)^2;f=frict(v0,d,rho,mju);fe=frict(ve,de,rho,mju);re=v0*d*rho/mju dpe=0.5*fe*le/de*rho*ve^2;dp=0.5*f*l/d*rho*v0^2; n=401;h=l/(n-1);v(1:n)=v0;p(1)=pe+dp+dpe; for i=2:n p(i)=p(i-1)-f*rho*v0^2*h/(2*d); end dt=h/a;tmax=3;itmax=tmax/dt; for it=1:itmax t=it*dt; for i=2:n-1 pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1); fa=frict(va,d,rho,mju);fb=frict(vb,d,rho,mju); ea=fa*h/(2*a*d)*va*abs(va); eb=fb*h/(2*a*d)*vb*abs(vb); pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+eb-ea); vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-eb-ea); end vc(1)=pump(it*dt);vb=v(2);pb=p(2); fb=frict(vb,d,rho,mju); pc(1)=pb+rho*a*(vc(1)-vb)+fb*h*rho/(2*d)*vb*abs(vb); va=v(n-1);vb=v(n);pa=p(n-1); fb=frict(vb,d,rho,mju); fa=frict(va,d,rho,mju); pc(n)=(pc(n-1)/h+pe/le-fb*rho*vb*abs(vb)/(2*d*de)*(de-d))/(1/h+1/le); vc(n)=va+1/(rho*a)*(pa-pc(n)-fa*h*rho/(2*d)*va*abs(va)); vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n); p=pc;v=vc; end function vrel=pump(t) if t<.1 vrel=0.1; elseif t<1 vrel=t; else vrel=1; end function f=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 pokud je někde chyba, tak asi tady

38. MOC elastická a tuhá trubice Le=1 Re=5250 pmax = 7785 Pa v0=0.35 m/s vmax=0.7 m/s a=8 m/s D=0.015 m Le=1 m L=0.18 Le=0.1 Re=5250 pmax = 3848 Pa v [m/s] 0.7 0.35 0.06 t [s] 0.03 NAP6 Lineární změna průtoku na vstupu (rampová funkce). Laminární i turbulentní režim. 1 0.18 0.18 0.1

39. MOC elastická a tuhá trubice NAP6 Rychlosti a tlak na konci elastické trubice (snížení frekvence při prodloužení trubky, zdá se mi to logické, zvětšila se setrvačná hmota a nezměnila elasticita) (L=0.18, Le=0.01, pmax=2946 Pa)(L=0.18, Le=0.1, pmax=3843 Pa). 0.1 0.01 0.18 0.18

40. Co je třeba si zapamatovat NAP6 Přednáška byla věnována klasifikaci parciálních diferenciálních rovnic druhého řádu, se zvláštním zřetelem na hyberbolický typ. Zapamatujte si alespoň to, co je na následující folii

41. Co je třeba si zapamatovat NAP6 Typ rovnice určují pouze koeficienty druhých derivací ty určují i rovnice tzv. charakteristik Charakteristiky jsou reálné, když b2-4ac>0 pak je rovnice hyperbolická (vlnová rovnice, kmitání) Charakteristiky je jen jedna, když b2-4ac=0 pak je rovnice parabolická (třeba rovnice vedení tepla) Reálné charakteristiky neexistují, když b2-4ac<0 pak je rovnice eliptická (Poissonova rovnice, statika, pružnost)