NUMERICK� ANAL�ZA PROCESu by 8ae1Pwz

VIEWS: 60 PAGES: 56

									  NAP5         NUMERICKÁ ANALÝZA
               PROCESů



                                       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
NAP5
       Modely                 ODE okrajový problém
  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
NAP5
       ODE výměníky tepla
NAP5
       ODE výměníky tepla
  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,
         dT1      k                                                      řešitelný stejně jako dřív,
               (T1  T2 )              T 1’                          třeba metodou Runge Kutta
         dx      W1                                              x
         dT2    k
                  (T1  T2 )            T 2’
         dx W2
                                                          L
   Tentýž případ, ale PROTIPROUD                                  okrajový problém, chybí
                                                               počáteční podmínka pro teplotu
         dT1    k                                                 druhého proudu pro x=0
               (T1  T2 )              T 1’
          dx    W1                                               x
         dT2    k
                 (T1  T2 )                                                          T 2’
         dx     W2
  (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é)
NAP5
       ODE výměníky tepla
  Rovnice teplotních profilů jednotlivých proudů lze zapsat v maticovém tvaru
             d                                                                speciální případ pro
                [T ]  [[ A]]  [T ]                                           souproud a pouze
                                                                                  dva kanály
             dx
                          T1                  k / W1 k / W1 
                  [T ]   
                         T         [[ A]]                     
                          2                  k /W      k / W2 
                                                      2          
  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í

           [T ]  [[U ]]  [ Z ]
  kde [[U]] je matice vlastních vektorů matice [[A]] splňující podmínku
          [[ A]]  [[U ]]  [[U ]]  [[]]              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).
NAP5
        ODE výměníky tepla
       Dosazením transformace do soustavy ODE
                              d
                                 [[U ]]  [ Z ]  [[ A]]  [[U ]]  [ Z ]        Čtvercová matice. Když je
                                                                                   [[U]] matice vlastních
                              dx                                                 vektorů, bude diagonální.
                               d
                                 [ Z ]  [[U ]]1  [[ A]]  [[U ]]  [ Z ]
                              dx
       získáme soustavu rovnic pro transformované teploty (vektor [Z] )
                               d
                                  [ Z ]  [[]]  [ Z ]
                               dx
                                             ODE pro Z jsou nezávislé, protože
                                                [[]] je diagonální matice
       Řešení ODE
                               Zi ( x)  di ei x         i  1, 2,...
                               [ Z ]  e[[  ]] x [d ]
       Návrat k původním teplotám
                                [T ( x)]  [[U ]]  [[e x ]]  [ d ]            koeficienty dj jsou dány
                                                                                 okrajovými podmínkami

                                Ti ( x)   U ij e j d j
                                                          x

                                              j
NAP5
         ODE výměníky tepla
 Konkrétní případ dvoukanálového souproudého výměníku
                          T1 ( x)  U11d1e1 x + U12 d 2e2 x
                          T2 ( x)  U 21d1e1 x + U 22 d 2e 2 x
  Vlastní problém pro matici [[A]] 2 x 2 se dá řešit analyticky:
                                                                                                 0        0      
                  k / W1 k / W1               1 W2                                                          
       [[ A]]  
                 k /W                [[U ]]  
                                                1  W 
                                                                                        [[]]          1   1 
                           k / W2                  1                                         0  k (W + W ) 
                         2
                                                                                                          1   2 


 Okrajové podmínky (v tomto případě spíš počáteční podmínky) pro x=0

   T1 '  1 W2  d1                                                      W T '+W2T2 '                T 'T '
   
   T '  1  W  d 
                                       a řešení soustavy              d1  1 1                   d2  1 2
                                                                                                          .
                                                                                                        W1 + W2
   2         1  2                                                       W1 + W2
    Ověření v MATLABu (k=1, w1=1, w2=2)         zdá se, že to vyšlo jinak? ne, jen se
                                             zaměnilo pořadí vlastních vektorů a ty se
    a=[-1 1;0.5 -0.5];                       normovaly na jednotkovou Eukleidovskou
    [u,l]=eig(a)                              normu. Na výsledek to nemá žádný vliv.
                   u=
                     -0.8944 -0.7071
                      0.4472 -0.7071
                   l=
                     -1.5000     0
                         0    0
NAP5
       ODE výměníky tepla
  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
                                                                                          d
                                                                                            [T ]  [[ A]]  [T ]
                                                                                         dx
   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]
                                                [T ( x)]  [[U ]]  [[e x ]]  [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).
