NUMERICK� ANAL�ZA PROCESu - PowerPoint

Document Sample
NUMERICK� ANAL�ZA PROCESu - PowerPoint Powered By Docstoc
					  NAP6         NUMERICKÁ ANALÝZA
               PROCESů




                                       Parciální diferenciální rovnice (PDE), klasifikace na
                                       hyperbolické, parabolické a eliptické


                                       Hyperbolické PDE (kmitání táhel, nosníků, ráz v potrubí)
                                       MOC-metoda charakteristik (stlačitelné proudění)




Rudolf Žitný, Ústav procesní a
zpracovatelské techniky ČVUT FS 2010
NAP6
       PDE   parciální diferenciální rovnice




                                               L.Wagner
NAP6
         PDE                parciální diferenciální rovnice
Když je pro popis systému nutný větší počet nezávisle proměnných (třeba prostorové souřadnice a čas)
musíme řešit místo obyčejných parciální diferenciální rovnice (PDE). Většinou jsou to diferenciální rovnice
s maximálně druhými derivacemi závisle proměnné  (výjimkou je např. biharmonická rovnice deformace
membrán se čtvrtými derivacemi):

Hyperbolické rovnice (Typické pro popis kmitání a vln, třeba elastických vln při rázové deformaci,
a pro nadzvukové proudění. Klíčovou roli hraje konečná rychlost šíření tlakové vlny.)

          2  2
               2  f                     Příklad: vlnová rovnice kmitání pružné tyče, ráz v potrubí
         t 2
               x
Parabolické rovnice (Typické pro evoluční problémy, třeba vývoj teplotního nebo koncentračního
profilu v čase, ale třeba i vývoj profilů v mezní vrstvě v závislosti na vzdálenosti od vstupu. Můžeme to
chápat i jako mezní případ hyperbolických rovnic s nekonečnou rychlostí šíření poruchy tlaku.)

            2 
              2  f                      Příklad: Fourierova rovnice vedení tepla
          t  x
Eliptické rovnice (Typické pro stacionární problémy popisu rozložení teplot, koncentrací,
deformací. Často jde o problém popisovaný parabolickými rovnicemi, ale až po ustálení, tj. pro
nekonečně dlouhý čas. I většina pružnostních problémů je popisována eliptickými rovnicemi. )

          2  2
                     f                   Příklad: Poissonova rovnice
         y 2
                x 2
NAP6
       PDE                    parciální diferenciální rovnice
   Typ PDE je dán pouze koeficienty druhých derivací. Každou PDE druhého řádu můžeme napsat ve
   tvaru       2
                     2
                               
                                2
             a        b      c 2  f
                 x 2
                         xy   y
    kde a,b,c,f mohou být libovolné funkce x,y i hledaného řešení , včetně jeho prvních derivací.
    Koeficienty a,b,c určují rovnice tzv. charakteristik, čar závislostí y(x), které splňují rovnici
                                      dy b  b 2  4ac
                                         
                                      dx     2a
    b2-4ac>0 hyperbolická rovnice (existuje dvojice reálných kořenů, tudíž dvojice charakteristik)
    b2-4ac=0 parabolická rovnice (existuje jen jeden kořen a jedna charakteristika)
    b2-4ac<0 eliptická rovnice (kořeny jsou imaginární a žádná reálná charakteristika neexistuje)

        hyperbolická PDE              y         parabolická PDE           y          eliptická PDE
                                                                                                                y
                       ( x, y)                    y  const   ( x, y)                              ( x, y)
                             y  x
              yx




                  x                                      x                                 x
   znázornění oblasti vlivu (modrá oblast pod charakteristikami) a okrajových podmínek (červeně) na hodnotu řešení (x,y)
NAP6
       PDE Hyperbolické PDE


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

                       2u x  2u x
                    A 2  EA 2  f x
                       t    x                                                                vidíte, že je to opravdu
                                                                                           hyperbolická rovnice, znaménka
                                                                                            druhých derivací jsou opačná


  Metoda vážených residuí založená na aproximaci ux (t , x)   uxj (t ) N j ( x)

                                     N u                                         N u
                                                                                                                        j
           b
                                2
                                              j xj                           2
                                                                                            j xj

            Ni ( EA                                  fx   A                                    )dx  0
                                     j                                            j

           a
                                     x   2
                                                                                  t   2
               Ni v roli váhové funkce          Nj v roli bázové funkce


  Integrace per partes aplikovaná na první člen
             Ni N j                          2uxj b
          b                    b

       [ EA x x dx]uxj  [  ANi N j dx] t 2   Ni f x dx
       j                     j
         a                     a                      a


                                                                                  matice hmot [[M]]                         budicí síly [[b]]
                   matice tuhosti [[K]]
