Numerick anal za proces
Download
1 / 41

NUMERICKÁ ANALÝZA PROCESů - PowerPoint PPT Presentation


  • 79 Views
  • Uploaded on

NUMERICKÁ ANALÝZA PROCESů . NAP 6. Parci ální diferenciální rovnice (PDE) , klasifikace na hyperbolické, parabolické a eliptické Hyperbolic ké PDE (kmitání táhel, nosníků, ráz v potrubí) MOC-metoda charakteristik (stlačitelné proudění).

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about ' NUMERICKÁ ANALÝZA PROCESů ' - vondra


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Numerick anal za proces
NUMERICKÁ ANALÝZA PROCESů

NAP6

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

Hyperbolické PDE (kmitání táhel, nosníků, ráz v potrubí)

MOC-metoda charakteristik (stlačitelné proudění)

Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010


Pd e parci ln diferenci ln rovnice
PDE parciální diferenciální rovnice

NAP6

L.Wagner


Pde parci ln diferenci ln rovnice
PDE parciální diferenciální rovnice

NAP6

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

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

  • Parabolickérovnice(Typické pro evoluční problémy, třeba vývoj teplotního nebo koncentračního profilu v čase, ale třeba i vývoj profilů v mezní vrstvě v závislosti na vzdálenosti od vstupu. Můžeme to chápat i jako mezní případ hyperbolických rovnic s nekonečnou rychlostí šíření poruchy tlaku.)

  • Eliptickérovnice(Typické pro stacionární problémy popisu rozložení teplot, koncentrací, deformací. Často jde o problém popisovaný parabolickými rovnicemi, ale až po ustálení, tj. pro nekonečně dlouhý čas. I většina pružnostních problémů je popisována eliptickými rovnicemi. )

Příklad: vlnová rovnice kmitání pružné tyče, ráz v potrubí

Příklad: Fourierova rovnice vedení tepla

Příklad: Poissonova rovnice


Pde parci ln diferenci ln rovnice1
PDE parciální diferenciální rovnice

NAP6

Typ PDE je dán pouze koeficienty druhých derivací. Každou PDE druhého řádu můžeme napsat ve tvaru

kde a,b,c,f mohou být libovolné funkce x,y i hledaného řešení , včetně jeho prvních derivací. Koeficienty a,b,curčují rovnice tzv. charakteristik, čar závislostí y(x), které splňují rovnici

b2-4ac>0 hyperbolická rovnice (existuje dvojice reálných kořenů, tudíž dvojice charakteristik)

b2-4ac=0 parabolická rovnice (existuje jen jeden kořen a jedna charakteristika)

b2-4ac<0 eliptická rovnice (kořeny jsou imaginární a žádná reálná charakteristika neexistuje)

hyperbolická PDE parabolická PDE eliptická PDE

y

y

y

x

x

x

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


Pde hyperbolick pde
PDE Hyperbolické PDE

NAP6

Všude tam, kde se projeví setrvačné síly a stlačitelnost, popisují problém hyperbolické parciální diferenciální rovnice, charakterizované reálnými charakteristikami (jejichž vlastnosti se někdy využívají při řešení hyperbolických PDE metodou charakteristik).

Tato metoda ale není jediná, zhusta se používá metoda konečných prvků, zvláště při analýze kmitání pružných konstrukcí (nosníků, skořepin,…). Lze využít výsledky získané v předchozí přednášce.


Pde kmit n t hla

fx(x), ux(x)

PDE kmitání táhla

NAP6

Rovnici popisující kmitání táhla (elastické tyče) získáme doplněním rovnice rovnováhy o zrychlující sílu

vidíte, že je to opravdu hyperbolická rovnice, znaménka druhých derivací jsou opačná

Metoda vážených residuí založená na aproximaci

Nj v roli bázové funkce

Ni v roli váhové funkce

Integrace per partes aplikovaná na první člen

matice hmot[[M]]

budicí síly [[b]]

matice tuhosti [[K]]


Pde kmit n t hla1
PDE kmitání táhla

NAP6

Uvažujme volné kmitání, bez budících sil

což je soustava jen obyčejných diferenciálních rovnic druhého řádu. Ale s konstantními koeficienty, takže existuje analytické řešení

jehož dosazením do předchozí rovnice dostaneme soustavu algebraických rovnic pro koeficienty amplitud kmitů

