NUMERICKÁ ANALÝZA PROCESů

1. NUMERICKÁ ANALÝZA PROCESů NAP5 Modely popisované obyčejnými diferenciálními rovnicemi – okrajové (ale důležité!) problémy. Metoda vážených residuí a metoda konečných prvků Výměníky tepla (využití analytického řešení) Potrubní sítě (tlaky a průtoky metodou konečných prvků) Táhla, nosníky a rotačně symetrické skořepiny Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. Modely ODE okrajový problém NAP5 Některé systémy jsou popisovány stejnými ODE jako dříve (tj. soustavou rovnic s prvními derivacemi), jenomže výchozí stav není definován,protože ne všechny podmínky jsou specifikovány pro stejnou hodnotu nezávisle proměnné (nehovoříme o počátečních, ale okrajových podmínkách). Zatímco počáteční problém se zpravidla týká nestacionárních dějů (nezávisle proměnnou je čas) je pro okrajový problém typické hledání ustáleného stavu (nezávisle proměnnou je prostorová souřadnice). Příklady: Teplotní profily podél trubek nebo deskových kanálů výměníků tepla Průběhy tlaků a průtoků v potrubní síti (ve stacionárním stavu) Deformace zatížené prutové konstrukce, nebo rotačních skořepin

3. ODE výměníky tepla NAP5

4. ODE výměníky tepla NAP5 Uvažujme nejjednodušší možný příklad výměníku tepla se dvěma paralelními proudy (souřadnice x je měřena od jednoho konce výměníku). Je to opakování příkladu, který byl řešený v kurzu „Heat processes“. Entalpická bilance při zanedbání axiálního vedení tepla SOUPROUD to je počáteční problém, řešitelný stejně jako dřív, třeba metodou Runge Kutta T1’ x T2’ L Tentýž případ, ale PROTIPROUD okrajový problém, chybí počáteční podmínka pro teplotu druhého proudu pro x=0 T1’ x T2’ (k je součinitel prostupu tepla vztažený na jednotku délky kanálu, Wi jsou tepelné kapacity proudů, uvažujeme je v obou případech kladné)

5. ODE výměníky tepla NAP5 Rovnice teplotních profilů jednotlivých proudů lze zapsat v maticovém tvaru speciální případ pro souproud a pouze dva kanály Když je matice [[A]] konstantní (když se nemění k) je systém lineární a soustava propojených ODE se dá převést na soustavu izolovaných ODE transformací kde [[U]] je matice vlastních vektorů matice [[A]] splňující podmínku aneb: když matici vynásobím vlastním vektorem, dostanu tentýž vlastní vektor (jen násobený konstantou) je diagonální matice vlastních hodnot1 2… (kolik rovnic tolik vlastních hodnot a vlastních vektorů). Výpočet matic [[]] a [[U]] je např. v MATLABu zajištěno jedinou funkcí, příkazem[U,L]=eig(A).

6. ODE výměníky tepla NAP5 Dosazením transformace do soustavy ODE Čtvercová matice. Když je [[U]] matice vlastních vektorů, bude diagonální. získáme soustavu rovnic pro transformované teploty(vektor [Z]) ODE pro Z jsou nezávislé, protože[[]] je diagonální matice Řešení ODE Návrat k původním teplotám koeficienty dj jsou dány okrajovými podmínkami

7. æ ö 0 0 æ ö 1 W ç ÷ ç ÷ 2 = L = 1 1 [[ U ]] [[ ]] ç ÷ ç ÷ - + 0 k ( ) - 1 W ç ÷ è ø W W 1 è ø 1 2 æ ö æ ö æ ö T ' 1 W d + - W T ' W T ' T ' T ' ç ÷ ç ÷ ç ÷ 1 2 1 = = = d 1 1 2 2 d 1 2 ç ÷ ç ÷ ç ÷ 1 2 - + + T ' 1 W d W W W W è ø è ø è ø 2 1 2 1 2 1 2 ODE výměníky tepla NAP5 Konkrétní případ dvoukanálového souproudého výměníku Vlastní problém pro matici [[A]] 2 x 2 se dá řešit analyticky:  Okrajové podmínky (v tomto případě spíš počáteční podmínky) pro x=0 a řešení soustavy . zdá se, že to vyšlo jinak? ne, jen se zaměnilo pořadí vlastních vektorů a ty se normovaly na jednotkovou Eukleidovskou normu. Na výsledek to nemá žádný vliv. Ověření v MATLABu (k=1, w1=1, w2=2) a=[-1 1;0.5 -0.5]; [u,l]=eig(a) u = -0.8944 -0.7071 0.4472 -0.7071 l = -1.5000 0 0 0

8. ODE výměníky tepla NAP5 Obecná metodika výpočtů teplotních profilů ve výměnících tepla a v sítích výměníků Výměníky tepla mají často více než jen 2 proudy (budeme uvažovat N proudů) a každý proud (kanál) je žádoucí rozdělit na subkanály (třeba proto, že jsou vstupní a výstupní hrdla jednotlivých proudů na různých místech, nebo že se podél kanálu mění součinitel prostupu teplači směr proudění). Uvažujme, že všech N proudů je rozděleno na M subkanálů (M>N), poloha vstupů těchto subkanálů je dána vektorem [x’] a poloha výstupůvektorem [x’’]. Subkanálu i (i=1,2,…,M) odpovídá teplota Ti(x), tepelná kapacita proudu Wi (konstantní), a součinitelé prostupu tepla kij se všemi ostatními subkanály (také konstantní). Tato data stačí k tomu, aby bylo možné napsat entalpické bilance všech M subkanálů a soustavu M diferenciálních rovnic Tato soustava se vyřeší dříve uvedeným postupem (nalezením vlastních hodnot a vlastních vektorů [[A]] a řešením separovaných rovnic).Výsledkem je vektor teplotních profilů v M subkanálech, vyjádřený M koeficienty zatím neznámého vektoru [d] Z této rovnice lze stanovit vstupní i výstupní teploty všech M subkanálů (dosazením souřadnic [x’],[x’’]). Některé vstupy/výstupy subkanálů jsou i vstupy/výstupy kanálů (tj.celého výměníku), některé vyjadřují jen vnitřní propojení. Tyto vazby (uspořádání) subkanálů vyjadřují binární matice [[G]]MxM jejíž řádky odpovídají vstupům a sloupce výstupům (Gij=1 znamená, že vstup i-tého subkanálu je výstupem j-tého subkanálu), matice vstupů výměníku [[G’]]MxN(G’ij=1 znamená, že vstup i-tého subkanálu je vstup j-tého proudu) a matice výstupů [[G’’]]NxM(G’’ij=1 znamená, že výstup i-tého proudu je i výstupem j-tého subkanálu).