NAP6
       PDE kmitání táhla
  Uvažujme volné kmitání, bez budících sil
                                      2u xj
                  K ij u xj  M ij             0
                                      t   2

  což je soustava jen obyčejných diferenciálních rovnic druhého řádu. Ale s
  konstantními koeficienty, takže existuje analytické řešení u xj (t )  uaj sin(t   )
  jehož dosazením do předchozí rovnice dostaneme soustavu algebraických rovnic
  pro koeficienty amplitud kmitů
                                               ( Kij   2 M ij )uaj  0
  Protože je soustava homogenní, bude mít netriviální (nenulové) řešení jen
  když je matice soustavy singulární, tj. když det(Kij   2 M ij )  0
  Čtvercové matice [[K]] i [[M]] jsou dané a mají N řádků. Determinant je de facto
  polynom N-tého stupně proměnné 2 a algebraická rovnice má N-kořenů, N-
  vlastních frekvencí. Tyto frekvence  spolu s vektory amplitud lze nalézt
  řešením vlastního problému
                                               [[ M ]] 1[[ K ]][ uk ]  k2 [uk ]
             v MATLABU lze tento problém řešit funkcí [u,omega]=eig(m\k)
NAP6
       PDE kmitání táhla
  Pro lineární bázové funkce popisu amplitud v obecném elementu

         N1 ( x)  1  x / L        N 2 ( x)  x / L
  snadno stanovíme matici tuhosti elementu (stejná jako dřív)
                   AE  1  1
          [[K ]]           
                   L  1 1 
                            
                                               račte si povšimnout, že součet
  stejně jako matici hmot                     všech prvků se rovná hmotnosti
                                                      celého elementu
                    AL  2 1 
         [[M ]]        
                 6  1 2
                        
  Matice hmot se často nahrazuje nahrazuje diagonalizovanou maticí, která
  odpovídá rozdělení hmotnosti táhla do jeho uzlů
                      AL  3 0 
         [[M d ]] 
                 6 
                     
                      0 3
                          
                          
   Výhoda diagonalizované matice hmot spočívá v tom, že inverze [[M]] je triviální.
NAP6
       PDE kmitání táhla MATLAB
  Příklad: Soustava ocelových táhel (test pro jedno táhlo)
  function [kl,ml,ig]=kloc(ie,x,con,a)           x=[0 1];con=[1 2];a=[1e-4];
  E=200e9;R=7000;                                nu=length(x);ne=length(a);n=nu;
  ae=a(ie);                                      k=zeros(n,n);
  i1=con(ie,1);i2=con(ie,2);                     m=zeros(n,n);
  ig=[i1 i2];                                    for e=1:ne
  le=abs(x(i1)-x(i2));                           [kl,ml,ig]=kloc(e,x,con,a);
  kl=E*ae/le*[1 -1;-1 1];                        k(ig(1:2),ig(1:2))=k(ig(1:2),ig(1:2))+kl(1:2,1:2);
  ml=R*ae*le/6*[2 1;1 2];                        m(ig(1:2),ig(1:2))=m(ig(1:2),ig(1:2))+ml(1:2,1:2);
                                                 end
                                                 [u,omega]=eig(m\k);

                                                                                                K 22     3E
   V tomto případě můžeme snadno zkontrolovat výsledek                                             
                                                                                                M 22     L2
   (pro diagonalizovanou matici hmot vychází
   vlastní frekvence menší        2E    )
                                    d 
                                           L2
NAP6
       PDE kmitání nosníků
  Možná důležitější je případ příčného kmitání nosníků, např. trubek výměníku
  tepla, podepřených v několika místech přepážkami a na koncích vetknutých
  do trubkovnice. Použijte výsledky předchozí přednášky, kde spojité příčné
  zatížení nosníku nahradíte setrvačnými silami
                          2u z  4u z
                       A 2  EI 4  0
                          t    x
               2 Ni  N j                       2uxj
           b          2          b

       [ EI x2 x2 dx]uxj  [  ANi N j dx] t 2  0
       j                       j
         a                       a



                                                          matice hmot [[M]]
                  matice tuhosti [[K]]



   Dál už je to stejné jako u táhla, použijete jenom dříve definované matice
   tuhosti a hmot odvozené pro nosníkový prvek.

        A nemusíte si dělat starosti s pravou stranou, protože místa podepření trubek i jejich vetknutí do
        trubkovnice jsou silné okrajové podmínky.
NAP6
       PDE Ráz v potrubí

        Nestacionární tok stlačitelné tekutiny v
        potrubních sítích, fungování
        kardiovaskulárního oběhu, to jsou
        typické příklady hyperbolických rovnic.
        Zvláštní případ reprezentuje
        hydraulický ráz, ke kterému dojde při
        náhlém uzavření ventilu v potrubí…
        projeví se náhlým vzrůstem tlaku
        (rázovou vlnou), který se šíří rychlostí
        zvuku v tekutině a stěně potrubí.
        Víte jak tu rychlost spočítat? Víte, jaký
        maximální tlak může vzniknout?

                                                                                 Barbara
                                                                                 Wagner
  M. Rohani, M.H. Afshar Simulation of transient flow caused by pump failure: Point-Implicit Method of Characteristics
  Annals of Nuclear Energy, Volume 37, Issue 12, December 2010, Pages 1742-1750