NAP5
       ODE výměníky tepla
   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

         [T ( x' )]M  [[G ' ]]MxN  [T ' ]N + [[G ]]MxM  [T ( x' ' )]M
         [[U ]]MxM  [[e x ' ]]MxM  [d ]M  [[G ' ]]MxN  [T ' ]N + [[G ]]MxM  [[U ]]MxM  [[e x '' ]]MxM  [d ]M
         ([[V ' ]]MxM  [[G ]]MxM  [[V ' ' ]]MxM )  [d ]M  [[G ' ]]MxN  [T ' ]N

                              U11e1x1 ' U12e2 x1 ' ... U1M eM x1 ' 
                                  1 x2 '       2 x 2 '             M x2 '
                                                                              
                             U e          U 22e          ... U 2 M e         
                  [[V ' ]]   21                                             
                                ...                                          
                             U e1xM '                              M x M ' 
                              M1              ...            U MM e          
   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

       [T '']N  [[G '']]NxM  [[V '']]MxM  ([[V ']]MxM  [[G ]]MxM  [[V '']]MxM ) 1  [[G ']]MxN  [T ']N

   Xing Luo, Meiling Li, Wilfried Roetzel: A general solution for one-dimensional multistream heat exchangers
   and their networks. International Journal of Heat and Mass Transfer 45 (2002) 2695–2705
NAP5
       ODE výměníky tepla
   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.
NAP5
         ODE výměníky tepla
  Ukázkový příklad dvoukanálového výměníku                                         T1’=90
                                                                          90
  délky L=1, souproud i protiproud. Řešení                                85

  metodou Runge Kutta v MATLABu (ode45)                                   80

                                                                          75

                                            definice ODE pro protiproud   70
               definice ODE pro
                                            se liší jen znaménkem druhé
                   souproud                                               65
                                                         rovnice
 function dtdx=dts(x,t)           function dtdx=dtp(x,t)                  60

 dtdx=zeros(2,1);                 dtdx=zeros(2,1);                        55
 w1=1; w2=2; k=1;                 w1=1; w2=2; k=1;                        50
 dtdx(1)=-k/w1*(t(1)-t(2));       dtdx(1)=-k/w1*(t(1)-t(2));
 dtdx(2)=k/w2*(t(1)-t(2));        dtdx(2)=-k/w2*(t(1)-t(2));
                                                                          45       T2’=40
                                                                          40
                                                                               0      0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1
     řešení pro souproud
    je počáteční problém
                                                                          90
                                  >> sol=ode45(@dts,[0 1],[90 40]);
                                  >> plot(sol.x,sol.y)                    80
                                                                                   T 1’
    řešení pro protiproud
                                                                          70
      metodou cik-cak
                                  s=ode45(@dtp,[0,1],[90,50]);
                                  for i=1:10
                                                                          60
                                  n=length(s.x);
                                  s=ode45(@dtp,[1,0],[s.y(1,n),40]);
                                                                          50
                                  n=length(s.x);
                                  s=ode45(@dtp,[0,1],[90,s.y(2,n)]);
                                  end                                     40                                                                    T 2’
                                  plot(s.x,s.y)
                                                                          30
                                                                               0     0.1    0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1
NAP5
        ODE sušárny, extraktory…
       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.
NAP5
       Potrubní sítě
 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)
NAP5
       Potrubní sítě
       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
                         v 1 v 2      p
                       ( +        )     2 v +  g x
                         t 2 x        x

       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
           v 1 v 2    p    m                       2  DA2                    hmotnostní průtok
         ( +        )           +  gx         k
           t 2 x      x k (m, x)                     m

                              Třecí ztráty              dArcy Weissbachův součinitel
                                                        tření, v lamináru např. 64/Re


       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…).
NAP5
       Potrubní sítě MKP, MVR
       Uvažujme speciální případ: stacionární průtok, konstantní průřez. Předchozí
       rovnici derivujme       p       k
                                 (k        )  gx
                            x        x            x
       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

                           p          k
                            (k )  g x       res( x)
                          x x         x
                                                            Kdyby bylo p(x) přesné řešení,
                                                             byla by funkce res(x) nulová
NAP5
        Potrubní sítě MKP, MVR
   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)

                              res( x)w( x)dx  0
   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
                             p ( x)   p j N j ( x)
                                       j 1

   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).

                                  N j          k 
                 L    N

                  ( p j
                 0    j 1       x
                                    (k
                                       x
                                          ) g x
                                                  x
                                                      ) wi dx  0   i  1, 2,..., N

       Výsledkem je soustava N lineárních algebraických rovnic pro N neznámých
                  N          L
                                      N j              L
                                                            k
                 
                 j 1
                      p j (  wi
                                 x
                                    (k
                                        x
                                            )dx)   wi g x
                                                             x
                                                                dx
                            0                      0
NAP5
       Metoda vážených residuí
  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é funkce wi(x))
   (např. metoda hraničních prvků)


  Metody konečných prvků (spojité váhové funkce)
                                                                              xi



  Metody konečných objemů (nespojité, po částech konstantní funkce)
   (metoda kontrolních objemů)
                               wi  h( x  xi 1 / 2 )  h( x  xi +1 / 2 )
                                                                              xi

  Kolokační metody (nulové reziduum v uzlech. Váhové funkce jsou delta funkce)
   (metoda konečných diferencí)
                                     wi   ( x  xi )
                                                                               xi
NAP5
       Konečné prvky a objemy
  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).




                xi    wi(x)                        xi    wi(x)
                          MKP váhové funkce                      MKO váhové funkce nejsou
               x         jsou spojité a existuje
                          jejich první derivace
                                                   x               spojité na hranicích
                                                                   kontrolních objemů


  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).
NAP5
       Konečné prvky a objemy


       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í.
NAP5
       Metoda konečných prvků
   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

                               Ni N j                 k
                    N      L                  L
                   p j (  k          dx)   Ni g x      dx
                   j 1    0
                                x x         0
                                                        x

   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.
NAP5
       Metoda konečných prvků
   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ý
                                                                        Ni
                                     1 1
                                     ke                                             Nj
                           [[ K ]]e                    
                                     Le  1            1                          i                             j
                                                                                                  Le
   Sčítání matic elementů do globální matice soustavy (a
                                                                                                              Podprogram
   vektoru pravé strany) se provádí algoritmem sestavení. V                             K=zeros(n,n);        výpočtu Ke
   cyklu přes všech Ne elementů se volá procedura výpočtu                               b=zeros(n,1);
                                                                                        for e=1:ne
   lokální matice pro daný průtok a geometrii a tato matice                              [Ke,be]=local(e);         Matice
   se přičte k matici globální. Korespondence mezi lokálními                             for i=1,2
                                                                                                                 konektivity

   indexy (1,2) a globálními indexy musí být definována v                                   ig=c(e,i);
   matici konektivity c, jejíž řádky odpovídají elementům a                                 b(ig)=b(ig)+be(i);
                                                                                            for j=1,2
   ve dvou sloupcích jsou indexy i,j uzlů každého elementu.                                     jg=c(e,j) ;
   Tím je vlastně definována topologie sítě.                                                   K(ig,jg)=K(ig,jg)+Ke(i,j);
                                                                                            end
   Pozn.: Korespondence mezi lokálními a globálními indexy bývá komplikovanější,         end
   protože uzlovým bodům většinou odpovídá více než jeden parametr (třeba posuvy)
                                                                                        end
   a elementy mívají různé počty uzlů (trojúhelníky,…). Matice konektivity pak má
   proměnnou délku řádků.
NAP5
       Metoda konečných prvků
   Fyzikální interpretace: Když vynásobíme lokální matici jednoho elementu [[K]]e
   vektorem vypočtených tlaků dostaneme hmotnostní průtok elementem (e),
   přesněji průtok do uzlu i a do uzlu j
                     K ii
                         e
                             K ij   pi   bie   m i 
                                e
                                                      e
                     e               e    e 
                    K
                     ji
                               e  
                             K jj   p j   b j   m 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.
NAP5
       Metoda konečných prvků
   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ř.
                                 Pásová                          Vroubená matice
                                 matice                             (sky line)




   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
NAP5
       Metoda konečných prvků
   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ě
                        [[ K ( p ( k 1) )]][ p ( k ) ]  [b]
   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.
   N                                 N    N
                                              Kim ( p( k 1) ) ( k 1)
   Kij ( p(k 1) ) p(jk 1)  bi +  ( p
  j 1                              j 1 m 1
                                                               pm + Kij ( p ( k 1) ))( p (jk )  p (jk 1) )  0
                                                      j

           residuum                                 residuum / p
NAP5
        Metoda konečných prvků
   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á).