9. doporučeno ke čtení ODE výměníky tepla NAP5 N koeficientů [d] se stanoví z okrajových podmínek (předepsané vstupní teploty v N proudech) a dále z požadavku, že výstupní teploty subkanálů jsou současně vstupní teploty subkanálů navazujících (to je definováno maticí [[G]]), což je chybějících M-N podmínek. Vstupní teploty všech subkanálů jsou tedy buď vstupní teploty celého výměníku nebo výstupní teploty navazujícího subkanálu, tedy Předchozí vztah je soustava M lineárních algebraických rovnic pro koeficienty [d]. Z nich lze pak např. vyjádřit výstupní teploty výměníku Xing Luo, Meiling Li, Wilfried Roetzel: A general solution for one-dimensional multistream heatexchangers and their networks. International Journal of Heat and Mass Transfer 45 (2002) 2695–2705

10. ODE výměníky tepla NAP5 Předchozí metoda je vhodná tehdy, když jsou ODE lineární a když tedy existuje analytické řešení. Pokud koeficienty prostupu tepla nejsou konstantní (ale jsou nezávislé na teplotě) je možné kanály vyměníku rozdělit na menší subkanály a každému přiřadit střední hodnotu k (viz předchozí text). Nu, a pokud by koeficienty prostupu na teplotách závisely, bylo by nutné celý proces iterativně opakovat (a v každé iteraci řešit vlastní problém pro matici [[A]] závisející na teplotě). Nebo je možné použít numerickou integraci (opakované řešení počátečního problému): Chybějící počáteční podmínky (pro x=0) se odhadnou a postupně zpřesňují opakovanou integrací (třeba metodou Runge Kutta) tak, aby byly co nejlépe splněny zbývající okrajové podmínky (pro x=L). Tento postup se nazývá metoda střelby a je to vlastně optimalizační problém, jehož parametry jsou chybějící počáteční podmínky optimalizované (v MATLABu třeba metodou fminsearch) tak, aby byl minimalizována norma rozdílu mezi okrajovými podmínkami a predikcí řešení pro x=L. Můžete zkusit i jiný postup spočívající v odhadu chybějících počátečních podmínek, provedení integrace od x=0 do x=L, a pak řešení v obráceném směru od x=L do x=0, ale s výchozími hodnotami (x=L) nahrazenými předepsanými okrajovými podmínkami. Celý postup se musí několikrát opakovat (cik cak metoda) a ne vždy konverguje.

11. ODE výměníky tepla NAP5 T1’=90 Ukázkový příklad dvoukanálového výměníku délky L=1, souproud i protiproud. Řešení metodou Runge Kutta v MATLABu (ode45) definice ODE pro protiproud se liší jen znaménkem druhé rovnice definice ODE pro souproud function dtdx=dts(x,t) dtdx=zeros(2,1); w1=1; w2=2; k=1; dtdx(1)=-k/w1*(t(1)-t(2)); dtdx(2)=k/w2*(t(1)-t(2)); function dtdx=dtp(x,t) dtdx=zeros(2,1); w1=1; w2=2; k=1; dtdx(1)=-k/w1*(t(1)-t(2)); dtdx(2)=-k/w2*(t(1)-t(2)); T2’=40 řešení pro souproud je počáteční problém >> sol=ode45(@dts,[0 1],[90 40]); >> plot(sol.x,sol.y) T1’ řešení pro protiproud metodou cik-cak s=ode45(@dtp,[0,1],[90,50]); for i=1:10 n=length(s.x); s=ode45(@dtp,[1,0],[s.y(1,n),40]); n=length(s.x); s=ode45(@dtp,[0,1],[90,s.y(2,n)]); end plot(s.x,s.y) T2’

12. ODE sušárny, extraktory… NAP5 Průtočné aparáty jako extraktory, reaktory nebo sušárny (souproudé, protiproudé) se řeší téměř stejně jako výměníky tepla. Jen k bilančním rovnicím přenosu tepla se přidávají analogické rovnice přenosu hmoty. Takže se opět používají metody numerické integrace (Runge Kutta), metoda střelby apod. Příklad: Bruce D.M., Giner S.A.: Mathematicall Modelling of Grain Drying in Counter-flow Beds. J.Agric.Engng.Res. (1993),pp.143-161 Protiproudá sušárna zrní. Řešení Euler, Runge Kutta, metoda střelby.

13. Potrubní sítě NAP5 Potrubní síť je tvořena přímými úseky, koleny a odbočkami. V širším slova smyslu do ní můžeme zahrnout i aparáty (čerpadla, akumulátory, ventily). Cílem je stanovení průběhu tlaků podél potrubí a rozdělení průtoků do odboček. Odlišný přístup vyžadují případy: • nestacionární tok (efekt zrychlení) • stlačitelné medium nebo stěny (ráz)

14. Potrubní sítě NAP5 Základním elementem je trubice (libovolného průřezu, který se může ve směru osy měnit). Charakteristikou proudění v určitém místě této trubice je tlak p(x) a charakteristická rychlost u. Poměrně obecně (i pro stlačitelné proudění) můžeme napsat Navierovu Stokesovu rovnici pro rychlost v ose Druhý člen pravé strany vyjadřuje třecí ztráty, které jsou funkcí průtoku, viskozity a geometrie kanálu. Tyto faktory zahrneme do koeficientu k: D ekvivalentní průměr, A je plocha průřezu, m hmotnostní průtok dArcy Weissbachův součinitel tření, v lamináru např. 64/Re Třecí ztráty Integrací této rovnice po délce trubky získáme Bernoulliho rovnici pro nestacionární proudění stlačitelné tekutiny (použijeme ji až později…).

