Vytvoření skriptů pro webové rozhraní předmětu Analýza a simulace technologických procesů

Vytvoření skriptů pro webové rozhraní předmětu Analýza a simulace technologických procesů M-files for the Internet Interface Used in the Subject Analysis and Simulation of Technological Processes.

Petr Tomášek

Bakalářská práce 2007

*** nascannované zadání str. 1 ***

*** nascannované zadání str. 2 ***

UTB ve Zlíně, Fakulta aplikované informatiky, 2007

4

ABSTRAKT Tato bakalářská práce se zabývá analýzou a simulací následujících technologických procesů: průtočný výměník tepla s promícháváním, zásobníky na kapalinu zapojené za sebou, trychtýřový zásobník s nekonstantním průřezem a kulový zásobník s nekonstantním průřezem. U těchto modelů se simulují statické a dynamické charakteristiky pomocí M-file souborů v programu Matlab. Cílem je propojení těchto souborů přes webové rozhraní za účelem podpory výuky analýza a simulace technologických procesů.

Klíčová slova: dynamická charakteristika, statická charakteristika, bilance, metoda prosté iterace, metoda Runge-Kutta, diferenciální rovnice

ABSTRACT This thesis deals with the analysis and simulation of the following technological processes: flow heat exchanger with stirring, two water tanks system connected one after another, conical tank with a non-constant diameter, and the spherical tank with non-constant diameter. Those modes are simulated by static and dynamic characteristics using M-files in the Matlab application. The aim is to make those files compatible using the Internet interface to improve and support the methods when teaching the subjects Analysis and Simulation of Technological Processes.

Keywords: dynamic characteristic, static characteristic, value balance, the method of simple iteration, the method Runge-Kutta, differential equation


5

Poděkování: Touto cestou bych chtěl srdečně poděkovat svému vedoucímu bakalářské práce Ing. Jiřímu Vojtěškovi za jeho trpělivost a ochotu ke spolupráci, za jehož odborného vedení byla práce dovedena do konečné podoby. Také bych chtěl na tomto místě poděkovat Ing. Františku Gazdošovi, Ph.D., jehož semináře Analýza a simulace technologických procesů mi rovněž pomohly pochopit některé záležitosti této problematiky.

Prohlašuji, že jsem na bakalářské práci pracoval samostatně a použitou literaturu jsem citoval. V případě publikace výsledků, je-li to uvolněno na základě licenční smlouvy, budu uveden jako spoluautor.

Ve …………………….

Zlíně Podpis diplomanta


6

OBSAH ÚVOD .............................................................................................................................. 8 I

TEORETICKÁ ČÁST ........................................................................................... 9

1.1

PRŮTOČNÝ VÝMĚNÍK TEPLA S PROMÍCHÁVÁNÍM .................................... 10

1.1.1 ANALÝZA PROCESU .......................................................................................... 10 1.1.2 TEPELNÉ BILANCE ............................................................................................ 10 1.1.3 ODCHYLKOVÝ MODEL ...................................................................................... 12 1.2

ZÁSOBNÍKY NA KAPALINU ZAPOJENÉ ZA SEBOU ...................................... 14

1.2.1 ANALÝZA PROCESU ........................................................................................... 14 1.2.2 BILANCE SYSTÉMU ............................................................................................ 14 1.2.3 ODCHYLKOVÝ MODEL ....................................................................................... 15 1.3

TRYCHTÝŘOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM.......................................................................................................... 18

1.3.1 ANALÝZA PROCESU ........................................................................................... 18 1.3.2 BILANCE SYSTÉMU ............................................................................................ 18 1.4

KULOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM.......................................................................................................... 21

1.4.1 ANALÝZA PROCESU ........................................................................................... 21 1.4.2 BILANCE SYSTÉMU ............................................................................................ 21 1.5

POUŽITÉ METODY PRO VÝPOČET STATICKÝCH A DYNAMICKÝCH CHARAKTERISTIK V PROGRAMU MATLAB.................................................. 24

1.5.1 METODA PROSTÉ ITERACE .................................................................................. 24 1.5.2 METODA RUNGE-KUTTA .................................................................................... 25 II

PRAKTICKÁ ČÁST............................................................................................ 28

2.1

PRŮTOČNÝ VÝMĚNÍK TEPLA S PROMÍCHÁVÁNÍM .................................... 29

2.1.1 STATICKÁ CHARAKTERISTIKA ........................................................................... 29 2.1.2 DYNAMICKÁ CHARAKTERISTIKA....................................................................... 32 2.2

ZÁSOBNÍKY NA KAPALINU ZAPOJENÉ ZA SEBOU ...................................... 36

2.2.1 STATICKÁ CHARAKTERISTIKA............................................................................. 36 2.2.2 DYNAMICKÁ CHARAKTERISTIKA ........................................................................ 39 2.3

TRYCHTÝŘOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM.......................................................................................................... 43

2.3.1 STATICKÁ CHARAKTERISTIKA............................................................................. 43 2.3.2 DYNAMICKÁ CHARAKTERISTIKA ........................................................................ 45 2.4

KULOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM..... 47


7

2.4.1 STATICKÁ CHARAKTERISTIKA............................................................................. 47 2.4.2 DYNAMICKÁ CHARAKTERISTIKA ........................................................................ 49 ZÁVĚR........................................................................................................................... 51 SUMMARY.................................................................................................................... 53 SEZNAM POUŽITÉ LITERATURY.............................................................................. 55 SEZNAM POUŽITÝCH SYMBOLŮ A ZKRATEK....................................................... 56 SEZNAM OBRÁZKŮ .................................................................................................... 57 SEZNAM TABULEK..................................................................................................... 58 SEZNAM PŘÍLOH ......................................................................................................... 59


8

ÚVOD Tato bakalářská práce se zabývá analýzou a simulací technologických procesů probíhajících v průtočném výměníku tepla s promícháváním, zásobnících na kapalinu zapojených za sebou, trychtýřovém zásobníku s nekonstantním průřezem a kulovém zásobníku s nekonstantním průřezem. Tato práce je rozdělena na teoretickou a praktickou část. V teoretické části se provádí analýza procesů v jednotlivých modelech, při níž se specifikují děje, které v procesu probíhají a určí se veličiny popisující daný model. Výsledkem je pak teoretický model, na jehož základě postupně sestavujeme matematický model procesu, ve kterém se využijí matematické rovnice vyjadřující známé zákony a vztahy. Z tohoto matematického popisu dále můžeme sestavit řešení modelu pro simulační program, kterým se zabývá praktická část. V praktické části se provádí simulace statických a dynamických vlastností jednotlivých modelů na základě matematického popisu z teoretické části. Pro simulaci statické charakteristiky se využívá metoda prosté iterace (MPI) a pro simulaci dynamické charakteristiky zase metoda Runge-Kutta. Na základě těchto metod je sestaven výpočetní program v jazyku Matlab. Vytvořené programy jsou propojené přes webové rozhraní Matlab web serveru (MWS). Tyto programy simulují chování daných systémů a umožňují následné vyhodnocení výsledků na základě zadaných hodnot.


I. TEORETICKÁ ČÁST

.

9


10

1.1 PRŮTOČNÝ VÝMĚNÍK TEPLA S PROMÍCHÁVÁNÍM 1.1.1 Analýza procesu Jedná se o výměník tvořený nádobou s pláštěm za předpokladu dokonalého promíchávání. Do nádoby natéká výchozí chlazená kapalina o objemovém průtoku q a teplotě Tv a do pláště chladící kapalina o objemovém průtoku qc a teplotě Tcv.. Pro zjednodušení modelu zanedbáme tepelnou kapacitu stěny oddělující obě kapaliny. Ostatní veličiny, které proces popisují jako koeficient přechodu tepla a průtoky, objemy, hustoty a měrná tepla obou kapalin budou konstantní. Systém průtočného výměníku tepla s chlazením v plášti je zobrazen na obr. 1., kde q jsou označeny jako objemové průtoky kapaliny v [m3min-1], V objemy v [m3], F přestupná plocha v [m2], α koeficient přechodu tepla v [kJ m-2K-1min-1], T teploty v [K] a index (⋅)c označuje chladicí kapalinu a bez indexu označujeme chlazenou kapalinu.

Obr. 1. Průtočný výměník tepla s chlazením v plášti.

1.1.2 Tepelné bilance Základní bilanční rovnice má tvar: VSTUP + ZDROJ = VÝSTUP + AKUMULACE


11

Slovní vyjádření bilance chlazené kapaliny

Na základě bilancí můžeme charakterizovat stav tohoto procesu hodnotami veličin, které jej popisují. Matematické vztahy mezi těmito veličinami jsou v závislosti na čase. Z tepelné bilance chlazené kapaliny získáme následující lineární diferenciální rovnici q ρ c p Tv = q ρ c p T + Fα (T − Tc ) + Vρ c p

1. řádu

dT dt

(1.1)

Tuto diferenciální rovnici můžeme dále upravit na tvar

dT q q Fα = ⋅ Tv − ⋅ T − (T − Tc ) dt V V Vρ c p

(1.2)

Pro počáteční podmínky T (0) = T s , kde ρ je hustota chlazené kapaliny v [kg m-3] a cp měrné teplo chlazené kapaliny v [kJ kg-1K-1]. Slovní vyjádření bilance chladící kapaliny