NAP6
       PDE Ráz v potrubí
  Formulace problému:
  Potrubím s proměnným průřezem A(t,x), protéká tekutina střední rychlostí v(t,x).
  Tekutina je stlačitelná a souvislost její hustoty a tlaku je charakterizována
  objemovým modulem stlačitelnosti K [Pa]          p                              K p
                                                                               
                                                                 Śouvisí s rychlostí
                                                                                                                              a
                                                                                          2

                                                                                                                         
                                                                           
                                                                                           zvuku v tekutině a
                                                                                   K
  Je to vlastně stejný problém jako dříve uvažované potrubní sítě, kdy cílem bylo
  stanovit průběhy tlaku a průtoky. Tentokrát však tlaky i rychlosti závisí i na čase,
  p(t,x), v(t,x).
  Tyto veličiny popisuje stejně jako ve stacionárním případě dvojice rovnic, rovnice
  kontinuity a bilance hybnosti (Bernoulliho rovnice). Ve zjednodušeném tvaru


   Rovnici kontinuity v této poněkud                                p     * v
   zvláštní podobě vysvětlím na                                         K      0
   následující folii
                                                                    t       x
                  K *  a2  
                                     K      takto se zohlední i
                                                                      v p
                                  1
                                     K A   pružnost potrubí
                                            (závislost průřezu na      ez  0
                                     A p   tlaku)                     t x       Tuto rovnici poznáváte: Bernoulliho
                                                                                                 rovnice kde je zanedbán
                                                                                             konvektivní člen (kinetická energie)
NAP6
        PDE Ráz v potrubí
       Rovnice kontinuity pro případ, kdy průřezová plocha trubky A(t,x) se
       mění po délce i v čase (elastická, nebo elasticko-plastická trubka, která
       se „nafukuje“ působením vnitřního přetlaku)
                         hmotnostní bilance

                                                          
                         elementu konstantní                                              akumulace


                                                   ( vA)  ( A)
                         délky dx
               ,A,v
                                                x         t
                         mění se hustota,                                                    souvislost
                         průřez i střední                                                  hustoty a tlaku
                         rychlost .


               dx                                                   1  1 p          1  1 p
                                            přítok-odtok                     ,            
                                                                     t K t           x K x
       Rozderivováním levé i pravé strany a využitím definice objemového modulu
       stlačitelnosti získáme zcela obecnou rovnici kontinuity
                                                                                              v K DA
       p    v     p K A    K A                                                      Dp
                                                                                            K       0
          K     v(       )      0                                                   Dt   x A Dt
       t    x     x A x    A t                                                      (to je totéž, jen zapsané v
                                                                                         materiálových derivacích)




                                                   Tento člen je v předchozí rovnici
                                                    zanedbán (dělá se to často, viz
                                                  např.Khamlichi, Wave motion 1995)
NAP6
        PDE Ráz v potrubí
       To, že je výše uvedený systém rovnic hyperbolický, vyplyne z eliminace tlaku
       nebo rychlosti (stačí derivovat rovnici kontinuity dle času, Bernoulliho rovnici
       dle souřadnice x, a odečíst). Výsledné rovnice pro tlak p(t,x)
               2 p    2  p
                          2
                    a         0,
               t 2      x 2
       nebo rychlost v(t,x)
                                                        K*
               v      2  v
                                        Ŕychlost
                                                   a
                2         2
                    a        0                        
                                        zvuku

               t 2
                         x 2


       jsou stejné a dokonce stejné jako hyperbolická rovnice kmitání táhla. Je
       tedy možné použít i stejnou metodu řešení (metodu vážených residuí).
       Objeví se ovšem malý technický problém: I když jsou rovnice stejné, mají
       různá řešení, protože jsou různé okrajové podmínky, např. časový průběh
       průtoku na jednom konci a časový průběh tlaku na druhém konci. V
       případě rázu v potrubí je to třeba skoková změna průtoku na jednom konci
       (zavření ventilu) a konstantní tlak na druhém konci (zásobník). Metodou
       vážených reziduí je tedy nutné řešit soustavu obou PDE.
NAP6
        PDE




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


Zábranský
NAP6
        Metoda charakteristik MOC
 Vraťme se k výchozí soustavě rovnice kontinuity a rovnice Bernoulliho
    p      v           v p             spoustu věcí jsme zanedbali, například unášivou rychlost, takže následující
       a 
         2
                0,        ez  0         analýza je správná jen tehdy, když rychlost zvuku je mnohem vyšší než
    t      x           t x             rychlost proudění. V řadě uváděných příkladů tento předpoklad splněn není.


  a tyto rovnice (první z nich vynásobme libovolnou nenulovou konstantou ) sečteme
             p 1 p        v      v
        (          )   (   a 2 )  ez  0
             t  x        t      x
  Zvolme libovolnou křivku v rovině t-x (parametrem křivky x() může být přímo čas t).
  Úplné diferenciály tlaku a rychlosti podél této křivky jsou      t
    dp p dx p dv v dx v                                                  x()
                    ,                                                    t=
     dt t dt x dt t dt x
  Jestliže bude křivka x(t) splňovat rovnice (současně)                                                           x
                dx 1               dx
                                      a2
                dt                dt
   získáme finální rovnici vyjádřenou pomocí úplných diferenciálů podél křivky x(t)
                  dp  dv                      (a tuto křivku nazveme charakteristikou diferenciální rovnice)
                     ez  0
                  dt  dt
