Numerick anal za proces
This presentation is the property of its rightful owner.
Sponsored Links
1 / 37

NUMERICKÁ ANALÝZA PROCESů PowerPoint PPT Presentation


  • 68 Views
  • Uploaded on
  • Presentation posted in: General

NUMERICKÁ ANALÝZA PROCESů. NAP 8. Metoda sítí s aplikací na parabolick é a eliptické rovnice Metodicky de facto totéž, co v minulé přednášce jen jiné příklady: Neustálený rychlostní profil při oscilačním toku v trubce (parabolická rovnice pro cylindrický s.s.)

Download Presentation

NUMERICKÁ ANALÝZA PROCESů

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ů

NAP8

Metoda sítí s aplikací na parabolické a eliptické rovnice

Metodicky de facto totéž, co v minulé přednášce jen jiné příklady:

Neustálený rychlostní profil při oscilačním toku v trubce (parabolická rovnice pro cylindrický s.s.)

Teplotní pole (parabolická rovnice) a elektrické pole (eliptická rovnice) při přímém ohmickém ohřevu.

Identifikace parametrů modelu popisovaného metodou konečných diferencí použitím algoritmu Nelder Mead (fminsearch)

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


Oscila n tok v trubce

Oscilační tok v trubce

NAP8

Při řešení neustáleného toku stlačitelné tekutiny v trubce nebyl ověřován radiální rychlostní profil. Tiše jsme předpokládali, že rychlostní profil odpovídá ustálenému toku (a že se třeba v laminárním režimu vyvine parabolický profil). Při rychle se měnícím axiálním gradientu tlaku to tak není. Mírou odchylky je frekvence tlakových pulzací a Womersleyho číslo

Mordechai Ardon


Oscila n tok v trubce1

Oscilační tok v trubce

NAP8

Prototypem oscilační toku je například tok krve v cévním systému. Pulzace tlaku vyvolává srdce (systola-diastola, opakující se s frekvencí cca 1Hz) a v cévách s větším průměrem (aorta) dochází k tomu, že krev u stěny teče opačným směrem než někde uprostřed. Je tomu tak proto, že největší rychlost a tedy i kinetická energie je největší u osy trubky, takže setrvačnost kapaliny překoná i opačně působící gradient tlaku při jeho náhlé změně. Naproti tomu, u stěny, kde je rychlost menší, nemusí gradient tlaku překonávat tak velkou hybnost kapaliny, a směr toku bude víceméně poslušně sledovat směr gradientu tlaku.

Když zanedbáme stlačitelnost kapaliny a tuhou trubku, popisuje celý problém Navierova Stokesova rovnice pro složku rychlosti kapaliny ve směru osy trubky, zapsaná v cylindrickém souřadném systému

V nekonečně dlouhé trubici musí být příčná složka rychlosti nulová a z NS rovnice v radiálním směru plyne, že gradient tlaku dp/dr=0, jinými slovy, tlak p nezávisí na poloměru. Axiální gradient tlaku pak může být pouze funkcí času t.


Oscila n tok v trubce2

Oscilační tok v trubce

NAP8

Pokud je gradient tlaku konstantní, vyvine se po jisté době parabolický rychlostní profil

Pokud je gradient tlaku harmonicky oscilující funkcí času, budou i rychlosti proudění harmonickou funkcí času, se stejnou frekvencí , jen fázově posunuté. Radiální rychlostní profil je pro tento speciální případ vyjádřen Besselovou funkcí J0 s komplexním argumentem

kde r je bezrozměrný poloměr,  je úhlová frekvence sinusových oscilací gradientu tlaku a Wo je Womersleyho číslo

Když je Wo<<1 můžeme předpokládat, že rychlostní profil je parabolický a docela přesně sleduje časové (pomalé) změny gradientu tlaku. Pak ani není nutné použít numerické řešení a třeba i změny průtoku lze počítat z běžné formy Bernoulliho rovnice.

Toto analytické řešení uvádím spíš jen pro zajímavost (v praxi se využívá při modelování toku krve v artériích), v dalším textu se zaměříme na numerické řešení, které je použitelné pro zcela libovolný časový průběh axiálního gradientu tlaku.


Oscila n tok v trubce3

r <0,1> r

1 2 i-1 i i+1 N

Oscilační tok v trubce

NAP8

Navier Stokesovu rovnici upravme zavedením bezrozměrného poloměru a času,

Diferenční přepis této rovnice na ekvidistantní síti uzlových bodů (pro jeden časový krok )

A to je soustava lineárních algebraických rovnic s tridiagonální maticí soustavy (vektor bi je diagonála, di je vektor pravé strany)

Na ose symetrie

Na stěně