NAP5
       Potrubní sítě - model cirkulace
       Příklad aplikace MKP – ale v trochu jiné podobě (není použita metoda
       vážených residuí, takže žádné integrály, žádné algoritmy sestavení, jen
       opakování Bernoulliho rovnice a rovnice kontinuity)




                                                Řešení toku kapaliny v umělém
                                                cirkulačním obvodu krve (fyzikální
                                                model laboratoře kardiovaskulární
                                                mechaniky).
NAP5
        Potrubní sítě - uzly
   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]             19          18                           Počet neznámých 2N-Nb,
         16                                                                             kde N je počet uzlů a Nb je počet
                                                    28
                15                                                                      koncových uzlů (v každém je
         27                                              17                             předepsaný tlak nebo průtok jako
              14                                                              22        okrajová podmínka)
                                                                20     29

                     9                         7
                      25                       24
  12                                                     6                  21
   26                    8                         5
                10
        11                                                           uvědomte si, že v určitém místě potrubní
                                                                    sítě není možné už z principu věci zadávat
  13                                                                současně tlak a současně průtok. To platí i
                                                                     pro předchozí Galerkinovu metodu řešení
                              4            3
                                  23

                                       2            n Předepsaný tlak p, počítaný průtok Q
                                                    m Počítaný tlak i průtok
                                                    k Uzel pro definici geometrie větvení
                                                          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
              -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                                 X [m]
NAP5
       Potrubní sítě - elementy
   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:                                              počet elementů
                                      2N-Nb =                           2Mp+3Mv                        větvení

                                                                počet elementů
                                  počet uzlů
                                                                    trubka
                                                                                      d       x2,y2,z2 p2,Q2
                                                   počet
                                               koncových uzlů




                                                                                  x1,y1,z1 p1,Q1

                                                                                      x2,y2,z2 p2,Q2
                    V1(2,3,4,23,d1,d2,d2)                                                            x3,y3,z3 p3,Q3

                   P1(1,2,d1)                                                    d1
                                                                                          x4,y4,z4

                                                                                 x1,y1,z1 p1,Q1
NAP5
        Potrubní sítě - Bernoulli
  Základem předchozího řešení jsou Bernoulliho rovnice, vzniklé integrací úvodní
  pohybové rovnice po proudnici mezi uzly 1 a 2

                       dv  2 2        p  p1
                 L12     + (v2  v1 ) + 2     + g x ( x2  x1 ) + ez12  0
                       dt 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ů,
            d                                  128   LQ         38 8Q2                           Hagenbachův faktor, výpočet
                            L          ez12       t 4 + (1.2 + ) 2 4                        lokálních ztrát, turbulentní viskozita
                                                   d          Re  d

            Q2
                                               8Q12               300Q1 128         QL           QL
            d2                  Q3     ez12    2 4 (k + 1) 12 +
                                                              t
                                                                         +    (  t 1 1 4 1 + t 2 2 4 2 )
                       L2                       d1                d13            d1           d2
       d1
                                                                p1  p2                 150
                  L1                                    12             (k + 1) 12 +
                                                                                    t

                                                                 v1 2
                                                                                        Re1
                                                                   2
       Q1
NAP5
       Potrubní sítě MATLAB
       MATLAB vzorové řešení: mainpq.m


       Soubory vstupních dat:

       xyz.txt   xyz                  (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í)
NAP5
       Táhla, nosníky, skořepiny