Protože je soustava homogenní, bude mít netriviální (nenulové) řešení jen když je matice soustavy singulární, tj. když

Čtvercové matice [[K]] i [[M]] jsou dané a mají N řádků. Determinant je de facto polynom N-tého stupně proměnné 2a algebraická rovnice má N-kořenů, N-vlastních frekvencí. Tyto frekvence  spolu s vektory amplitud lze nalézt řešením vlastního problému

v MATLABU lze tento problém řešit funkcí [u,omega]=eig(m\k)


Pde kmit n t hla2
PDE kmitání táhla

NAP6

Pro lineární bázové funkce popisu amplitud v obecném elementu

snadno stanovíme matici tuhosti elementu (stejná jako dřív)

račte si povšimnout, že součet všech prvků se rovná hmotnosti celého elementu

stejně jako matici hmot

Matice hmot se často nahrazuje nahrazuje diagonalizovanou maticí, která odpovídá rozdělení hmotnosti táhla do jeho uzlů

Výhoda diagonalizované matice hmot spočívá v tom, že inverze [[M]] je triviální.


Pde kmit n t hla matlab
PDE kmitání táhla MATLAB

NAP6

Příklad: Soustava ocelových táhel (test pro jedno táhlo)

x=[0 1];con=[1 2];a=[1e-4];

nu=length(x);ne=length(a);n=nu;

k=zeros(n,n);

m=zeros(n,n);

for e=1:ne

[kl,ml,ig]=kloc(e,x,con,a);

k(ig(1:2),ig(1:2))=k(ig(1:2),ig(1:2))+kl(1:2,1:2);

m(ig(1:2),ig(1:2))=m(ig(1:2),ig(1:2))+ml(1:2,1:2);

end

[u,omega]=eig(m\k);

function [kl,ml,ig]=kloc(ie,x,con,a)

E=200e9;R=7000;

ae=a(ie);

i1=con(ie,1);i2=con(ie,2);

ig=[i1 i2];

le=abs(x(i1)-x(i2));

kl=E*ae/le*[1 -1;-1 1];

ml=R*ae*le/6*[2 1;1 2];

V tomto případě můžeme snadno zkontrolovat výsledek

(pro diagonalizovanou matici hmot vychází vlastní frekvence menší )


Pde kmit n nosn k
PDE kmitání nosníků

NAP6

Možná důležitější je případ příčného kmitání nosníků, např. trubek výměníku tepla, podepřených v několika místech přepážkami a na koncích vetknutých do trubkovnice. Použijte výsledky předchozí přednášky, kde spojité příčné zatížení nosníku nahradíte setrvačnými silami

matice hmot[[M]]

matice tuhosti [[K]]

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

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


Pde r z v potrub
PDE Ráz v potrubí

NAP6

Nestacionární tok stlačitelné tekutiny v potrubních sítích, fungování kardiovaskulárního oběhu, to jsou typické příklady hyperbolických rovnic.

Zvláštní případ reprezentuje hydraulický ráz, ke kterému dojde při náhlém uzavření ventilu v potrubí… projeví se náhlým vzrůstem tlaku (rázovou vlnou), který se šíří rychlostí zvuku v tekutině a stěně potrubí.

Víte jak tu rychlost spočítat? Víte, jaký maximální tlak může vzniknout?

Barbara Wagner

M. Rohani, M.H. AfsharSimulation of transient flow caused by pump failure: Point-Implicit Method ofCharacteristics

Annals of Nuclear Energy, Volume 37, Issue 12, December 2010, Pages 1742-1750


Pde r z v potrub1
PDE Ráz v potrubí

NAP6

Formulace problému:

Potrubím s proměnným průřezem A(t,x), protéká tekutina střední rychlostí v(t,x).Tekutina je stlačitelná a souvislost její hustoty a tlaku je charakterizována objemovým modulem stlačitelnosti K [Pa]

Śouvisí s rychlostí zvuku v tekutině a

Je to vlastně stejný problém jako dříve uvažované potrubní sítě, kdy cílem bylo stanovit průběhy tlaku a průtoky. Tentokrát však tlaky i rychlosti závisí i na čase, p(t,x), v(t,x).

Tyto veličiny popisuje stejně jako ve stacionárním případě dvojice rovnic, rovnice kontinuity a bilance hybnosti (Bernoulliho rovnice). Ve zjednodušeném tvaru

Rovnici kontinuity v této poněkud zvláštní podobě vysvětlím na následující folii