Zkuste si to odvodit: na ose je problém s dělením poloměrem r1 – musíte použít limitní přechod (třeba l’Hopitalovo pravidlo)


Oscila n tok v trubce4

u>0

u<0

Oscilační tok v trubce

NAP8

Tato tridiagonální soustava algebraických rovnic se dá velmi efektivně řešit faktorizací, algoritmizovanou např. v MATLABovské funkci tridag

q()-průtok

function u = tridag(a,b,c,r,n)

bet=b(1);

u(1)=r(1)/bet;

for j=2:n

gam(j)=c(j-1)/bet;

bet=b(j)-a(j)*gam(j);

u(j)=(r(j)-a(j)*u(j-1))/bet;

end

for j=n-1:-1:1

u(j)=u(j)-gam(j+1)*u(j+1);

end

Úplný výpis MATLABovského programu bude uveden až za okamžik. Zde je jen ukázka typického časového průběhu průtoku, který jak patrno není ve fázi se zadaným harmonickým průběhem gradientu tlaku.

Druhý obrázek (radiální rychlostní profil) dokumentuje existenci zpětného toku v ose trubky – pro harmonicky oscilující gradient tlaku by to měla být Besselova funkce J0.


Bernoulliho rovnice

Bernoulliho rovnice

NAP8

Z praktického hlediska má asi největší význam posouzení toho, do jaké míry ovlivní změna rychlostního profilu Bernoulliho rovnici. Otázka: Mohu použít následující tvar Bernoulliho rovnice pro výpočet časově proměnné střední rychlosti (nebo průtoku)?

Rovnici si odvoďte z rovnováhy sil působících na váleček kapaliny o poloměru R a tloušťce dx.

Jistá nesrovnalost takto napsané Bernoulliho rovnice je v tom, že zrychlení se odvozuje z předpokladu konstantního radiálního profilu rychlosti, zatímco třecí ztráty z představy parabolického profilu rychlosti.

Tatáž Bernoulliho rovnice, ale zapsaná s bezrozměrným časem a gradientem tlaku

ověřte, že to je analytické řešení (odvodí se metodou variace konstant z řešení homogenní rovnice)

Odchylka predikovaných rychlostí nebo průtoků závisí na tom, jak velké jsou oscilace tlaku; srovnání s přesným (numerickým) řešením ukazuje, že amplituda průtoků je dle Bernoulliho rovnice, založené na tlakových ztrátách odpovídajících parabolickému rychlostního profilu, zhruba o 10% vyšší (parabolický profil má menší tlakovou ztrátu). Lepší přesnosti se dá dosáhnout aproximací rychlostního profilu polynomem 4 stupně, viz následující text...


P ibli n e en ns rovnice

Přibližné řešení NS rovnice

NAP8

Vyjděme z Navier Stokesovy rovnice ve tvaru použitém pro numerické řešení (s bezrozměrným časem a radiální souřadnicí)

a hledejme přibližné řešení ve tvaru, který zajistí automatické splnění okrajových podmínek v ose a na stěně

(1)

Tento rychlostní profil by měl v každém časovém kroku splnit počáteční podmínky (pro čas =0), vyjádřené např. dvojicí hodnot c1(0), c2(0), a také Navier Stokesovu rovnici

a to pro libovolné hodnoty r a pro libovolnou zadanou funkci P()


P ibli n e en ns rovnice1

Přibližné řešení NS rovnice

NAP8

Porovnáním koeficientů u mocnin poloměru r (r0, r2, r4) získáme soustavu obyčejných diferenciálních rovnic

Teď ovšem máme problém, protože jsme získali soustavu tří diferenciálních rovnic pro pouze dvě neznámé c1,c2. A ty tři rovnice jsou lineárně nezávislé, takže přesné řešení (c1(),c2()) asi existovat nebude, můžeme jen hledat přijatelný kompromis.

  • Lze zkusit například tyto varianty:

  • c2=0 (ignorujeme druhou a třetí rovnici), řešení

  • ignorujeme třetí rovnici a řešíme soustavu 2 rovnic metodou vlastních čísel

  • eliminujeme c2 z prvních dvou rovnic předpokládajíce přibližnou platnost třetí rovnice (přibližně konstantní c2 v řešeném časovém kroku)

  • V následujícím textu rozebereme variantu c), protože je nejjednodušší:

Eliminace c2() z první rovnice

a jejím řešením je

pro konstantní gradient tlaku P

Funkce c2() plyne z algebraické rovnice vzniklé sečtením předchozích 3 rovnic.

Limitní případ pro  má řešení c1()=-P/4 a c2()=0 což správně odpovídá plně vyvinutému parabolickému rychlostnímu profilu.


V en residua p ibli n e en

Vážená residua–přibližné řešení

NAP8