NAP5
        Táhla, nosníky, skořepiny
       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ů.
       Naprosto 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.
NAP5
        Táhla, nosníky, skořepiny
       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ů
NAP5
        Táhla
       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 x modul pružnosti)
       a vztlakový člen může být nahrazen spojitým vnějším zatížením táhla fx
                        2u x                        fx(x), ux(x)
                     EA 2 + f x  0
                       x
       Metoda vážených residuí tedy vede na stejnou soustavu rovnic pro posuvy uzlů
                                           Kij
                                                                    bi
                                Ni N j
                     N     L                     L
                   uxj (  EA          dx)   Ni f x dx
                   j 1    0
                                 x x         0                         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 i matice tuhosti jednoho elementu je totožná                          L
                                                                                 1   du
                                                                          W   ( EA( x ) 2  f xu x )dx
                                                                                 2   dx
                                     EAe  1 1
                                                                              0

                                                                         a hledali její minimum.
                          [[ K ]]e            
                                      Le  1 1 
NAP5
       Táhla
   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)
       y
                          u                         ux1         Fx1   c    0
                                                     
                uy
                                                                        
                                u 1   c s 0 0   u y1        Fy1   s      
                                                                                 0   F 1 
                                                                           
                     ux
                                                                                       
                               u 2   0 0 c s   ux 2 
                                                  
                                                                   Fx 2   0   c   F 2 
                                                                
                                                                   F  0         
                          x
                                                    u                  
                                                     y2          y2         s
NAP5
       Táhla
   Rovnice táhla v lokálním souřadném systému
                  EAe  1 1  u 1   F 1 
                                        
                   Le  1 1   u 2   F 2 

   se po dosazení transformačních vztahů změní na
                    c 0                       ux1   Fx1 
                                                          
               EAe  s 0   1 1 c s 0 0   u y1   Fy1 
                                                   
                Le  0 c   1 1  0 0 c s   ux 2   Fx 2 
                                               
                                               u   F       
                    0 s                       y2   y2 
   Pronásobením matic získáme finální podobu matice tuhosti táhla pro souřadný
   systém x,y           2
                         c    cs c cs 
                                    2

                                           2           >> t=[c^2 c*s;c*s s^2];
                     EA cs      s 2
                                     cs s 
           [[ K ]]e  e  2                              >> ke=[t -t;-t t]
                      Le  c cs c 2    cs 
                         
                          cs  s 2 cs      
                                         s2 
                                                                   všímněte si jak elegantně se
                                                                       v MATLABu dá napsat
                                                                      matice tuhosti elementu
NAP5
        Táhla                                      to je celý program,
                                                        včetně dat,
                                                    sestavení i řešení



       Příklad                                       5
                                                                         x=[0 1 0 1 2]; y=[0 0 1 1 2];
           y                                                             con=[1 3; 2 3; 2 4; 3 4; 3 5; 4 5];
                                e5                                       a=[4e-4 4e-4 4e-4 4e-4 4e-4 4e-4];
                                          e6                             nu=length(x)
           3            e4                                               ne=length(a)
                                      4                                  n=2*nu;                                        jediným příkazem se dá
                                                                                                                       příčíst lokalní matice 4x4
                                                                         k=zeros(n,n);b=zeros(n,1);
                                                                                                                       ke globální matici tuhosti
            e1          e2           e3                                  for e=1:ne
                                                                           [kl,ig]=kloc(e,x,y,con,a);
                                                     x                     k(ig(1:4),ig(1:4))=k(ig(1:4),ig(1:4))+kl(1:4,1:4);
                                                                         end
            1                    2                                       for i=1:4
                                                                            k(i,:)=0;k(i,i)=1;              nulové posuvy v uzlech 1,2 a
                                                                         end                               zatížení -100N v uzlu 5 (ošklivé
       vektor korespondence                                                                              řešení, platí jen pro danou topologii)
                                                                         b(n)=-100;
       lokálních a globálních
               indexů                     matice konektivity
                                                                         p=k\b;

   function [kl,ig]=kloc(ie,x,y,con,a)                                                            Gaussova eliminace.
   E=200e9;                                                                                      Vektor p jsou výsledná
                                                                                                        posunutí
   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
NAP5
        Nosníky (trubky)
   Průhybová rovnice momentové rovnováhy nosníku je diferenciální
   rovnice čtvrtého řádu
                                                                                                      Fx
                                                                                 z
                      4u z  2u z                                                                         x
                   EI 4 + Fx 2 + kuz  f z
                     x     x                                                         fz

   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

                               u                        u
   residua je
               b
                        4
                                   zj
                                         2
                                         Nj                  zj   Nj
                Ni ( EI
                               j
                                              + Fx       j
                                                                        + k  uzj N j  f z )dx  0
               a
                               x   4
                                                         x   2
                                                                             j

       Dvojím použitím integrace per partes snížíme řád derivací
                             2 Ni  N j Ni N j
                      b             2                              b

                   [ ( EI x2 x2  Fx x x + kNi N j )dx]uzj   Ni f z dx
                   j a                                             a

   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ý.