takto se zohlední i pružnost potrubí (závislost průřezu na tlaku)

Tuto rovnici poznáváte: Bernoulliho rovnice kde je zanedbán konvektivní člen (kinetická energie)


Pde r z v potrub2

,A,v

dx

PDE Ráz v potrubí

NAP6

Rovnice kontinuity pro případ, kdy průřezová plocha trubkyA(t,x) se mění po délce i v čase (elastická, nebo elasticko-plastická trubka, která se „nafukuje“ působením vnitřního přetlaku)

hmotnostní bilance elementu konstantní délky dx

akumulace

mění se hustota, průřez i střední rychlost .

souvislost hustoty a tlaku

přítok-odtok

Rozderivováním levé i pravé strany a využitím definice objemového modulu stlačitelnosti získáme zcela obecnou rovnici kontinuity

(to je totéž, jen zapsané v materiálových derivacích)

Tento člen je v předchozí rovnici zanedbán (dělá se to často, viz např.Khamlichi, Wave motion 1995)


Pde r z v potrub3
PDE Ráz v potrubí

NAP6

To, že je výše uvedený systém rovnic hyperbolický, vyplyne z eliminace tlaku nebo rychlosti (stačí derivovat rovnici kontinuity dle času, Bernoulliho rovnici dle souřadnice x, a odečíst). Výsledné rovnice pro tlak p(t,x)

nebo rychlost v(t,x)

Rychlost zvuku

jsou stejné a dokonce stejné jako hyperbolická rovnice kmitání táhla. Je tedy možné použít i stejnou metodu řešení (metodu vážených residuí).

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


PDE

NAP6

Pro hyperbolické rovnice se ale často využívá jiná metoda, která využívá specifickou vlastnost hyperbolických rovnic a tou je existence dvou reálných charakteristik.

Nazývá se metoda charakteristik (MOC).

Doporučená literatura Wylie, Streeter: Fluid Transients. McGraw Hill, 1978

Zábranský


Metoda charakteristik moc
Metoda charakteristik MOC

NAP6

Co vlastně chceme řešit?

Proudění stlačitelné tekutiny v trubce nebo nestlačitelné kapaliny v elastické trubce (která se může „nafukovat“ při zvýšení tlaku) – výsledný efekt stlačitelnosti je úplně stejný. Ve hře jsou dva vlivy: setrvačnost tekutiny a tlak vyvolaný jejím stlačením. Cílem je stanovení střední rychlosti v (tudíž i průtoku q) a tlaku p, přičemž obě tyto veličiny (v,p) se mění v čase i po délce trubky.

Budeme uvažovat různé okrajové podmínky na koncích trubky: buď zadávaný (a časově proměnný) průtok, nebo zadávaný (a časově proměnný) tlak.


Metoda charakteristik moc1
Metoda charakteristik MOC

t

x() t=

x

NAP6

Vraťme se k výchozí soustavě rovnice kontinuity a rovnice Bernoulliho

spoustu věcí jsme zanedbali, například unášivou rychlost, takže následující analýza je správná jen tehdy, když rychlost zvuku je mnohem vyšší než rychlost proudění. V řadě uváděných příkladů tento předpoklad splněn není.

a tyto rovnice (první z nich vynásobme libovolnou nenulovou konstantou ) sečteme

Zvolme libovolnou křivku v rovině t-x (parametrem křivky x() může být přímo čas t). Úplné diferenciály tlaku a rychlosti podél této křivky jsou

Jestliže bude křivka x(t=) splňovat rovnice (současně)

získáme finální rovnici vyjádřenou pomocí úplných diferenciálů podél křivky x(t=)

(a tuto křivku nazveme charakteristikou diferenciální rovnice)


Metoda charakteristik

C

x2=-at

x1=at

h/a

B

A

h

Metoda charakteristik

NAP6

Předchozí rovnice má dvě řešení

a jim odpovídají dvě charakteristické křivky a diferenciální rovnice, které je třeba na nich integrovat

Třecí ztráty vyjádřené součinitelem třecích ztrát f

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

např.

přibližnost integrace spočívá jen v zanedbání proměnnosti třecích ztrát podél charakteristiky


Metoda charakteristik

NAP6

Řešení této soustavy vyjádříme v explicitním tvaru

a můžeme ho ihned využít pro numerické řešení rychlostí a tlaků na časové hladině C ze známých hodnot v bodech A,B na staré časové hladině.