Na těchto šedivých stránkách je uvedeno alternativní řešení vypočtu koeficientů c1(t), c2(t) metodou vážených reziduí.

Řešení metodou vážených reziduí je možná lepší (konzistentnější) než předchozí heuristické řešení.


V en residua p ibli n e en1

Vážená residua–přibližné řešení

NAP8

Dosazením získáme residuum NS rovnice

Funkce c1 c2 musí splňovat podmínku nulového váženého rezidua alespoň pro váhové funkce w1(r)=1 a w2(r)=r, tj. rovnice

V maticovém zápisu

A to už je standardní formulace počátečního problému pro 2 rovnice


V en residua p ibli n e en2

Vážená residua–přibližné řešení

NAP8

Z praktického hlediska je možná lepší řešit přímo průtoky

Po dosazení původních rozměrových veličin a doplnění zrychlení, které působí proti směru osy x, získáme soustavu dvou obyčejných diferenciálních rovnic pro průtoky

Pro konstantní gradient tlaku lze nalézt řešení evoluce vývoje průtoků v analytickém tvaru (řešením vlastních vektorů a vlastních čísel) jako lineární kombinaci exponenciálních funkcí.


V en residua p ibli n e en3

Vážená residua–přibližné řešení

NAP8

Speciální případ kdy gradient tlaku je konstantní a řešení je ustálené

Je zřejmé, že řešení se redukuje na plně vyvinutý parabolický profil, stěnová složka mizí, a řešení splňuje Hagen Poisseuille zákon pro laminární tok v trubce s kruhovým průřezem.


Bernoulliho rovnice1

Bernoulliho rovnice

NAP8

Z rychlostního profilu získáme integrací průtok nebo střední objemovou rychlost

Tento výsledek je možné srovnat s analogickým řešením Bernoulliho rovnice

kde exponent m=8 odpovídá konstantnímu, zatímco m=4 parabolickému radiálnímu rychlostnímu profilu v nestacionárním členu (du/dt).

Porovnání výsledků přesného numerického řešení, Bernoulliho rovnice a aproximace radiálního rychlostního profilu polynomem 4-tého stupně je naprogramováno jako MATLABovský program…


Bernoulliho rovnice2

rychlostní profily pro stejnou

střední objemovou rychlost=1 m/s

bezrozměrný poloměr

Bernoulliho rovnice

NAP8

Na počátku každého časového kroku (v čase =0) je obecně třeba zadat kompletní radiální rychlostní profil jako počáteční podmínku. Za předpokladu, že rychlostní profil je parabolický, stačí zadat jen průtok (nebo střední objemovou rychlost), v případě s rychlostním profilem 4 stupně musí být podmínky dvě, např. c1(0) a c2(0), nebo střední objemová rychlost a relativní velikost stěnového průtoku uw, z nichž lze koeficienty c1,c2 snadno spočítat


Oscila n tok v trubce5

Oscilační tok v trubce

NAP8

Zadaný průběh gradientu tlaku , počítají se časové průběhy průtoků qfd (Finite Dif.), qb (Bernoulli), qode (Ordinary Dif.Eq.), včetně matic rychlostních profilů

for i=1:nt

%finite differences

for ir=1:nr-1

d(ir)=-p(i)+u(ir)/dt;

end

u(1:nr)=tridag(a,b,c,d,nr);

q=0;

for ir=1:nr

q=q+u(ir)*r(ir);

end

q=2*3.141*dr*q;

qfd(i)=q;

ufd(i,1:nr)=u(1:nr);

%bernoulli

umean0=(umean0+p(i)/8)*exp(-8*dt)-p(i)/8;

qb(i)=umean0*3.141;

for ir=1:nr

ub(i,ir)=2*umean0*(1-r(ir)^2);

end

%ODE

c1(i)=(c10+p(i)/4)*exp(-16*dt/3)-p(i)/4;

c2(i)=(-p(i)-4*c1(i))/12;

qode(i)=3.141/2*(c1(i)+c2(i)/3);

c10=c1(i);

for ir=1:nr

uode(i,ir)=c1(i)*(1-r(ir)^2)+c2(i)*(r(ir)^2-r(ir)^4);

end

end

nt=1501;

nr=21;

a=zeros(nr);b=zeros(nr);c=zeros(nr);d=zeros(nr);u=zeros(nr);

ufd=zeros(nt,nr);uode=zeros(nt,nr);ub=zeros(nt,nr);

qfd=zeros(nt);qb=zeros(nt);qode(nt)=0;t=zeros(nt);

tmax=3;dt=tmax/(nt-1);t=linspace(0,tmax,nt);

% defince casoveho prubehu gradientu tlaku

p(1:nt)=-1-3*sin(10*t);