NAP5
        Nosníky (trubky)
    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                 2
                                               uz
                                         EI         M ( x)
                                              x 2
   kde druhá derivace průhybové čáry je mírou zakřivení (1/d 2u/dx2 je poloměr kruhového oblouku).

   Ohybový moment může být způsoben nejenom příčnou silou (třeba reakcí M=F z(b-x)), ale i axiální silou Fx
                                                      2M F        2u z               4u z   2u z
               M F ( x)  (u z (b)  u z ( x)) Fx            Fx 2               EI 4 + Fx        0
                                                       x 2       x                  x      x 2
   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 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í
                            b x
                                                       M                     2M                  4u z
                                                                  b
                 M ( x)     
                             0
                                    f z ( x +  )d 
                                                        x
                                                              f z ( )d 
                                                               x
                                                                              x 2
                                                                                    f z ( x)  EI 4  f z ( x)
                                                                                                  x
                b
           x                                               Tento přechod úplně triviální není,
                                                          jde o derivaci integrálu s proměnnou
                M(x) M(x)                                       horní mezí i integrandem

                            fz
NAP5
         Nosníky (trubky)
   Vraťme se ještě jednou k aplikaci integrace per partes na rovnici vážených reziduí (s obecnou váhovou funkcí w(x))

                                             4u z    2u z
                                     b

                                      w(EI x4 + Fx x2 + kuz  f z )dx  0
                                     a

       Dvakrát opakovaná integrace na první člen se čtvrtou derivací a jedna integrace per partes pro druhý člen dává
               2 w  2u z w uz                           3uz     uz       w  2uz b
        b

         (EI x2 x2  Fx x x + w(kuz  f z ))dx + [w(EI x3 + Fx x )  EI x x2 ]a  0
        a
       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)                    F (b)
                                                                                                                           z



              2 w  2u z      w uz                             w b
        b                                                                                         M(b)

         x x2
        a
          (EI 2            Fx
                               x x
                                      + w(kuz  f z ))dx  [wFz +
                                                                  x
                                                                     M ]a  0
                                                                                                                                 M(b)
                                                                                                      Fz(a)

       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í rovnice a 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.
NAP5
            Nosníky (trubky)
       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
                        2u                              w uz          u L
                                                    L                                 L
                     EI 2z  M ( x)
                       x
                                                    x x
                                                    0
                                                      EI        dx  [wEI z ]0 =   wMdx
                                                                         x        0

        Vetknutý nosník zatížený konstantním momentem                                 Tatáž úloha,ale kubické bázové funkce. V každém
        M, použity lineární bázové funkce. V uzlech je                                uzlu je dvojice parametrů, posunutí a natočení.
        jediný parametr, posunutí uz. Nulové natočení ve                              Nulové natočení lze předepsat, ale řešení s velkým
        vetknutí (x=0) ani nejde předepsat.                                           počtem elementů tu podmínku vlastně ignoruje.
       5                                                                                  5

       4             přesné řešení
                               Mx 2                 U řešení uvedených na obrázcích       4        20 elementů
                   uz ( x) 
       3

       2                       2 EI                       je tento člen zanedbán,
                                                       předpokládáme, že platí tzv.       3                                          1 element
       1
                                                      přirozené okrajové podmínky.
       0                                                  Není to příčina selhání?        2                            2 elementy
       -1

       -2
                                                                                          1               3 elementy
                                       3 elementy
       -3
                                                                                          0
                20 elementů
       -4                                                                                                                Přesné řešení
       -5                                                                                 -1
            0       0.2          0.4   0.6       0.8         1                                 0    0.2          0.4        0.6      0.8         1
NAP5
        Nosníky (trubky)
                                       uz L
       Zanedbání členu        [ wEI      ]0          znamená, že předpokládáme nulové duz/dx v koncovém bodě L, a rovnice
                                       x
       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.
          5                                                    5                                               Řešení uvažující člen
                                                                                                                           uz L
                   Lineární bázové a váhové funkce                                                                [ wEI      ]0
          4                                                    4        Kubické bázové a váhové funkce                     x
          3                                                    3
                                                                          Libovolný počet
          2
                                                                             elementů
                     3 elementy           20 elementů          2


          1                                                    1


          0                                                    0


          -1                                                   -1
               0   0.2     0.4      0.6         0.8       1         0      0.2      0.4     0.6     0.8   1


       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ů.
NAP5
       Nosníky (trubky)
       Finální formulace Galerkinovy metody vážených residuí pro Bernoulliovský
       nosník na pružném podkladě
                                                              Fz

                                                          M
                                                                   Fx
                                        z
                                                                        x
                                               fz



              2 Ni  N j      Ni N j                                            N
       b             2                                      b

  
  j
    [  ( EI
              x 2 x 2
                           Fx
                                x x
                                        + kNi N j )dx]uzj   Ni f z dx + [ Ni Fz + i M ]b
                                                                                    x
                                                                                         a
      a                                                     a




                       Matice tuhosti                         Spojité
                                                                                    Osamělé síly
                       nosníku [[K]]                          zatížení
                                                                                     a momenty



            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í!
NAP5
          Nosníky (trubky)
      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
                                                                                                    Bázová funkce mající ve všech
                                 Bázová funkce rovná 0 ve všech
                                                                                                    uzlech nulové první derivace a v
                                 uzlech a v jednom uzlu má první
                                                                                                        jednom uzlu se rovná 1
                                        derivaci rovnou 1                     N zi ( x)

                      N di (x)

                           i           i+1
                                                  ( x  xi 1 )2 ( x  xi )     i          i+1
                                                                                                      ( x  xi 1 )2 (3xi  xi 1  2 x)
                                      N di ( x)                                          N zi ( x) 
 Rovnice průhybové
 čáry v nosníku 1-2
                                                      ( xi  xi 1 )2                                            ( xi  xi 1 )3

                      u z ( x)  u z1 N z1 ( x) + 1 N d 1 ( x) + u z 2 N z 2 ( x) +  2 N d 2 ( x)
      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
                                         [[ K ]]  [[ K M ]]  Fx [[ K F ]] + [[ M ]]
       [[KM]] je maticí ohybové a [[KF]] geometrické tuhosti a [[M]] je maticí hmot.