Metoda charakteristik1
Metoda charakteristik

NAP6

Alternativní (a možná názornější) formulace je analogií postupu, který jsme použili při řešení soustavy diferenciálních rovnic výměníků tepla, viz Whitham: Linear and nonlinear waves, Wiley 1974. Soustavu dvou diferenciálních rovnic pro tlak a rychlost zapíšeme v maticovém tvaru

a použijeme transformaci vektoru W takovou, že se ve výsledných rovnicích objeví vždy jen jedna neznámá

pokud je [[Q]] maticí vlastních vektorů matice [[A]] bude tato matice diagonální

Pro náš konkrétní případ jsou vlastní vektory a vlastní hodnoty matice [[A]]

funkcí Z1 získáme integrací rovnice

podél charakteristikydx/dt=-a

podél charakteristikydx/dt=a

funkcí Z2 získáme integrací rovnice


Moc hydraul r z matlab
MOC hydraul.ráz MATLAB

NAP6

Následující příklad ukazuje řešení konkrétního problému stlačitelného proudění v trubce, kde na vstupu je zadaný tlak (např. konstantní p0) a na výstupu je zadaný časový průběh průtoku, např. rychle (téměř skokově) klesající průtok, což je situace odpovídající rychlému zavírání ventilu.

Otázkou je především to, jaké zvýšení tlaku uzavírání ventilu způsobí (hydraulický ráz v potrubí).


Moc hydraul r z matlab1
MOC hydraul.ráz MATLAB

C

x2=-at

x1=at

h/a

B

A

h

NAP6

Vraťme se k předchozímu výsledku integrace dle dvojice charakteristik

který popisuje rychlosti a tlaky ve vnitřních uzlech sítě. Do okrajových uzlů se lze dostat jen integrací dle jedné charakteristiky. Druhou potřebnou rovnicí je okrajová podmínka, např. předepsaný průběh tlaků na vstupu a předepsaná (třeba nulová) rychlost na výstupu pro případ hydraulického rázu.

Zavírání ventilu lze v MATLABU popsat jako funkci času, např.

function vrel=valve(t)

if t<.1

vrel=1;

else

vrel=exp(-10*(t-.1));

end

C

v bodě C je tlak daný okrajovou podmínkou, dopočítá se jen rychlost integrací přes charakteristiku vlevo

p=p0

v(t)

B

A


Moc hydraul r z matlab2

Průběhy tlaku (pro různé hustoty sítě), vždy je dobré takto testovat chybu aproximace

n=301

n=101

čas (do 3 s)

MOC hydraul.ráz MATLAB

NAP6

Trubice L=1m, D=0.01 m, rychlost zvuku a=1 m/s, ustálená dopředná rychlost v=0.6325 m/s.

l=1;d=0.01;rho=1000;f=0.1;a=1;p0=2e3;

v0=(p0*2*d/(l*rho*f))^0.5

n=101;h=l/(n-1);v(1:n)=v0;p(1)=p0;

for i=2:n

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

end

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

for it=1:itmax

t=it*dt;

for i=2:n-1

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

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

vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-fhr*(vb*abs(vb)+va*abs(va)));

end

pc(1)=p0; vb=v(2);pb=p(2);

vc(1)=vb+(pc(1)-pb)/(a*rho)-fhr*vb*abs(vb);

vc(n)=v0*valve(t); va=v(n-1);pa=p(n-1);

pc(n)=pa-rho*a*(vc(n)-va)+f*h*rho/(2*d)*va*abs(va);

vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n);

p=pc;v=vc;

end

pmax=max(max(pres))/p0

x=linspace(0,1,n);

time=linspace(0,tmax,itmax);

contourf(x,time,pres,30)

V čase t=1.1 už dorazila zpětná vlna na vstup

V čase t=0.1 začíná zavírat výstupní ventil


Moc hydraul r z matlab3

v dobré takto testovat chybu aproximaceáleček se po nárazu zastaví a zkrátí (a pak se zase odrazí, ale o to teď nejde)

MOC hydraul.ráz MATLAB

NAP6

Vliv rychlosti zavírání ventilu na maximální přetlak v uzávěru (není příliš výrazný)

vrel=exp(-50*(t-.1)); velmi rychlé zavření

v0 =0.6325 m/s a=1 m/s=1000 kg/m2