r=linspace(0,1,nr);dr=1/(nr-1);dr2=dr^2;

%pocatecni rychlostni profil stabilizovany

for i=1:nr

u(i)=-p(1)/4*(1-r(i)^2);

end

umean0=-p(1)/8;

c10=-p(1)/4;

for i=2:nr-1

a(i)=-(r(i)+r(i-1))/(2*r(i)*dr2);

c(i)=-(r(i+1)+r(i))/(2*r(i)*dr2);

b(i)=1/dt-a(i)-c(i);

end

b(1)=1/dt+4/dr2;c(1)=-4/dr2;d(1)=-p(1)+u(1)/dt;

a(nr)=0;b(nr)=1;d(nr)=0;u(nr)=0;


Oscila n tok v trubce6

Oscilační tok v trubce

NAP8

Výsledky jsou znázorněny v grafech

figure(1)

plot(t(1:nt),qfd(1:nt),'r',t(1:nt),qb(1:nt),'b',t(1:nt),qode(1:nt),'g')

figure(2)

hold off

for it=nt-100:50:nt

plot(r(1:nr),uode(it,1:nr),'g',r(1:nr),ufd(it,1:nr),'r',r(1:nr),ub(it,1:nr),'b')

hold on

end

q-diference

q-Bernoulli

q-aproximace


P m o hm ick oh ev

Přímý Ohmický ohřev

NAP8

Vnitřní (objemové) zdroje tepla vyvolané průchodem elektrického proudu nebo absorpcí mikrovlnného záření

Ardon


O hm ick oh ev masa

Napájení 50 Hz, konstantní napětí

Rovinné nerezové elektrody

plátek masa

Teplotní sondy

REFLEX R232

Nortech TS-RX4

MULTIPLEXER AGILENT 20CH 34901A

TECTRA - LMG 95

Měření výkonu

GPIB

Ohmický ohřev masa

NAP8

Na obrázku je schema aparatury pro experimentální přímý ohmický ohřev plátku masa (sekaná), který je sevřen mezi nerezovými elektrodami. Napájení 24V, frekvence 50Hz. Měří se elektrický výkon, teploty masa (optická vlákna) a teploty elektrod (termočlánky).

Cílem je stanovení elektrické vodivosti masa a její závislost na teplotě.

Kdyby byl ohřev naprosto rovnoměrný (teploty v každém místě stejné) stačil by jednoduchý integrální model a vyhodnocení elektrické vodivosti přímo z výkonu a teploty (viz dále). Nerovnoměrný ohřev je ale nutné modelovat řešením FK rovnice s proměnným objemovým zdrojem tepla.

Je to podobný případ jako sušení zrna (podobné rovnice a metody řešení). Při elektrickém ohřevu ale musíme kromě parabolické rovnice vedení tepla (FK rovnice s akumulačním a zdrojovým členem) současně řešit eliptickou rovnici rozložení elektrického potenciálu.


O hm ick oh ev masa integr ln model

E

Ohmický ohřev masa integrální model

NAP8

Pokud by byla teplota masa všude stejná a intenzita elektrického pole E konstantní, bylo by možné napsat jednoduchou integrální bilanci

Zanedbání tepelné kapacity elektrod

měrná elektrická vodivost rostoucí lineárně s teplotou

S-povrch, V-objem vzorku

jejímž řešením je exponenciální nárůst teploty

Když je R>0 bude rychlost ohřevu klesat, pro R<0 bude exponenciálně růst (akcelerace ohřevu je ovlivněna vysokou teplotní citlivostí elektrické vodivosti , tj. vysokou hodnotou koeficientu 1).

Následující numerické řešení je vlastně jen zpřesnění tohoto integrálního modelu a jeho analytického řešení.


O hm ick oh ev masa numerick model

Ohmický ohřevmasa numerický model

NAP8

Schiele


Metoda s t ohm ick oh ev

he

kontaktní vrstva W × W

T0

z

m-maso

W

hk

2H

Metoda sítí – ohmický ohřev

NAP8

Geometrie: ohřívaný vzorek masa má tloušťku 2H a příčné rozměry W x W

Experimenty ukazují, že přímý ohmický ohřev je ovlivněn i elektrickým a termickým odporem kontaktní vrstvy, a že elektrické vlastnosti masa i této kontaktní vrstvy závisí nejenom na teplotě, ale i na přítlačné síle elektrod.