NAP6
          Metoda charakteristik
                             dx 1
   Předchozí rovnice              a 2 má dvě řešení
                     1       dt 
         1, 2  
                     a
   a jim odpovídají dvě charakteristické křivky a diferenciální rovnice, které je
   třeba na nich integrovat                      Třecí zpráty vyjádřené
                                                          součinitelem třecích ztrát f
        dx1
            a           1 dp    dv f  v | v |
                                             0                                       C
        dt               a dt    dt     2D                                                   x2=-at
                                                                        x1=at
                                                                                                           h/a
        dx2                1 dp    dv f  v | v |
             a                              0                    A                       B
         dt                a dt    dt    2D                                   h

   Přibližná integrace těchto rovnic ve směru charakteristik vede na
   soustavu dvou algebraických rovnic pro neznámé pC, vC
                                                                                                      tC

                                       fh
                                                                                              např.        dp
        1                                                                                              dt dt  p     pA
          ( pC  pA )   (vC  vA )      vA | vA | 0
                                                                                                                 C
                                                                                                      tA

        a                              2aD
        1                             fh
        ( pC  pB )   (vC  vB )      vB | vB | 0
        a                             2aD                              přibližnost integrace spočívá jen v
                                                                        zanedbání proměnnosti třecích
                                                                            ztrát podél charakteristiky
NAP6
          Metoda charakteristik
       Řešení této soustavy vyjádříme v explicitním tvaru
              1 p A  pB               fh
          vC  (          v A  vB      (v A | v A | vB | vB |))
              2    a                 2aD
                 a pA  pB                   fh
          pC      (         (vA  vB )      (vB | vB | vA | vA |))
                 2    a                     2aD
   a můžeme ho ihned využít pro numerické řešení rychlostí a tlaků na časové
   hladině C ze známých hodnot v bodech A,B na staré časové hladině.
NAP6
             Metoda charakteristik
Alternativní (a možná názornější) formulace je analogií postupu, který jsme použili při řešení soustavy
diferenciálních rovnic výměníků tepla. Soustavu dvou diferenciálních rovnic pro tlak a rychlost zapíšeme v
maticovém tvaru
   [W ]          [W ]                                 p                   0             a2            0 
          [[ A]]        [ B]                  [W ]   
                                                       v          [[ A]]  
                                                                             1 / 
                                                                                                    [ B]  
                                                                                                             e 
                                                                                                                 
    t             x                                                                      0            z
 a použijeme transformaci vektoru W takovou, že se ve výsledných rovnicích objeví vždy jen jedna neznámá

   [W ]  [[Q]][Z ]
         [ Z ]               [ Z ]          [ Z ]                      [ Z ]
   [[Q]]         [[ A]][[Q]]         [ B]          [[Q]]1[[ A]][[Q]]         [[Q]]1[ B]
           t                  x              t                          x
                                                                                                                            0 
                                                                                              [[Q]]1[[ A]][[Q]]  [[]]   1
                                         pokud je [[Q]] maticí vlastních vektorů
                                                                                                                                  
                                                                                                                             0 2 
                                         matice [[A]] bude tato matice diagonální
                                                                                                                           
Pro náš konkrétní případ jsou vlastní vektory a vlastní hodnoty matice [[A]]
            1       1                  1 a                    aez                1                    a 0 
   [[Q]]   1       1   [[Q]]1  1                                                              [[]]  
                       
                                                 [[Q]]1[ B] 
                                      2 1  a 
                                                                                     
                                                                                    1                       0  a
                                                                                                                    
            a     a                                          2                                             

W1  p  Z1  Z 2                                                         Z1    Z    a e z
                          funkcí Z1 získáme integrací rovnice                  a 1           podél charakteristiky dx/dt=-a
                                                                           t    x      2
            1
W2  v       ( Z1  Z 2 ) funkcí Z získáme integrací rovnice             Z 2   Z   aez podél charakteristiky dx/dt=a
           a                      2                                           a 2 
                                                                           t     x    2
NAP6
                MOC hydraul.ráz MATLAB
         Vraťme se k předchozímu výsledku integrace dle dvojice charakteristik
                    1 p A  pB                  fh
                vC   (          v A  vB          (v A | v A | vB | vB |))
                    2    a                   2aD
                    a p  pB                      fh
                pC  ( A         (vA  vB )        (vB | vB | vA | vA |))
                         2           a                           2aD
        který popisuje rychlosti a tlaky ve vnitřních uzlech sítě. Do okrajových uzlů se
        lze dostat jen integrací dle jedné charakteristiky. Druhou potřebnou rovnicí je
        okrajová podmínka, např. předepsaný průběh tlaků na vstupu a předepsaná
        (třeba nulová) rychlost na výstupu pro případ hydraulického rázu.

                                                                                             Zavírání ventilu lze v MATLABU
   p=p0                                                                                      popsat jako funkci času, např.
                C
                                                                                  v(t)                           function vrel=valve(t)
                                                                                                                 if t<.1
                                                                                                                        vrel=1;
                                                                                                                 else
                        B                                       A                                                       vrel=exp(-10*(t-.1));
                pC  p B    fh                                                  fh                              end
   vC  v B                  vB | vB |          pC  p A  a (vC  v A )        vA | vA |
                  a       2aD                                                  2D

 1                             fh                             1                              fh
 ( pC  pB )   (vC  vB )      vB | vB | 0                  ( pC  pA )   (vC  vA )      vA | vA | 0
 a                             2aD                             a                              2aD