pmax =2180 Pa

n=101

Přibližná analýza maximálního tlaku: kinetická energie „válečku“ tekutiny pohybujícího se rychlostí v se po nárazu promění v deformační energii

vrel=exp(-10*(t-.1)); rychlé zavření cca 0.1 s

v

čas (do 6 s)

L

vrel=exp(-0.1*(t-.1)); pomalé zavírání cca 10s

To odpovídá Žukovského rovnici, viz další folie

rychlost zvuku

Pro náš případ tedy pmax=632 Pa, což je jen cca 30% skutečné hodnoty. Ta analýza je jen orientační, nestlačuje se rovnoměrně celá trubice.


Moc hydraul r z matlab4
MOC hydraul.ráz MATLAB dobré takto testovat chybu aproximace

NAP6

Souvislost mezi rychlostí zvuku, rychlostí sloupce kapaliny a přetlakem p vyjadřuje několik století stará rovnice hydraulického rázu (Young 1808, v modernější podobě pak Žukovského rovnice 1898 pro ráz v elastické trubce)

Rovnice je tak známá, že se už většinou ani neuvádí odvození.

dv rychlost pístu

arychlost čela rázové vlny

Odvození vychází z hmotnostní a hybnostní bilance kontrolního objemu, který se pohybuje konstantní rychlostí rázové vlny a směrem vpravo

kontrolní objem pohybující se rychlostí zvuku vpravo

Předpoklad konstantní rychlosti až do místa rázové vlny

bilance hmoty

v

tok hybnosti je hmotnostní průtok krát rychlost

rychlost tekutiny vzhledem k rozhraní, které se pohybuje rychlostí zvuku a

p,

bilance hybnosti

Předpoklad skokového poklesu hustoty i tlaku v zoně rázové vlny

Po dosazení rovnice bilance hmotnosti plyne rovnice hydraulického rázu


Moc hydraul r z matlab5
MOC hydraul.ráz MATLAB dobré takto testovat chybu aproximace

NAP6

Tyto černé stránky můžete bez obav přeskočit. Jsou jen ukázkou toho, jak numerický model zrealističtit (výpočtem realistického součinitele třecích ztrát) a současně demonstrovat skutečnost, že pro numerické výpočty se MATLAB zase tak moc nehodí, je to totiž jen interpret a tudíž pomalý.

Následující program byl beze změny přepsán do FORTRANU (místo MATLABovských cyklů for i=2:n-1, se napíše do i=2,n-1, místo rozhodovacích příkazů typu if re<2100se napíše if(re.lt.2100)then atd.).

Výpočet následující úlohy zavírání ventilu trval na mém počítači 2s ve FORTRANu a 125 s v MATLABu (cca 60 krát pomalejší).


Moc hydraul r z matlab6
MOC hydraul.ráz MATLAB dobré takto testovat chybu aproximace

NAP6

Program uvažující závislost součinitele třecích ztrát na Re, a identifikující laminární/turbulentní režim toku. Viz funkce frict.

Blasiův vztah pro třecí ztráty

function f=frict(v,d,rho,mju)

re=abs(v)*d*rho/mju;

if re<2100

f=64/re;

else

f=0.316/re^0.25;

end

for it=1:itmax

t=it*dt;

for i=2:n-1

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

fa=frict(va,d,rho,mju);

fb=frict(vb,d,rho,mju);

ea=fa*h/(2*a*d)*va*abs(va);

eb=fb*h/(2*a*d)*vb*abs(vb);

pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+eb-ea);

vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-eb-ea);

end

pc(1)=p0;

vb=v(2);pb=p(2);

fb=frict(vb,d,rho,mju);

eb=fb*h/(2*a*d)*vb*abs(vb);

vc(1)=vb+(pc(1)-pb)/(a*rho)-eb;

vc(n)=v0*valve(t);

va=v(n-1);pa=p(n-1);

fa=frict(va,d,rho,mju);

pc(n)=pa-rho*a*(vc(n)-va)+fa*h*rho/(2*d)*va*abs(va);

vres(it,1:n)=vc(1:n);

pres(it,1:n)=pc(1:n);

p=pc;v=vc;

end

konstantní tlak na vstupu

p0=1e3;l=1;d=0.01;rho=1000;mju=0.001;a=4;

vlam=p0*d^2/(32*mju*l);

re=vlam*d*rho/mju

if re<2100

v0=vlam