15. Potrubní sítě MKP, MVR NAP5 Uvažujme speciální případ: stacionární průtok, konstantní průřez. Předchozí rovnici derivujme V této rovnici se vyskytuje pouze tlak a hmotnostní průtok (který je konstantní – rovnice zachování hmoty) byl eliminován. Levá strana vyjadřuje třecí ztráty, pravá strana vztlak (když je hustota závislá na teplotě a teplota podél kanálu se mění). Koeficient k je v laminárním režimu toku konstantní, v turbulentním režimu s rostoucím Reynoldsovým číslem klesá. Tuto obyčejnou diferenciální rovnici využijeme pro vysvětlení principů metody vážených residuí (MVR) a její speciální varianty, metody konečných prvků (MKP). Co je to residuum diferenciální rovnice? To co zbude na pravé straně, když všechny členy převedeme nalevo Kdyby bylo p(x) přesné řešení, byla by funkce res(x) nulová

16. Potrubní sítě MKP, MVR NAP5 Metoda Konečných Prvků Metoda Vážených Residuí Metoda vážených residuí vybírá z možných variant řešení p(x) takovou, která anuluje integrály residua násobeného zvolenou váhovou (testovací) funkcí w(x) Hledané řešení se aproximuje lineárním modelem, lineární kombinací bázových funkcí Nj(x) (mohou to být třeba po elementech definované polynomy) N neznámých parametrů pj vyžaduje alespoň N rovnic. Proto je třeba zvolit alespoň N různých váhových funkcí wi(x). Výsledkem je soustava N lineárních algebraických rovnic pro N neznámých

17. Metoda vážených residuí xi xi NAP5 • Váhové funkce wi(x,y,z) jsou navrhovány předem, víceméně nezávisle na hledaném řešení. Mohou být různého typu a tomu odpovídají různé numerické metody řešení • Spektrální methody (analytické váhové funkcewi(x)) • Metody konečných prvků (spojité váhové funkce) • Metody konečných objemů(nespojité, po částech konstantní funkce) • Kolokační metody (nulové reziduum v uzlech. Váhové funkce jsou delta funkce) (Jako bázové a váhové funkce se používají ortogonální polynomy nebo goniometrické funkce. Analytické váhové funkce používá např. i metoda hraničních prvků BEM.) (např. Galerkinova metoda používající stejné váhové a bázové funkce. MKP je nejčastější metoda řešení pružnostních problémů, asi znáte software ANSYS) (též metoda kontrolních objemů. Nejčastěji používaná metoda v mechanice tekutin, tj. proudění, rozšířený software FLUENT) xi (též metoda konečných diferencí nebo metoda sítí. Vhodné pro jednoduché geometrie, snadné programování.)

18. Konečné prvky a objemy wi(x) xi x wi(x) xi MKP váhové funkce jsou spojité a existuje jejich první derivace MKO váhové funkce nejsou spojité na hranicích kontrolních objemů x NAP5 Konstrukce (potrubní síť) je rozdělena na konečné elementy s uzlovými body na hranicích elementů (element je úsek potrubí mezi dvěma uzly). Každému uzlu (i) je přiřazena jedna váhová funkce wi(x), jedna bázová funkce Ni(x) a jedna počítaná hodnota (tlak pi). V metodě konečných prvků je definován průběh řešení uvnitř elementu interpolací uzlových hodnot na hranici. U konečných objemů definuje průběh řešení jen jediný vnitřní uzel, takže uvnitř konečného objemu to nemůže být nic jiného než konstanta. Pozn.: Kompartment modely jsou vlastně prototypem metody konečných objemů. Cílem bylo stanovit koncentrace a teploty, stejné v celém kompartmentu. U metody konečných elementů je cílem stanovit hodnoty na hranicích elementů, tj. parametry proudů, které kompartmenty propojují. V MKO je už z principu zachována bilance kontrolního objemu, v MKP to přesně platit nemusí (a to je výhoda MKO).

19. Konečné prvky a objemy NAP5 Terminologie FEM (Finite Element Method). Upřímně řečeno nevím, jaký je rozdíl mezi konečným prvkem a konečným elementem. Někdo pod pojmem konečný element chápe jen geometrii, třeba trojúhelník, a teprve, když se k němu přidruží bázové funkce nazve ho konečný prvek (nebo obráceně ☺). FVM (Finite Volume Method). Stejně tak nevím jaký je rozdíl mezi konečným a kontrolním objemem. Možná, že konečný objem je speciální případ kontrolního objemu, používaný v numerice a snaží se vystihnout optickou podobnost mezi konečnými prvky a konečnými objemy. FD (Finite Differences). A také nevím jaký je rozdíl mezi metodou sítí a metodou konečných diferencí.

20. Metoda konečných prvků NAP5 Bázové funkce Ni(x) se zpravidla konstruují tak, že jsou rovny 1 v uzlu i a ve všech ostatních uzlech (1,…,N) jsou nulové. Bázová funkce Ni(x) je tedy nenulová jen v těch elementech, které obsahují uzel číslo i. Jako váhové funkce lze použít bázové funkce (použijeme např. po elementech lineární průběhy tlaku i váhových funkcí). Ztotožnění váhových a bázových funkcí (wi(x)=Ni(x)) je přízračným rysem Galerkinovy metody. Integrand integrálu residua je součin váhové funkce a druhé derivace bázové funkce. Protože je v MKP váhová funkce spojitá lze ji derivovat a upravit integrál metodou per partes Kij bi Koeficienty matice soustavy algebraických rovnic Kij i vektoru pravé strany bi jsou integrály přes celou potrubní síť (L je souhrnná délka všech větví). Tyto integrály se počítají jako součet integrálů přes jednotlivé konečné elementy.