Pro jednoduchost budeme vliv kontaktních odporů i efekt působení tlaku zanedbávat.

  • Předpoklady:

  • Kontaktní odpor je zanedbatelný, takže teplota masa je stejná jako teplota elektrody v místě kontaktu. Stejný je i elektrický potenciál (zanedbáme elektrickou dvouvrstvu, polarizační efekty, elektrolýzu).

  • Maso je homogenní, a termofyzikální vlastnosti jsou konstantní. Teplotně závislá je pouze měrná elektrická vodivost.

  • Povrch masa je chlazen přirozenou konvekcí vzduchem (konstantní teplota vzduchu Ta) a kontaktem s elektrodou (elektroda je rovněž chlazena přirozenou konvekci vzduchem). Elektrická vodivost elektrod (ocel) je tak velká, že potenciál U je v celé elektrodě konstantní a žádný elektrický výkon se v ní negeneruje.


Metoda s t fk rovnice

Metoda sítí – FK rovnice

NAP8

Teplotní pole v mase popisuje FK rovnice se zdrojem tepla, který je úměrný kvadrátu intenzity elektrického pole (což je gradient potenciálu U).

 je měrná elektrická vodivost masa, která závisí lineárně na teplotě

Tuto rovnici je třeba doplnit počátečními podmínkami (konstantní teplota T0 v čase 0) a okrajovými podmínkami (přestup tepla  na povrchu vzorku a teplota elektrod na ploškách z=0 a z=2H).

Rozložení teplot na povrchu elektrod je třeba stanovit současným řešením FK rovnice v elektrodách, ale tentokrát bez objemového zdroje tepla.


Metoda s t fk rovnice1

x

z

W

y

Metoda sítí – FK rovnice

NAP8

3D řešení s 3D sítí uzlových bodů je časově náročné (cílem je identifikace parametrů 0 a 1 a optimalizační algoritmus bude tudíž vyžadovat opakované řešení této diferenciální rovnice). Je tudíž žádoucí formulaci zjednodušit převedením na 1D problém. Zbavit se derivací x,y v příčném směru lze integrací FK rovnice přes celý průřez (teplota T(t,x,y,z) je po integraci nahrazena střední teplotou v místě z T(t,z)). Zanedbáme teplotní závislost tepelné vodivosti  a složky intenzity elektrického pole v příčném směru (na povrchu, který je elektricky izolovaný, i v ose musí být U/x i U/y nulové).

Tento člen vznikne integrací druhých derivací ve směrech x a y a použitím okrajových podmínek typu

Poznámka: Techniku redukce dimenze problému integrací 3D rovnic přes méně důležité rozměry jste použili v předmětu Tepelné procesy nebo Přenos hybnosti, tepla a hmoty např. při stanovení teplotního profilu podél chladicího žebra nebo při převodu Prandtlových rovnic mezní vrstvy na Karmánovy integrální rovnice.

W


Metoda s t fk rovnice2

Metoda sítí – FK rovnice

NAP8

Objemový zdroj tepla je vyjádřen intenzitou elektrického pole (dominantní složkou U/z), která se podél souřadnice z mění (mění-li se teplota a tedy i elektrická vodivost). Alternativně lze tento člen vyjádřit pomocí hustoty elektrického proudu I [A/m2], která se dle souřadnice z nemění, závisí jen a pouze na čase

V takto formulované FK rovnici se vyskytuje pouze teplota (a ne elektrický potenciál). Nelinearitu zdrojového členu lze odstranit linearizací (ověřte Taylorovým rozvojem vzhledem k teplotě T0):

Hustotu elektrického proudu I (i když je konstantní) je třeba stanovit v každém časovém kroku řešením Laplaceovy rovnice elektrického potenciálu


Metoda s t program

Metoda sítí – program

NAP8

Rámcové schéma numerického řešení (vývojový diagram)

Hlavní program

Podprogram MODEL

Počáteční podmínky (teploty a

potenciál ve všech uzlech)

Čtení geometrie a měřených dat

(časy, výkony, teploty)

Cyklus časových kroků

Použití algoritmu Nelder Mead pro

minimalizaci součtu čtverců

odchylek numerické predikce

(podprogram MODEL)

a měřených teplot a výkonů.

Optimalizované parametry jsou 0 1

Cyklus iterací

Sestavení a řešení soustavy rovnic

pro teplotní pole

Sestavení a řešení soustavy rovnic

pro elektrický potenciál

Výpočet hustoty proudu I

Konec iterací

Konec cyklu časových kroků


Metoda s t diskretizace fk

H

h

1 2 3 W P E N-1 N

Metoda sítí – diskretizaceFK

NAP8

Jednoduchá geometrie je dobrý předpoklad pro použití metody sítí nebo metody kontrolních objemů. Vzhledem k symetrii stačí řešit teplotní a potenciálové pole v intervalu z:<0,H>, který pokryjeme sítí N ekvidistantních uzlů. Každému uzlu přiřadíme teplotu TP a potenciál UP.