NAP6
           MOC hydraul.ráz MATLAB
 Trubice L=1m, D=0.01 m, rychlost zvuku a=1 m/s, ustálená dopředná rychlost v=0.6325 m/s.
                                                                                     Průběhy tlaku (pro různé hustoty sítě), vždy je
  l=1;d=0.01;rho=1000;f=0.1;a=1;p0=2e3;
                                                                                     dobré takto testovat chybu aproximace
  v0=(p0*2*d/(l*rho*f))^0.5                                                          3
  n=101;h=l/(n-1);v(1:n)=v0;p(1)=p0;                                                                                     2000

  for i=2:n                                                                        2.5                                   1800

     p(i)=p(i-1)-f*rho*v0^2*h/(2*d);                                                                                     1600

  end
  dt=h/a;tmax=3;itmax=tmax/dt;fhr=f*h/(2*a*d);
                                                                                     2
                                                                                                      n=301              1400

                                                                                                                         1200

  for it=1:itmax                                                                   1.5
                                                                                                                         1000

     t=it*dt;                                                                                                            800
                                                                                     1
     for i=2:n-1                                                                                                         600

         pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1);                                  0.5                                   400

         pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+fhr*(vb*abs(vb)-va*abs(va)));                                                  200

         vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-fhr*(vb*abs(vb)+va*abs(va)));              0
                                                                                       0   0.2    0.4   0.6   0.8   1
                                                                                     3
     end                                                                                                                 2000
     pc(1)=p0; vb=v(2);pb=p(2);                                                                                          1800
                                                                                   2.5
     vc(1)=vb+(pc(1)-pb)/(a*rho)-fhr*vb*abs(vb);                                                                         1600
     vc(n)=v0*valve(t); va=v(n-1);pa=p(n-1);                                         2                                   1400
     pc(n)=pa-rho*a*(vc(n)-va)+f*h*rho/(2*d)*va*abs(va);              čas (do 3 s)                                       1200
     vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n);                                   1.5
                                                                                                                         1000
     p=pc;v=vc;                                                                                       n=101              800
  end                                                                                1
                                                                                                                         600
  pmax=max(max(pres))/p0                                                                                                 400
                                                                                   0.5
  x=linspace(0,1,n);                                                                                                     200
  time=linspace(0,tmax,itmax);                                                       0
                                                                                       0   0.2    0.4   0.6   0.8   1
  contourf(x,time,pres,30)
  NAP6
                       MOC hydraul.ráz MATLAB
      Vliv rychlosti zavírání ventilu na maximální přetlak v uzávěru (není příliš výrazný)
               6
                                                         2000
                       vrel=exp(-50*(t-.1));
               5
                       velmi rychlé zavření
                                                         1800

                                                         1600
                                                                v0 =0.6325 m/s          a=1 m/s                 =1000 kg/m2
n=101          4                                         1400
                                                                pmax =2180 Pa
                                                         1200
               3
                                                         1000



                                                                Přibližná analýza maximálního tlaku: kinetická
                                                         800
               2
                                                         600

               1                                         400

                                                         200
                                                                energie „válečku“ tekutiny pohybujícího se rychlostí v
               0
               6
                 0      0.2    0.4   0.6       0.8   1

                                                         2000
                                                                se po nárazu promění v deformační energii
                       vrel=exp(-10*(t-.1));             1800
               5
                       rychlé zavření cca 0.1 s          1600
                                                                                         váleček se po nárazu