21. Metoda konečných prvků Nj Ni j i Le NAP5 V obecném konečném elementu jsou nenulové jen dvě bázové (a současně váhové) funkce odpovídající dvěma uzlům i,j. Příspěvky jednoho elementu se tedy uplatní jen u 4 prvků matice soustavy (Kii, Kjj Kij Kji). Protože jsou derivace Ni(x) v elementu konstantní je výpočet integrálu příspěvku elementu snadný Sčítání matic elementů do globální matice soustavy (a vektoru pravé strany) se provádí algoritmem sestavení. V cyklu přes všech Ne elementů se volá procedura výpočtu lokální matice pro daný průtok a geometrii a tato matice se přičte k matici globální. Korespondence mezi lokálními indexy (1,2) a globálními indexy musí být definována v matici konektivity c, jejíž řádky odpovídají elementům a ve dvou sloupcích jsou indexy i,j uzlů každého elementu. Tím je vlastně definována topologie sítě. Podprogram výpočtu Ke K=zeros(n,n); b=zeros(n,1); for e=1:ne [Ke,be]=local(e); for i=1:2 ig=c(e,i); b(ig)=b(ig)+be(i); for j=1:2 jg=c(e,j); K(ig,jg)=K(ig,jg)+Ke(i,j); end end end Matice konektivity Pozn.: Korespondence mezi lokálními a globálními indexy bývá komplikovanější, protože uzlovým bodům většinou odpovídá více než jeden parametr (třeba posuvy) a elementy mívají různé počty uzlů (trojúhelníky,…). Matice konektivity pak má proměnnou délku řádků.

22. Metoda konečných prvků NAP5 Fyzikální interpretace: Když vynásobíme lokální matici jednoho elementu [[K]]evektorem vypočtených tlaků dostaneme hmotnostní průtok elementem (e), přesněji průtok do uzlu i a do uzlu j Sestavení lokálních matic tedy vyjadřuje požadavek na to, aby byl součet všech orientovaných průtoků v každém uzlu nulový. • Vygenerovaná matice soustavy [[K]] je v této fázi řešení ještě singulární, protože nejsou zahrnuty okrajové podmínky (a rozložení tlaků není jednoznačné). Okrajové podmínky je třeba předepsat ve všech koncových uzlech sítě: • Koncové tlaky (příslušný řádek matice soustavy se nahradí jedničkou na diagonále a na pravou stranu se dosadí předepsaný tlak). • Pokud je v koncovém uzlu předepsán průtok, dosadí se hodnota průtoku přímo jako pravá strana bi.

23. Metoda konečných prvků Vroubená matice (sky line) Pásová matice NAP5 Pro řešení soustav lineárních algebraických rovnic se používají buď finitní přímé metody (příkladem je Gaussova eliminační metoda, frontální metoda u MKP) nebo metody iterační (např. Gaussova Seidelova nadrelaxační metoda, nebo metoda sdružených gradientů). Matice soustavy je řídká (v jednom řádku, který odpovídá uzlu v němž se spojují jen dva elementy, jsou jen 3 nenulové prvky) a proto se používají i speciální metody úsporného ukládání řídkých matic, např. Nejde jen o úsporu paměti, ale o rychlost výpočtu. Např. pro Gaussovu eliminaci s plnou a nesymetrickou maticí soustavy N x N je třeba N3/3 operací typu násobení. U pásové matice se počet operací sníží na N.Nb2/2, kde Nb je šířka pásu. U relativně malých soustav (méně rozsáhlých potrubních sítí) stačí standardní Gaussova eliminace a v MATLABu napsat řešení soustavy [[K]][p]=[b] příkazem p=K\b

24. Metoda konečných prvků NAP5 Matice soustavy rovnic pro uzlové hodnoty tlaků je konstantní pouze při laminárním toku. V obecném případě turbulentního proudění, nebo proudění nenewtonských kapalin je funkcí Reynoldsova čísla, respektive počítané tlakové ztráty (k(p)). Soustava algebraických rovnic je pak nelineární a musí se řešit iteračně kde k je index iterace (každá iterace vyžaduje sestavení nové matice [[K]]). Pro řešení soustavy nelineárních rovnic se místo výše uvedené metody substitucí někdy používá Newtonova Raphsonova metoda, která je přímočarým zobecněním Newtonovy metody tečen hledání kořene transcendentní rovnice. Metoda je založena na Taylorově rozvoji rezidua, což vede opět na soustavu lineárních rovnic pro k-tou iteraci tlaku, navíc je však třeba počítat i matici derivací průtokových součinitelů vzhledem k tlaku. residuum  residuum / p

25. Metoda konečných prvků NAP5 Historická poznámka: Galerkinova metoda a MKP V ruské odborné literatuře se místo Galerkinovy metody hovoří o Bubnovově metodě. s poukazem na jeho práci publikovanou v roce 1914 (Bubnov I.G.: Stroitelnaja mechanika korablja. II. SPB-Tipografija morskogo ministerstva – Bubnov byl konstruktér válečných lodí). Galerkin publikoval své práce až v třicátých letech Galerkin B.G. K voprosu ob issledovajii naprjazenij i deformacij v uprugom izotropnom tele, Doklady Akademii Nauk SSSR, 1930, Ser.A, No.14, s.353-358. Metoda konečných prvků jako specifická aplikace Galerkinovy metody však vzniká až mnohem později, v souvislosti s analýzou havárií britského proudového letadla Comet. Turner M.J. et al: Stiffness and deflection analysis of complex structures. J. of Aeronautical Sciences, 23, pp.805-823 (1956). Historická poznámka: FEMINA FEMINA je MKP program orientovaný především na 1D elementy (program byl vyvinut na ústavu procesní a zpracovatelské techniky a je stále dostupný na webových stránkách). Používá výše uvedený algoritmus Galerkinovy metody pro stanovení rozložení tlaků v potrubní síti s elementy typu trubka, ventil, výměník tepla. Ambicí programu byla integrace výpočtu proudění (se zvláštním zřetelem na nenewtonské, např. tixotropní kapaliny), přenosu tepla i hmoty (axiální disperze, reakce) a současně i pevnostní výpočty. Nevýhodou programu je omezení na kvazistacionární problémy (neumožňuje řešit třeba problém rázu v potrubí) a především to, že už není dále vyvíjen (na víceprocesorových počítačích někdy selhává).

26. Potrubní sítě - model cirkulace NAP5 Tok kapaliny v umělém cirkulačním obvodu krve (fyzikální model v laboratoři kardiovaskulární mechaniky). Ukázka toho, že MKP se dá formulovat i jinak než Galerkinovsky (není použita metoda vážených residuí, tudíž žádné integrály, žádné algoritmy sestavení, jen opakování Bernoulliho rovnice a rovnice kontinuity). V předchozím řešení Galerkinovou metodou byla v každém uzlu jediná neznámá (tlak p) a průtoky Q v konečných elementech se počítaly až ex post.