V metodě sítí požadujeme, aby v každém uzlu byla splněna příslušná diferenciální rovnice (počet podmínek je potom stejný jako počet neznámých). První a druhé derivace v diferenciálních rovnicích se nahradí konečnými diferencemi (numerickou aproximací derivací v bodě). Pro uzly P=2,3,…,N tudíž píšeme algebraické rovnice pro triádu sousedních teplot TW TP TE na nové hladině času


Metoda s t diskretizace fk1

Metoda sítí – diskretizaceFK

NAP8

Sdružením členů u neznámých TW TP TE vyjádříme koeficienty soustavy rovnic s tridiagonální maticí

bP

aP

cP

dP

Tuto rovnici píšeme i pro uzel na ose, tj. pro P=N, s tím, že TE=TW=TN-1, čímž aplikujeme i okrajovou podmínku symetrie v ose (T/z=0). Okrajovou podmínku na elektrodě (P=1) je ovšem nutné vyjádřit jinak, ne metodou sítí (konečných diferencí), ale kontrolních objemů.

Metodu konečných diferencí lze použít v takových uzlech, kde je víceméně hladký průběh řešení, a kde nejsou výrazné nespojitosti koeficientů řešené diferenciální rovnice. Pokud to tak není, je vhodnější použít metodu kontrolních objemů, která vyšetřovanou oblast rozdělí na menší podoblasti (kontrolní objemy) z nichž každý obklopuje „svůj“ uzlový bod. Místo požadavku na splnění diferenciální rovnice v izolovaném uzlu se požaduje splnění diferenciální rovnice v „průměru“, ve smyslu integrální bilance kontrolního objemu.


Metoda s t diskretizace fk2

1 2 3

he

Metoda sítí – diskretizace FK

NAP8

Kontrolní objem, který přiřadíme uzlu číslo 1 zahrnuje nalevo elektrodu, v pravé části maso. Z pravé části vtéká do kontrolní objemu entalpický tok T/z, z levé části odtéká tok povrchem elektrody. Jistá komplikace je v tom, že povrch elektrody Ae je větší než je průřez ohřívaného masa A=W2 a kontrolní objem má i svou tepelnou kapacitu, v níž převažuje tepelná kapacita elektrody. Velmi přibližnou integrální bilanci tohoto kontrolního objemu můžeme vyjádřit takto

tepelná kapacita elektrody vztažená k průřezu masa

tepelné ztráty povrchem elektrody

Nahrazením derivací diferencemi získáme analogon (integrální bilanci kontrolního objemu)

d1

c1

b1


Metoda s t diskretizace fk3

Metoda sítí – diskretizace FK

NAP8

Okrajová podmínka na elektrodě bude asi největším zdrojem nepřesností numerické predikce, protože je založena na předpokladech, které zcela určitě nejsou splněny: například teplota elektrody není stejná, protože součinitel teplotní vodivosti není nekonečně velký (u nerezové oceli je docela malý cca 15 W/m/K). Asi by tedy nebylo správné uvažovat celkovou tepelnou kapacitu elektrody, ale jen její část (redukcí povrchu elektrody Ae). Nejistoty se týkají i součinitele přenosu tepla z povrchu elektrody do okolního vzduchu – jde o přirozenou konvekci. Nejistoty související s odhadem  se ovšem týkají i chlazení povrchu masa,  na horizontálním a na vertikálním povrchu určitě nebude stejné. V následujícím programu bude použita korelace, kterou navrhl Churchill


Metoda s t diskretizace fk4

Metoda sítí – diskretizace FK

NAP8

Diskretizace rovnice pro elektrický potenciál i okrajové podmínky jsou mnohem jednodušší

WPE

kde

Na elektrodě předepíšeme přímo polovinu napájecího napětí (střední hodnotu u napětí střídavého) a v ose symetrie nulu


Program matlab

Program MATLAB

NAP8

Feininger


Ohmick oh ev matlab

Ohmický ohřev MATLAB

NAP8

Funkce fmeat, která popisuje numerickou predikci teplot i výkonů

function [vtime,vtemean,vtmmean,vpower] = fmeat ...

(h,w,u,he,ae,rkm0,rkm1,rlamm,n,tmax,ntime,ta,tinit)

%

% h - tloustka vzorku [m]

% w - prurez vzorku W x W [m]

% u - napeti [V] (rozdil potencialu na elektrodach)

% he - tloustka elektrody

% ae - pomer plochy elektrody a prurezu vzorku

% rkm0,rkm1 - koeficienty el.vodivosti masa kappa=rkm0+rkm1*T

% rlamm - tepelna vodivost masa

% N -celkovy pocet uzlu (uzel N na ose symetrie)

% tmax -max.cas

% ntime-pocet casovych kroku

% Ta - teplota vzduchu (okoli) ve stupnich Celsia

% Tinit(N) - pocatecni rozlozeni teplot