Z tepelné bilance chladící kapaliny získáme druhou lineární diferenciální rovnici 1. řádu q c ρ c c pc Tcv + Fα (T − Tc ) = q c ρ c c pc Tc + Vc ρ c c pc

dTc dt

(1.3)

Opět upravený tvar diferenciální rovnice chladící kapaliny, který má tvar dTc q c q Fα = ⋅Tcv + ⋅ (T − Tc ) − c ⋅Tc d t Vc Vc ρ c c pc Vc

(1.4)


12

Pro počáteční podmínky Tc (0) = Tcs , kde ρc je hustota chladící kapaliny v [kg m-3] a cpc měrné teplo chladící kapaliny v [kJ kg-1K-1]. Můžeme konstatovat, že model popsaný rovnicemi (1.1), (1.3) je lineární, protože vztahy mezi vstupními a výstupními veličinami jsou lineární. Jedná se také o model systému se soustředěnými parametry, protože v procesu dochází k dokonalému promíchávání obou kapalin a tím pádem obě sledované stavové veličiny jsou nezávislé na souřadnicích.

1.1.3 Odchylkový model Odchylkový tvar získáme následujícím způsobem. Nejdříve model uvedeme do ustáleného stavu. Ustálený stav získáme anulováním časových derivací v rovnicích (1.2), (1.4) a dále nahradíme vstupní a výstupní veličiny jejich ustálenými hodnotami. Poté si vyjádříme dvě stavové veličiny T s a Tcs z rovnic ustáleného stavu. Následně změníme vstupní a výstupní hodnoty a jejich odchylky od ustáleného stavu původních veličin nám dají odchylkový model ve tvaru x1 (t ) = ∆T (t ) = T (t ) − T s , x 2 (t ) = ∆Tc (t ) = Tc (t ) − Tcs

(1.5)

u1 (t ) = ∆Tv (t ) = Tv (t ) − Tvs , u 2 (t ) = ∆Tcv (t ) = Tcv (t ) − Tcvs

(1.6)

Dále po dosazení odchylkových rovnic (1.5), (1.6) do rovnic (1.1), (1.3) dostaneme odchylkový tvar Vρ c p Vc ρ cc pc

dx1 + (q ρ c p + F α) x1 − F α x2 = q ρc p u1 dt

dx2 + ( qcρ c c pc + F α ) x2 − F α x1 = qcρ cc pc u 2 dt

(1.7) (1.8)

s nulovými počátečními podmínkami x1(0) = x2(0) = 0. Rovnice ve tvaru (1.7), (1.8) upravíme tak, že osamostatníme derivace výstupních veličin a zjednodušíme model zavedením konstant a11 , a12 , a 21 , a 22 , b11 , b22 . Model pak dostaneme ve tvaru dx1 = a11x1 + a12 x2 + b11 u1 dt

(1.9)

dx 2 = a 21 x1 + a 22 x 2 + b22 u 2 dt

(1.10)


kde

q Fα a11 = −  +  V Vρ c p 

b22 =

qc . Vc

 ,  

a12 =

Fα , V ρc p

a21 =

Fα , Vcρc c pc

13 q Fα  a22 = −  c + ,  Vc Vcρcc pc   

b11 =

q , V

Odchylkový tvar rovnic (1.9), (1.10) nám slouží k počítání dynamické charakteristiky modelu s nulovými počátečními podmínkami.


14

1.2 ZÁSOBNÍKY NA KAPALINU ZAPOJENÉ ZA SEBOU 1.2.1 Analýza procesu Matematický model válcových zásobníků na kapalinu je určený pro sledování výšky hladiny kapaliny v něm v závislosti na změnách přítoku. Do obou zásobníků přitékají kapaliny q1v a q2v a na odtoku z nádrže vytékají přes ventil. Tyto 2 propojené zásobníky jsou válcového tvaru a jsou zapojené za sebou. Oba zásobníky jsou ve stejné výšce, jsou otevřené a průřezy zásobníků jsou stejné. Systém zásobníků na kapalinu je zobrazen na obr. 2., kde q jsou objemové průtoky kapaliny v [m3s-1], h výšky hladin v [m] a F průřezy zásobníků v [m2] bez ohledu na indexy.

Obr. 2. Zásobníky na kapalinu zapojené za sebou

1.2.2 Bilance systému Základní bilanční rovnice má tvar: VSTUP + ZDROJ = VÝSTUP + AKUMULACE Slovní vyjádření bilance pro oba zásobníky


15

Při bilancování celého systému je bilancovanou veličinou množství kapaliny, které můžeme vyjádřit objemem V1 pro první zásobník a V2 pro druhý zásobník. Potom získáme bilanční rovnice ve tvaru dvou diferenciálních rovnic q1v = q1 +

dV1 dV , q 2v + q1 = q 2 + 2 dt dt

(1.11)

Změny objemů V1 a V2 můžeme také vyjádřit jako změnu výšek hladin v jednotlivých zásobnících v závislosti na průřezech zásobníků. V našem případě jsou průřezy zásobníků stejné. Tento vztah je vyjádřen tvarem dV = F dh a rovnici (1.11) tedy můžeme přepsat na tvar F1

dh dh1 + q1 = q1v , F2 2 + q 2 − q1 = q 2v dt dt

(1.12)

Pro počáteční podmínky h1 (0) = h1s , h2 (0) = h2s , což znamená že na počátku odpovídá výška hladiny ustálenému stavu. Jelikož průtoky kapalin vytékají ze zásobníků přes ventily, musíme brát v úvahu závislost hydrostatických tlaků působících na tyto ventily. Jedná se o druhou odmocninu z rozdílu tlaků kapaliny před a za ventilem. Pro náš případ jsou hydrostatické tlaky úměrné výškám hladin v zásobnících, takže průtoky můžeme vyjádřit dvěma nelineárními rovnicemi q1 = k1 h1 − h2 , q 2 = k 2 h2

(1.13)

kde k1, k2 jsou konstanty ventilu. Tento matematický model můžeme popsat dvěmi diferenciálními rovnicemi (1.12) a dvěmi nelineárními rovnicemi (1.13). Na základě těchto dvou nelineárních rovnic můžeme tvrdit, že se jedná o nelineární model.

1.2.3 Odchylkový model Pro získání odchylkového tvaru nejdříve převedeme model do ustáleného stavu. Z ustáleného stavu potom můžeme určit počáteční podmínky, které nám následně pomohou k sestavení linearizovaného modelu a který bude současně modelem odchylkovým. Pro ustálený stav platí, že průtok q1s = q1sv a

q 2s = q1sv + q 2s v

dosazení těchto vztahů do rovnic (1.13) obdržíme dvě rovnice

a výška hladiny h1s > h2s . Po


q1sv = k1 h1s − h2s , q1sv + q 2s v = k 2 h2s

16

(1.14)

Nyní si z rovnic (1.14) vyjádříme ustálené hodnoty výšek hladin h1s a h2s . Potom z výšek hladin a vstupních průtoků provedeme odchylky od počátečního ustáleného stavu x j (t ) = ∆ h j (t ) = h j (t ) − h sj , u j (t ) = ∆ q jv (t ) = q jv (t ) − q sjv , pro j = 1, 2

(1.15)

Dále pak převedeme rovnice (1.12) do odchylkového tvaru F1

d ∆h1 d ∆h2 + ∆ q1 = ∆ q1v , F2 + ∆ q 2 − ∆ q1 = ∆ q 2 v dt dt

(1.16)

Následně dochází k linearizaci modelu, když pomocí Taylorova rozvoje nahradíme odchylky průtoků lineárními členy v okolí pracovního bodu ( h1s , h2s ). Pro tyto odchylky postupně dostaneme s

s

k1 h1s − h2s  ∂q   ∂q  k1 ∆ q1 ≈  1  ∆ h1 +  1  ∆h2 = ( ∆h1 − ∆h2 ) = ( ∆h1 − ∆h2 ) = 2 ( h1s − h2s )  ∂ h1   ∂ h2  2 h1s − h2s

 ∂q ∆ q 2 ≈  2  ∂ h2

q1s

( ∆h1 − ∆h2 ) = K 1 ( ∆h1 − ∆h2 )

(1.17)

k 2 h2s  q 2s k2  ∆ h2 = h ∆h2 = ∆ = ∆h2 = K 2 ∆h2 2 2 h2s 2 h2s  2 h2s

(1.18)

=

2 (h1s − h2s ) s

Vyjádříme-li koeficienty K1 a K2, které jsou závislé na poloze pracovního bodu a rovněž na počátečním ustáleném stavu dostaneme rovnice K1 =

q1s q 2s , K = 2 2 (h1s − h2s ) 2 h2s

(1.19)

Při další úpravě dosadíme rovnice (1.17) a (1.18) spolu s rovnicemi (1.15) do rovnic (1.16) a tím už dostaneme linearizovaný odchylkový model systému ve tvaru dx1 + K 1 ( x1 − x 2 ) = u1 dt

(1.20)

dx 2 + K 2 x 2 − K 1 ( x1 − x 2 ) = u 2 dt

(1.21)

F1 F2


17