27. Potrubní sítě - uzly NAP5 Parametry uzlů jsou u této varianty tlak p i objemový průtok Q. To ale znamená, že uzly nemohou být přímo v místě větvení, kde není průtok definován! Y [m] Počet neznámých 2N-Nb,kde N je počet uzlů a Nb je počet koncových uzlů (v každém je předepsaný tlak nebo průtok jako okrajová podmínka) 19 18 16 28 15 17 27 14 22 20 29 7 9 12 24 25 21 6 8 26 5 10 uvědomte si, že v určitém místě potrubní sítě není možné už z principu věci zadávat současně tlak a současně průtok. To platí i pro předchozí Galerkinovu metodu řešení 11 13 4 3 23 Předepsaný tlak p, počítaný průtok Q n 2 Počítaný tlak i průtok m Uzel pro definici geometrie větvení k nesouvisí s uzlovými parametry, jeho souřadnice jen přesněji definují geometrii větvení a umožní zpřesnit výpočet ztrátových koeficientů. 1 X [m] -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

28. x2,y2,z2 p2,Q2 x3,y3,z3 p3,Q3 x2,y2,z2 p2,Q2 d d1 x4,y4,z4 x1,y1,z1 p1,Q1 x1,y1,z1 p1,Q1 V1(2,3,4,23,d1,d2,d2) P1(1,2,d1) Potrubní sítě - elementy NAP5 Dva typy konečných elementů: dvouzlový element trubka generuje dvě rovnice (Bernoulliho rovnice proudnice mezi uzly 1,2 a rovnice kontinuity) a tříuzlový element větvení tři rovnice (Bernoulliho rovnice proudnic 1,2 a 1,3 plus rovnice kontinuity). Rovnice soustavy nevznikají sčítáním příspěvků elementů, elementy je generují přímo ve finálním tvaru. Počet neznámých parametrů souhlasí s počtem takto generovaných rovnic: 2N-Nb = 2Mp+3Mv počet elementů větvení počet elementů trubka počet uzlů počet koncových uzlů

29. d L Q2 Q3 d2 L2 d1 L1 Q1 Potrubní sítě - Bernoulli NAP5 Základem předchozího řešení jsou Bernoulliho rovnice, vzniklé integrací úvodní pohybové rovnice po proudnici mezi uzly 1 a 2 Poslední člen ez12 je ztrátová energie [J/kg] závislá na průtoku Q, zdánlivé viskozitě, místních ztrátách, viz E.Fried, E.I.Idelcik: Flow resistance, a design guide for engineers, E.I.Idelcik: Handbook of hydraulic resistances. to byste měli znát z hydromechanických procesů, Hagenbachův faktor, výpočet lokálních ztrát, turbulentní viskozita

30. Potrubní sítě MATLAB NAP5 MATLAB vzorové řešení: mainpq.m Soubory vstupních dat: xyz.txt x y z (souřadnice uzlů) cb.txt i pq (okrajové podmínky i-uzel, +tlak, -průtok jako parametry pq) cp.txt i j d (matice konektivity pro trubky, průměr trubky) cv.txt i j k l d1 d2 d3(matice konektivity pro větvení, průměry trubek, úhel větvení)

31. Táhla, nosníky, skořepiny NAP5 V této kapitole uvidíte, že metody výpočtu tlaků a průtoků v potrubní síti mají mnoho společného s výpočtem sil a posuvů příhradových konstrukcí.

32. Táhla, nosníky, skořepiny MKP NAP5 Jaký je rozdíl mezi táhlem a nosníkem? Obé jsou tyčky, jenže táhla v příhradové konstrukci jsou spojena klouby (na táhlo nepůsobí žádné momenty), zatímco nosníky jsou ve spojích svařeny a deformují se tudíž i působením ohybových a torzních momentů. Stejné je to s rozdílem mezi skořepinou a membránou. Membrána je analogií táhla (nepůsobí v ní žádné ohybové momenty, ve skořepině ano). Inženýrská teorie nosníků je založena na předpokladu, že průřez nosníku bude i v deformovaném stavu rovinný. Nejčastěji se přijímá i předpoklad, že průřez nosníku zůstane i v zatíženém stavu kolmý na (deformovanou) střednici nosníku. To je Bernoulliho nosník. Timošenkův nosník uvažuje i to, že rovina průřezu se vlivem smykových sil natočí i vzhledem ke kolmici osy nosníku. Bernoulliho teorie je adekvátním popisem deformace štíhlých, zatímco Timošenkova teorie spíše „tlustých“ nosníků. Stejná analogie platí u skořepin. Kirchhoffova teorie desek odpovídá Bernoulliho nosníku (průřez skořepiny zůstává kolmý k deformované střednici), zatímco Mindlinova skořepina (Mindlinova deska) odpovídá Timošenkovu nosníku.

33. Táhla, nosníky, skořepiny MKP NAP5 Všechny následující příklady mají něco společného: Řešení bude založeno na deformační metodě, kdy výsledkem řešení soustavy rovnic jsou posuvy a natočení uzlů (výpočty napětí jsou dopočítány z posuvů až ex post, tento postprocessing nechávám na vás). Deformační metoda odpovídá Lagrangeovu principu minima celkové potenciální energie. Deformační metoda není jediná možná a jako primární počítanou veličinu bylo možné vzít průběhy napětí nebo momentů, což by odpovídalo Castiglianovu principu doplňkové energie (ale není to běžné a ani výhodné). Základním výsledkem každého příkladu jsou matice elementů, implementované jako MATLABovské funkce, které si můžete zkopírovat a použít pro řešení konkrétních problémů

34. Táhla MKP Kij bi fx(x), ux(x) NAP5 Příhradová konstrukce složená z táhel je z hlediska řešení totožná s metodou výpočtu tlaků v potrubní síti. Místo tlaků je posunutí uzlů táhla, místo součinitele třecího odporu je koeficient tuhosti (průřez táhla * modul pružnosti) a vztlakový člen může být nahrazen spojitým vnějším zatížením táhla fx Metoda vážených residuí tedy vede na stejnou soustavu rovnic pro posuvy uzlů Poznámka: FUNKCIONÁLOVé by se k témuž výsledku dostali asi jednodušeji, spočítali by totiž celkovou potenciální energii táhla a hledali její minimum. a i matice tuhosti jednoho elementu je totožná