čas (do 6 s)
               4                                         1400

                                                         1200
                                                                                      v zase odrazí, ale o to teď nejde)
                                                                                         zastaví a zkrátí (a pak se
               3
                                                         1000


               2
                                                         800

                                                         600
                                                                               L
               1                                         400
                                                                        1 2                        1         1 2
                                                         200
                                                                 Ek       v0 LA        Ep         pLA     p LA
               0
               6
                 0       0.2   0.4    0.6      0.8   1

                                                         1800
                                                                        2                          2        2K
                       vrel=exp(-0.1*(t-.1));

                                                                          pmax   Kv0   av0
               5                                         1600
                       pomalé zavírání cca 10s
                                                         1400
               4
                                                         1200

               3                                         1000

                                                         800                                                  rychlost zvuku
               2
                                                         600



                                                                  Pro náš případ tedy pmax=632 Pa, což je jen cca 30% skutečné hodnoty. Ta
                                                         400
               1


                                                                  analýza je jen orientační, nestlačuje se rovnoměrně celá trubice.
                                                         200

               0
                   0     0.2   0.4    0.6      0.8   1
NAP6
              MOC hydraul.ráz MATLAB
         Souvislost mezi rychlostí zvuku, rychlostí sloupce kapaliny a přetlakem p
         vyjadřuje několik století stará rovnice hydraulického rázu (Young 1808, v
         modernější podobě pak Žukovského rovnice 1898 pro ráz v elastické trubce)
         Rovnice je tak známá, že se už většinou ani neuvádí odvození.

dv rychlost pístu             a rychlost čela rázové vlny      dp  adv
                                                            Odvození vychází z hmotnostní a hybnostní bilance
                                                            kontrolního objemu, který se pohybuje konstantní
                    Předpoklad konstantní                   rychlostí rázové vlny a směrem vpravo
                     rychlosti až do místa

                                                                (   d )(a  dv)  a  ad  dv
                          rázové vlny


     v                                                                                         tok hybnosti je
                                                             rychlost tekutiny vzhledem
                                                                k rozhraní, které se          hmotnostní průtok
                                                             pohybuje rychlostí zvuku a         krát rychlost

   p,
                                                               (   d  )(a  dv)2  p  dp   a 2  p
                    Předpoklad skokového                        a 2 d   2a 2 d   dp  0  dp  a 2 d 
                    poklesu hustoty i tlaku
                      v zoně rázové vlny
                                                                                          Po dosazení rovnice bilance
                                                                                            hmotnosti plyne rovnice
                                                                                              hydraulického rázu
NAP6
       MOC hydraul.ráz MATLAB

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

                                    Blasiův vztah
function f=frict(v,d,rho,mju)
re=abs(v)*d*rho/mju;
                                    pro třecí ztráty    for it=1:itmax
if re<2100                                                 t=it*dt;
    f=64/re;                                               for i=2:n-1
else                            konstantní tlak                pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1);
                                                               fa=frict(va,d,rho,mju);
    f=0.316/re^0.25;              na vstupu                    fb=frict(vb,d,rho,mju);
end
                                                               ea=fa*h/(2*a*d)*va*abs(va);
                                                               eb=fb*h/(2*a*d)*vb*abs(vb);
p0=1e3;l=1;d=0.01;rho=1000;mju=0.001;a=4;                      pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+eb-ea);
vlam=p0*d^2/(32*mju*l);                                        vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-eb-ea);
re=vlam*d*rho/mju                                          end
if re<2100                                                 pc(1)=p0;
    v0=vlam                                                vb=v(2);pb=p(2);
    f=64/re;                                               fb=frict(vb,d,rho,mju);
else                                                                                               zadaná rychlost
                                                           eb=fb*h/(2*a*d)*vb*abs(vb);
    v0=(p0*d^(5/4)/(0.158*mju^0.25*rho^0.75*l))^(4/7)      vc(1)=vb+(pc(1)-pb)/(a*rho)-eb;            na výstupu
    re=v0*d*rho/mju;                                       vc(n)=v0*valve(t);
    f=0.316/re^0.25                                        va=v(n-1);pa=p(n-1);
end                                                        fa=frict(va,d,rho,mju);
n=401;h=l/(n-1);v(1:n)=v0;p(1)=p0;                         pc(n)=pa-rho*a*(vc(n)-va)+fa*h*rho/(2*d)*va*abs(va);
for i=2:n                                                  vres(it,1:n)=vc(1:n);
    p(i)=p(i-1)-f*rho*v0^2*h/(2*d);                        pres(it,1:n)=pc(1:n);
end                                                        p=pc;v=vc;
dt=h/a;tmax=3;itmax=tmax/dt;                            end
NAP6
                 MOC hydraul.ráz MATLAB
 Voda, a=4 m/s, L=1m, D=0.01 m, p0=1 kPa, v0=1.14 m/s, Re=31000
                 Rozložení tlaku                                                                          Vývoj rychlosti
        3                                                                                         3

                                                                                                                                              0.6
                                                                         3000
       2.5                                                                                      2.5
                                                                                                                                              0.4

                                                                         2000
        2                                                                                         2                                           0.2


                                                                         1000                                                                 0
       1.5                                                                                      1.5

                                                                                                                                              -0.2
                                                                         0
        1                                                                                         1
                                                                                                                                              -0.4
                                                                         -1000
       0.5                                                                                      0.5
                                                                                                                                              -0.6

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


                                                                         průběhy tlaku na zavíracím ventilu
                                   5000


                                   4000




                                                                                                      To, že se pulzace zrychlily, je
                                   3000


                                   2000


                                   1000
                                                                                                      způsobeno vyšší zvolenou
                                                                                                      rychlostí zvuku (tužší elastická
                                      0


                                   -1000


                                   -2000
                                                                                                      trubice) a=4 m/s. Průlet tlakové
                                                                                                      vlny pak odpovídá času 0.25 s.
                                   -3000
                                           0   0.5     1   1.5       2    2.5    3