Počáteční podmínky budou nulové x1(0) = x2(0) = 0. Rovnice ve tvaru (1.20), (1.21) upravíme a zjednodušíme zavedením konstant

a11 , a12 , a 21 , a 22 , b11 , b22 . Model pak

dostaneme ve tvaru

kde a11 = −

dx1 = a11 x1 + a12 x2 + b11 u1 dt

(1.22)

dx2 = a21 x1 + a22 x2 + b22 u2 dt

(1.23)

K1 K K K + K2 1 1 , a12 = 1 , a21 = 1 , a22 = − 1 , b11 = , b22 = . F1 F1 F2 F2 F1 F2

Tento linearizováný odchylkový model obsahuje v konstantách aij již zmíněné koeficienty K 1 a K2 a tím pádem je závislý na počátečním ustáleném stavu systému a na poloze pracovního bodu. S

nulovými počátečními podmínkami použijeme tento odchylkový tvar rovnic (1.22), (1.23) k počítání dynamické charakteristiky modelu.


18

1.3 TRYCHTÝŘOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM 1.3.1 Analýza procesu V modelu trychtýřového zásobníku se opět budeme zabývat změnou výšky hladiny v závislosti na změně přítoku. Do zásobníku nám přitéká kapalina qv, která vytéká přes ventil. Tento zásobník je nekonstantního průřezu. Systém trychtýřového zásobníku je zobrazen na obr. 3., kde q jsou objemové průtoky kapaliny v [m3min-1] bez ohledu na index, h je výška hladiny v [m], H je výška zásobníku v [m], d je průměr zásobníku ve výšce hladiny v [m] a D je průměr zásobníku v [m].

Obr. 3. Trychtýřový zásobník s nekonstantním průřezem

1.3.2 Bilance systému Základní bilanční rovnice má tvar: VSTUP + ZDROJ = VÝSTUP + AKUMULACE


19

Slovní vyjádření bilance pro trychtýřový zásobník

Bilancovanou veličinou je množství kapaliny, které vyjádříme objemem V a systém popíšeme diferenciální rovnicí ve tvaru qv = q +

dV dt

(1.24)

Tuto diferenciální rovnici lze přepsat na tvar, ve kterém změnu objemu vyjádříme pomocí vztahu dV = F dh a pak dostáváme rovnici qv = q + F ⋅

dh dt

(1.25)

Odtok ze zásobníku vytéká přes ventil. Na ventil působí hydrostatické tlaky (viz. zásobníky na kapalinu), které jsou úměrné výšce hladiny v zásobníku. Proto pro tento případ platí nelineární vztah q = k h , který po dosazení do rovnice (1.25) nám umožní tuto rovnici přepsat na tvar qv = k h + F ⋅

dh dt

(1.26)

Jelikož je průřez F kruhového tvaru, můžeme ho také vyjádřit jako F=

π ⋅ d2 4

(1.27)

a Pro pravoúhelný trojúhelník obecně platí vztah tgα = , který můžeme využít v našem b modelu. Pravoúhlý trojúhelník lze vidět na obr. 3., kde je pro nás strana trojúhelníku a=d/2 a b=h. Tento vztah potom upravíme pro náš model a dostaneme ho ve tvaru tgα =

d/2 ⇒ d / 2 = tgα ⋅ h h

(1.28)


20

a zároveň platí tgα =

D/ 2 D = H 2H

(1.29)

Potom po dosazení do rovnice (1.28) dostáváme d 2 D2 2 d D = tgα ⋅ h = ⋅h ⇒ = ⋅h 2 2H 4 4H 2

(1.30)

Dále pro úpravu plochy průřezu dosadíme rovnici (1.30) v umocněném tvaru do rovnice (1.27) a získáme následující vztah F=

π ⋅ d 2 π ⋅ D2 2 = ⋅h 4 4H 2

(1.31)

Nyní můžeme přepsat bilanční rovnici (1.26) na tvar qv = k h + F ⋅

dh π ⋅ D 2 2 dh =k h+ ⋅h ⋅ dt dt 4H 2

(1.32)

Následně osamostatníme na levou stranu derivaci výšky hladiny h podle času t z rovnice (1.32) a dostaneme

(

)

qv − k h ⋅ 4 H 2 dh q v − k h = = dt π ⋅ D 2 2 π ⋅ D 2 ⋅ h2 ⋅ h 4H 2

(1.33)

Tato výsledná rovnice se použije pro počítání dynamické charakteristiky, kde počáteční podmínkou bude h(0) = h s .


21

1.4 KULOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM 1.4.1 Analýza procesu Model kulového zásobníku je podobný trychtýřovému zásobníku, liší se pouze tvarem a postupem při odvozování bilanční rovnice. Do kulového zásobníku přitéká kapalina qv, u které změníme přítok a budeme zjišťovat změnu výšky hladiny v zásobníku. Přítok kapaliny qv vytéká přes ventil. Zásobník má nekonstantní průřez. Systém kulového zásobníku je zobrazen na obr. 4., kde q jsou objemové průtoky kapaliny v [m3min-1] bez ohledu na index, h je výška hladiny v [m] a D je průměr zásobníku v [m].

Obr. 4. Kulový zásobník s nekonstantním průřezem

1.4.2 Bilance systému Základní bilanční rovnice má tvar: VSTUP + ZDROJ = VÝSTUP + AKUMULACE


22

Slovní vyjádření bilance pro kulový zásobník

Diferenciální rovnice popisující daný děj systému je stejná jako u předešlého modelu. Tvar rovnice je qv = q +

dV dt

(1.34)

Po stejných úpravách jako u trychtýřového zásobníku můžeme rovnici (1.34) přepsat na tvar qv = k h + F ⋅

dh dt

(1.35)

Opět využijeme pro průřez obsah kruhu a vyjádříme ho ve tvaru F=

π ⋅ d2 4

(1.36)

Dále využijeme Pythagorovy věty pro pravoúhlý trojúhelníka a podle obr. 4. platí 2

 D  D  x +h−  =    2 2

2

2

(1.37)

Po úpravě této rovnice dostáváme x 2 = D⋅ h − h 2 = h ⋅ (D − h) A jestliže x 2 =

(1.38)

d2 , pak můžeme dosadit rovnici (1.38) do (1.36) a získáme 4 F =π ⋅ (D − h) ⋅ h

(1.39)

Bilanční rovnice (1.35) nám potom dá tvar qv = k h + F ⋅

dh dh = k h + π ⋅ (D − h ) ⋅ h ⋅ dt dt

(1.40)


23

Po osamostatnění derivace na levou stranu, dostaneme výslednou diferenciální rovnici pro počítání dynamiky s počáteční podmínkou h(0) = h s . Rovnice má tvar qv − k h dh = dt π ⋅ (D − h ) ⋅ h

(1.41)

Při tvorbě analýzy procesů a odvozování matematických modelu v teoretické části parafrázuji použitou literaturu [4].


24

1.5 POUŽITÉ METODY PRO VÝPOČET STATICKÝCH A DYNAMICKÝCH CHARAKTERISTIK V PROGRAMU MATLAB Pro výpočet statických charakteristik se nejčastěji využívá metody prosté iterace, co se dynamických metod týče, pak metody Runge-Kutta. Tyto metody jsou pro modely popisované v této práci plně dostačující. Obě jsou zahrnuty v matematických funkcích programu Matlab, jehož integrované prostředí je z tohoto důvodu využito pro řešení a simulaci modelů.

1.5.1 Metoda prosté iterace Pro postupy využívané při simulaci ustáleného stavu systému je tedy využito metody prosté iterace (MPI). Tato metoda se používá pro řešení soustav nelineárních rovnic. Pro systémy lineárních rovnic se prakticky nepoužívá. MPI tedy můžeme použít pro numerické řešení rovnic, které mohou být ve tvaru f ( x) = 0

(1.42)

Iterační úlohy jsou založeny na řešení rovnic (1.42) převedených do ekvivalentního tvaru x = g( x)

(1.43)

Tato rovnice slouží na nalezení pevných bodu funkce g. Pro ekvivalentnost úloh platí, jestliže ξ je pevný bod funkce g, pak ξ je kořen funkce f a naopak. Dále platí věta, jestliže g ∈C [a,b], g: [a,b]→ [a,b] pak funkce g má v intervalu [a,b] pevný bod. A jestliže g navíc splňuje Lipshitzovu podmínku s konstantou q ve vztahu g( x) − g( y) ≤ q x − y , kde 0 ≤ q < 1

(1.44)

pro ∀x, y ∈[a, b] pak g má v intervalu jediný pevný bod. Za předpokladu, že g ∈C [a,b], g: [a,b]→ [a,b] se pak zvolí počáteční aproximace x 0 ∈[a, b] a poté se generuje posloupnost {x k }k =0 následujícím způsobem ∞

x k +1 = g(x k ) , k=0,1,2…

(1.45)


25

Pak funkci g nazýváme iterační funkcí a uvedenou metodu iterační metodou nebo metodou prosté iterace. Uvedená metoda patří mezi jednokrokové metody. Obecně pak pro

j-

krokové metody platí vztah

(

)

x k +1 = g x k , x k −1 ,...,x k − j +1 , k=0,1,2…

(1.46)