35. Táhla y u uy ux  x NAP5 Tímto způsobem bychom zatím mohli řešit jen triviální problém natahování seriově řazených táhel (pružinek) v jediném směru x. Posunutí uzlů prostorově uspořádaných táhel je ovšem vektor se složkami ux,uy (u rovinné), uz (trojrozměrné) konstrukce. Lokální matice tedy musí být transformovány do globálního souřadného systému, který je společný pro všechny elementy. Jestliže má například táhlo lokální matici tuhosti s rozměrem 2 x 2, je třeba ji transformovat na matici 4 x 4 pro rovinné konstrukce (v každém uzlu jsou dvě složky posunutí ux, uy). Transformace z lokálního do globálního souřadného systému se týká jen rotace (jediná rotace prvku u rovinné konstrukce, resp. až tři rotace u prostorové soustavy). Lokální souřadný systém (osa ) Globální systém x,y (c=cos, s=sin)

36. Táhla NAP5 Rovnice táhla v lokálním souřadném systému (kde je jediná prostorová souřadnice ve směru osy) se po dosazení transformačních vztahů změní na Pronásobením matic získáme finální podobu matice tuhosti táhla pro souřadný systém x,y >> t=[c^2 c*s;c*s s^2]; >> ke=[t -t;-t t] všímněte si jak elegantně se v MATLABu dá napsat matice tuhosti elementu

37. Táhla 5 y e5 e6 3 e4 4 e3 e2 e1 x 2 1 to je celý program, včetně dat, sestavení i řešení NAP5 Příklad x=[0 1 0 1 2]; y=[0 0 1 1 2]; con=[1 3; 2 3; 2 4; 3 4; 3 5; 4 5]; a=[4e-4 4e-4 4e-4 4e-4 4e-4 4e-4]; nu=length(x) ne=length(a) n=2*nu; k=zeros(n,n);b=zeros(n,1); for e=1:ne [kl,ig]=kloc(e,x,y,con,a); k(ig(1:4),ig(1:4))=k(ig(1:4),ig(1:4))+kl(1:4,1:4); end for i=1:4 k(i,:)=0;k(i,i)=1; end b(n)=-100; p=k\b; jediným příkazem se dá příčíst lokalní matice 4x4 ke globální matici tuhosti nulové posuvy v uzlech 1,2 a zatížení -100N v uzlu 5 (ošklivé řešení, platí jen pro danou topologii) vektor korespondence lokálních a globálních indexů matice konektivity Gaussova eliminace. Vektor p jsou výsledná posunutí function [kl,ig]=kloc(ie,x,y,con,a) E=200e9; ae=a(ie); i1=con(ie,1);i2=con(ie,2); ig=[2*i1-1 2*i1 2*i2-1 2*i2]; le=((x(i1)-x(i2))^2+(y(i1)-y(i2))^2)^0.5; c=(x(i2)-x(i1))/le;s=(y(i2)-y(i1))/le; td=E*ae/le*[c^2 c*s;c*s s^2]; kl=[td -td;-td td]; výpočet matice tuhosti elementu

38. Nosníky (trubky) Fx z x fz NAP5 Průhybová rovnice momentové rovnováhy nosníku je diferenciální rovnice čtvrtého řádu kde I je moment setrvačnosti průřezu, Fx je axiálně působící síla, k je součinitel tuhosti elastického podkladu a fz(x) spojité příčné zatížení. Odpovídající integrál residua je Dvojím použitím integrace per partes snížíme řád derivací přičemž byly zatím vynechány členy na hranici (koncích nosníku), které vznikají integraci per partes – to je konzistentní s předpokladem, že na koncích platí buď silné okrajové podmínky (vetknutí), nebo je okraj nezatížený.

39. Nosníky (trubky) b  x M(x) M(x) fz NAP5 Pokud vás nezajímají detaily, tak tyto šedivé stránky přeskočte Uvedená diferenciální rovnice průhybové čáry je založena na Bernoulliho modelu. Tedy: Průřez nosníku zůstává kolmý ke střednici, osová napětí jsou na střednici nulová a po průřezu se mění lineárně. Tomu odpovídá ohybový moment M a diferenciální rovnice kde druhá derivace průhybové čáry je mírou zakřivení (1/d2u/dx2 je poloměr kruhového oblouku). Ohybový moment může být způsoben nejenom příčnou silou (třeba reakcí M=Fz(b-x)), ale i axiální silou Fx Tato diferenciální rovnice je základem pro vyšetřování stability nosníků zatížených axiální silou (Eulerova stabilita). Ohybový moment M(x) je ovšem vyvolán především silami, které působí kolmo (ve směru z). U osamělých sil je jejich moment pouze lineární funkcí souřadnice, takže druhá derivace momentu je nulová a v diferenciální rovnici průhybové čáry se tento příspěvek vůbec neobjevil (zadaný okrajový moment nebo jeho první derivace se ovšem promítnou do vektoru zatížení, viz integrace per partes, popsaná dále). Spojité zatížení tlakem nebo reakcí pružného podepření nosníku ve směru z ale vyvolá ohybový moment, jehož druhá derivace nulová není Tento přechod úplně triviální není, jde o derivaci integrálu s proměnnou horní mezí i integrandem

40. Nosníky (trubky) Fz(b) M(b) M(b) Fz(a) NAP5 Vraťme se ještě jednou k aplikaci integrace per partes na rovnici vážených reziduí (s obecnou váhovou funkcí w(x)) Dvakrát opakovaná integrace na první člen se čtvrtou derivací a jedna integrace per partes pro druhý člen dává Kromě integrálu se objeví nové členy příspěvků okrajových podmínek na koncích nosníku. Vyjádříme je prostřednictvím momentů a osamělých sil (na základě právě odvozených relací mezi momenty a průhybovou čarou) V předchozím odvozování Galerkinovy metody jsme tyto členy zanedbali. To odpovídá situaci, kdy na konci nosníku jsou nulové momenty M i osamělé síly Fz nebo jsou tam fixovány průhyby a první derivace průhybu jako okrajové podmínky vetknutí či kloubového uchycení (tzv. silné okrajové podmínky). V uzlech, kde jsou předepsány silné okrajové podmínky totiž neurčují průhyby a jejich natočení diferenciální rovnicea váhové funkce w i jejich první derivace dw/dx mohou být nulové. Zmíněný člen je nenulový jen tehdy, když je na volném konci předepsáno zatížení momentem M nebo silou Fz – v tomto případě tvoří pravou stranu řešené soustavy rovnic.