f=64/re;

else

v0=(p0*d^(5/4)/(0.158*mju^0.25*rho^0.75*l))^(4/7)

re=v0*d*rho/mju;

f=0.316/re^0.25

end

n=401;h=l/(n-1);v(1:n)=v0;p(1)=p0;

for i=2:n

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

end

dt=h/a;tmax=3;itmax=tmax/dt;

zadaná rychlost na výstupu


Moc hydraul r z matlab7
MOC hydraul.ráz MATLAB dobré takto testovat chybu aproximace

NAP6

Voda, a=4 m/s, L=1m, D=0.01 m, p0=1 kPa, v0=1.14 m/s, Re=31000

Vývoj rychlosti

Rozložení tlaku

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

To, že se pulzace zrychlily, je způsobeno vyšší zvolenou rychlostí zvuku (tužší elastická trubice) a=4 m/s. Průlet tlakové vlny pak odpovídá času 0.25 s.


Moc hydraul r z matlab8
MOC hydraul.ráz MATLAB dobré takto testovat chybu aproximace

NAP6

Téměř stejným způsobem lze řešit případ, kdy je na vstupu zadaný průtok (objemové čerpadlo) a na výstupu je konstantní tlak (např. nulový)

function vrel=pump(t)

vrel=0.5*sin(3*t);

v(t)

p=0

B

A

vc(1)=pump(it*dt);vb=v(2);pb=p(2);

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

va=v(n-1);pa=p(n-1); pc(n)=0;

vc(n)=va-(pc(n)-pa)/(rho*a)-fhr*va*abs(va);


MOC elasti dobré takto testovat chybu aproximacecká a tuhá trubice

NAP6

Řešení, které je asi úplně špatně


Moc elastick a tuh trubice
MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

NAP6

Uvažujme případ toku nestlačitelné tekutiny elastickou trubicí (např. cévou) na níž je napojena trubice (třeba se stejným průměrem), ale dokonale tuhá. Tato zdánlivě nevinná kombinace okrajových podmínek představuje pro numeriku problém. V tuhé trubce a pro nestlačitelnou kapalinu je rychlost zvuku nekonečná, charakteristiky jsou rovnoběžné s osou trubky a časový krok vychází nekonečně malý.

Tuhá trubice. Rychlost zvuku je nekonečně velká.

Elastická trubice. Rychlost zvuku a je dána její elasticitou.

v(t)

pe(t)

L

Le

skoková změna impedance (odporu) vyvolá odraz vln, které interferují s vlnami, které postupují směrem toku vpravo


Moc elastick a tuh trubice1
MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

Le

NAP6

Nekorektní (?) okrajová podmínka na konci (podmínka s první derivací)

D

C

A

e

Vraťme se k původní Bernoulliho rovnici, která by měla platit i ve vyústění elastické trubice

a odečtěme od ní Bernoulliho rovnici pro tuhou trubici


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

NAP6

Program v MATLABu

le=0.001;de=0.1;pe=0;l=1;d=0.01;rho=1000;f=0.1;a=1;p0=0;v0=0;

n=401;h=l/(n-1);v(1:n)=v0;p(1)=p0;

for i=2:n

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

end

dt=h/a;tmax=3;itmax=tmax/dt;

fhr=f*h/(2*a*d);

for it=1:itmax

t=it*dt;

for i=2:n-1

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

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

vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-fhr*(vb*abs(vb)+va*abs(va)));

end

vc(1)=pump(it*dt);vb=v(2);pb=p(2);

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

va=v(n-1);vb=v(n);pa=p(n-1);

pc(n)=(pc(n-1)/h+pe/le-f*rho*vb*abs(vb)/(2*d*de)*(de-d))/(1/h+1/le);

vc(n)=va+1/(rho*a)*(pa-pc(n)-f*h*rho/(2*d)*va*abs(va));

vres(it,1:n)=vc(1:n); pres(it,1:n)=pc(1:n);

p=pc;v=vc;

end

pmax=max(max(pres))

x=linspace(0,1,n);

time=linspace(0,tmax,itmax);

contourf(x,time,pres,30)


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

function vrel=pump(t)

vrel=0.5*sin(3*t);

n=401,Le=2, pmax=803

n=401,Le=1, pmax=709

n=401,Le=0, pmax=669

n=401,Le=0.1, pmax=669

n=401,Le=0.5, pmax=669