%

% Vystupy:

% vtime(ntime) - casy simulace

% vtemean(ntime),... - stredni teplot v elektrode, a v mase

% vpower(ntime),... - dissipovane vykony [W]

%

% numericka simulace

% prostorove kroky a fixni hodnoty

dz=double(h/(n-1));

dt=double(tmax/(ntime-1));

rhom=980.;

cpm=2700.;

rlame=15.0;

rhoe=7800.0;

cpe=460.0;

% pocatecni profil teplot a vodivosti

for i=1:n

tv0(i)=tinit(i);

end

% pocatecni profil rozlozeni napeti resenim tridiagonalni soustavy

av(1:n)=0;

bv(1:n)=1;

cv(1:n)=0;

dv(n)=0;

dv(1)=u/2;

for i=2:n-1

tm2=(tv0(i-1)+tv0(i))/2;

tp2=(tv0(i+1)+tv0(i))/2;

dv(i)=0;

rm2=(rkm0+rkm1*tm2);

rp2=(rkm0+rkm1*tp2);

av(i)=-rm2;

cv(i)=-rp2;

bv(i)=rm2+rp2;

end

uv0=tridag(av,bv,cv,dv,n);

% Prekopirovani pocatecniho rozlozeni teplot a napeti

uv=uv0;

tv=tv0;

% hustota elektrickeho proudu [A/m^2]

cur=(uv(3)-uv(2))/dz*(rkm0+rkm1*(tv(3)+tv(2))/2);

av-vektor pod diagonálou, bv-vektor diagonálních prvků a cv vektor prvků nad diagonálou matice soustavy

řešení soustavy n-rovnic s tridiagonální maticí soustavy


Ohmick oh ev matlab1

Ohmický ohřev MATLAB

NAP8

… pokračování funkce fmeat

% Nove rozlozeni potencialu

av(1:n)=0;

bv(1:n)=1;

cv(1:n)=0;

dv(1)=u/2;

dv(n)=0;

for i=2:n-1

tm2=(tv(i-1)+tv(i))/2;

tp2=(tv(i+1)+tv(i))/2;

dv(i)=0;

rm2=(rkm0+rkm1*tm2);

rp2=(rkm0+rkm1*tp2);

av(i)=-rm2;

cv(i)=-rp2;

bv(i)=rm2+rp2;

end

uv=tridag(av,bv,cv,dv,n);

cur=(uv(2)-uv(3))/dz*(rkm0+rkm1*(tv(3)+tv(2))/2);

% konec cyklu iterace na hladine IT

end

% vyhodnoceni strednich teplot v casovem kroku itime

tmm=0;

% vyhodnoceni strednich vykonu (kontakt, maso / Watty)

for i=1:n

tmm=tmm+tv(i);

end

vtemean(it)=tv(1);

vtmmean(it)=tmm/n;

vpower(it)=cur*u*w^2;

% konec cyklu casovych kroku - nova casova hladina

tv0=tv;

uv0=uv;

end

end

elegantnější by bylo omezit počet iterací testováním rozdílu vypočtených proudových hustot (while…). Mně se zdá, že 2 iterace úplně stačí…

% Cyklus casovych kroku

for it=1:ntime

vtime(it)=double(dt*(it-1));

% cyklus iteraci na jedne hladine casu

for iterat=1:2

% Teplotni profil

for i=2:n-1

% maso

alpha=(0.02+1.3*(w^3*(tv(i)-ta))^0.25)/w;

av(i)=-rlamm/dz^2;

cv(i)=-rlamm/dz^2;

bv(i)=rhom*cpm/dt+2*rlamm/dz^2+4/w*alpha+cur^2*rkm1/(rkm0+rkm1*tv(i))^2;

dv(i)=rhom*cpm*tv0(i)/dt+4/w*alpha*ta+cur^2*(rkm0+2*rkm1*tv(i))/(rkm0+rkm1*tv(i))^2;

end

% Rovnice v ose

alpha=(0.02+1.3*(w^3*(tv(n)-ta))^0.25)/w;

av(n)=-2*rlamm/dz^2;

cv(n)=0;

bv(n)=rhom*cpm/dt+2*rlamm/dz^2+4/w*alpha+cur^2*rkm1/(rkm0+rkm1*tv(n))^2;

dv(n)=rhom*cpm*tv0(n)/dt+4/w*alpha*ta+cur^2*(rkm0+2*rkm1*tv(n))/(rkm0+rkm1*tv(n))^2;

% Rovnice elektrody - korekce na odpor (zda se nevyznamna)

alphae=(0.02+1.3*(w^3*(tv(1)-ta))^0.25)/w;

alpha=1/(1/alphae+he/rlame);

cv(1)=-rlamm/dz;