41. Nosníky (trubky) NAP5 Silné a slabé okrajové podmínky Výsledkem sestavení matic tuhosti elementů je singulární matice a soustava rovnic není jednoznačně řešitelná. Je třeba odstranit stupně volnosti pohybu celé soustavy jako tuhého tělesa (posuvy a natočení celku) dodáním silných okrajových podmínek, což jsou předepsané posuvy/ natočení v určitých uzlech (nebo tlaky či teploty u potrubních sítích). Implementace silných okrajových podmínek je snadná: řádek matice tuhosti odpovídající silné okrajové podmínce se vynuluje, na diagonálu se dosadí jednička a do vektoru pravé strany předepsané posunutí nebo natočení (nebo tlak či teplota). Každé přibližné řešení pak alespoň ve fixovaných uzlech bude přesně splňovat předepsané posuvy a natočení. Je ale otázkou zda tyto předepsané podmínky bude splňovat i limita numerických řešení pro zjemňující se síť elementů. Ne vždy, alespoň ne tehdy, když diferenciální rovnice prostě neumožňuje předepsat některý parametr jako silnou okrajovou podmínku. Příkladem je rovnice nosníku zapsaná s druhou derivací průhybovéčáry M Vetknutý nosník zatížený konstantním momentem M, použity lineární bázové funkce. V uzlech je jediný parametr, posunutí uz. Nulové natočení ve vetknutí (x=0) ani nejde předepsat. Tatáž úloha,ale kubické bázové funkce. V každém uzlu je dvojice parametrů, posunutí a natočení. Nulové natočení lze předepsat, ale řešení s velkým počtem elementů tu podmínku vlastně ignoruje. přesné řešení U řešení uvedených na obrázcích je tento člen zanedbán, předpokládáme, že platí tzv. přirozené okrajové podmínky. Není to příčina selhání? 20 elementů 1 element 2 elementy 3 elementy 3 elementy 20 elementů Přesné řešení

42. Nosníky (trubky) Kubické bázové a váhové funkce Libovolný počet elementů NAP5 znamená, že předpokládáme nulové duz/dx v koncovém bodě L, a rovnice Zanedbání členu metody vážených residuí se ji opravdu snaží splnit (je to patrné z předchozích obrázků, kde aproximace průhybové čáry lineárními i kubickými polynomy tuto nechtěnou okrajovou podmínku respektují). Když tento okrajový člen zahrneme do matice tuhosti, výsledek se zlepší (u kubických polynomů bude dokonce řešení přesné), ale okrajovou podmínku nulového natočení v místě vetknutí u lineárních polynomů stejně nesplníme. Je to dáno tím, že u diferenciální rovnice druhého řádu lze jako silnou okrajovou podmínku uplatnit pouze posunutí. První derivace posunutí (natočení) můžeme jako silnou okrajovou podmínku předepsat jen u rovnic čtvrtého řádu. Řešení uvažující člen Lineární bázové a váhové funkce 20 elementů 3 elementy Okrajové podmínky, které a priori zajišťuje každá aproximace řešení se nazývají silné, zatímco okrajové podmínky přenesené do integrální formule MWR aplikací integrace per partes jsou slabé. Pro diferenciální rovnice 2s-tého řádu jsou slabé okrajové podmínky vyjádřeny minimálně s-tými normálovými derivacemi, okrajové podmínky s nižšími derivacemi jsou silné. U diferenciálních rovnic druhého řádu jsou silné okrajové podmínky předepsané hodnoty a slabé podmínky zahrnují první derivace, zatímco u rovnice čtvrtého řádu (třeba ohyb desek) jsou silné okrajové podmínky hodnoty i první derivace (posuvy a natočení, tedy kinematické podmínky), slabé podmínky se týkají až druhých nebo třetích derivací (zatěžující podmínky, druhá derivace odpovídá momentům M=EId2u/dx2, třetí derivace posouvajícím silám F=-EI d3u/dx3). Splnění slabých okrajových podmínek je věcí správné volby integrální formule u metody vážených residuí a nelze si je vynutit fixováním uzlových parametrů.

43. Nosníky (trubky) Fx z Fz x fz M Spojité zatížení Matice tuhosti nosníku [[K]] Osamělé síly a momenty NAP5 Finální formulace Galerkinovy metody vážených residuí pro Bernoulliovský nosník na pružném podkladě Pozn.: Nezapomínejte, že je to stále jen silně zjednodušený model nosníku, fungující pouze pro malé deformace, štíhlé nosníky (neuvažuje smykové síly jako Timošenkův model) a dokonce ani nemá axiální tuhost. Model kalkuluje jen s ohybovou tuhostí!

44. Nosníky (trubky) Bázová funkce mající ve všech uzlech nulové první derivace a v jednom uzlu se rovná 1 Bázová funkce rovná 0 ve všech uzlech a v jednom uzlu má první derivaci rovnou 1 i i+1 i i+1 NAP5 Podstatný rozdíl oproti problému s táhly spočívá v tom, že v integrandu jsou nyní druhé derivace bázových funkcí a není tedy možné použít například lineární aproximační polynomy – bázové funkce je třeba konstruovat tak, aby měly spojité první derivace, používají se tzv.Hermiteovy kubické polynomy Rovnice průhybové čáry v nosníku 1-2 Výsledkem dosazení Hermiteových polynomů do integrálů je opět soustava lineárních algebraických rovnic, jejíž matice je součtem tří matic [[KM]] je maticí ohybové a [[KF]] geometrické tuhosti a [[M]] je  maticí hmot.

45. Nosníky (trubky) NAP5 Pro matice jednoho elementu o délce L získáme integrací tyto výsledky zkuste prosím prointegrovat alespoň jeden prvek některé matice uzlové parametry(příčný posun a natočení neutrální osy) Matice tuhosti Matice geometrické tuhosti Matice hmot