Za již zmíněných předpokladů, kde g ∈C [a,b], g: [a,b]→ [a,b] pro počáteční konfiguraci x 0 ∈[a, b] je posloupnost {x k }k =0 , x k = g(x k −1 ) konvergentní a platí, že ∞

lim x k = ξ k →∞

(1.47)

kde ξ je pevný bod funkce g. Nechť pro funkci g jsou splněny předchozí předpoklady a

{ }

pro posloupnost x k

∞ k =0

, x 0 ∈[a, b] , x k = g(x k −1 ) dále platí x k −ξ ≤

qk 0 1 x − x , ∀k ≥ 1 1− q

(1.48)

Při popisu metody prosté iterace vycházím ze zdroje [1] uvedeného v seznamu použité literatury.

1.5.2 Metoda Runge-Kutta Metoda Runge-Kutta patří mezi krokové metody a bere do úvahy i členy vyšších řádů. Potřebné derivace funkce f(t,x) počítá složitější diferenční metodou pomocí dalších pomocných bodů mezi sousedními uzly v síti. Obecně lze metody Runge-Kutta zapsat následovně i −1   ki = f t + αi h, xn + h∑βij k j  j =1  

(1.49)

p

xn+1 = xn + h∑wi ki

(1.50)

i =1

Koeficienty u těchto metod jsou vypočteny tak, aby metoda řádu p odpovídala Taylorovu polynomu funkce x(t) stejného řádu.


26

Metod Runge-Kutta je více. Klasická metoda 4. řádu je součástí funkce ode45 a je vyjádřena z rovnic (1.49), (1.50) a vypadá následovně k1 = f (t n , xn )

(1.51)

h h  k2 = f  t n + , xn + k1  2 2 

(1.52)

h h  k3 = f t n + , xn + k2  k4 = f (t n + h, xn + k3 h) 2 2 

(1.53)

k4 = f (t n + h, xn + k3 h)

(1.54)

xn+1 = xn +

h (k1 + 2k2 + 2k3 + k4 ) 6

(1.55)

Pomocné hodnoty ki představují derivace stavu systému ve speciálních bodech na začátku, na konci a uprostřed intervalu t n , t n+1 . V případě fázového prostoru vyšší dimenze m (m>1) je stav systému popsán vektorem x a dynamický systém vektorovou funkcí f(t,x). Algoritmus Runge-Kutta je pak třeba psát ve vektorové formě, takže pomocné hodnoty ki budou m-složkové vektory. Pak soustava rovnic vyjádřená ze vztahu (1.51) vypadá následovně k11 = f1 (t n , xn1, , xn2 ,...,xnm )

(1.56)

k12 = f 2 (t n , xn1, , xn2 ,...,xnm )

(1.57)

… k1m = f m (t n , xn1, , xn2 ,...,xnm )

(1.58)

Metodu Runge-Kutta lze použít pro řešeni obyčejných diferenciálních rovnic (ODR). Díky tomu je tato metoda součástí matematické funkce Matlabu ode45, která je použita při řešení dynamických charakteristik. Funkce ode45 využívá metody Runge-Kutta pro numerické řešení ODR, které obsahují derivace jen podle jedné proměnné. Ode45 lze použít pro řešení ODR jak prvního tak i vyššího řádu. Syntaxe funkce ode45 má následující tvar: [T, Y] = ode45(‘název funkce’, časový interval, počáteční podmínky)

UTB ve Zlíně, Fakulta aplikované informatiky, 2007 kde název funkce

27

je název M-file souboru popisujícího soustavu diferenciálních rovnic, v němž je obsažena funkce se syntaxí dy = název funkce(t, y), ve které je dy sloupcový vektor, představující hodnoty derivací pro vektor hodnot y v čase t

časový interval

je vektor o dvou prvcích – počáteční a koncový čas řešení

počáteční podmínky je vektor počátečních podmínek y0, pro který platí y(t0) = y0 T

je sloupcový vektor obsahující časové okamžiky řešení

Y

je matice vlastního řešení

Při definici funkce Runge-Kutta jsem čerpal ze zdroje [2] a při popisu funkce ode45 ze zdroje [3], které jsou uvedené v seznamu použité literatury.


II. PRAKTICKÁ ČÁST

28


29

2.1 PRŮTOČNÝ VÝMĚNÍK TEPLA S PROMÍCHÁVÁNÍM 2.1.1 Statická charakteristika Statická charakteristika popisuje model v ustáleném stavu, který je charakterizován rovnicemi q ρ c p Tvs = q ρ c p T s + Fα (T s − Tcs )

(2.1)

q c ρ c c pc T cvs + Fα (T s − T cs ) = q c ρ c c pc T cs

(2.2)

Jde o rovnice původního modelu u nichž byla anulována derivace podle času. Po vyjádření stavových veličin T s , Tcs jsou tyto rovnice ve tvaru Ts =

Tcs =

q ρ c p Tvs + FαTcs

(2.3)

q ρ c p + Fα q c ρ c c pc Tcvs + FαT s

(2.4)

q c ρ c c pc + Fα

Ve statické charakteristice se počítá závislost těchto ustálených teplot T s a Tcs na změnách intervalu průtoku q a qc. Při změně intervalu průtoku q se průtok qc nemění a počítá se průběh teplot T s a Tcs pro ustálený stav v daném intervalu průtoku q. Pro změnu intervalu průtoku qc to platí obdobně. V následujících tabulkách jsou uvedeny použité hodnoty pro počítání statické charakteristiky. Tab. 1. Tabulka zadaných proměnných pro statickou charakteristiku Zadané proměnné Začátek intervalu pro různý průtok q

Hodnoty 0 [m3 min-1]

Krok intervalu pro různý průtok q

0,0005 [m3 min-1]

Konec intervalu pro různý průtok q

1 [m3 min-1]

Začátek intervalu pro různý průtok qc

0 [m3 min-1]

Krok intervalu pro různý průtok qc

0,0005 [m3 min-1]

Konec intervalu pro různý průtok qc

1 [m3 min-1]


30

Tab. 2. Tabulka zadaných veličin Konstanty

Hodnoty

Objemový průtok chlazené kapaliny q

0,08 [m3 min-1]

Objemový průtok chladící kapaliny qc

0,03 [m3 min-1]

Objem chlazené kapaliny V

1,2 [m3]

Objem chladící kapaliny Vc

0,64 [m3]

Hustota chlazené kapaliny ρ

985 [kg m-3]

Hustota chladící kapaliny ρc

998 [kg m-3]

Měrné teplo chlazené kapaliny cp

4,05 [kJ kg-1 K-1]

Měrné teplo chladící kapaliny cpc

4,18 [kJ kg-1 K-1]

Přestupná plocha F Koeficient přechodu tepla α

5,5 [m2] 43,5 [kJ m-2 min-1K]

Teplota chlazené kapaliny na vstupu Tvs

323 [K]

Teplota chladící kapaliny na vstupu Tcvs

293 [K]

Vyhodnocením jsou 4 grafy simulující statické chování průtočného výměníku tepla na základě zadaných vstupních hodnot. Grafy jsou zobrazeny na obr. 5.


31

Obr. 5. Průběhy ustálených teplot T s a Tcs v závislosti na změnách průtoku q a qc

Na obr. 5. je vidět průběh teplot chlazené a chladící kapaliny v závislosti na změně průtoku. Reakce na změnu průtoku q u teploty chlazené kapaliny T s se projevila zvýšením a následným ustálením teploty. Teplota chladící kapaliny Tcs reagovala také zvýšením a následným ustálením své teploty. Po odezvě změny průtoku q je ustálená teplota chlazené kapaliny vyšší než ustálená teplota chladící kapaliny. Odezva na změnu průtoku qc se projevila u teploty chlazené i chladící kapaliny snížením a následným ustálením teplot. Při této změně průtoku je ustálená teplota chlazené kapaliny vyšší než ustálená teplota chladící kapaliny. Z těchto teplot v ustáleném stavu se pak mohou určit polohy pracovních bodů. Program statické charakteristiky je v příloze P I pod názvem Průtočný výměník – statická charakteristika. Program je propojený přes Matlab web server (MWS), kde se přes webové rozhraní zadávají vstupní hodnoty a naopak vypočtené hodnoty se pošlou opět přes MWS na internetové stránky. Na obr. 6. je zobrazen vzhled tabulky pro zadávání vstupních hodnot z webové stránky.


32

Obr. 6. Tabulka pro zadávání vstupních hodnot s tabulkou konstant

2.1.2 Dynamická charakteristika Dynamická charakteristika se řeší z odchylkového modelu jako funkce času při nulových počátečních podmínkách. Rovnice mají tvar dx1 = a11x1 + a12 x2 + b11 u1 dt

(2.5)

dx 2 = a 21 x1 + a 22 x 2 + b22 u 2 dt

(2.6)

V dynamice sledujeme jak se změní průběh teploty chlazené kapaliny T a chladící kapaliny Tc při skokové změně průtoku q a qc za daný časový průběh simulace. V následující tabulce jsou uvedeny použité hodnoty pro počítání dynamické charakteristiky.


33

Tab. 3. Tabulka zadaných proměnných pro dynamickou charakteristiku Zadané proměnné

Hodnoty

Doba simulace Tsim