NAP6
         MOC hydraul.ráz MATLAB
  Téměř stejným způsobem lze řešit případ, kdy je na vstupu zadaný průtok
  (objemové čerpadlo) a na výstupu je konstantní tlak (např. nulový)
  function vrel=pump(t)
  vrel=0.5*sin(3*t);

                      v(t)

                                                                                                    p=0

                                                B                                A
                                                                 fh                              pC  p A    fh
                                   pC  p B  a (vC  v B )        vB | vB |       vC  v A                  vA | vA |
                                                                 2D                                 a       2aD


  vc(1)=pump(it*dt);vb=v(2);pb=p(2);
  pc(1)=pb+rho*a*(vc(1)-vb)+f*h*rho/(2*d)*vb*abs(vb);

  va=v(n-1);pa=p(n-1); pc(n)=0;
  vc(n)=va-(pc(n)-pa)/(rho*a)-fhr*va*abs(va);
NAP6
       MOC elastická a tuhá trubice




Řešení, které je asi úplně špatně
NAP6
        MOC elastická a tuhá trubice
        Tato zdánlivě nevinná kombinace okrajových podmínek představuje pro
        numeriku problém. V tuhé trubce a pro nestlačitelnou kapalinu je rychlost
        zvuku nekonečná, charakteristiky jsou rovnoběžné s osou trubky a
        časový krok vychází nekonečně malý.



        Elastická trubice. Rychlost            Tuhá trubice. Rychlost zvuku je
        zvuku a je dána její elasticitou.      nekonečně velká.


 v(t)                                                                                       pe(t)


                         L                                                 Le

                                                   skoková změna impedance (odporu)
                                                    vyvolá odraz vln, které interferují s
                                                   vlnami, které postupují směrem toku
                                                                  vpravo
NAP6
        MOC elastická a tuhá trubice
  Nekorektní (?) okrajová podmínka na konci (podmínka s první derivací)

                                                       fh                        v p  pe f  v | v |
                         pC  p B  a (vC  v B )        vA | vA |                                 0
                                                       2D                         t   Le     2 De
                                                                           D C
                                                                                               Le

                                                                            A                                        e
                                                 1               fh                   pC  pD pC  pe f  v | v |
                                     vC  vA       ( p A  pC      vA | vA |)                                  ( De  D)  0
                                                 a              2D                       h      Le     2 DDe


       Vraťme se k původní Bernoulliho rovnici, která by měla platit i ve vyústění
       elastické trubice
                           v p f  v | v |
                                           0
                           t x       2D
        a odečtěme od ní Bernoulliho rovnici pro tuhou trubici
                            p p  pe f  v | v |
                                                ( De  D)  0
                            x   Le    2 DDe
NAP6
       MOC elastická a tuhá trubice
  Program v MATLABu
  le=0.001;de=0.1;pe=0;l=1;d=0.01;rho=1000;f=0.1;a=1;p0=0;v0=0;
  n=401;h=l/(n-1);v(1:n)=v0;p(1)=p0;
  for i=2:n
     p(i)=p(i-1)-f*rho*v0^2*h/(2*d);
  end
  dt=h/a;tmax=3;itmax=tmax/dt;
  fhr=f*h/(2*a*d);
  for it=1:itmax
     t=it*dt;
     for i=2:n-1
         pa=p(i-1);pb=p(i+1);va=v(i-1);vb=v(i+1);
         pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+fhr*(vb*abs(vb)-va*abs(va)));
         vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-fhr*(vb*abs(vb)+va*abs(va)));
     end
     vc(1)=pump(it*dt);vb=v(2);pb=p(2);
     pc(1)=pb+rho*a*(vc(1)-vb)+f*h*rho/(2*d)*vb*abs(vb);

     va=v(n-1);vb=v(n);pa=p(n-1);
     pc(n)=(pc(n-1)/h+pe/le-f*rho*vb*abs(vb)/(2*d*de)*(de-d))/(1/h+1/le);
     vc(n)=va+1/(rho*a)*(pa-pc(n)-f*h*rho/(2*d)*va*abs(va));
     vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n);
     p=pc;v=vc;
  end
  pmax=max(max(pres))
  x=linspace(0,1,n);
  time=linspace(0,tmax,itmax);
  contourf(x,time,pres,30)
          NAP6
                               MOC elastická a tuhá trubice
                 Tlakové profily v elastické trubici (a=1 m/s, L=1 m, D=De=0.01 m)
                     function vrel=pump(t)                                   3
                                                                                                                         600     3
                     vrel=0.5*sin(3*t);                                                                                                                                    600
          0.5
                                                                            2.5                                                 2.5
          0.4                                                                                                            400
                                                                                                                                                                           400
          0.3
                                                                             2                                                   2
          0.2
                                                                                                                         200                                               200
          0.1

            0                                                               1.5
                                                                                                                                1.5
                                                                                                                         0                                                 0
          -0.1

          -0.2                                                               1
                                                                                                                                 1
          -0.3                                                                                                           -200                                              -200

          -0.4
                                                                            0.5
                                                                                                                                0.5
          -0.5                                                                                                                                                             -400
                 0      0.5     1     1.5         2    2.5   3                                                           -400
                                                                                            n=401,Le=0, pmax=669                                n=401,Le=1, pmax=709
                                                                             0
                                                                                  0   0.2    0.4    0.6   0.8      1             0
                                                                                                                                      0   0.2     0.4   0.6    0.8     1

 3                                                                            3
                                                                     600                                                  600 3

                                                                                                                                                                           600