NAP5
               Nosníky (trubky)
  Pro matice jednoho elementu o délce L získáme integrací tyto výsledky

                                         uzlové parametry
                      u zi            (příčný posun a natočení
                      
                      
                                              neutrální osy)

                      i
                                                                                                         12 6 L        12  6L 
                                                                                                                       6 L 2 L2 
                                                                                     2 Ni  N j
                                                                                            2
       1                                                                        L
                                                                                                     EI      4 L2                
                                                                   [[ K M ]]   EI              dx  3
   0.8
                                                                                     x 2 x 2       L                  12 6 L 
                                                                                                                                 
                                                                                0
                                                                   Matice tuhosti
   0.6
                                                  N1                                                     sym                4 L2 
                                                  N2
   0.4                                            N3



                                                                                                          36 3L      36 3L 
                                                  N4




                                                                                                                     3L  L2 
   0.2


                                                                                  N i N j                    4 L2
                                                                               L
                                                                                                     1                        
       0
                                                                   [[ K F ]]                 dx 
                                                                                   x x            30 L              36 3L 
           0                    0.5                    1




                                                                                                                              
   -0.2
                                                                               0
                                                                    Matice geometrické tuhosti
                                                                                                          sym            4 L2 
                N1  1  3( x / L) 2 + 2( x / L)3 
                                                                                              156 22 L 54 13L 
                N 2  x(1  2 x / L + ( x / L) 2 ) 
                                                                               L                     4 L2 13L 3L2 
                N3  3( x / L) 2  2( x / L)3                    [[ M ]]   kNi N j dx 
                                                                                            kL                      
                N  x (  x / L + ( x / L) 2 ) 
                                                                                          420           156 22 L 
                                                                                                                  
                                                                                0
                     4
                                                                    Matice hmot
                                                                                                 sym           4 L2 
NAP5
       Nosníky (trubky) MATLAB
  Procedura generování lokálních matic elementu nosník (ei-momenty setrvačnosti,k-tuhost podkladu)
                                                                                12 6 L        12  6L 
                                                                                    4 L2      6 L 2 L2 
                                                                            EI
       function [kstif,kgeom,mass,ig]=kloc(ie,x,y,con,ei,k)      [[ K M ]]  3                          
                                                                            L                  12 6 L 
       EI=ei(ie);E=200e9;                                                                               
       i1=con(ie,1);i2=con(ie,2);                                               sym                4 L2 
       ig=[2*i1-1 2*i1 2*i2-1 2*i2];
       l=((x(i1)-x(i2))^2+(y(i1)-y(i2))^2)^0.5;                                    36 3L      36 3L 
       kstif=EI*E/l^3*[12 6*l -12 6*l;6*l 4*l^2 -6*l 2*l^2;                            4 L2   3L  L2 
                                                                              1                        
               -12 -6*l 12 -6*l;6*l 2*l^2 -6*l 4*l^2];           [[ K F ]] 
                                                                             30 L              36 3L 
       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];                            sym            4 L2 
       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];                   156 22 L 54 13L 
                                                                                    4 L2 13L 3L2 
                                                                           kL                      
                     uzlové                                      [[ M ]] 
        u zi      parametry                                              420           156 22 L 
        
                                                                                                
        i                                                                     sym           4 L2 

  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.
NAP5
        Příhradová konstrukce 1/3
  Procedura generování matice tuhosti ocelového elementu nosník + táhlo
                                                           uy   z
   function [kstif,ig]=klocp(ie,x,y,con,ei,ea)                                                AL2                      AL2              
   EI=ei(ie);EA=ea(ie); E=200e9;                                                              I    0  0                        0   0 
   i1=con(ie,1);i2=con(ie,2);                            ux                                                             I               
   ig=[3*i1-2 3*i1-1 3*i1 3*i2-2 3*i2-1 3*i2];                                                    12 6 L               0      12 6 L 
   l=((x(i1)-x(i2))^2+(y(i1)-y(i2))^2)^0.5;                                               EI 
                                                                                             
                                                                                                      4 L2              0      6 L 2 L2 
                                                                                                                                         
                                                                                 [[ K ]]  3
   c=(x(i2)-x(i1))/l;s=(y(i2)-y(i1))/l;                                                   L                        AL2                  
   r=[c s 0;-s c 0;0 0 1];z=zeros(3,3);                                                                                         0   0 
                                                                                                                    I                   
   Q=[r z;z r];
                                                                                                                               12 6 L 
   km=EI*E/l^3*[EA*l^2/EI 0 0 -EA*l^2/EI 0              0;
                                                                                              sym                                  4 L2 
                 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;
                                                                                                 to je jen předchozí matice [[KM]] rozšířená o první a
                 0        6*l 2*l^2     0      -6*l 4*l^2];                                      čtvrtý řádek a sloupec, kam je dosazena matice
   kstif=Q'*km*Q;                                                                                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
                           C       S    0    0    0  0
                                                                                    uzlu i (výsledek výpočtu)                 zobecněná zatížení v uzlu
                            S     C 0       0 0 0
   uy      z                                                                                                                i (pravá strana soustavy)
                            0      0    1    0 0 0
                        Q                        ,               C  cos , S  sin
                                                                                                          u xi                      Fxi 
                            0      0    0    C S 0                                                      
                            0                                                                                                        
                ux                  0    0    S C 0                                                     u yi 
                                                                                                                                      Fyi 
                                                                                                                                M 
                            0
                                   0    0    0 0 1                                                     i                         i