150 [s]

Skoková změna průtoku chlazené kapaliny q

2 [m3 min-1]

Skoková změna průtoku chladící kapaliny qc

3 [m3 min-1]

Ostatní zadané veličiny jsou stejné jako v tabulce na obr. 2. Výstupem jsou 2 grafy dynamického chování průtočného výměníku, které jsou zobrazeny na obr. 7. a obr. 8. 2.5

2

1.5 ] K[ T 1

0.5

0

0

50

100

150

t[min]

Obr. 7. Průběh teploty chlazené kapaliny v závislosti na čase


34

2.5

2

1.5 ] K[ c T

1

0.5

0

0

50

100

150

t[min]

Obr. 8. Průběh teploty chladící kapaliny v závislosti na čase

Na obr. 7. můžeme pozorovat časový průběh teploty chlazené kapaliny při skokové změně průtoku q. Graf na obr. 8. ukazuje průběh teploty chladící kapaliny při skokové změně průtoku qc. U obou grafů se dynamické vlastnosti projevily přechodovou charakteristikou 1. řádu, která udává závislost skokové změny na čase. Z přechodové charakteristiky se pak může určit časová konstanta a zesílení soustavy. Po odezvě skokové změny a vyhodnocení přechodové charakteristiky lze vidět, že teplota chladící kapaliny Tc je větší než teplota chlazené kapaliny T. Program dynamické charakteristiky je uveden v příloze P II pod názvem Průtočný výměník – dynamická charakteristika a tvoří ho 3 M-file soubory (vymenik.m, vymenik_hlavni.m, vymenik_konstanty.m), které spolu navzájem komunikují. Dynamika se v programu počítá pomocí funkce ode45 pro řešení diferenciálních rovnic. Vzhled webové stránky programu propojeného přes MWS je zobrazen na obr. 9.



35


36

2.2 ZÁSOBNÍKY NA KAPALINU ZAPOJENÉ ZA SEBOU 2.2.1 Statická charakteristika Ve statické charakteristice je uveden model v počátečním ustálením stavu, pro který platí vztahy q1s = q1sv a q 2s = q1sv + q 2s v . Za těchto podmínek je ustálený stav ve tvaru q1sv = k1 h1s − h2s , q1sv + q 2s v = k 2 h2s

(2.7)

kde k1, k2 jsou konstanty ventilu. Po vyjádření ustálených výšek hladin z rovnic (2.7) dostaneme dvě rovnice, které nám umožní určení počátečních podmínek a také možnost vypočítat polohu pracovního bodu ( h1s , h2s ). Rovnice mají tvar 2

 qs + qs   qs  h1s =  1v  + h2 s , h2 s =  2 v 1v   k1   k2 

2

(2.8)

s počátečními podmínkami h1 (0) = h1s , h2 (0) = h2s , což znamená že na počátku odpovídá výška hladiny ustálenému stavu. Ve statické charakteristice počítáme změnu výšky hladin v zásobnících. Tyto změny budeme sledovat v závislosti na změnách přítoku q1v a q 2 v . To znamená, že při změně přítoku q1v bude přítok q 2 v konstantní a spočítá se ustálená výška hladiny h1s v prvním zásobníku a h2s v druhém zásobníku. Totéž se provede při změně přítoku q 2 v . V následujících tabulkách jsou uvedeny použité hodnoty pro počítání statické charakteristiky.

Tab.4. Tabulka zadaných proměnných pro statickou charakteristiku


Zadané proměnné Začátek intervalu pro různý průtok q1v

37

Hodnoty 0 [m3 s-1]

Krok intervalu pro různý průtok q1v

0,1 [m3 s-1]

Konec intervalu pro různý průtok q1v

1 [m3 s-1]

Začátek intervalu pro různý průtok q2v

0 [m3 s-1]

Krok intervalu pro různý průtok q2v

0,1 [m3 s-1]

Konec intervalu pro různý průtok q2v

1 [m3 s-1]


Hodnoty

Objemový průtok kapaliny v prvním zásobníku q1sv

0,5 [m3 s-1]

Objemový průtok kapaliny v druhém zásobníku q 2sv

0,3 [m3 s-1]

Výška hladiny v prvním zásobníku h1s

1,5 [m]

Výška hladiny v druhém zásobníku h2s

1 [m]

Průřez prvního zásobníku F1

2 [m2]

Průřez druhého zásobníku F2

2 [m2]

Vyhodnocením jsou 4 grafy, které charakterizují statické chování zásobníku na kapalinu za daných podmínek. Grafy jsou zobrazeny na obr. 10.


38

Obr. 10. Průběhy výšek hladin h1s a h2s v závislosti na změnách průtoku q1v a q2v

Na obr. 10. můžeme pozorovat reakci na změnu jednotlivých průtoku v podobě průběhu výšek hladin v jednotlivých zásobnících. Odezva na změnu průtoku q1v v prvním zásobníku se projevila zvýšením hladiny v obou zásobících, přičemž platí h1s > h2s . Při změně průtoku q2v se výšky hladin opět zvýší a h1s bude zase větší než h2s . Po ustálení výšek hladin se mohou určit polohy pracovních bodu h1s , h2s v okolí kterých se provádí linearizace. Program statické charakteristiky je v příloze P III pod názvem Zásobníky na kapalinu – statická charakteristika. Program je propojený přes MWS na stejném principu jako u předešlého modelu. Na obr. 11. je zobrazen vzhled tabulky pro zadávání vstupních hodnot z webové stránky.


39


2.2.2 Dynamická charakteristika Dynamickou charakteristiku budeme řešit podobně jako u modelu průtočného výměníku tepla. Akorát u tohoto modelu se nejdříve musí nelineární model linearizovat.. Proto se dynamika bude počítat z linearizovaného odchylkového modelu jako funkce času. Počáteční podmínky budou nulové. Rovnice jsou ve tvaru dx1 = a11 x1 + a12 x2 + b11 u1 dt

(2.9)

dx2 = a21 x1 + a22 x2 + b22 u2 dt

(2.10)

V dynamice pozorujeme změnu výšek hladin v jednotlivých zásobnících při skokových změnách přítoků do jednotlivých nádrží zásobníků v závislosti na časovém průběhu simulace. V následující tabulce jsou uvedeny použité hodnoty pro počítání dynamické charakteristiky.


40

Tab. 6. Tabulka zadaných proměnných pro dynamickou charakteristiku Zadané proměnné

Hodnoty

Doba simulace Tsim

200 [s]

Skoková změna průtoku q1v

0,09 [m3 s-1]

Skoková změna průtoku q2v

0,05 [m3 s-1]

Ostatní zadané veličiny jsou stejné jako v tabulce na obr. 5. Výstupem jsou 2 grafy, které charakterizují dynamické chování zásobníku na kapalinu za daných podmínek. Grafy jsou zobrazeny na obr. 12. a obr. 13. 12

10

8

] m [ 1 h

6

4

2

0

0

50

100 t[s]

Obr. 12. Průběh výšky hladiny v prvním zásobníku v závislosti na čase

150


41

10 9 8 7 6 ] m [ 2 h

5 4 3 2 1 0

0

50

100

150

t[s]

Obr. 13. Průběh výšky hladiny v druhém zásobníku v závislosti na čase

Na obr. 12. můžeme pozorovat časový průběh výšky hladiny h1 při skokové změně průtoku q1v. Na obr. 13. lze vidět časový průběh výšky hladiny h2 při skokové změně průtoku q2v. Po odezvě skokových změn je výška hladiny v prvním zásobníku vyšší než výška hladiny v druhém zásobníku. Průběh změn výšek hladin se projevil přechodovými charakteristikami 1. řádu. Z těchto přechodových charakteristik dále můžeme určit časové konstanty a zesílení soustavy pro jednotlivé průběhy výšek hladin. Program dynamické charakteristiky je uveden v příloze P IV pod názvem Zásobníky na kapalinu – dynamická charakteristika a tvoří ho 3 M-file soubory (zasobnik.m, zasobnik_hlavni.m, zasobnik_konstanty.m). Vzhled webové stránky programu propojeného přes MWS je zobrazen na obr. 14.



42


43

2.3 TRYCHTÝŘOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM 2.3.1 Statická charakteristika Statickou charakteristiku vyjadřuje rovnice počátečního ustáleného stavu ve tvaru q vs = k h s

(2.11)

Z této rovnice můžeme vypočítat ustálenou výšku hladiny  qs h =  v  k s

  

2

(2.12)

Ve statické charakteristice budeme řešit ustálenou výšku hladiny v zásobníku v závislosti na změně přítoku qv. V následujících tabulkách jsou uvedeny použité hodnoty pro počítání statické charakteristiky. Tab. 7. Tabulka zadaných proměnných pro statickou charakteristiku Zadané proměnné

Hodnoty

Začátek intervalu pro různý průtok qv

0 [m3 min-1]

Konec intervalu pro různý průtok qv

1 [m3 min-1]

Konstanta ventilu

2


Hodnoty

Průměr zásobníku D

2 [m]

Výška zásobníku H

3 [m]

Výška hladiny h s

2 [m]

Objemový průtok kapaliny q vs

0,6 [m3 min-1]

Vyhodnocením je zobrazený graf na obr. 15., který simuluje statické chování trychtýřového zásobníku.


44

Obr. 15. Průběh výšky hladiny v závislosti na změně přítoku

Na obr. 15. je vidět průběh výšky hladiny v závislosti na změně přítoku. Po odezvě na změnu přítoku a zvolení zadané konstanty ventilu bude výška hladiny stoupat než se dostane do rovnovážného stavu. Program statické charakteristiky je v příloze P V pod názvem Trychtýřový zásobník – statická charakteristika. Na obr. 16. je zobrazen vzhled tabulky pro zadávání vstupních hodnot z webové stránky.



45

2.3.2 Dynamická charakteristika Dynamickou charakteristiku budeme počítat z diferenciální rovnice popisující daný děj. Počáteční podmínky budou h(0 ) = h s . Rovnice má tvar

(

)

dh q v − k h ⋅ 4 H 2 = dt π ⋅ D 2 ⋅ h2

(2.13)

kde k je konstanta ventilu, která je vyjádřena z ustáleného stavu a má tvar k=

q vs

(2.14)

hs

V dynamice se zabýváme změnou výšky hladiny zásobníku při skokové změně přítoku do zásobníku v závislosti na časovém průběhu simulace. V následující tabulce jsou uvedeny použité hodnoty pro počítání dynamické charakteristiky. Tab. 9. Tabulka zadaných proměnných pro dynamickou charakteristiku Zadané proměnné Doba simulace Tsim Skoková změna průtoku qv

Hodnoty 150 [min] 0,7 [m3 min-1]

Ostatní zadané veličiny jsou stejné jako v tabulce na obr. 8. Výstupem je graf simulující dynamické chování trychtýřového zásobníku zobrazený na obr. 17.


46

2.8 2.7 2.6 2.5 ] m [ h

2.4 2.3 2.2 2.1 2

0

50

100

150

t[min]

Obr. 17. Průběh výšky hladiny v závislosti na čase

Na obr. 17. můžeme vidět časový průběh výšky hladiny v zásobníku. Z reakce na skokovou změně přítoku do zásobníku, lze vidět že výška hladiny v zásobníku stoupne. Průběh výšky hladiny se projevil přechodovou charakteristikou, která opět může vypovědět o zesílení soustavy či o časové konstantě. Program dynamické charakteristiky je uveden v příloze P VI pod názvem Trychtýřový zásobník – dynamická charakteristika. Vzhled webové stránky programu propojeného přes MWS je zobrazen na obr. 18.



47

2.4 KULOVÝ ZÁSOBNÍK NA KAPALINU S NEKONSTANTNÍM PRŮŘEZEM 2.4.1 Statická charakteristika Statickou charakteristiku řešíme stejně jako u modelu trychtýřového zásobníku. Rovnice počátečního ustáleného stavu je q vs = k h s

(2.15)

Ustálená výška hladiny vyjádřená z předešlé rovnice má tvar  qs h =  v  k s

  

2

(2.16)

V tomto zásobníku ve statické charakteristice budeme opět počítat změnu výšky hladiny v ustáleném stavu v závislosti na změně přítoku qv. V následujících tabulkách jsou uvedeny použité hodnoty pro počítání statické charakteristiky. Tab. 10. Tabulka zadaných proměnných pro statickou charakteristiku Zadané proměnné

Hodnoty

Začátek intervalu pro různý průtok qv

0 [m3 min-1]

Konec intervalu pro různý průtok qv

1 [m3 min-1]

Konstanta ventilu

2


Hodnoty

Průměr zásobníku D

5 [m]

Výška hladiny h s

3 [m]

Objemový průtok kapaliny q vs

0,6 [m3 min-1]

Vyhodnocením je graf zobrazený na obr. 19., který charakterizuje statické chování kulového zásobníku na základě zadaných hodnot.


48

Obr. 19. Průběh výšky hladiny v závislosti na změně přítoku Na obr. 19. můžeme pozorovat průběh výšky hladiny v závislosti na změně přítoku. Chování systému je podobné jak u trychtýřového zásobníku. Při změně přítoku a zadané konstanty bude výška hladiny stoupat než se dostane do ustáleného stavu. Program statické charakteristiky je v příloze P VII pod názvem Kulový zásobník – statická charakteristika. Na obr. 20. je zobrazen vzhled tabulky pro zadávání vstupních hodnot z webové stránky



49

2.4.2 Dynamická charakteristika Dynamickou charakteristiku budeme počítat z diferenciální rovnice popisující daný děj. S počátečními podmínkami h(0 ) = h s . Rovnice má tvar qv − k h dh = dt π ⋅ (D − h ) ⋅ h

(2.17)

kde k je konstanta ventilu, která je vyjádřena z ustáleného stavu a má tvar k=

q vs

(2.18)

hs

V dynamice se zabýváme stejným problémem jako u trychtýřového zásobníku. V následující tabulce jsou uvedeny použité hodnoty pro počítání dynamické charakteristiky. Tab. 12. Tabulka zadaných proměnných pro dynamickou charakteristiku Zadané proměnné Doba simulace Tsim Skoková změna průtoku qv

Hodnoty 800 [min] 0,75 [m3 min-1]

Ostatní zadané veličiny jsou stejné jako v tabulce na obr. 11. Vyhodnocením je graf zobrazený na obr. 21., který charakterizuje dynamické chování kulového zásobníku.


50

4.8 4.6 4.4 4.2 4 ] m [ h

3.8 3.6 3.4 3.2 3

0

100

200

300

400 t[min]

500

600

700

800

Obr. 21. Průběh výšky hladiny v závislosti na čase Na obr. 21. můžeme vidět časový průběh výšky hladiny v zásobníku. Reakce na skokovou změnu přítoku se opět projevila zvýšením výšky hladiny v zásobníku v podobě přechodové charakteristiky. Program dynamické charakteristiky je uveden v příloze P VIII pod názvem Kulový zásobník – dynamická charakteristika. Vzhled webové stránky programu propojeného přes MWS je zobrazen na obr. 22.

Obr. 22. Tabulka pro zadávání vstupních hodnot s tabulkou konstant V praktické části jsem používal literaturu [5], [6].


51

ZÁVĚR Hlavním cílem této bakalářské práce bylo vytvoření M-file souborů v programu Matlab, které budou schopny odsimulovat statické a dynamické charakteristiky technologických modelů s využitím MWS. Propojení programů přes MWS zajišťuje možnost chování modelů prakticky vyzkoušet a pozorovat přes webové rozhraní. Zde uživatel zadá vstupní hodnoty, ty jsou odeslány ke zpracování programem a vypočtené hodnoty se následně pošlou opět přes MWS na internetové stránky uživateli. M-file soubory bylo tedy potřeba vytvořit takovým způsobem, aby program byl schopen zpracovat a vyhodnotit jakékoliv libovolné zadání vstupních hodnot. Proto analýza, která programování předcházela, musela být provedena na zcela obecné úrovni. Po odvození matematického modelu průtočného výměníku tepla můžeme konstatovat, že se jedná o lineární model se soustředěnými parametry. Při simulaci statické charakteristiky se model převedl do ustáleného stavu a dynamická charakteristika se simulovala z odchylkového tvaru. Simulování statického chování vnitřní proměnné za daných podmínek se projevilo změnou teplot obou kapalin s přechodem do ustáleného stavu. Při změně průtoku chlazené kapaliny se teploty zvýšily a při změně průtoku chladící kapaliny obě teploty zase klesly. Z ustáleného stavu těchto teplot se dále může určit poloha pracovního bodu. Při simulaci dynamického chování při nulových počátečních podmínkách se průběh teplot projevil přechodovou charakteristikou, která může vypovědět o časové konstantě a o zesílení soustavy. U matematického modelu zásobníku na kapalinu zapojených za sebou se ukázalo, že se jedná o nelineární model. Tento model bylo potřeba zlinearizoval a vyjádřit odchylkový tvar potřebný pro simulaci dynamiky. Při simulaci statické charakteristiky výšky hladin stoupnou a ustálí se na nové hodnotě rovnovážného stavu, ze kterého se může určit poloha pracovního bodu v okolí kterého tento model linearizujeme. U simulace dynamického chování při nulových počátečních podmínkách se výška hladin projevila přechodovou charakteristikou, ze které se opět může zhodnotit časová konstanta a zesílení soustavy. U trychtýřového i kulového zásobníku se rovněž jedná o nelineární modely. Při simulaci statického chování u těchto modelů se výšky hladin ustálí na vyšší hodnotě, než ve které byly před změnou přítoku. Při pozorování dynamického chování za počáteční podmínky


52

odpovídající ustálené výšce hladiny se změna výšek hladin v trychtýřovém i kulovém zásobníku projevila přechodovou charakteristikou. Modely se chovají při simulacích adekvátně podle zadaných vstupních veličin a výsledky jsou reálné. Mohou tak uživateli sloužit pro lepší pochopení problematiky statických a dynamických charakteristik.


53

SUMMARY The main aim of this thesis was to create M-files using the Matlab application that would be able to simulate static and dynamic characteristics of the technological models applying MWS. The interconnection of the programs via MWS enables the user to examine the models in practice and watch it through the Internet interface. There, the user might enter chosen input values that are then sent to be processed by the program and the computed values are consequently sent back to the user’s web pages via MWS. Thus, it was necessary to formulate the M-files in such a way that the program would be able to process and compute an arbitrary selection of any input values. Therefore, the analysis that preceded programming had to be performed at a general level. After the mathematic model of the flow heat exchanger is defined, we may observe that it represents a linear model with condensed parameters. When simulating the static characteristic, the model transformed into the equilibrium state and the dynamic characteristic was simulated from the deviation state. The simulation of the static behaviour of the internal variable was manifested by change of temperatures of both liquids with transition to the equilibrium state. This state allows to determine the position of the working point. When simulating the dynamic behaviour using null initial conditions the temperature value progress was demonstrated by transitional characteristic that can document the time constant and consolidation of the system. The mathematic model of the system with two water tanks is linear. It was necessary to linearise this model and to define the deviation needed for the dynamism simulation. When simulating the static characteristic, the levels of the liquid rise and stabilise at a new position of the equilibrium state that might be used to determine the position of the working point located in the area of the actual model linearisation. When simulating the dynamic behaviour in the null initial conditions the liquid levels were demonstrated by transition characteristic that might be used again for the determination of the time constant and the consolidation of the system. Both the spherical and conical tanks represent non-linear models. When simulating the static behaviour, the levels of the liquids arrive at a new equilibrium state where their levels elevate. Monitoring the dynamic behaviour when the initial conditions were set to a level


54

that corresponded to the stabilised liquid level, the difference between the liquid levels in both the spherical and conical tanks was demonstrated by the transitional characteristic. The models behave adequately in simulations accordingly to the input values and the results are real. Thus, they may serve as a tool for better comprehension of the issue of static and dynamic characteristics.


55

SEZNAM POUŽITÉ LITERATURY [1] TOBOLA, Petr. Matematické základy – Numerické řešení soustavy rovnic [online], poslední aktualizace: 2007-05-23 [cit. 2007-05-23]. Dostupné z WWW: < http://www.fi.muni.cz/~ptx/PA010/Slides/PA010_11.pdf > [2] MACUR, Jiří . Simulace dynamických systémů v jazyce Java [online], poslední aktualizace: 2006-05-01 [cit. 2007-05-15].Dostupné z WWW: [3] Západočeská univerzita v Plzni, citace: před_matl13.doc [online], poslední aktualizace: 2004-10-20 [cit. 2007-05-15].Dostupné z WWW: < http://home.zcu.cz/~ruz/zs/pred_matl3.doc> [4] Noskievič, P. Modelování a identifikace systémů. MONTANEX a.s., Ostrava, 1999, ISBN 80-7225-030-2. [5] BAKOŠOVÁ Monika, FIKAR Miroslav, ČIRKA Luboš. Základy automatizácie. ISBN 80-227-1831-9 [6] KARBAN, Pavel. Výpočty a simulace v programech MATLAB a Simulink. Vydal: Computer Press, a.s. ISBN 978-80-251-1448-3


SEZNAM POUŽITÝCH SYMBOLŮ A ZKRATEK MWS

Matlab web server.

ODR

Obyčejné diferenciální rovnice.

MPI

Metoda prosté iterace

56


57

SEZNAM OBRÁZKŮ Obr. 1. Průtočný výměník tepla s chlazením v plášti…………………………………… 10 Obr. 2. Zásobníky na kapalinu zapojené za sebou…………………………………….... 14 Obr. 3. Trychtýřový zásobník s nekonstantním průřezem………………………….…... 18 Obr. 4. Kulový zásobník s nekonstantním průřezem………………………………….... 21 Obr. 5. Průběhy ustálených teplot T s a Tcs v závislosti na průtocích q a qc…………….

31

Obr. 6. Tabulka pro zadávání vstupních hodnot s tabulkou konstant…………………..

32

Obr. 7. Průběh ustálené teploty chlazené kapaliny v závislosti na čase………………...

33

Obr. 8. Průběh ustálené teploty chladící kapaliny v závislosti na čase…………………

34

Obr. 9. Tabulka pro zadávání vstupních hodnot s tabulkou konstant…………………..

35

Obr. 10. Průběhy ustálených výšek hladin h1s a h2s v závislosti na průtocích q1v a q2v……………………………………………………………….……....

38

Obr. 11. Tabulka pro zadávání vstupních hodnot s tabulkou konstant…………………

39

Obr. 12. Průběh ustálené výšky hladiny v prvním zásobníku v závislosti na čase……...

40

Obr. 13. Průběh ustálené výšky hladiny v druhém zásobníku v závislosti na čase……..

41


42

Obr. 15. Průběh výšky hladiny v závislosti na přítoku…………………………………

44


44

Obr. 17. Průběh ustálené výšky hladiny v závislosti na čase…………………………..

46


46

Obr. 19. Průběh výšky hladiny v závislosti na průtoku………………………………..

48


48

Obr. 21. Průběh ustálené výšky hladiny v závislosti na čase…………………………..

50


50


58

SEZNAM TABULEK Tab. 1. Tabulka zadaných proměnných pro statickou charakteristiku…………………... 29 Tab. 2. Tabulka zadaných veličin………………………………………………………..

30

Tab. 3. Tabulka zadaných proměnných pro dynamickou charakteristiku……………….

33

Tab.4. Tabulka zadaných proměnných pro statickou charakteristiku…………………… 37 Tab. 5. Tabulka zadaných veličin………………………………………………………..

37


40

Tab. 7. Tabulka zadaných proměnných pro statickou charakteristiku…………………..

43

Tab. 8. Tabulka zadaných veličin…………………………………….………………….

43


45

Tab. 10. Tabulka zadaných proměnných pro statickou charakteristiku.………………… 47 Tab. 11. Tabulka zadaných veličin………………………………………………………

47

Tab. 12. Tabulka zadaných proměnných pro dynamickou charakteristiku……………… 49


SEZNAM PŘÍLOH PI:

Průtočný výměník – statická charakteristika

PII:

Průtočný výměník – dynamická charakteristika

PIII:

Zásobníky na kapalinu – statická charakteristika

PIV:

Zásobníky na kapalinu – dynamická charakteristika

PV:

Trychtýřový zásobník – statická charakteristika

PVI:

Trychtýřový zásobník – dynamická charakteristika

PVII: Kulový zásobník – statická charakteristika PVIII: Kulový zásobník – dynamická charakteristika

59

PŘÍLOHA P I: PRŮTOČNÝ VÝMĚNÍK – STATICKÁ CHARAKTERISTIKA. vymenik_statika.m clc clear all close all V=1.2;

% m^3

Vc=0.64;

% m^3

ro=985;

%kg/m^3

roc=998;

%kg/m^3

cp=4.05;

%kJ/kg*K

cpc=4.18;

%kJ/kg*K

F=5.5;

%m^2

a=43.5;

%kJ/m^2*min*K

Tv_s=323;

%K

Tcv_s=293;

%K

q_z = input(' zadejte zacatek intervalu pro ruzny prutok q: '); q_kr = input(' zadejte krok intervalu pro ruzny prutok q: '); q_k = input(' zadejte konec intervalu pro ruzny prutok q: '); qc_z = input(' zadejte zacatek intervalu pro ruzny prutok qc: '); qc_kr = input(' zadejte krok intervalu pro ruzny prutok qc: '); qc_k = input(' zadejte konec intervalu pro ruzny prutok qc: ');

q = q_z:q_kr:q_k; qc = 0.03; for i = 1:length(q) T_s(1,i)=(qc*roc*cpc*Tcv_s*F*a+q(i)*ro*cp*Tv_s*F*a+qc*roc*cpc*q(i)*ro*cp*Tv_s)/(q(i)*ro*cp*F*a+q c*roc*cpc*q(i)*ro*cp+qc*roc*cpc*F*a); end subplot (2,2,2); plot(q,T_s(1,:)); xlabel('q[m^3/min]'); ylabel('T_s[K]'); for i = 1:length(q) Tc_s(1,i)=(q(i)*ro*cp*qc*roc*cpc*Tcv_s+qc*roc*cpc*Tcv_s*F*a+q(i)*ro*cp*Tv_s*F*a)/(q(i)*ro*cp*F*a +qc*roc*cpc*q(i)*ro*cp+qc*roc*cpc*F*a); end subplot (2,2,1); plot(q,Tc_s(1,:)); xlabel('q[m^3/min]');

ylabel('Tc_s[K]');

q = 0.08; qc = q_z:q_kr:q_k; for i = 1:length(qc) T_s(1,i)=(qc(i)*roc*cpc*Tcv_s*F*a+q*ro*cp*Tv_s*F*a+qc(i)*roc*cpc*q*ro*cp*Tv_s)/(q*ro*cp*F*a+qc(i )*roc*cpc*q*ro*cp+qc(i)*roc*cpc*F*a); end subplot (2,2,4); plot(qc,T_s(1,:)); xlabel('qc[m^3/min]'); ylabel('T_s[K]'); for i = 1:length(qc) Tc_s(1,i)=(q*ro*cp*qc(i)*roc*cpc*Tcv_s+qc(i)*roc*cpc*Tcv_s*F*a+q*ro*cp*Tv_s*F*a)/(q*ro*cp*F*a+q c(i)*roc*cpc*q*ro*cp+qc(i)*roc*cpc*F*a); end subplot (2,2,3); plot(qc,Tc_s(1,:)); xlabel('qc[m^3/min]'); ylabel('Tc_s[K]');

PŘÍLOHA P II: PRŮTOČNÝ VÝMĚNÍK – DYNAMICKÁ CHARAKTERISTIKA. vymenik.m function dx = vymenik(t,x) global u A B dx = zeros(2,1); dx(1) = A(1,1)*x(1) + A(1,2)*x(2) + B(1,1)*u(1); dx(2) = A(2,1)*x(1) + A(2,2)*x(2) + B(2,2)*u(2); vymenik_hlavni.m clear all close all clc Tsim = input(' Zadejte casovy prubeh simulace: '); u(1) = input(' Zadejte skokovou zmenu prutoku q: '); u(2) = input(' Zadejte skokovou zmenu prutoku qc: '); global u A B vymenik_konstanty; A(1,1)=-(q/V + (F*a)/(V*ro*cp)); A(1,2)=(F*a)/(V*ro*cp); A(2,1)=(F*a)/(Vc*roc*cpc); A(2,2)=-(qc/Vc + (F*a)/(Vc*roc*cpc)); B(1,1)=q/V; B(2,2)=qc/Vc; [T,Y] = ode45(@vymenik,[0 Tsim],[0 0]); plot(T,Y(:,1)) xlabel('t[min]') ylabel('T[K]') figure plot(T,Y(:,2),'r--') xlabel('t[min]') ylabel('Tc[K]') vymenik_konstanty.m q=0.08;

% m^3/min

qc=0.03;

% m^3/min

V=1.2;

% m^3

Vc=0.64;

% m^3

ro=985;

%kg/m^3

roc=998;

%kg/m^3

cp=4.05;

%kJ/kg*K

cpc=4.18;

%kJ/kg*K

F=5.5;

%m^2

a=43.5;

%kJ/m^2*min*K

Tv_s=323;

%K

Tcv_s=293;

%K

PŘÍLOHA P III: ZÁSOBNÍKY NA KAPALINU – STATICKÁ CHARAKTERISTIKA. zasobniky_statika.m clear all close all clc q1v_s = 0.5;

% m^3/s

q2v_s = 0.3;

% m^3/s

h1_s=1.5;

%m

h2_s=1;

%m

q1v_z = input(' zadejte zacatek intervalu pro ruzny prutok q1v: '); q1v_kr = input(' zadejte krok intervalu pro ruzny prutok q1v: '); q1v_k = input(' zadejte konec intervalu pro ruzny prutok q1v: '); q2v_z = input(' zadejte zacatek intervalu pro ruzny prutok q2v: '); q2v_kr = input(' zadejte krok intervalu pro ruzny prutok q2v: '); q2v_k = input(' zadejte konec intervalu pro ruzny prutok q2v: ');

q1v = q1v_z:q1v_kr:q1v_k; q2v_s = 0.3; k2 = (q1v_s+q2v_s)/sqrt(h2_s); for i = 1:length(q1v) h2s(1,i) = ((q1v(i)+q2v_s)/k2)^2; end subplot (2,2,2); plot(q1v,h2s(1,:)); xlabel('q1v[m^3/s]'); ylabel('h2^s[m]'); k1 = (q1v_s)/sqrt(h1_s-h2_s); for i = 1:length(q1v) h1s(1,i) = h2s(1,i)+(q1v(i)/k1)^2; end subplot (2,2,1); plot(q1v, h1s(1,:)); xlabel('q1v[m^3/s]'); ylabel('h1^s[m]');

q1v_s = 0.5; q2v = q2v_z:q2v_kr:q2v_k;

for i = 1:length(q2v) h2s(1,i) = ((q1v_s+q2v(i))/k2)^2; end subplot (2,2,4); plot(q2v,h2s(1,:)); xlabel('q2v[m^3/s]'); ylabel('h2^s[m]'); for i = 1:length(q2v) h1s(1,i) = h2s(1,i)+(q1v_s(1)/k1)^2; end subplot (2,2,3); plot(q2v, h1s(1,:)); xlabel('q2v[m^3/s]'); ylabel('h1^s[m]');

PŘÍLOHA P IV: ZÁSOBNÍKY NA KAPALINU – DYNAMICKÁ CHARAKTERISTIKA. zasobniky.m function dx = zasobniky(t,x) global u A B dx = zeros(2,1); dx(1) = A(1,1)*x(1) + A(1,2)*x(2) + B(1,1)*u(1); dx(2) = A(2,1)*x(1) + A(2,2)*x(2) + B(2,2)*u(2); zasobniky_hlavni.m clear all close all clc Tsim = input(' Zadejte cas simulace: '); u(1) = input(' Zadejte skokovou zmenu prutoku q1v: '); u(2) = input(' Zadejte skokovou zmenu prutoku q2v: ');

global u A B zasobniky_konstanty; A(1,1)=-K1/F1; A(1,2)=K1/F1; A(2,1)=K1/F2; A(2,2)=-((K1+K2)/F2); B(1,1)=1/F1; B(2,2)=1/F2; [T,Y] = ode45(@zasobniky,[0 Tsim],[0 0]); plot(T,Y(:,1)) title('graf c.1:Zasobniky na kapalinu') xlabel('t[s]') ylabel('h1[m]') figure plot(T,Y(:,2),'r--') title('graf c.2:Zasobniky na kapalinu') xlabel('t[s]') ylabel('h2[m]') zasobniky_konstanty.m q1v_s=0.5;

% m^3/s

q2v_s=0.3;

% m^3/s

h1_s=1.5;

%m

h2_s=1;

%m

F1=2;

% m2

F2=2;

% m2

K1=q1v_s/(2*(h1_s-h2_s)); K2=q2v_s/(2*h2_s);

PŘÍLOHA P V: TRYCHTÝŘOVÝ ZÁSOBNÍK – STATICKÁ CHARAKTERISTIKA. trychtyrovy_zas_statika.m clear all close all clc D=2;

%m

H=3;

%m

h_s=2; qv_s=0.6;

%m % m^3/min

qv_z = input(' zadejte zacatek intervalu pro ruzny prutok qv: '); qv_k = input(' zadejte konec intervalu pro ruzny prutok qv: '); k = input(' zadejte konstantu ventilu: ');

qv = qv_z:0.1:qv_k; k=qv_s/sqrt(h_s); for i = 1:length(qv) h(1,i) = ((qv(i))/k)^2; end plot(qv, h(1,:)); xlabel('qv[m^3/min]'); ylabel('h[m]');

PŘÍLOHA P VI: TRYCHTÝŘOVÝ ZÁSOBNÍK – DYNAMICKÁ CHARAKTERISTIKA. trychtyrovy_zasobnik.m function dh = trychtyrovy_zasobnik(t,h,qv) dh = zeros(1,1); D=2;

%m

H=3;

%m

h_s=2; qv_s=0.6;

%m % m^3/min

k=(qv_s/sqrt(h_s)); dh = (qv-k*sqrt(h))/((pi*D^2*h^2)/(4*H^2)); trychtyrovy_zasobnik_hlavni.m clear all close all clc qv = input(' Zadejte skokovou zmenu prutoku qv: '); Tsim = input(' Zadejte cas simulace Tsim: '); [T,h] = ode45(@trychtyrovy_zasobnik,[0 Tsim],2, [], qv); plot(T,h) title('Trychtyrovy zasobnik na kapalinu') xlabel('t[min]') ylabel('h[m]')

PŘÍLOHA P VII: KULOVÝ ZÁSOBNÍK – STATICKÁ CHARAKTERISTIKA. kulovy_zas_statika.m clear all close all clc D=5;

%m

h_s=3; qv_s=0.6;

%m % m^3/min

qv_z = input(' zadejte zacatek intervalu pro ruzny prutok qv: '); qv_k = input(' zadejte konec intervalu pro ruzny prutok qv: '); k = input(' zadejte konstantu ventilu: '); qv = qv_z:0.1:qv_k; k=qv_s/sqrt(h_s); for i = 1:length(qv) h(1,i) = ((qv(i))/k)^2; end plot(qv, h(1,:)); xlabel('qv[m^3/min]'); ylabel('h[m]');

PŘÍLOHA P VIII: KULOVÝ ZÁSOBNÍK – DYNAMICKÁ CHARAKTERISTIKA. kulovy_zasobnik.m function dh = kulovy_zasobnik(t,h,qv) dh = zeros(1,1); D=5; h_s=3; qv_s=0.6;

%m %m % m^3/min

k=qv_s/sqrt(h_s); dh = (qv-(k*sqrt(h)))/(pi*(D-h)*h); kulovy_zasobnik_hlavni.m clear all close all clc qv = input(' Zadejte skokovou zmenu prutoku qv: '); Tsim = input(' Zadejte cas simulace Tsim: '); [T,Y] = ode45(@kulovy_zasobnik,[0 Tsim],3, [], qv); plot(T,Y) title('Kulovy zasobnik na kapalinu') xlabel('t[min]') ylabel('h[m]')

Vytvoření skriptů pro webové rozhraní předmětu Analýza a simulace technologických procesů

Recommend Documents