2.5                                                                         2.5
                                                                     400                                                  400 2.5

                                                                                                                                                                           400
 2                                                                            2                                           200 2
                                                                     200

                                                                                                                                                                           200
1.5                                                                         1.5                                                 1.5
                                                                                                                          0
                                                                     0
                                                                                                                                                                           0
 1                                                                            1                                           -200 1
                                                                     -200
                                                                                                                                                                           -200
0.5                                                                         0.5                                           -4000.5

                              n=401,Le=0.1, pmax=669                 -400
                                                                                            n=401,Le=0.5, pmax=669                              n=401,Le=2, pmax=803       -400

 0                                                                            0                                           -600 0
      0               0.2       0.4         0.6       0.8        1                0   0.2     0.4   0.6    0.8       1                0   0.2     0.4   0.6    0.8     1
NAP6
       MOC elastická a tuhá trubice




  Experimenty prováděné v naší laboratoři neoperovaly s harmonickým
  průběhem průtoku, ale jen s jeho dočasným zvýšením. I skutečná rychlost
  zvuku byla vyšší. A o to větší jsou problémy s numerickým řešením.
NAP6
           MOC elastická a tuhá trubice
Elastická trubice (AORTA, rychlost šíření pulzní vlny=rychlost zvuku cca 8 m/s) a
na ni napojená tuhá trubice (plexisklo). Puls průtoku na vstupu generuje
elektronický řízená SUPERPUMP (Varimex). Výstup do uklidňovací nádrže.




Ukázka experimentů z naší laboratoře (pulzní vlna způsobuje deformační vlnu aorty, a ta
je monitorována dvojicí vysokorychlostních kamer, DIC=Digital Image Correlation). Pro
řešení byla použita metoda sítí Lax Wendroff a MUSCL (Monotone Upstream Scheme for
Conservation Laws, programováno ve Fortranu). Problém je v tom, že metoda nefunguje,
když je na elastickou trubici napojená dokonale tuhá trubice (v modelech bylo nutné
modelovat tuto nástavnou trubici také jako elastickou, jen s vysokou tuhostí). Matematici
říkají, že je to tím, že v této hyperbolické rovnici nelze použít okrajovou podmínku s
derivací. Pokud mají pravdu, je následující řešení chybné…
NAP6
         MOC elastická a tuhá trubice
Lineární změna průtoku na vstupu (rampová funkce). Laminární i turbulentní režim.

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

  v0=0.35 m/s vmax=0.7 m/s                         0.18
                                                              Le=1 Re=5250 pmax = 7785 Pa
                                                                                                              6000
                      D=0.015 m         a=8 m/s    0.16

                                                   0.14                                                       4000

                                                   0.12
                                                                                                              2000
                                                    0.1
                                                                 0.18                  1                      0
                      L=0.18          Le=1 m       0.08

                                                   0.06
                                                                                                              -2000

                                                   0.04
                                                                                                              -4000
                                                   0.02
         v [m/s]
                                                     0                                                        -6000
                                                          0      0.05   0.1     0.15       0.2   0.25   0.3

  0.7                                                         Le=0.1 Re=5250 pmax = 3848 Pa
                                                   0.18
                                                                                                              3000
                                                   0.16

  0.35                                             0.14                                                       2000

                                                   0.12
                                                                                                              1000
                                                    0.1
                                                                 0.18         0.1
                                                                                                              0
                   0.03        0.06   t [s]        0.08

                                                   0.06                                                       -1000


                                                   0.04
                                                                                                              -2000

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

                         1                                                                        1

                       0.9                                                                       0.9

                       0.8                                                                       0.8
             v [m/s]




                       0.7                                                                       0.7

                       0.6                         0.18             0.01                         0.6                               0.18                0.1
                       0.5                                                                       0.5

                       0.4                                                                       0.4


                             0   0.1   0.2   0.3     0.4   0.5     0.6   0.7   0.8   0.9   1           0   0.1   0.2   0.3   0.4     0.5   0.6   0.7     0.8   0.9   1
                       400                                 t [s]                               2500

                       300                                                                     2000

                                                                                               1500
                       200
                                                                                               1000
   p [Pa] (L=0.18)




                       100
                                                                                                500
                         0
                                                                                                  0
                     -100
                                                                                                -500
                     -200
                                                                                               -1000
                     -300                                                                      -1500
                     -400                                                                      -2000
                             0   0.1   0.2   0.3     0.4   0.5     0.6   0.7   0.8   0.9   1           0   0.1   0.2   0.3   0.4     0.5   0.6   0.7     0.8   0.9   1
                                                           t [s]
NAP6
        Co je třeba si zapamatovat


       Přednáška byla věnována klasifikaci parciálních diferenciálních rovnic
       druhého řádu, se zvláštním zřetelem na hyberbolický typ. Zapamatujte si
       alespoň to, co je na následující folii
NAP6
       Co je třeba si zapamatovat
                                                             2  2    2
   Typ rovnice určují pouze koeficienty druhých derivací a 2  b      c 2  f
                                                            x   xy   y
   ty určují i rovnice tzv. charakteristik dy b  b 2  4ac
                                              
                                           dx     2a

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:9
posted:5/27/2012
language:
pages:39