NAP5
         Příhradová konstrukce 2/3
  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        I        ( D2  D14 )
                                                                                                       4

  Čtvercový průřez nosníku 1x1 cm.                                                             64
                                                                                   3
  Plocha 1e-4  m2, moment                    Kvadratický moment průřezu pro I  bh
  setrvačnosti I=8.3333e-10 m4               obdélníkový profil                 12
                                                                                                                   to je celý program,
            y                                                                                                           včetně dat,
                                                                                                                    sestavení i řešení


           2          e2             3                     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)
                                 Fy=-100 N
                                                           n=3*nu;
                                                           k=zeros(n,n);b=zeros(n,1);
                                                           for e=1:ne
            e1                                               [kl,ig]=klocp(e,x,y,con,ei,ea);
                                 x                           k(ig(1:6),ig(1:6))=k(ig(1:6),ig(1:6))+kl(1:6,1:6);
            1                                              end
                                                           for i=1:3
                                                              k(i,:)=0;k(i,i)=1;
                                                           end                               nulové posuvy i natočení v uzlu
                                                           b(n-1)=-100;                        1 a zatížení -100N v uzlu 3
                                                           p=k\b;
NAP5
        Příhradová konstrukce 3/3
  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))
                                     kl je lokální matice tuhosti a vektor
       f=                           ig vybere 6 odpovídajících posuvů z
                                        vektoru globálních posunutí p.
             0            Fx2
        100.0000          Fy2
        100.0000         M2     rameno síly 1m, síla 100 N

             0            Fx3
       -100.0000          Fy3
         0.0000           M3
NAP5
       Rotačně symetrické nádoby
   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)
NAP5
          Rotačně symetrické nádoby
      Element reprezentující kuželový prstenec má dva uzly a v každém 3 uzlové
      parametry: posuv v axiální směru ux, radiálním směru uy a natočení .
                                                                                  N
                                   r=y                                  2                                                  u x1 
                                                                                                                           
                                               ur1                          M                                             u y1 
                                                     rz1=                                                                 
                                                                                                                  [u ]   1 
                                               1                                                                           ux2 
                                                          ux1                                                             u 
                                                                                                                           y2 
                                                                r                                                          
                                                                                                                           2
                                                                                                x


        Matice tuhosti elementu má tedy rozměr 6x6 a počítá se jako součin
                                                                                                          [[ K ]]e  rL[[ B ]]T  [[ D ]]  [[ B ]]
                                                                                   matice derivací

  matice elastických                  1            0          0      0          bázových funkcí
                                                                                                       c  s     0      c         s    0 
                                      
      konstant                             1         0          0       0                                   L                      L       
                                                   h2        h2                                       0          0      0              0 
                                Eh  0     0                            0                                  2r                      2r      
                                                 12        12                                       1 0   0     1 0             0     1 
               [[ D]]5 x 5
                               1  2                                                 [[B ]]5 x 6   
                                                    h2      h2                                       L 0   0
                                                                                                                    sL
                                                                                                                       0 0
                                                                                                                                          sL 
                                      0   0                            0                                         2r                    2r 
                                                  12        12                                                    L                     L
                                                                   5(1   )                           s c         s c                  
                                      0   0         0          0                                                  2                     2 
                                                                     12 
NAP5
       Rotačně symetrické nádoby
   Cílem této partie není odvozování, ani bližší vysvětlování těchto matic, ale
   jen praktický návod, jak sestrojit matici elementu i vektor pravé strany
       function [kstif,bp,ig]=kloc(ie,x,y,con,he,pe)
       % kstif(6x6) matice tuhosti, bp(6) zatížení tlakem, ig(6) indexy globální soustavy               c  s 0       c       s     0 
       % ie index elementu,x(),y() globalni souradnice, con(:,2) konektivita                                 L                L        
                                                                                                        0        0     0             0 
       % he() tloušťky stěn, pe() vnitřní přetlaky                                                          2r                2r       
                                                                                                      1 0     0 1
                                                                                         [[B ]]5 x 6  
       E=200e9;mi=0.28;h=he(ie);p=pe(ie);                                                                               0      0      1 
       i1=con(ie,1);i2=con(ie,2);                                                                     L 0       sL                  sL 
                                                                                                              0         0      0
       ig=[3*i1-2 3*i1-1 3*i1 3*i2-2 3*i2-1 3*i2];                                                              2r                  2r 
                                                                                                                 L                   L
       r1=y(i1);r2=y(i2);                                                                               s c           s c            
       l=((x(i1)-x(i2))^2+(r1-r2)^2)^0.5;                                                                        2                   2 
       r=(r1+r2)/2;
       c=(x(i2)-x(i1))/l;s=(r2-r1)/l;
       b=1/l*[-c -s 0 c s             0;                                                            1      0    0            0      
                                                                                                                                     
                0 l/(2*r) 0 0 l/(2*r) 0;                                                             1     0    0             0 
                0 0 -1 0 0            1;                                                                   h2    h2                 
                0 0 s*l/(2*r) 0 0 s*l/(2*r);                                                Eh      0 0                        0 
               -s -c l/2 s -c l/2];                                         [[ D]]5 x 5                   12   12                   
                                                                                          1   2
                                                                                                           h  2
                                                                                                                  h2                  
       d=E*h/(1-mi^2)*[1 mi 0 0 0;                                                                  0 0                        0 
                           mi 1 0 0 0;                                                                     12   12                   
                           0     0 h^2/12 mi*h^2/12 0;                                                                     5(1   ) 
                                                                                                    0 0     0    0                   
                           0     0 mi*h^2/12 h^2/12 0;                                                                       12 
                           0     0    0 0 5*(1-mi)/12];
       kstif=r*l*b'*d*b;
       pl=p*l/8;
       bp=[-pl*s*(3*r1+r2) pl*c*(3*r1+r2) 0 -pl*s*(3*r2+r1) pl*c*(3*r2+r1) 0];
                                                                                               lokální vektor pravé
                                                                                             strany (zatížení tlakem)