46. Nosníky (trubky) MATLAB uzlové parametry NAP5 Procedura generování lokálních matic elementu nosník(ei-momenty setrvačnosti,k-tuhost podkladu) function [kstif,kgeom,mass,ig]=kloc(ie,x,y,con,ei,k) EI=ei(ie);E=200e9; i1=con(ie,1);i2=con(ie,2); ig=[2*i1-1 2*i1 2*i2-1 2*i2]; l=((x(i1)-x(i2))^2+(y(i1)-y(i2))^2)^0.5; kstif=EI*E/l^3*[12 6*l -12 6*l;6*l 4*l^2 -6*l 2*l^2; -12 -6*l 12 -6*l;6*l 2*l^2 -6*l 4*l^2]; kgeom=1/(30*l)*[36 3*l -36 3*l;3*l 4*l^2 -3*l -l^2; -36 -3*l 36 -3*l;3*l -l^2 -3*l 4*l^2]; mass=k*l/420*[156 22*l 54 -13*l;22*l 4*l^2 13*l -3*l^2; 54 13*l 156 -22*l;-13*l -3*l^2 -22*l 4*l^2]; Tyto matice můžeme použít pro sestavení několika elementů, ale orientovaných jen ve směru osy x (matici hmot využijeme při analýze kmitání nosníků, matici geometrické tuhosti pro vyšetřování stability). V každém uzlu jsou dva počítané parametry, průhyb uz a natočení průhybu . Nemůžeme však použít transformaci souřadnic pro natočené nosníky v příhradové konstrukci stejně jako u táhel, protože tento nosníkový element má nulovou tuhost ve směru x. Doplnění matice tuhosti nosníku o tuhost prvku typu táhlo je uvedeno v následujícím příkladu.

47. Příhradová konstrukce 1/3 uy z  uy ux z ux NAP5 Procedura generování matice tuhosti ocelového elementu nosník + táhlo function [kstif,ig]=klocp(ie,x,y,con,ei,ea) EI=ei(ie);EA=ea(ie); E=200e9; i1=con(ie,1);i2=con(ie,2); ig=[3*i1-2 3*i1-1 3*i1 3*i2-2 3*i2-1 3*i2]; l=((x(i1)-x(i2))^2+(y(i1)-y(i2))^2)^0.5; c=(x(i2)-x(i1))/l;s=(y(i2)-y(i1))/l; r=[c s 0;-s c 0;0 0 1];z=zeros(3,3); Q=[r z;z r]; km=EI*E/l^3*[EA*l^2/EI 0 0 -EA*l^2/EI 0 0; 0 12 6*l 0 -12 6*l; 0 6*l 4*l^2 0 -6*l 2*l^2; -EA*l^2/EI 0 0 EA*l^2/EI 0 0; 0 -12 -6*l 0 12 -6*l; 0 6*l 2*l^2 0 -6*l 4*l^2]; kstif=Q'*km*Q; to je jen předchozí matice [[KM]] rozšířená o první a čtvrtý řádek a sloupec, kam je dosazena matice tuhosti táhla. transformace matice tuhosti do globálního s.s. Transformace matice tuhosti z lokálního do globálního systému [[K]]=[[Q]]’[[Ke]][[Q]] zobecněná posunutí v uzlu i (výsledek výpočtu) zobecněná zatížení v uzlu i (pravá strana soustavy)

48. Příhradová konstrukce 2/3 y e2 3 2 e1 x 1 NAP5 Souřadnice uzlů, matice konektivity nosníků, sestavení, zadání okrajových podmínek a řešení soustavy je prakticky stejné jako u soustavy táhel. Kvadratický moment průřezu pro trubku Čtvercový průřez nosníku 1x1 cm. Plocha 1e-4 m2, moment setrvačnosti I=8.3333e-10 m4 Kvadratický moment průřezu pro obdélníkový profil to je celý program, včetně dat, sestavení i řešení x=[0 0 1]; y=[0 1 1]; con=[1 2; 2 3]; ei=[8.3333e-10 8.3333e-10];ea=[1e-4 1e-4]; nu=length(x) ne=length(ea) n=3*nu; k=zeros(n,n);b=zeros(n,1); for e=1:ne [kl,ig]=klocp(e,x,y,con,ei,ea); k(ig(1:6),ig(1:6))=k(ig(1:6),ig(1:6))+kl(1:6,1:6); end for i=1:3 k(i,:)=0;k(i,i)=1; end b(n-1)=-100; p=k\b; Fy=-100 N nulové posuvy i natočení v uzlu 1 a zatížení -100N v uzlu 3

49. Příhradová konstrukce 3/3 NAP5 Výsledkem řešení je vektor uzlových parametrů p(1:9), tzv. zobecněné posuvy p = -0.0000 ux1 0.0000 uy1 0 1 0.3000 ux2 -0.0000 uy2 -0.6000 2 0.3000 ux3 -0.8000 uy3 -0.9000 3 Vnitřní síly a momenty se dopočítají jako součin matice tuhosti vybraného elementu krát vektor posunutí (uzlů tohoto elementu). Například pro element číslo 2 f=kl*p(ig(1:6)) f = 0 Fx2 100.0000 Fy2 100.0000 M2 rameno síly 1m, síla 100 N 0 Fx3 -100.0000 Fy3 0.0000 M3 kl je lokální matice tuhosti a vektor ig vybere 6 odpovídajících posuvů z vektoru globálních posunutí p.

50. Rotačně symetrické nádoby NAP5 V rotačně symetrických nádobách zatížených např. vnitřním tlakem existují normálová membránová napětí (vzorečky typu pr/s), ale v místech přechodů, kde se mění křivost nebo tloušťka, se objeví i napětí ohybová a smyková (vlny napětí s dosahem typicky (rs)). Nádobu je pak třeba chápat jako skořepinu. MKP systémy, např. ANSYS, COSMOS,… vždy obsahují prostorové skořepinové prvky (zakřivené trojúhelníky, či čtyřúhelníky), ale pro aplikačně důležité rotační skořepiny se symetrickým zatížením stačí jednorozměrné dvouzlové elementy. Asi nejúspěšnější element tohoto typu navrhl už asi před 30 lety ing.Vykutil (VÚT Brno), je implementován ve FEMINě a dá se snadno napsat i v MATLABu. Přesný popis elementu uvádí Schneider P., Vykutil J.:Aplikovaná metoda konečných prvků, skripta VÚT Brno, 1997. ukázky grafických výstupů programu FEMINA (Vykutilovský element)