NAP6

Tlakové profily v elastické trubici (a=1 m/s, L=1 m, D=De=0.01 m)


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

NAP6

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


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

NAP6

Elastická trubice (AORTA, rychlost šíření pulzní vlny=rychlost zvuku cca 8 m/s) a na ni napojená tuhá trubice (plexisklo). Puls průtoku na vstupu generuje elektronický řízená SUPERPUMP (Varimex). Výstup do uklidňovací nádrže.

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


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

NAP6

Lineární změna průtoku na vstupu (rampová funkce). Laminární i turbulentní režim.

le=0.1;de=0.015;pe=0;mju=0.001;l=.18;d=0.015;rho=1000;a=8;

v0=pump(0);ve=v0*(d/de)^2;f=frict(v0,d,rho,mju);fe=frict(ve,de,rho,mju);re=v0*d*rho/mju

dpe=0.5*fe*le/de*rho*ve^2;dp=0.5*f*l/d*rho*v0^2;

n=401;h=l/(n-1);v(1:n)=v0;p(1)=pe+dp+dpe;

for i=2:n

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

end

dt=h/a;tmax=3;itmax=tmax/dt;

for it=1:itmax

t=it*dt;

for i=2:n-1

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

fa=frict(va,d,rho,mju);fb=frict(vb,d,rho,mju);

ea=fa*h/(2*a*d)*va*abs(va); eb=fb*h/(2*a*d)*vb*abs(vb);

pc(i)=a/2*((pa+pb)/a+rho*(va-vb)+eb-ea);

vc(i)=0.5*((pa-pb)/(rho*a)+va+vb-eb-ea);

end

vc(1)=pump(it*dt);vb=v(2);pb=p(2); fb=frict(vb,d,rho,mju);

pc(1)=pb+rho*a*(vc(1)-vb)+fb*h*rho/(2*d)*vb*abs(vb);

va=v(n-1);vb=v(n);pa=p(n-1); fb=frict(vb,d,rho,mju); fa=frict(va,d,rho,mju);

pc(n)=(pc(n-1)/h+pe/le-fb*rho*vb*abs(vb)/(2*d*de)*(de-d))/(1/h+1/le);

vc(n)=va+1/(rho*a)*(pa-pc(n)-fa*h*rho/(2*d)*va*abs(va));

vres(it,1:n)=vc(1:n);

pres(it,1:n)=pc(1:n);

p=pc;v=vc;

end

function vrel=pump(t)

if t<.1

vrel=0.1;

elseif t<1

vrel=t;

else

vrel=1;

end

function f=frict(v,d,rho,mju)

re=abs(v)*d*rho/mju;

if re<1

f=64;

elseif re<2100

f=64/re;

else

f=0.316/re^0.25;

end

pokud je někde chyba, tak asi tady


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

Le=1 Re=5250 pmax = 7785 Pa

v0=0.35 m/s vmax=0.7 m/s

a=8 m/s

D=0.015 m

Le=1 m

L=0.18

Le=0.1 Re=5250 pmax = 3848 Pa

v [m/s]

0.7

0.35

0.06 t [s]

0.03

NAP6

Lineární změna průtoku na vstupu (rampová funkce). Laminární i turbulentní režim.

1

0.18

0.18

0.1


MOC elastická a tuhá trubice dobré takto testovat chybu aproximace

NAP6

Rychlosti a tlak na konci elastické trubice (snížení frekvence při prodloužení trubky, zdá se mi to logické, zvětšila se setrvačná hmota a nezměnila elasticita)

(L=0.18, Le=0.01, pmax=2946 Pa)(L=0.18, Le=0.1, pmax=3843 Pa).

0.1

0.01

0.18

0.18


Co je třeba si zapamatovat dobré takto testovat chybu aproximace

NAP6

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


Co je třeba si zapamatovat dobré takto testovat chybu aproximace

NAP6

Typ rovnice určují pouze koeficienty druhých derivací ty určují i rovnice tzv. charakteristik

Charakteristiky jsou reálné, když b2-4ac>0 pak je rovnice hyperbolická

(vlnová rovnice, kmitání)

Charakteristiky je jen jedna, když b2-4ac=0 pak je rovnice parabolická

(třeba rovnice vedení tepla)

Reálné charakteristiky neexistují, když b2-4ac<0 pak je rovnice eliptická

(Poissonova rovnice, statika, pružnost)


ad