NAP5
         Rotačně symetrické nádoby
       Výsledkem řešení soustavy lineárních rovnic [[K]][u]=[b] je vektor
       zobecněných posuvů, pro element s uzly 1,2
                                                                                                                       N
                         u x1 
                                                                          r=y                               2
                         u y1                                                        ur1                        M
                                                                                          rz1=
                                                                                                      
                 [u ]   1                                                           1
                         ux2                                                                  ux1
                        u 
                         y2 
                                                                                                                                 x
                         2

       Ze zobecněných posuvů lze stanovit            Eh                     Eh u x 2  u x1       u  u r1          u + u r1
                                             N           (  +   )         [         cos + r 2     sin  +  r 2     ]
       5 zobecněných sil, které na element          1  2
                                                                           1  2
                                                                                     L                 L                2R
       působí:                                       Eh                     Eh        u  u x1       u  u r1          u + u r1
                                             N           (  +   )         [ ( x 2     cos + r 2     sin  ) + r 2     ]
                                                    1  2
                                                                           1  2
                                                                                           L              L                2R
       N [N/m] síla ve směru meridiánu
                                                     Eh3                        Eh3        1           + 1
       N [N/m] síla v obvodovém směru       M             (  +   )             [ 2     +  sin  2     ]
                                                  12(1   )
                                                          2
                                                                             12(1   )
                                                                                     2
                                                                                            L              2R
       M [N] moment meridiánový                     Eh3                        Eh3         1         + 1
                                             M             (  +   )             [ 2     + sin  2     ]
       M [N] moment obvodový                     12(1   )
                                                          2
                                                                             12(1   )
                                                                                     2
                                                                                             L            2R

       Q [N/m] příčná smyková síla                   5 Eh      u  u x1        u  u r1        + 1
                                             Q               [ x2      sin  + r 2     cos + 2     ]
                                                  12 (1 +  )      L                L            2
NAP5
       Co je třeba si pamatovat


       Přednáška byla věnována řešení okrajových problémů
       popisovaných obyčejnými diferenciálními rovnicemi. Je dost obsáhlá
       a poměrně důležitá, je třeba si toho zapamatovat trochu víc než
       obvykle…
NAP5
       Co je třeba si pamatovat
   Jak využít řešení vlastního problému pro transformaci soustavy
   provázaných diferenciálních rovnic na soustavu rovnic vždy jen pro
   jedinou závisle proměnnou          [[ A]]  [[U ]]  [[U ]]  [[]]
   Princip metody vážených reziduí. Co je to bázová funkce.
   Co je to váhová funkce. Typy bázových a váhových funkcí                   res( x)w( x)dx  0
                                                                                         N
   v metodách spektrálních, konečných prvků, konečných                       p ( x)   p j N j ( x)
   objemů.                                                                               j 1


   Použití metody konečných prvků na výpočet                   p       k
                                                              (k )  g x
   rozložení tlaků v potrubní síti (nestačitelná kapalina, x x          x
   stacionární průtok)                              N     L
                                                            N N
                                                                          L
                                                                                                       k
                                                p j (  k                        dx)   Ni g x
                                                                        i      j
                                                                                                           dx
   Konečný prvek typu trubka a jeho              j 1               x x                              x
                                        1
                                                              0                              0
   matice průtoku              ke    1
                      [[ K ]]e          
                                Le  1 1 
   Jak se provádí sestavení matic jednotlivých          for e=1:ne
   elementů do globální matice celého systému?            [kl,ig]=kloc(e);
                                                          k(ig(1:2),ig(1:2))=k(ig(1:2),ig(1:2))+kl(1:2,1:2);
                                                        end
NAP5
       Co je třeba si pamatovat
   Jaký je rozdíl mezi táhlem a nosníkem, mezi membránou a skořepinou.
   Aplikace metody konečných prvků na soustavy táhel. Matice tuhosti táhla.

                                2u x                                               EAe  1 1
                             EA 2 + f x  0                              [[ K ]]e            
                               x                                                    Le  1 1 
   Transformace matice tuhosti a zatížení z lokálního souřadného systému         ux1 
   elementu  do globálního souřadného systému x,y                               
                                                            u 1   c s 0 0   u y1 
                                                                          
                                                            u 2   0 0 c s   ux 2 
                                                                                 
                                                                                u 
                                                                                 y2 
  Aplikace metody konečných prvků na nosníky. Typ aproximace průhybové čáry
  (kubické polynomy)     4u       2u
                              EI            z
                                                + Fx         z
                                                                 + kuz  f z
                                      x   4
                                                       x   2


  Matice ohybové a geometrické tuhosti, matice hmot
                     2 Ni  N j                               N N j
              L             2                                L                             L

   [[ K M ]]   EI              dx              [[ K F ]]   i       dx        [[ M ]]   kNi N j dx
               0
                     x 2 x 2                               0
                                                                x x                      0

								
To top