bv(1)=rhoe*cpe*he*ae/dt+rlamm/dz+alpha*ae;

av(1)=0;

dv(1)=tv0(1)*rhoe*cpe*he*ae/dt+alpha*ae*ta;

% Nova iterace teplot na hladine it

tv=tridag(av,bv,cv,dv,n);

s tou rovnicí teplot elektrody je možné si pohrát a zkoušet různé varianty…


Ohmick oh ev matlab2

Ohmický ohřev MATLAB

NAP8

… pro řešení soustavy rovnic s tridiagonální maticí soustavy byla použita pomocná procedura tridag

function u = tridag(a,b,c,r,n)

bet=b(1);

u(1)=r(1)/bet;

for j=2:n

gam(j)=c(j-1)/bet;

bet=b(j)-a(j)*gam(j);

u(j)=(r(j)-a(j)*u(j-1))/bet;

end

for j=n-1:-1:1

u(j)=u(j)-gam(j+1)*u(j+1);

end


Ohmick oh ev matlab3

Ohmický ohřev MATLAB

NAP8

Funkce fdev počítající odchylku predikce modelu a experimentálních dat

function sumsq =fdev(p,h,w,u,he,ae,rlamm,n,ta,wel,wmeat,wpower,ndata,timdata,teldata,tmdata,powerdata)

% vypocet odchylek modelu - teplot i vykonu pro konkretni experiment

% sumsq - soucet ctvercu odchylek

% vstupy

% p(1)=rkm0, p(2)=rkm1

% wel,wmeat,wpower - vahove koeficienty pro soucet ctvercu

% ndata pocet casovych kroku v datech

% timdata(1:ndata) casy (predpoklad ekvidistantni vzorkovani)

% teldata(1:ndata) teploty elektrod

% tmdata(1:ndata) teploty masa (nebo sondy uprostred)

% powerdata(1:ndata) vykon ve watech celkovy

tmax=timdata(ndata)-timdata(1);

dtdata=tmax/(ndata-1);

nmult=10;

dt=dtdata/nmult;

ntime=tmax/dt+1;

% provedeni simulace pro pocatecni podminky teplot

tinit(1)=teldata(1);

for i=2:n

tinit(i)=tmdata(1);

end

[vtime,vtemean,vtmmean,vpower] = fmeat(h,w,u,he,ae,p(1),p(2),rlamm,n,tmax,ntime,ta,tinit);

sumsq=0;

iv=0;

for i=1:nmult:ntime

iv=iv+1;

sumsq=sumsq+wel*(vtemean(i)-teldata(iv))^2;

sumsq=sumsq+wmeat*(vtmmean(i)-tmdata(iv))^2;

sumsq=sumsq+wpower*(vpower(i)-powerdata(iv))^2;

end

end

numerické řešení bude prováděno s 10x menším časovým krokem než je vzorkování dat. To je otázka přesnosti-třeba vyzkoušet.


Ohmick oh ev matlab4

Ohmický ohřev MATLAB

NAP8

Skript čtení dat a volání optimalizační funkce fminsearch aplikované na fdev

% Fixni data

h=0.01;

w=0.02;

u=24;

he=0.002;

ae=2;

rlamm=0.71;

% parametry site, a casovy krok

n=101;

tmax=50;

ntime=101;

% teploty

ta=30;

t0=30;

% experimentalni data t, te,tm,power

rdata=load('fmeat.load');

% vahy

wel=.1;

wmeat=.1;

wpower=1;

[pnew,sums,exitflag]=fminsearch(@(x)fdev(x,h,w,u,he,ae,rlamm,n,ta,wel,wmeat,wpower,length(rdata),rdata(:,1),rdata(:,2),rdata(:,3),rdata(:,4)),[0.2,0.03])

Vzorek dat souboru fmeat.load (čas, teploty elektrod a masa, příkony…)

0.000e+000 3.029e+001 3.075e+001 1.762e+001

5.000e-001 3.073e+001 3.140e+001 1.785e+001

1.000e+000 3.058e+001 3.198e+001 1.805e+001

1.500e+000 3.078e+001 3.157e+001 1.805e+001

2.000e+000 3.045e+001 3.203e+001 1.825e+001

2.500e+000 3.046e+001 3.331e+001 1.820e+001

3.000e+000 3.118e+001 3.337e+001 1.876e+001

3.500e+000 3.119e+001 3.334e+001 1.871e+001

….

menší váha je přiřazena teplotám – věříme víc predikovaným výkonům

použití anonymní funkce, viz. druhá přednáška NAP2

počáteční odhady

pnew je vektor identifikovaných parametrů masa 0 1

indikátor úspěšnosti fminsearch (překročení max. počtu kroků…)


  • Login