fyzikální principy modelování - zákony zachování a konstitutivní vztahy, empirické uzávěry; druhy modelů podle hloubky (počet souřadnic) a šířky

Matematické modelování oběhu motorů Základy modelování oběhu motoru

Jan Macek ČVUT v Praze Výzkumné centrum Josefa Božka

Základy modelování oběhu motoru

1.

Úvod.

2.

Základní prvek modulárního modelu s využitím MKO (FVM) - podklad formulace zákonů zachování s konvekcí, difusí a zdrojovými členy.

3.

Zákon zachování pro hmotnost látek: popis chemické stechiometrie a kinetiky.

4.

Zákon zachování hybnosti.

5.

Konstitutivní vztahy - rovnice stavu a jejich derivace.

6.

Zákon zachování energie.

7.

Formulace rovnice pro společný tlak vícezónového modelu.

8.

Princip inverzního algoritmu pro vyhodnocení průběhu hoření.

9.

Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. Matematické modelování oběhu motorů

2

Základy modelování oběhu motoru

10.

Numerické metody a specifické vlastnosti pravých stran zákonů zachování (stiffness).

11.

Dosazení do pravých stran ODR: - hoření - určení počátku spontánního hoření (zejména průtahu vznětu), semiempirický model s Vibeho funkcemi, určení jejich parametrů, detonační spalování - prostup tepla stěnami (přestupní součinitel, tepelné odpory stěn, výpočet teploty stěn) - průtoky ventily - průtok a účinnost turbiny, model turbodmychadla

12.

Mechanická účinnost motoru - výpočet ztrát.

13.

Příklady výstupů pro výfukový systém a turbodmychadlo (interpretace výsledků)

Matematické modelování oběhu motorů

3

Základy modelování oběhu motoru 14.

Soustava rovnic pro 1-D nestacionární proudění a její vlastnosti. Charakteristiky, řešení Cauchyho úlohy.

15.

Vznik nestacionární rázové vlny. Pojem slabého řešení, důležitopst MKO.

16.

Nestacionární sdílení tepla v polomasivu, typy okrajových úloh. Iterační řešení okrajové podmínky 3. druhu při proměnlivé teplotě i součiniteli přestupu tepla.

17.

Další příklady interpretace výstupů

•

adiabatický motor a motor bez chlazení

•

rozbor hoření při simulaci

•

rozbor naplnění válce, pomůcky pro odhad dynamických jevů a vnitřního chlazení

18.

Jak kalibrovat model při omezených znalostech z experimentů.


4

Základy modelování oběhu motoru: 1. Úvod – fyzikální principy modelování - zákony zachování a konstitutivní vztahy, empirické uzávěry; – druhy modelů podle hloubky (počet souřadnic) a šířky (okolí pracovní látky v motoru - rozsah agregátů motoru a respektování časových měřítek); – lagrangeovský a eulerovský (bilanční) přístup, zónové a vícerozměrné modely.


5

Základy modelování oběhu motoru: 1. Úvod • druhy modelů podle hloubky – počet souřadnic v prostoru 0-D, 1-D, 2.5-D a 3-D – ... a v čase - vždy nestacionární z hlediska vysokých frekvencí, při modelování dynamiky pevných těles i z hlediska nízkých frekvencí (přechodové režimy)

• ... a podle šířky – rozsah agregátů motoru: • • • • • • •

válec potrubí turbodmychadlo, chladič plnicího vzduchu zařízení pro tvorbu směsi klikový mechanismus regulátory mechanismus rozvodu ...


6

1. Úvod - druhy modelů: lagrangeovský a eulerovský (bilanční) přístup, zónové a vícerozměrné modely 0-D a 3-D CFD modely lze popsat obdobným FVM algoritmem s různou fyzikální hloubkou. Pro 1-D je typický zónový přístup se zahrnutím setrvačných sil. • Lagrange - vazba na látku, • Euler - vazba na materiální hranice dílů motoru.

Eulerovský model konečného objemu pro vazbu na 1-D (0-D) model Matematické modelování oběhu motorů

7

Základy modelování oběhu motoru: Úvod - druhy modelů Lagrangeovský multizónový model (vázaný na konkrétní část látky a její termodynamické i chemické změny) a jeho zobecnění

Základní nevýhoda: omezený počet sousedících zón, jinak deformace celé sítě


Lagrangeovský model konečného objemu rozšířený o propustnou stěnu umožňující interakci s okolím 8

2. Základ eulerovského modelu Bilance v konečném objemu (intuitivní odvození) Nejjednodušší názorný příklad: zákon zachování hmotnosti pro složku plynu (látku) x : • změna hmotnosti uvnitř objemu v čase = • přítok z okolí (znaménko a stav v okolí jsou důležité) + • - odtok do okolí (znaménko a stav v objemu!) + • + změna v důsledku zdrojů nebo propadů uvnitř kontrolního objemu (zde chemické reakce, měnící složení, ne však celkovou hmotnost) • • dm x dm xCH = σ px m p − σ x m o + dt dt Matematické modelování oběhu motorů

9

2. Struktura modulárního modelu s využitím MKO. Cíl: Univerzální rovnice pro 0 - 3-D modely. Prostředek: Různá formulace spojovacích členů Obecně se odvodí rovnice zachování pro extenzitní veličinu Φ (měrnou veličinu φ) v pohyblivém konečném objemu (KO). Leibnizův přístup je vhodný pro kontinuum (odvození pomocí Reynoldsova transportního teorému + Leibnizovy věty)

d d t

→ → → φρ δ V = ∫∫ φ ρ  w G − w F  . n δ A + ∫∫∫   Vj A j = ∂V j

Vazby mezi KO: - objemový tok >0...přítok do j - tok přenášené veličiny: jak nutno interpolovat? Matematické modelování oběhu motorů

→

→

∫∫ χ ρ ∇ φ . n δ A + Pφ

,j

− Dφ , j

Aj

→ → → A V j ,i =  w G , j ,i − w F , j ,i . n j, i j, i •  • Φ j ,i = V j ,i φ j ,i ρ j ,i •

10

2. Struktura modulárního modelu s využitím MKO (FVM) vhodná pro moduly ohraničené pevnými stěnami. Základní toky. Základní modul modelu: • hmotnostní toky  → → m = ρ V = ρ  wF . n  A   •

•

• hybnostní toky • →  → → → m w F = ρ  wF . n  A w F 



• působící síly • entalpické toky

[

T T T  → → m h = ρ  wF . n  A c p T (T − Tref ) + l T + ∆hCHref 0 ref   •

]

• upwind popis přenášené veličiny •  1 + signV  Φ j ,i = V j , i  φ i ρ i 2  • sdílení tepla  •

•

j ,i

•

1 + signV +φjρ j 2

γj,i

• sdílení práce Matematické modelování oběhu motorů

j ,i

   

δj,i 11

3. Zákon zachování hmotnosti pro látky směsi • • dm x dm xCH = σ px m p − σ x m o + dt dt

Zákon zachování hmotnosti: změna = přítok látky-odtok látky+zdroje látky (viz snímek 9) se aplikuje především na: – Průtok ventily a turbinou - viz dále – Vstřik paliva a tvoření směsi - mimo tento kurs – Chemické děje – zde se uplatní zdrojové členy založené na • rychlosti chemické transformace (empirické modely, rovnováha, kinetika) • stechiometrických vztazích (zachování hmotnosti prvků)


12

3. Zákon zachování hmotnosti pro látky Vektorové vyjádření lineárních kombinací (vlastností směsí)  mO 2  1 Vektorový zápis pro výpočet     m N 2  extenzitních veličin z intenzitních 1 T { s m} =   ; {1} =   ; {1} = {1 1 1 ...} (lineární kombinace platné pro mCO 2  1 směsi ideálních plynů).     M  M Používají se maticové součiny, s T aplikované obvykle jen na řádkové mr = {s m} {s r} = ∑ mi ri 1 transponované sloupcové - 1*s a  m1 sloupcové vektory s*1. Při součinu   řádkový*sloupcový je výsledkem m 2  skalární součin, tedy 1 hodnota.    M Násobení matic m*n a n*k: řádky *   mn sloupce, vždy stejnolehlé prvky {s σ } =  T  vynásobit a sečíst, vznikne matice {s 1} {s m} m*k. T mr = {s m} {s r} Matematické modelování oběhu motorů

13

3. Zákon zachování hmotnosti pro látky Zachování látek (tokem přenášená je hmotnostní koncentrace) : • přítok γ=0 nebo 1, V°>0, • odtok δ=0 nebo 1, V°<0 • zdroje (propady): chemická kinetika.

{s c} = {s ρ } = d {s m} j d t

{s m} Vi

{s m} j  d {s m}CH  {s m}i  = ∑V j , i . .γ j, i + .δ j, i  +  V  Vj d t i  i  Nj •

j

Zdrojový člen chemických reakcí: • Guldberg-Waage - Arrhenius pro hlavní reakční proměnnou; • stechiometrická transformace na ostatní složky reakce.


14

3. Zákon zachování hmotnosti pro látky Vektorové vyjádření vlastností směsí a chemických reakcí Výběr hlavní reakční proměnné - jedné z látek v reakci - je libovolný, avšak obvykle se volí mezi reagenty. Při její změně jsou změny množství (zde hmotnosti, ale lze použít i látková množství, tedy kmol) ostatních látek vázány stechiometrickými poměry. Příklad pro sumární reakci oxidace CO:

CO + 1 2 O2 → CO2 { s ∆m} = {s C} ∆mmain

28kg + 16kg → 44kg ∆mCO  − 28 28     ∆ m = − 16 28   ∆mCOreag O2  ∆m   44 28    CO2  

Vektor hlavních reakčních proměnných má počet složek r<s, pro jednu reakci je samozřejmě r=1. Pro sumární výsledek více reakcí je vektor hlavních proměnných třeba násobit transformační maticí ||C||, obsahující sloupce stechiometrických koeficientů příslušných reakcí. Matematické modelování oběhu motorů

15

3. Zákon zachování hmotnosti pro látky - stechiometrická transformační matice, hlavní složky reakcí Příklad transformační části matice ||C|| pro disociaci produktů hoření a pro sumární podchycení hoření oktanu; matice pokračuje oxidací N2 a lze ji libovolně rozvádět dále) R eaction &

H2

H 2O

CO

CO2

C 8H 18

O2

-16/2

+ 16/18

-16/28

+ 16/44

-400/114

H2

-2/2

2/18

H 2O

18/2

-18 /1 8

m ain com p on ent ⇒

162/114

CO

-28/28

28 /4 4

CO2

4 4/28

-44/44

1

-114/114

C 8H 18 N2 N⋅ NO O⋅


352/114

16

3. Zákon zachování hmotnosti pro látky Chemická kinetika • chemická kinetika

d {s m}CH

• stechiometrická transformace;

d t

j

=

s*r C .

d {r m}CH

j

d t

Zdrojový člen chemických reakcí Guldberg-Waage-Arrhenius •

příklad pro ternární pro nejobecnější ternární reakci (obvyklé jsou však binární). Koncentrace jsou vždy pro reagenty - tady není volba možná (na rozdíl od volby hlavní reakční proměnné na pravé straně – ta může být i produktem).

d { r m}CH d t

j

 = V . r K (T ).∏ C Rx R R 

A   − Tj −b    K . T j .e x y z . mX <0 . mY < 0 . m Z< 0   = r x + y + z −1   Vj    j


17

3. Zákon zachování hmotnosti pro látky rekapitulace Zachování látek: d {s m} j N • = ∑V • přítok γ=0 nebo 1, V°>0, j

• odtok δ=0 nebo 1, V°<0 • chemická kinetika.

d t

i

{s m} j  {s m}i  d {s m}CH . + .δ j, i  + j , i .  V γ j, i  V d t i j   d {s m}CH

Zdrojový člen chemických reakcí: • stechiometrická transformace;

d t

d {r m}CH

j

=

s*r C .

d {r m}CH

j

j

d t

  −b  K . T j .e x y z = r . . . m X < 0 m Y < 0 m Z< 0  x + y + z −1 V j    j A − Tj

j

d t

• Guldberg-Waage-Arrhenius.

Zdrojový člen chemických změn lze rozdělit pro podstatné – sig – a ostatní chemické reakce (důležité pro inverzní algoritmy): obvykle podstatné právě pro r sig=1

d {s m}CHj d t


=

s*1 C

d msig

18

d t

j

+

s*r -1 C

d {r −1 m}CH d t

j

Zákony zachování pro konečné objemy: 4. Zákon zachování hybnosti ve vektorové formě platné pro kartézské souřadnice (není zapotřebí pro 0-D) Zákon zachování hybnosti je použitelný pro minimálně 1-D, pokud se mají respektovat setrvačné síly uvnitř KO. Principiálně je nutný pro 2-D a více (i v případech, kdy jsou setrvačné síly malé). Pro křivočaré souřadnice nutno doplnit transformační členy zrychlení (odstředivé, část Coriolosova atd.). Plochy, na nichž působí tlak (síla zevně na tekutinu), jsou orientovány vnější normálou. Plochy Aτ , na nichž působí smyková napětí, nutno orientovat ve shodě s konstitutivní rovnicí pro vliv viskozity, případně pro další napětí v tekutině. Obvykle se použije Stokesova hypotéza o tenzoru vazkých napětí

doplněná pro turbulentní proudění o Boussinesqovu turbulentní viskozitu (vypočtenou z modelu turbulence) nebo Reynolds Stress Model. →

Nj  • d wj = ∑ V mj. d t i  

→ →   → mi  m j → +wj j ,i . wi δ j, i  − pi Ai , j n i , j − τ i, j . Aτ i, j  −  V γ j, i Vj i   

→ → d mj − wj . + ∑ F j ,k ,kapka→ plyn + ∑ F j , p , prekázka→ plyn + ... d t kj pj →


19

5. Konstitutivní rovnice pro konečné objemy Derivace stavové rovnice plynu Obecný technický tvar (jediná rovnice pro směs s empirickými součiniteli odvozenými z vlastností jednotlivých látek nelineární transformací). p = p (T , V , { s m}) ;

∂p dT j ∂p dV j ∂p dmi , j = + +∑ . . . ; dt ∂T dt ∂V dt ∂mi dt

dp j dp j

dT j dt

=

dt

−

∂p ∂V

.

dV j

s

−∑

dt ∂p

i

∂p ∂mi

∂T

.

dmi , j dt T p jV j = { s r } { s m }j Tj

T dV j dT j { s r } d{ s m }j 1  T  T j +  { s r } .{ s m } j . = − pj. dt Vj  dt dt dt Vj d pj d { s m }j dV j T dT j V j . dt − { s r } . dt .T j + p j . dt = T dt { s r } .{ s m } j

d pj

Ideální plyn (lineární kombinace vlastností složek)


20

  

5. Konstitutivní rovnice pro konečné objemy Derivace stavové rovnice pro vodní páru Z technického tvaru rovnice stavu pro jednosložkovou soustavu s proměnlivou hmotností plyne při daném měrném objemu V

d dp(T , v) ∂p dT ∂p m ∂p dT ∂p 1 dV ∂p V dm = + . = + − . . dt ∂T dt ∂v dt ∂T dt ∂v m dt ∂v m 2 dt dp ∂p 1 dV ∂p V dm − + dT dt ∂v m dt ∂v m 2 dt = ∂p dt ∂T

Rovnice IAPWS 97 (jediná rovnice pro Gibbsovu funkci – v bezrozměrném tvaru γ a v závislosti na normovaném tlaku a teplotě; derivacemi se pak odvozují další stavové veličiny) se pak musí derivovat jako implicitní funkce pro neznámý tlak v = f ( p (v, T ), T ) =

Tref RT p  ∂γ  ; τ = γπ ; γπ =   ; π = p ref p ref T  ∂π  T

RTTref R ∂f γπ − γ 2 πτ  RT  p p T p ref τγ πτ − γ π ∂p 1 ∂p ref ref ∂T = − = = ; = − = γ  ππ ∂f RT ∂v ∂f  p ref 2 ∂T T γ ππ  γ 2 ππ ∂p ∂p p ref −1


21

5. Konstitutivní rovnice pro konečné objemy Derivace stavové rovnice pro vodní páru Přímo a jednodušeji

V dv 1 dV V dm df (T , p) ∂f dT ∂f dp = m= − 2 = = . + . dt dt m dt m dt dt ∂T dt ∂p dt 1 dV V dm ∂f dp − − dT m dt m 2 dt ∂p dt = ∂f dt ∂T d

Rovnice IAPWS 97 (jediná rovnice pro Gibbsovu funkci – v bezrozměrném tvaru γ a v závislosti na normovaném tlaku a teplotě; derivacemi se pak odvozují další stavové veličiny) se pak musí derivovat jako implicitní funkce pro neznámý tlak v = f ( p(v, T ), T ) =

Tref RT p  ∂γ  ; τ = γπ ; γπ =   ; π = p ref p ref T  ∂π  T

RTTref ∂f RT ∂f R = γ ; = γ − γ πτ ππ π ∂p p ref 2 ∂T p ref p ref T 2


22

5. Konstitutivní rovnice pro konečné objemy Derivace stavové rovnice pro mokrou vodní páru V mokré páře je teplota syté páry funkcí tlaku. Parciální derivace stavu jsou jen podle tlaku a suchosti v − v′ h = h( p, x ) = xh′′( p ) + (1 − x )h′( p ) ; v = xv ′′( p ) + (1 − x )v ′( p ) ; x = v ′′ − v′ dh dx  dh′′ dh′  dp  = (h′′ − h′) +  x + (1 − x ) dt dt  dp dp  dt 1 dV V dm x dv ′′ + (1 − x ) dv ′ − dx m dt m 2 dt dp dp dp − = dt v ′′ − v′ v ′′ − v′ dt

m

 (h′′ − h ′)  dv ′′ dh (h ′′ − h ′)  dV V dm  dv ′   dh′′ dh′  dp  x  +  x  = − + (1 − x ) + (1 − x )   − m ′ ′ ′ dt v ′′ − v ′  dt m dt  v v dp dp dp dp −     dt   V (h′′ − h ′)  dv ′′ dv ′   dh′′ dh′  dp  x  +  x  = m + + (1 − x ) + (1 − x ) ′ ′ ′ m v v dp dp dp dp −     dt  m • ( h′′ − h ′)  dV V dm  dm n  •  = − − ∑ m pi h pi − moi h  − ∑ α Q Ai (Tst ,i - T )  +h v ′′ − v ′  dt m dt  dt  i i 


23

Zákony zachování pro konečné objemy 6. Zachování energie ve formě pro celkovou vnitřní energii nebo entalpii HK v zóně K

N

i

i

K

N

i

i

(

)

•

)

•

d (U + Emech ) = dQ − ∑ pi dVi + ∑ ui + emech ,i m i dt

(

d (H + Emech ) = dQ + ∑Vi dpi + ∑ hi + emech ,i m i dt hK = h + emech  w 2j   d H j + m j   N j  •  H .γ 2 H Kj .δ j ,i Ki j ,i  = + V j ,i . ∑ dt Vj i   Vi  + ∑Q + V j . •

k

k, j


  + α Q .Ai , j .(Ti - T j  

d pj dt

24

 ) + 

6. Zachování energie ve formě pro celkovou entalpii v zóně j a derivace termické stavové rovnice

Pro zónu se zanedbatelnou kinetickou energií po derivaci celkové entalpie v zóně

H = H (T , p , { s m }T ) dH dt N

=

j

∑ i

j

=

 • V j , i 

∂H ∂T

.

dT j

+

dt

 H Ki .γ .  Vi 

j ,i

+

∂H ∂p

.

H Kj .δ Vj

dp

s

j

dt j ,i


∑ i

∂mi

.

dm i , j dt

  + α Q .A i , j .(T i - T j   d pj

+ ∑ Q& k, j + V j . k

+

∂H

dt

25

=  ) + 

6. Zachování energie ve formě pro celkovou entalpii a derivace termodynamické i termické stavové rovnice ∂H     ∂H s s  dp j  ∂H dmi , j dT  ∂T ∂H ∂p dmi , j  . .∑ . + ∑ − − V j . + =  p p ∂ ∂ m dt dt m dt p ∂ ∂ ∂ i i i   i    ∂T    ∂T Nj  •   H Ki γ j ,i H Kj δ j ,i    ( ) = ∑ V j , i + + α Q Ai , j Ti - T j  +   V V i   i j    ∂H ∂p dV j • + ∑ Q + ∂T . . k, j ∂ p ∂ V dt k ∂T

Pro ideální plyn:

∂H T ∂ {h } T = {m } = {m } {c p } ∂T ∂T  {m }T {r }  ∂  T  T V { ∂p m } {r }   = = V ∂T ∂T ∂H = hi ∂m i ∂p ∂m i


26

 {m }T {r } ∂  T V =  ∂m i

    = ri T V

6. Zachování energie ve formě pro celkovou entalpii a derivace termodynamické i termické stavové rovnice pro vodní páru H = m h(T , p) ;

dH dm ∂h dT ∂h dp =m . + m . + h. dt ∂T dt ∂p dt dt

• dh n  • dp dm  m m = ∑ m pi h pi − moi h + ∑α Q Ai (Tst ,i - T ) + V −h dt dt dt  i i 

1 dV V dm ∂f dp − − Nj •  dp ∂h m dt m 2 dt ∂p dt  ∂h dm  •  m . +  m − V  = ∑ m pi h pi − moi h + ∑α Q .Ai .(Tst ,i - T) − h ∂f ∂T dt  i i   ∂p dt ∂T ∂h ∂f     m   Nj   • •  m ∂h − V − ∂T ∂p  dp = m h − m h + α .A .(T - T) −  h − ∂h 1 V  dm − ∂h 1 dV ∑ Q i st,i oi   pi pi  ∂p dt ∑ ∂f ∂T ∂f m  dt ∂T ∂f dt   i i      ∂T  ∂T  ∂T   π=

Tref p ; τ = p ref T −1

 ∂p  RT ∂p p ref τγ πτ − γ π = = γ  ; 2 ππ ∂v  p ref ∂ T T γ ππ  ∂h ∂ (RTref γ τ ) = RTref γ πτ ; ∂h = ∂ (RTref γ τ ) = − RTref γ ττ Tref2 = − Rτ 2γ ττ = ∂p ∂p p ref ∂T ∂T T


27

7. Rovnice pro společný tlak v 0-D objemu s jednou hlavní exotermickou reakcí (hoření paliva) Zadání hoření buď ve formě úbytku hlavní reakční proměnné (paliva) rozhodující reakce hoření sig (vstup paliva do reakce, možná dodatečná korekce na nedokonalé hoření se spalitelnými Tref produkty) nebo přímo vývinu tepla Qcorr (pak ∆hCH =0 ) V dp  = − κ-1 d t   − 

{s h } − κ T {s r }T  κ-1  T

{s h } − κ T κ-1 T

j {s r }

T

Nj •  dmsig • T C + + + Q V i  {s hK }i ∑ corr s*1 d t i  

{s m}i Vi

N { m }  κ dV   j •  {s m }i − p γ i + s δ i  + − ∑V i  V V 1 κ − d t  i  i  

{ c } { m}

κ =

p

s

{ c } { m} − { r } { m} T

s

p

T

s

s

;

T

{s h }T −

T

s

γ i + {s h }

{s h } =

{

s

κ T κ-1

{s m} V

   α + − ) A (T T δi  ∑ i i i +  i 

T  j {s r }  s*r-1 C

d

{r −1 m}CH



T ref ∆hCH + c p (T − T ref

)}

s

w uj 2 + w vj 2 + w wj 2 3  2 {s hK }= {s h }+ (K + k ){s 1} == {s h } +  + w ′ j {s 1} 2 2  


28

d t

8. Zachování energie pro ideální plyn ve formě pro tlak. Inverzní algoritmus. ODE řešená pro d p/d t , pokud jsou známy přeměny látek (simulace oběhu s daným přívodem tepla) nebo inverzně pro hlavní reakci, pokud znám průběh tlaku z měření. Možnosti dalších korekcí zejména u děleného spalovacího prostoru.

dm sig V d p  κ T  T +  {s h } − = T {s r }  s*1 C κ-1 dt  d t κ-1  Nj •    T {s m }i T {s m } = ∑ V i  {s h K }i γ i + {s h } δ i  + ∑ α i Ai (T i − T ) + Vi V i    i  • dV  − + Q corr −  p κ −1 d t 

κ

 − 

{s h }T

κ − T κ-1

{s h }T − T { } r j s

κ T κ-1

T  { } r  s*r -1 C j s



d

{r −1 m }CH d t

Nj •   {s m }i  {  s m} γi +  ∑V i  δ i  V   i  Vi 


29

+

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. Obecná forma systému Obyč. Dif. (E)rovnic Vektorový zápis pro neznámé složení o s látkách; úplná soustava obsahuje navíc hybnost a energii i derivace stavové rovnice s+3=y Obecná forma ODE soustavy po transformaci pomocí zákonů zachování do KO Matematické modelování oběhu motorů

 mO 2     m N 2  {s m} =   mCO2     M

 f 1    f 2  {s f } =    M    f n 

y* y A .

d { y f (t)} d t

d { y f (t)} = d t

=

y* y A

30

{ F [t, {

−1

y

y

{ [

f( t )}

T

]}

. y F t, {y f( t )}

T

]}

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH.

• OBEH je 0-D (zónový) model přeplňovaných i nepřeplňovaných vznětových i zážehových motorů s možností využití výkonové turbiny (turbocompound), regulace turbodmychadla a dvoustupňového přeplňování. • Počítá ustálený stav pro přímou i nepřímou úlohu o rovnováze turbodmychadla. • Iteruje některé další parametry oběhu.


31

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. Zákon zachování hmotnosti pro SacíPotr., VAlec, VýfukovéPotrubí • • dm x dm xCH = σ px m p − σ x m o + dt dt

a pro průtoky Kompresor, Sání, Výfuk, Turbinu Zákon zachování energie d ∑ mx .hx •

x

dt

•

= m p h0 p − mo h0 + ∑ x

dQxCH dQ dp − ∑ iCH + V dt dt dt i

dx dQxCH dQiCH = H um p Q ; = Siαi (T − Tsi ) dt dt dt Důsledky předpokladu 0-D a jedné zóny + vystačí se se zachováním hmotnosti a energie; - lze respektovat jen časové závislosti, ne tlakové vlny.


32

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. • Soustava se používá pro 3 seriově zapojené nádoby (SP, VA, VP), přičemž do výfukového potrubí se může superponovat přítok z více fázově posunutých stejných válců. • Pro řešitelnost 3*(s+1) rovnic o 3*(s+2) stavových proměnných (s složek, teplota, tlak ve 3 nádobách s předpokládaným homogenním stavem uvnitř) je zapotřebí dalších 3 rovnic – zde stavové rovnice pro ideální směs plynů ve 3 objemech, které kromě toho nutno derivovat pro vyjádření neznámé derivace tlaku. • Entalpie složek jsou funkcí pouze teploty. Změna objemu nádoby (válce) je funkcí pouze času (tedy úhlu otočení kliky a úhlové rychlosti, která též může být funkcí času).


33

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. • Rovnice dále obsahují neznámé proměnné uvnitř nádoby, mezi nádobami a na hranicích modelu (před první, za třetí nádobou, na stěnách nádob). • Uvnitř nádoby – neznámé toky tepla v důsledku sdílení (přestupu) do stěn lze vyjádřit jako funkci teploty plynu, teploty stěny („okrajová podmínka“), součinitele přestupu tepla (závisí po konečném dosazení do kriterální rovnice rychlostifunkci času, na teplotě a tlaku plynu), teplosměnného povrchu (nanejvýš funkce času); – neznámé rychlosti chemických reakcí, zejména výsledné rychlosti hoření paliva, spojené s uvolněním reakčního tepla – zde se aproximuje empirickou funkcí času, závisející na provozních parametrech motoru, případně na rychlosti tvoření nehomogenní směsi). Matematické modelování oběhu motorů

34

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. • Rovnice dále obsahují neznámé proměnné uvnitř nádoby, mezi nádobami a na hranicích modelu (před první, za třetí nádobou, na stěnách nádob). • Mezi nádobami – neznámé hmotnostní podíly a entalpie na přítoku do nádoby jsou dány stavem v předchozí nádobě proti proudu (upwind) nebo „okrajovou“ podmínkou na hranici modelu (např. atmosférou); – neznámé hmotnostní toky mezi nádobami se určí z (klidového) tlaku a teploty proti proudu ve spojení nádob, z tlaku po proudu (statického, lze však při předpokladu nulové rychlosti v nádobě klidové a statické veličiny splývají) a z efektivního průřezu spojovacího prvku, který obsahuje empirické parametry ztrát a kontrakce proudu a může být funkcí času, případně zmíněných stavů. Matematické modelování oběhu motorů

35

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. • Rovnice dále obsahují neznámé proměnné uvnitř nádoby, mezi nádobami a na hranicích modelu (před první, za třetí nádobou, na stěnách nádob). • Na hranicích modelu – atmosférický tlak a teplota v sání nepřeplňovaného motoru a na výstupu z výfuku; – teploty stěn závisejí na středním toku tepla a tepelném odporu stěny zřejmě nutná iterace; – tlak a teplota na výstupu z kompresoru závisejí na příkonu kompresoru, daném tlakovým poměrem, účinností, otáčkami a průtokem vzduchu; teplota na výstupu z chladiče plnicího vzduchu závisí kromě prostupního součinitele na teplotě za kompresorem – zřejmě nutná iterace; průtok kompresorem je pak dán jako střední nebo okamžitá hodnota tlakem v sacím potrubí; – příkon kompresoru u turbodmychadla je zajištěn výkonem turbiny (průtok, entalpický spád, účinnost), který lze spočítat jako okamžitou hodnotu, avšak pro předpoklad ustálených otáček turbodmychadla nutno integrovat a stanovit střední hodnotu po výpočtu celého oběhu; – průtok turbinoyu závisí na jejích otáčkách, avšak i ty lze stanovit až z výkonové rovnováhy turbodmychadla – další důvod pro iteraci.


36

9. Obecná formulace soustavy obyčejných DR pro 0-D model – příklad OBEH. Pro numerické řešení je nutno vzít v úvahu • • • •

•

vektorovou formulaci soustavy ODE dle úvodu kap. 9 nutno po dosazení z rovnice stavu vyřešit vůči derivacím; řešit numericky zdánlivě počáteční problém při použití iterace (přiblížení k neznámým okrajovým veličinám vypočteným v předchozím oběhu); pak jsou použitelné např. explicitní metody s volbou kroku typu Runge-Kutty; zdánlivě počáteční problém, ve skutečnosti - při hledání ustáleného stavu motoru – je to okrajová úloha pro nalezení periodického řešení. Je řešitelná opakováním výpočtu, pokud prostá iterace konverguje. Tedy: nutnost iterace mezi oběhy. Ukončení výpočtu po ustálení řešení. Jak posoudit? Přesné řešení neznáme, při každém opakování vznikne určitá odchylka, která nemá u všech proměnných monortónní průběh. Proto je navíc žádoucí kontrolovat přesnost splnění zákonů zachování ve středních hodnotách za celý oběh („rezidua“). Matematické modelování oběhu motorů

37

10. Numerické metody Obecné vlastnosti Z názoru i z integrální formulace MKO pro nezávislých proměnných (a ze zkušenosti s 0-D modely): Musí se brát zřetel především na časové derivace - metody, pro něž se v principu řeší počáteční metoda • explicitní (např. R.- K); • implicitní, zejména v úpravě pro stiff systémy; • časový krok shora omezen přesností předpovědi z vypočtené pravé strany (derivace) - potíže u stiff rovnic; • ... a zdola zaokrouhlovacími chybami. Pro hledání ustálených řešení (např. periodických) přechází počáteční úloha v okrajovou, ovšem řešenou iteračně postupným zpřesňováním počátečních podmínek („metoda střelby“). Pro PDE (viz dále) jde o Cauchyho úlohu s elegantním fyzikálním významem pro hyperbolické (zobecněné vlnové) PDE (viz dále).


38

10. Numerické metody Obecné vlastnosti Časové derivace závisejí na úhlu otočení kliky a úhlové rychlosti motoru. S ohledem na význam fázování vůči klikovému hřídeli např. derivace zdvihu, průtokové křivky ventilů - i vůči času (všechny průtoky p-t, tj. řízené pouze rozdílem tlaků) jsou problematické výpočty při nízkých otáčkách. Proto počítat charakteristiky vždy od vysokých otáček a využívat již zkonvergované počáteční podmínky, přitom zmenšovat úhlový krok zhruba úměrně otáčkám. Při nízkém zatížení se vyrovnávají tlaky na vstupu a výstupu nádob. Potíže s derivací závislosti průtoku na rozdílu tlaků, tj přibližně • mi = sign(pi − p )Aeff 2 ρ ( pi − p ) Všechny průtoky klesají, mezioběhová iterace konverguje pomalu. Proto potíže také při nízkém zatížení. Matematické modelování oběhu motorů

39

10. Numerické metody Vlastnosti pravých stran Příklad stiffness: Potíže s neomezeností derivace pravé strany obsahující průtok závisející na rozdílu tlaků. Pro příklad přibližně z Bernoulliho rovnice • Rovnice pro průtok do a z mi = sign(pi − p )Aeff 2 ρ ( pi − p ) nádoby spojené se zásobníky s neomezeným tlakem, ideální plyn, malé rychlosti, const. cv: • • dm = m1 + m2 dt dT dm cvm + c vT = dt dt •

= c p m1 (T1γ 1 + Tδ 1 ) + •

+ c p m2 (T2γ 2 + Tδ 2 ) Matematické modelování oběhu motorů

40

10. Numerické metody Závislost stability na kroku a na nespojitostech derivace pravých stran


41

10. Numerické metody Závislost stability na kroku


42

10. Numerické metody Obecné vlastnosti Bilance: Soustava ODE by měla i při numerické integraci zachovávat bilancované veličiny (hmotnost, energie, ...). Její zápis a numerické řešení pro měrné veličiny násobené hmotností (a v součinu derivované) obecně konzervativní být nemusí (není zaručeně v případě Parciálních Diferenciálních rovnic PDE). Kromě toho nemusí být konzervativní přibližné ošetření výjimečných stavů (např. zpětných průtoků s vrstvenou směsí plynů v potrubích), respektování středních izoentropických exponentů při reálných expanzích - průtok ventilem, turbinou atp. Proto se počítají bilance hmoty a energie a vyhodnocuje jejich nedodržení. Chyby do 1% u hmotnosti a do 2% u energie jsou přijatelné.


43

11. Dosazení do pravých stran ODR Průběh hoření • Zákon zachování hmotnosti a energie: – Chemické děje - změna hmotnosti látek a vývin tepla • Určení počátku spontánního hoření (průtah vznětu). • Vibeho funkce pro empirickou náhradu průběhu chemických reakcí. • Popis průběhu hoření a jeho zvláštnosti (dvou-

až třífázové hoření – kinetický a difusní plamen vznětového motoru, klepání zážehového motoru). Matematické modelování oběhu motorů

44

11. Dosazení do pravých stran ODR Průběh hoření Určení průtahu vznětu (také indukční doby pro mez klepání). Předpoklady výpočtu průtahu na základě množství zreagovaného paliva: o vznětu rozhoduje kritické množství připraveného paliva, stejné ve stacionárním i nestacionárním případě (poměr objemu a povrchu pro tepelnou bilanci zárodku plamene?); koncentrace paliva v heterogenní směsi je dána optimálním přebytkem vzduchu (součást konstanty vztahu pro průtah); palivo se během přípravy doplňuje z okolí, takže předplamennou reakcí v místě vznětu se koncentrace nemění; teplotní pole je pro zárodek plamene homogenní; závislost na teplotě se obvykle omezuje z důvodů měřitelnosti jen na exponenciální část.


45

11. Dosazení do pravých stran ODR Průběh hoření Určení průtahu vznětu (také indukční doby pro mez klepání). dmmain dc a cOb 2 e = V main = − KVcmain dt dt ∆mmain , krit ∆t p

a = K ′V 1− a cmain

  M pO2 O2 1− a  RT ′ = KV  λ Ot  

pOb 2 Tb

e

−

TA T

−

E RT

 mmain  M = − K ′ main  V  

 c 1 − a  O2 = K ′V  λ Ot 

     

a

 pO2   RT

a

 pOb − TA 2  e T =  Tb 

a

  b TA 1 − a p a + b − TA O2  pO2 e − T = K ′′ V e T a a+b  Tb λ T   T

T

A ∆mmain , krit λa a + b − (a + b ) TA a + b − (a + b ) T T pO2 e = K pT pO2 e ∆t p = K ′′V 1− a

dt p = ∆t p

⇒

∫ 0

TA a + b − (a + b ) T K pT pO2 e

dmmain , krit ⇒ ∆mmain , krit

dmmain , krit ∆mmain , krit = =1= ∆mmain , krit ∆mmain , krit

∆t p

∫ 0

b

dt p TA a + b − (a + b ) T K pT pO2 e


T

 − TA  e 

Předpoklad (další): a=1 - jinak by se musel sledovaný a proměnlivý objem zahrnout do integrace.

Výsledkem je integrální rovnice s neznámou mezí intergační proměnné (času), která se řeší současně s numerickou integrací ODE zkusmo, tj. postupným načítáním integrandu a porovnáváním výsledku s mezní hodnotou 1. Při prvním překročení se okamžik vznětu v rámci kroku zpětně interpoluje. 46

11. Dosazení do pravých stran ODR Průběh hoření Určení průtahu vznětu (také indukční doby pro mez klepání). Výsledek simulačního programu dějů ve válci a úplné chemie Časy do vznícení pro iso-oktan 1000

100

t [ms]

10

1 p= p= p= p=

0.1

20 bar 40 bar 20 bar, D-E 40 bar, D-E

0.01 0.6

0.8

1

1.2

1.4

1000/T [1/K]


1.6

1.8

Výsledek je možno shrnout do závislosti průtahu na teplotě, tlaku a pro homogenní směs i na koncentraci kyslíku. Pak dosadit tuto náhradu do integrálního vztahu pro průtah. Toto řešení je o několik řádů rychlejší než přímá integrace rovnic chemické kinetiky během simulace dějů v motoru.

47

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce pro normovaný průběh hoření xQ Předpokládá se, že palivo vstoupivší do reakce vyvine ihned všechno reakční teplo; palivo má hmotnost mpal nebo látkové množství Npal H uηH (t )m pal ,r (t ) m pal N pal QH = ≈ 1− 0 = 1− 0 xQ = 0 0 Předpokládá se rychlost H u m palηH , KH H u m palηH ,KH m pal N pal reakce úměrná rychlosti vzniku aktivních center řetězových reakcí, dané okamžitým množstvím paliva, které zbývá:

−

dt

! dN akt = = N akt . f akt (t ) = N akt .k .t m dt t

N pal N

Konec hoření KH se musí definovat (zde již pro 95% shořelého paliva), neboť reakce dospívá ke konci asymptoticky:

dN pal

0 pal

=

m pal m

0 pal

=

m 0pal − m pal ,r (t ) m

0 pal

−

=e

∫ f akt (t )dt t PH

t

∫ f akt (t )dt

−

xQ = 1 − e ln

K=

t PH

( ); 1 1− x KH m +1

∆t KH


= 1− e

− K . t m +1

xQ = 1 − e

; xKH = 1 − e

 1 −ln   1− x KH

  t −t PH    ∆t KH

48

   

− K . ∆t KH m +1

m +1

= 1− e

!

= 0.95

 t −t −3 PH  ∆t  KH

   

m +1

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce xQ= aIHR a jejímaxima časová derivace Rychlost hoření poloha rychlosti ROHR dxQ dt

( t = 3(m + 1)

 t −tPH  − t PH m − 3. ∆t KH e  m+1

)

   

∆t KH


m+1

1  m+ 1

 m ; t dxq max = t PH + ∆t KH    3(m + 1) 

49

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce - IHR a její časová derivace - ROHR VIBE FUNCTION 1.200

0.04500

0.04000

1.000

0.03500

IHR = x Q [1]

0.03000

Integral HR 1 Integral HR 2 Integral HR 3 ROHR 1 ROHR 2 ROHR 3

0.600 m=3

m = 0.4 0.400

m = 1.5

0.02500

0.02000

0.01500

0.200 0.01000 0.000 -20

0

20

40

60

1000.00500

80

-0.200

0.00000 Crank Angle [deg. CA]


50

ROHR = d xQ/d al [1/°]

0.800

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce Toto vyhovuje zhruba pro zážehové motory při 2,5>m>1,4 a při předpokladu nenulového „opoždění“ zážehu (mezi přeskokem jiskry a počátkem měřitelného vývinu tepla, např. pro xQ =0.01).

x 1% =

 t −t − 3. 1% PH  ∆t 1 − e  KH

   

m+1

1  m+ 1

1  1  = 0.01; t1% = t PH + ∆t KH  ln   3  1 − x1%  !

Lépe je kombinovat 2 - 3 Vibeho funkce (rychlé počáteční hoření u vznětových motorů nebo komůrkového zážehu, hlavní hoření SI nebo CI, dohořívání SI nebo CI). Dohořívání může začít zhruba v maximu rychlosti hoření, rychlé počáteční hoření se odhadne z množství paliva vstříknutého během průtahu vznětu nebo v komůrce SI. Matematické modelování oběhu motorů

51

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce – kombinace více průběhů

xQ =

Q p (α ) m p .H u .ηchem

m +1   α −α PH ,i  i   − 3.     ∆α s 0 − 95%,i   . xi ; = ∑ 1 − e   i     


∑ xi = 1 i

52

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce – přepočet parametrů Přepočet parametrů Vibeho zákona z referenčních (při známém, např. změřeném nebo odhadnutém stavu motoru) pro vznětové nebo zážehové motory (úhel zážehu nebo vznětu αPH , úhel hoření ∆αs a parametr(y) tvaru m) podle přebytku vzduchu, otáček motoru a stavu náplně válce na počátku hoření, případně i během něho (teplota, tlak, složení u zážehových motorů, průtah vznícení u vznětových motorů) atp. Např. pro vznětové motory lze úhel hoření odhadnout pro dané a klasické vstřikovací zařízení (ne CR) a pro daný tvar spalovacího prostoru jako ∆α s 0−95% ∆α s 0−95%

 λ = ∆α s ,ref . λ  ref  λ = ∆α s ,ref . λ  ref


   

−0.6

   

−0.6

 nM . n  M ,ref

0. 5

  pro nM ≥ nM ,lim  

pro nM < nM ,lim 53

11. Dosazení do pravých stran ODR Průběh hoření Vibeho funkce – přepočet parametrů Rostoucí otáčky motoru úhel hoření a průtah zážehu-vznětu prodlužují. •

U zážehových motorů má úhel hoření minimum v okolí stechiometrické směsi (dobré hodnoty kolem 40° kliky), u vznětových klesá s rostoucím přebytkem vzduchu.

•

Parametr tvaru pro zážehové motory m=1,5, pak však nutno zavést podle parametrů náplně válce “průtah” zážehu (ve skutečnosti směs hoří bezprostředně po přeskoku jiskry, ale bez viditelného vývinu tepla převyšujícího chyby měření). V OBEH se “průtah” počítá do shoření 1% paliva. Vždy existuje dohořívání cca 5-10% paliva s velkým ∆αs 0-95% a malou hodnotou m, často se však zanedbává (ovlivňuje však teplotu výfukových plynů a tím zpracovatelný výkon pro turbinu). Spaliny recirkulované z předešlého oběhu, nízký tlak a teplota hoření dále zpomalují.


54

Průběh hoření - přepočet parametrů • U vznětových motorů je zapotřebí počítat průtah vznícení a množství paliva vstříknutého během něho. Na něm závisí parametr m, lépe však počítat s kombinací dvou Vibeho zákonů, z nichž první respektuje spalování připravené směsi (x1 odpovídá poměrnému množství připraveného paliva) při nízkém m1 cca 0,1 a ∆αs 0-95%, 1 cca 10°. Zbytek na difusní hoření s m2 cca 1,5. Úhel difusního hoření klesá s rostoucím přebytkem vzduchu. • Přepočet úhlu hoření u vznětových motorů v závislosti na otáčkách se použije jen pro otáčky větší než zadané referenční, pro menší jen vliv přebytku vzduchu (odpovídá zhruba závislostem na turbulenci). • Obecně má na parametry oběhu primární vliv poloha počátku hoření a jeho úhel (délka), parametr tvaru má vliv menší.


55

Průběh hoření - přepočet parametrů • • • • • • • • • • •

U vznětových motorů průtah PRUTW(P,T)=1.668E-4*RNM*(P[MPa])**(-1.19)*EXP(4650./T) !...WOLFER !PRUTS(P,T)=2.381E-5*RNM*(P[MPa])**(.386)* EXP (4644/T) !Shipinski-Myers-Uyehara !PRUTI(PARC,T,V,ALKH)=6.E-3*RNM*(del tau fyz+ (.02651*(P[MPa])**(-.7)+.07345*(P)**(-1.8))*EXP(3930/T)) !Sitkei

U zážehových motorů FUN(RNM)=1.+400./RNM+8.E5/RNM/RNM GUN(RNM)=1.33-660./RNM FUL(RLA)=2.2*RLA*RLA-3.74*RLA+2.54 GUL(RLA)=2.*RLA*RLA-3.4*RLA+2.4 FUL87(RLA)=2.2*RLA*RLA-3.74*RLA+3.7 Matematické modelování oběhu motorů

56

Průběh hoření - přepočet parametrů •

U zážehových motorů GUNO=GUN(RNM)/GUN(RNM0H) FUNO=FUN(RNM)/FUN(RNM0H) IF(RNM.LT.800.) FUNO=1. IF(RNM.LT.800.) GUNO=1. CHIPOM=CHIRES/CHI0 IF(CHIPOM.LT.0.5) THEN CHIPOM=0.5 ELSE IF(CHIPOM.GT.2.) THEN CHIPOM=2. END IF DALZ=DALPP*(FUNO*(70.-AL1Z)/(70.-AL1Z0)*(2.16*TVA0/TVA-1.16)*(PVA/PVA0)**(-.17) *(.088*CHIPOM+.912)*FUL(RLA)/FUL(RLAM0)) DAS=DAS0*GUL(RLA)/GUL(RLAM0)*(1.33*TVA0/TVA-.33)*(PVA/PVA0) *(-.28) *(.237*CHIPOM+.763)*GUNO ALPH=AL1Z+DALZ RM=RM0 AVIBE=A0*(DAH/DAS)**(RM+1.) ALPH=ALPH-DAH*KDAL*(.0100503/AVIBE)**(1./(RM+1.))


57

Průběh hoření - přepočet parametrů • • • • • • • • • • • • • • • • •

30015 PVAPOM=PVA/PVA0 C PREPOCET PARAMETRU ZAKONA HORENI PODLE TEORIE LACINA VU 87 C PLATI PRO PH DEFINOVANY OD SHORENI 10% PALIVA ! FUNO= 0.45*RNM/RNM0H+0.55 GUNO=-0.26*RNM/RNM0H+1.26 FUAL1Z=0.45*AL1POM*AL1POM-0.6*AL1POM+1.15 GUAL1Z=(3.36*AL1POM-2.36)**(-0.13) FUPVA=0.31*(1./PVAPOM-0.99)**(-0.25)+0.9 IF(PVA.GT.PVA0)FUPVA=1.6/PVAPOM-0.6 ALPH=AL1Z+DALP0*FUNO*FUAL1Z*(1.16*TVA0/TVA-0.16)*PVAPOM 1**(-0.08)*(0.05*CHI/CHI0+0.95)*FUL87(RLA)/FUL87(RLAM0) DAS=DAS0*GUL(RLA)/GUL(RLAM0)*(1.33*TVA0/TVA-0.33)*FUPVA* 1(0.237*CHI/CHI0+0.763)*GUNO*GUAL1Z ALPH=ALPH-DAS*0.03512**(1./(RM+1.)) 'lambda mimo meze 1.1-1.3''T VA mimo meze 510-610K pro predpoved zak. hor. !' 'p VA mimo meze 0.3-0.6MPa pro predpoved zak. hor. !''CHI mimo meze .029-.053 pro predpoved zak. hor. !''n M mimo meze 600-1500/min pro predpoved zak. hor. 'al 1Z mimo meze -46;-27st. pro predpoved zak. hor. !'


58

Určení meze klepání zážehových motorů vysokého měrného výkonu Meze klepání, měrná spotřeba paliva a teplotní pole hlavy pro různé podmínky chlazení


59

Dosazení do pravých stran DR: Přestup tepla do stěn • Zákon zachování energie - odvod tepla do stěn – model stěn s místní teplotou prakticky nezávislou na čase (během jednoho oběhu) – rozdělení na píst - hlava - odkrytá část válce – dno hlavy s ventily a píst reprezentován konstantní plošně střední teplotou, u válce je střední teplota odkryté části funkcí úhlu kliky – teplosměnné průřezy - u válce funkce úhlu kliky, jinak konstantní – nutno zavést fyzikálně zdůvodněné střední hodnoty teploty plynu a součinitele přestupu tepla Matematické modelování oběhu motorů

60

Dosazení do pravých stran ODR Přestup tepla do stěn- teploty stěn • Zákon zachování energie - odvod tepla do stěn – model s tepelnými odpory – pro energetickou bilanci se odvod tepla charakterizuje KCH: • •

Q CH = K CH Q p τ per •

Q CH =

τ per

∫ Sα (T − T

s

∫ Sα .T dτ

)dτ

0

;

τ per

Ts =

0

•

Q CH .τ per

− τ per

τ per

∫ Sα dτ

∫ Sα dτ

0

0

τ per •

Q CH =

S

ρs

∫ Sα .(T − T

s

(Ts − TCH ) =

0


τ per

)dτ

(

= S .α . T − Ts 61

)

Dosazení do pravých stran ODR Teploty stěn - tepelné odpory

Tg

Sg

αg

T1

•

Q = S g α g (Tg − Ti ) =

ρ=

1

α

kS ref =

T2

λ S

nebo S ref

ρΣ

=

Sλ

δ

(T1 − T2 ) = S Σ k (Tg − TW )

δ λ S ref S gα g

=

+


αW SW

S ref S ref δ Sλ

1 = ∑ Ri i

1

ρi

∑S i

+

S ref SW α W 62

i

TW

Dosazení do pravých stran ODR Přestup tepla do stěn- teploty stěn Model vedení tepla a jeho přestupu na chlazené straně podrobnější náhrada vestavěná též do OBEH Odpory=1/přestup nebo vedení tepla na příslušné teplosměnné ploše. Možno doplnit kondenzátory pro nestacionární proces. αg

Tg

Tol

αo

Tw


63

Zjednodušené zadání chlazení válce z teploty horního okraje • •

Menis rozdeleni teploty podel vlozky Chces se vratit v teto casti vstupu

•

K o n e c

? (y/n/end/back) ? (y/n/end/back)

3.casti vstupu Teplotní funkce vložky

1

• • • •

Menis charakteristiku turbiny 0.9 Menis zavislosti na pomernem tlaku Menis zavislosti na rychlostnim pomeru 0.8 Prubezne otacky pri KSI 1.479596

? (y/n/end/back) ? (y/n/end/back) ? (y/n/end/back)

AAVLO AVL2

0.7

•

Chces se vratit v teto casti vstupu

•

K o n e c

? (y/n/end/back)

0.6

4.casti vstupu

0.5 0.4 0.3 0


50

100

64

150

200

Přestup tepla

Válec o vrtání D, zdvihu Z, zdvihovém objemu Vz; okamžitý tlak p a (prostorově střední) teplota T, střední rychlost pístu cs při otáčkách nM, obvodová rychlost víru ve válci (na 0,7 D) cobv. Stav na počátku komprese x značen PZ, tlak ve válci protáčeného motoru pbez spalování. Eichelbergův vztah (cca 1930) pro čtyřdobé i dvoudobé naftové motory (vhodný stále pro přeplňované naftové motory)

α = 2 ,1. 3 c s . p.T [ kcal.m −2 .h −1 .K −1 ; m.s −1 ; kp.cm −2 ; K ] =

2 ,1.1,163 9 ,80665.10 4

.3 c s . p.T [ W .m −2 .K −1 ; m.s −1 ; Pa ; K

nM .Z 30 = 2 ,51.10 − 3.3 nM .Z p.T [ W .m −2 .K −1 ; min −1 ; m ; Pa ; K ] = 7 ,799.10 − 3.......

; pro


cs =

65

Přestup tepla Woschni (1967) z analogie s turbulentně protékaným kanálem Nu = C . Re 0.8 . Pr 0.33 po dosazení empirických mocninných vztahů pro látkové vlastnosti a po opravě rychlosti směrodatné pro konvekci na ”zvýšenou turbulenci během spalování” (ve skutečnosti je v této korekci zahrnut i vliv radiace tepla ze svítivého plamene):  c α = 110. D −0.2 .p 0.8 . T −0.53 .  2.28 + 0 ,308. obv cs 

  V .T  . c s + 3,24.10 − 3. z PZ . p − pbez spalování  VPZ .pPZ  

(

)

0.8

[ kcal.m− 2 .h −1 .K −1 ; m ; kp.cm −2 ; K ; m.s −1 ; dm 3 ; kp.cm − 2 ; K ] cobv 110.1,163 −0.2 0.8 −0.53   = + D p T . . . . 2 . 28 0 , 308 .  0.8  cs  9 ,80665.10 4

(

)

 V .T  . c s + 3,24.10 − 3. z PZ . p − pbez spalování VPZ .pPZ 

(

[ W .m −2 .K −1 ; m; Pa ; K ; m.s −1 ; dm 3 ; Pa ; K ] nM .Z 30 = 2 ,51.10 − 3.3 nM .Z p.T [ W .m −2 .K −1 ; min −1 ; m ; Pa ; K ] = 1.299.10 −2.......

; pro

cs =


66

  

)

0.8

Přestup tepla Woschni a Zapf – výměna náplně α = 110. D−0.2 .p 0.8 . T −0.53 . (6,18. c s )0.8 [kcal.m −2 .h −1 .K −1 ; m ; kp.cm −2 ; K ; m.s −1 ] =

110.1,163

(9,80665.10 )

4 0.8

.D −0.2 .p 0.8 . T −0.53 . (6 ,18. c s )0.8 [W .m −2 .K −1 ; m; Pa ; K ; m.s −1 ; dm 3 ; Pa ; K ] ;

nM .Z : 30 = 3 ,58.10 − 3.D −0.2 .p 0.8 . T −0.53 . (nM .Z )0.8 [W .m −2 .K −1 ; m ; Pa ; K ; min −1 ]

pro c s =

Pro přídavné rozvíření použita konstanta 3,7.10-3 v posledním vztahu. Pro motory spalující připravenou směs oba vzorce dávají nižší hodnoty než odpovídá měření. Snad jde o vliv izolující vrstvy porézního grafitu na teplosměnném povrchu motoru s difusním plamenem. Existují další použitelné vztahy (Bargende- 2 zónový SI, Morel-Jennings, Annand, Sitkei, Pflaum - komůrkové motory,...) Matematické modelování oběhu motorů

67

Přestup tepla na straně plynu - Woschni Nepřeplňovaný zážehový motor - přestup tepla 2000 min-1 2500

1800

Střední teplota ve válci [K]

2000

1400 1200

1500 1000 800 1000 600 400

500

Teplota Součinitel přestupu tepla

200

0 -150

-100

-50

0 0

50

100

Úhel kliky [st.]


68

150

Součinitel přestupu tepla Woschni [W.m-2.K-1]

1600

Simulace a analýza hoření pomocí zónového modelu válce zážehového motoru

1.

Úvod

2. Základní r ovnice zónových modelů

3. P roblém ř ešení soustavy O DR pro zónové modely

4. M odel hoření vazebné vztahy a geometrie

5. P říklady výsledků

6. Závěry

Geometrické údaje (plochy, objemy) zón v závislosti na poloze geometricky definovaného (např. hemisférického) plamene a poloze pístu Automaticky generovaná geometrická data pomocí CAD programu Úhel kliky Plocha fronty plamene

(poloha pístu

Vzdálenost od svíčky

Spalovací prostor 4 ventilové hlavy


69

Simulace a analýza hoření pomocí zónového modelu válce zážehového motoru Úvod

2500 ABV ABP

2000

ostatní díly

ABF ABH

1500 1000 500

úhel natočení klikového hřídele a [o]

Úhel kliky [°]


70

160.0

140.0

120.0

100.0

80.0

60.0

40.0

20.0

0.0

0 -20.0

6. Závěry

-40.0

3000

výfukový ventil

-60.0


-80.0

-100.0


-120.0

-140.0


-160.0

-180.0


[W.m .K ]

Součinitele přestupu tepla

1.

součinitel pro přestup tepla do vložky válce, pístu, výfukového ventilu a hlavy válců -2 -1 [W.K-1.m-2]

Součinitele přestupu tepla podle Bargendeho pro různé části pracovního prostoru

Energetická bilance Ve výstupech programu OBEH (porovnatelná s měřením, tč. kalibrace pro naftové motory John Deere). Závislost na přesně provedených měřeních v širším rozsahu otáček a zatížení vždy existuje, pokud mají výsledky mít relativní chybu lepší než 5-10% (podle velikosti tepelného či entalpického toku). ENERGETICKA BILANCE QP PH = 0.493612E+05 Q ETAH=-0.176030E+04 QP PL = 0.593587E+03 QP PAL= 0.000000E+00 QO PI = 0.191645E+05 QO PHV= 0.103009E+05 QO VP = 0.251782E+04 QO TVS= 0.770221E+03 QO TVY= 0.000000E+00 QO OLT= 0.000000E+00 QO VYF= 0.185508E+05 QO CHL=-0.230152E-04

DELTAQ=

0.82 %

Pomerne hodnoty QP PH = 100.00 % QO PI = 38.82 % QO TVY= 0.00 % QO OLT= 0.00 %

DELTAM=

0.10 %

Q ETAH= QO PHV=

-3.57 % 20.87 %

QP PL = QO VP =

QO VYF=

37.58 %

QO CHL=


1.20 % 5.10 % 0.00 %

QP PAL= QO TVS= Q ZBYT=

71

0.00 % 1.56 % -3.94 %

Využití modelu motoru pro stanovení okrajových podmínek; mechanické ztráty motoru Teplotní namáhání: prostup tepla stěnami

Vyvinuté prostředky umožňují • Přímý přenos středních hodnot okrajových podmínek 3. druhu do MKP výpočtů teplotně zatížených dílů. • Inverzní algoritmus pro výpočet sumárních teplotních odporů dílů pro zadání do modelů oběhu (OBEH, GT Power) pomocí výpočtu středních povrchových teplot v rozsahu jednotlivých konečných odporů. • Řešení interakce pohyblivého pístu a válce z hlediska přestupu tepla z plynu, z pístu a vývinu třecího tepla na kroužcích a plášti pístu. • Nestacionární nízkofrekvenční děje při přechodových režimech (síť odporů doplněná o kondenzátory a převedená tím ze soustavy lineárních rovnic na soustavu ODE). Lze použít i pro potrubní systémy s katalyzátory. Matematické modelování oběhu motorů

72

Dosazení do pravých stran ODR Průtoky škrcením mezi objemy • Zákon zachování hmotnosti: – Průtok ventily a turbinou


73

Průtok pro proudění v subsonické a transsonické oblasti.

Průtok škrcením mezi moduly modelu m& = ρ 2 .A2 .w2 max ; w2 max =

〈wa

2

;

2

ρ2=

 p2    r.T  2

〈 ρ 01 .

ε∗ κ −1  ∗ κ  1 − η.1 − ε   

ι κ −1    p2  κ   κ − 1  κ −1  p 2  1 ∗   < w2 = 2.c p .(T01 − T2 ) = 2.c p .T01 .η .1 −    ,η = 1 + ζ . ε = 1 − η.(κ + 1)  p p01  subsonic 01  N      


74

Průtok kanálem s ventilem - možnost kombinace globálního 0-D (1-D) a 3-D detailního modelu Průtok s odtržením a kontrakcí, výstupní ztráta kinetické energie: popis pomocí průtokového součinitele*účinná plocha pro rychlost stanovenou z vratného procesu plně vyhovuje pro Re>Rekrit;

Izočáry velikosti rychlosti při stacionárním proudění soustavou sací kanál-ventil do válce motoru.


možné problémy pro M>1 – avšak u výfukového ventilu vcelku vyhovuje. 75

Předpověď vírového čísla kanálu - dodatečné zadání pro kumulaci momentu hybnosti Schéma experimentu Aerodynamický experiment a výpočet v systému FLUENT pro zjištění průtokových a vírových vlastností kanálu - závislost na uspořádání sítě konečných objemů. Možnost předpovědi výsledků změn geometrie kanálu. POZOR – pro přímý výpočet i pro vyhodnocení experimentu musí být použit stejný model, méně (co do podrobností průtoku) bývá více!


76

Průtokové vlastnosti - ventilový rozvod Ventil zadán pomocí poměrného průtokového průřezu (geometrie+kontrakce+ztr áta kinetické energie) v závislosti na zdvihu;

Poměrný průtokový průřez sacího ventilu 0.9 0.8 0.7 0.6

mu*sigma

0.5

zdvihu na úhlu kliky.

0.4 0.3

Zdvih sacího ventilu

0.025

0.2 0.1

0.02

0.8 0.7 0.6

0

0.05

Průtokový součinitel sacího ventilu 0.1 0.15 0.2

0.25

Zdvih ventilu [m]

0

0.9

zdvih sest bok vule

poměrný zdvih h/d ref

0.3

0.35

0.4

0.015 FIS sest bok

0.01

0.5 0.4 0.3 0.2 0.1 0

0.005

0 0


50

100

150 úhel kliky [°KH]

77

200

250

300

Potrubní okolí motoru - tlakové ztráty • Ovlivnění řešení podle integrálních hodnot během jednotlivých oběhů • tlakové ztráty v potrubích s ustáleným průtokem – odhad z proudění nestlačitelné tekutiny, • má výhodu v možnosti sumace místních ztrát – ale jen pokud se vzájemně neovlivňují

1 1 ς m& 2 2 ∆p = ς . . ρ .w = . 2 2 2 Aref ρ stredni


78

Model turbodmychadla: průtok jako podmínka na konci výfukového potrubí Průtoková charakteristika turbiny při optimální účinnosti K 27 ..21

• 6.0 5.5

10.21

m red [kg.s-1.K^0.5.bar-1]

5.0

13.21

4.5

14.21

4.0

17.21

3.5

19.21 21.21

3.0

Power Turbine

2.5

Turbocharger

2.0

Turbocharger

1.5 1.0 0.5 0.0 1.0

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

2.0

2.1

pi T


79

2.2

2.3

2.4

2.5

•

Model turbiny: charakteristika průtoku a výkon okamžité hodnotyúčinnosti, pro určení průtoku, účinnostiturbiny a výkonu; pozor - redukované hodnoty definovány jinak než u kompresoru

•

•

mT

µT =

pT 1 ψ (π T ) rTT 1

Aref

mTred = f π T , nT 1 ψ (π T ) r

(

= Aref

TT 1

)

• •

mTred =

• mT TT 1 ; mT = µT Aref pT 1

pT 1 ψ (π T ) rTT 1 κ

κ +1

2 κ +1  − 2.κ  − κ 2  κ −1 κ + 1  κ −1 ψ (π T ) = . π T − π T κ  pro π T ≤  , jinak ψ * = κ .    κ − 1  2  κ +1   

cs = 2.∆iTs

nTred

1 −κ T   = 2.c pT .TT 1 1 − π T κ T   

n = TD TT 1

;

; x=

ηTs = f (π T , x ) ;


π .DT 2 .nTD 60.cs

=

π 60 2.c pT

•

PTs = mT .∆iTs

DT 2

1 − πT •

;

80

nTred 1 −κ T

κT

PT = mT .∆iTs .ηTs

;

Model turbodmychadla: průtok a účinnost turbiny Úplná charakteristika turbiny při optimální účinnosti •


81

Hltnostní součinitel turbin K 27 v různém provedení (různý referenční průřez) při optimální účinnosti

•


82

Model turbodmychadla: průtok jako podmínka na konci výfukového potrubí, výkon a účinnost do integrálních hodnot

1-D model turbiny Matematické modelování oběhu motorů

83

Podklady pro vstupy - charakteristika turbiny bez redukce na rychlostní poměr při optimální účinnosti

•


84

Vstupy - redukovaná charakteristika turbiny na rychlostní poměr při optimální účinnosti

•


85

Vstupy - redukovaná charakteristika turbiny na rychlostní poměr při optimální účinnosti - výhodnější vztažení hltnostního součinitele na redukované otáčky •


86

Model turbodmychadla: průtok jako podmínka na konci výfukového potrubí, výkon a účinnost • model turbodmychadla: turbina - průtok a účinnost: zadání do integrálních hodnot

Zadání

Zadání ηmax

xopt = f x (π T ) nTred ,opt = f n (π T ) = =

60.xopt

π .DT 2

Zadání xopt

nT ( xopt ) TT 1

=

1 −κ T   2c pT 1 − π T κ T   

 x = fη , x =opt. (π T ). gη   xopt 

   

 x µT = f µ , x =opt . (nTred ). g µ   xopt 

   

η sT

Zadání µnom Zadání


87

Zadání charakteristiky turbiny • •

• • • • • • 1.2 • 1 • 0.8 • 0.6

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

Menis charakteristiku turbiny 1.3 Menis zavislosti na pomernem tlaku Posledni hodnota 1.2 tlak.pomeru ma byt > 5 ! 1.1 Pocet= 7 1 Pom.tlak Max.ise.ucin. Opt.rych.pomer 1.00 0.8500 0.5700 0.9 1.50 0.8500 0.5700 0.8 2.00 0.8600 0.5700 0.7 2.50 0.8400 0.5700 1.3 3.00 0.8300 0.5700 0.6 1.2 3.50 0.8000 0.5700 0.5 10.00 0.7000 1.1 0.5700 1 1.5 2 2.5 XPIT 1 Pocet (0..BEZE ZMENY) = 0.9 Chces zadat soucinitel hltnosti Pocet= 7 0.8 0.5 Max.ise.ucin. 1 mu T,BOpt.rych.pomer 1.5 Pom.tlak eta T,B 0.7

? (y/n/end/back) ? (y/n/end/back)

etaTmax Jmen.reakce xopt r 0.0000 j 0.0000 mu T 0.0000 0.0000 0.0000 0.0000 0.0000 3 3.5

Souc.hltnosti 1.0200 1.0900 1.1450 1.1710 1.1770 1.1800 1.1800 4 etaTmax xopt rj ? (y/n/end/back) mu T

2 Prut.souc.

0.6 Menis zavislosti na rychlostnim pomeru Udej posl.hodnotu vztazneho 0.5 rychl.pomeru > 1E4 !

•Matematické

1

modelování XKSI oběhu motorů

1.5

2

? (y/n/end/back) 2.5 88 XPIT

3

3.5

4

Model turbiny: výkon a účinnost do integrálních hodnot • turbina - střední hodnoty výkonu pro srovnání s příkonem kompresoru a určení plnicího tlaku

cs = 2.∆iTs

1 −κ T   = 2.c pT .TT 1 1 − π T κ T    720

•

PTs = mT .

cs2 2

;

PT =

cs2 2

∫ Pdα

.ηTs

;

P =

; x=

πDT 2nTD

720 •

∫ mT ∆iTsηTs dα

ηT =

0 720 •

∫ mT ∆iTs dα

0

=

PT

___

PTs


___

60 2

PTs •

m 89

0

720

Model turbodmychadla: kompresor • Ovlivnění řešení podle integrálních hodnot během jednotlivých oběhů – model turbodmychadla: kompresor - účinnost - odečte se přednostně ručně, součinitel skluzu pro výpočet otáček se vyhodnotí pomocným programem Č Z Strakon ice a.s. T u rbo D ivision

ČZ Strakon ice a.s. T u rbo D ivision

3.5

TK 1 pKref =m . . & LK pK 1 TKref 4.0

m & Kred

ΠD

C12

3.0

Standard conditions:

ΠD

p 0 = 981 mbar

3.5

p 0 = 981 /T mbar n Tdre d =n Td *(sqrtT 0 1) n

=n *(sqrtT /T 1 )

Tdre d Td 0 m Dre d = m D *sqrt(T 1 /T 0 )*p 0 /p

nKred = nTD

160 000

T =-1 293 K

0 n Tdre d = [min ]

n Tdre d = 180 000

2.5

K36

T 0 = 293 KStandard conditions:

TKref

Π D =p 2/p 1

n Tdre d = 90 000

n Tdre d = [min -1]

m Dre d =m D *sqrt(T 1/T 0 )*p 0 /p 1 ) Π D= p 2/p 1

3.0

TK 1

80 000 2.5

140 000 2.0

70 000

2.0

120 000 70 %

65 %

78 %

79 %

60 000

75 % 70 %

1.5

100 000

60 %

80 000

60 %

1.5

50 %

40 000 30 000

60 000 1.0 0.00

65 %

50 000

0.05

0.10

1.0

0.15

0


m Dre d

0.1

0.20

0.2

0.3

90

0.4

0.5

m Dred

0.6

Model turbodmychadla: kompresor

Č Z Strakon ice a.s. T u rbo D ivision

4.0

ΠD

TK 1 pKref Standard conditions: . & integrálních LK . Kred = m Ovlivnění řešením&podle p = 981 mbar pK 1 TKref

•

3.5

0

hodnot během

T 0 = 293 K

jednotlivých oběhů n Tdre d =n Td *(sqrtT 0 /T 1) -1

n Tdre d = [min ]

nKred = nTD

TKref

n Tdre d = 90 000

– model turbodmychadla:TKkompresor - účinnost 1 odečte se přednostně ručně, součinitel skluzu pro výpočet otáček se vyhodnotí pomocným programem m Dre d =m D *sqrt(T 1 /T 0)*p 0 /p 1 )

Π D =p 2/p 1

3.0

80 000

2.5

70 000

K36

2.0

78 %

79 %

60 000

75 % 70 % 65 %

50 000

60 %

1.5

40 000 30 000

1.0 0

0.1

0.2


0.3

0.4

91

0.5

m D re d

0.6

Model turbodmychadla: kompresor • Ovlivnění řešení podle integrálních hodnot během jednotlivých oběhů – model turbodmychadla: kompresor - otáčky a stlačení jako integrální (střední) hodnoty

∆ iK =

PT .ηmTD •

κ K −1   κ   K  p  − 1 c pK .TK 1 . K 2  p  K 1   =

mK

∆iK = u2 .ct 2 =

ηKs

π .nTD .D2 D

ct 2

60 ct 2 = µ K .(u2 − wr 2 .tgβ 2 )


Jako 1. přiblížení pro β2=0 a TK1=Tref se podle charakteristiky kompresoru vyhodnotí pro určení otáček

µK =

92

∆i K  π .nTDred .D2 D    60  

2

Potrubní okolí motoru - chladič plnicího vzduchu a tlakové ztráty • Ovlivnění řešení podle integrálních hodnot během jednotlivých oběhů • chladič plnicího vzduchu ηCH =

TK 2 − TCHL 2 TK 2 − TCHW 1

• tlakové ztráty v potrubích s ustáleným průtokem

1 1 ς m& 2 2 ∆p = ς . . ρ .w = . 2 2 2 Aref ρ stredni


93


94


95


96


97

Model turbodmychadla: optimalizace • Přímá a nepřímá úloha – model daného turbodmychadla: výpočet otáček a stlačení při daném režimu motoru; hledání vhodného turbodmychadla iterací - potřebný průřez turbiny (pro požadovaný vyšší plnicí tlak snížit) a potřebný vnější průměr oběžného kola (měnit tak, aby střední rychlostní poměr byl blízký optimálnímu); – přímé hledání potřebného průřezu podle přebytku / nedostatku výkonu (1. přiblížení podle potřebné účinnosti kompresoru při daném průřezu turbiny např. metodou regula falsi); průměr turbiny nutno iterovat, jak zmíněno nahoře.

Příklady průběhu tlaku s různými výfukovými systémy a počty válců při srovnatelném plnicím tlaku u motoru ČKD 27.5B8 jsou uvedeny dále. Matematické modelování oběhu motorů

98

Mechanické ztráty - model Poloempirický model mechanických ztrát: • Kinetostatický výpočet excentrického klikového mechanismu s úplnými účinky všech dílů. • Stribeckovy křivky pro třecí dvojice. • Přímá návaznost na výpočet tlaku v OBEH. • Dynamický výpočet nerovnoměrnosti chodu tuhého mechanismu se zpětnou vazbou na setrvačné účinky. Matematické modelování oběhu motorů

99

Mechanické ztráty - model • Kalibrace podle experimentů s postupným odstraňováním zdrojů ztrát a se simulací termodynamiky procesu v protáčeném válci. • Nenáročný na geometrické vstupní údaje. • Předpokládá rozložení tlaku na kroužcích podle zobecněných termodynamických výpočtů. Matematické modelování oběhu motorů

100

Mechanické ztráty - Stribeckovy křivky

STRIBECK CURVES | 4 cylinder 3017 min-1 | 34.5 kW (brake) | 90 % COEFFICENT OF FRICTION [1]

0.200 0.150 0.100 0.050 0.000 -0.050 -0.100 -0.150 -0.200 -0.250

-0.200

-0.150

-0.100

-0.050

Sommerfeld Number (1 st ring) [1]


0.000

0.050

0.100

mu R 1 ring [1]

101

0.150

0.200

mu R 2 ring [1]

Mechanické ztráty - průběh tlaku na pístu a kroužcích

PRESSURES AT RINGS | ŠKODA 3EA111 5000 min-1 | 28.8 kW (brake) | 78 % 5.0000

4.5000

4.0000

PRESSURE [ MPa ]

3.5000

3.0000

2.5000

2.0000

1.5000

1.0000

0.5000

-310

0.0000 -210 -110 CRANK ANGLE [ deg. ]

-10

90


190

290 CYLINDER

390 1st ring

102

490 2 ring

590

Mechanické ztráty - kalibrace modelu na protáčeném motoru s korekcí na termodynamické ztráty ve válci MOTORING OF ŠKODA 781.135 ENGINE Engine Speed [min-1] 0

500

1000

1500

2000

2500

3000

3500

4000

4500

5000

10

Brake torque (simulated) [N.m] 5

Torques

[N.m]

Brake torque measured 1 [N.m] Brake torque measured 2 [N.m]

0

Luboil pump torque [N.m] -5

Valve gear torque w/t HL [N.m] Thermodynamic loss torque [N.m]

-10

Valve gear torque with HL [N.m] -15

Torque of cooling pump and alternator [N.m] -20


103

Mechanické ztráty - ověření modelu Motoring of ŠKODA 781.15 (5 Main Bearings) 0 0

1000

1500

2000

2500

3000

-5

Motoring Torque [N.m]

Beze změny nastavení prametrů modelu: výpočet a měření na motoru 1500 cm3 s pětkrát uloženou klikou

500

-10

-15

-20

Brake torque (simulated) [N.m] Brake torque measured 1 [N.m] Brake torque measured 2 [N.m]

-25

Engine Speed [min-1]


104

3500

4000

4500

Mechanické ztráty - model rozvodu OHV s hydraulickými vymezovači vůle a bez nich SROVNÁNÍ MĚŘENÍ A VÝPOČTU ŠKODA 781.135

Hnací moment na klikovém hřídeli [N.m]

7 6 5 4 3 Měření bez HZ

2

Výpočet 781.135 bez HZ Měření s HZ

1

Výpočet 781.135 s HZ 0 1000

1500

2000

2500

3000

3500

4000

Otáčky kliky [min-1]


105

4500

5000

Mechanické ztráty - síly tření na PK FRICTION FORCES AT PISTON | 4 cylinder 3017 min-1 | 34.5 kW (brake) | 90 % 150.000

FORCE [N]

100.000 50.000 0.000 -50.000 -100.000 -150.000 -200

-100

0

100

200

300

400

500

CRANK ANGLE [ deg. ] Friction: 1 ring

2 ring


3 ring

106

Scrapper ring

Mechanické ztráty - výkon ztracený třením

LOST POWER [W]

FRICTION LOSS POWER | 4 cylinder 3017 min-1 | 34.5 kW (brake) | 90 % 1000 900 800 700 600 500 400 300 200 100 0 -100 -200

pístní čep

PK

plášť pístu hlavní ložiska

-100

0

100

200

300

400

500

CRANK ANGLE [ deg. ] Rings

Pi skirt


Pi pin

Conrod/Crank

107

Crank main

Mechanické ztráty - síly na pístu FORCES AT PISTON AND CR SMALL END | 4 cylinder 3017 min-1 | 34.5 kW (brake) | 90 %

1500 1000 FORCE [N]

500 0 -500 -1000 -1500 -2000 -2500 -200

-100

0

100

200

300

400

500

CRANK ANGLE [ deg. ] Inertia P

Normal P


Friction P

Perpend. to CR

108

Normal tiltin

Mechanické ztráty - momenty na klikovém hřídeli pro 1V a všechny válce. Moment tření a moment na pohon ventilů (stupnice vlevo)

TORQUE [N.m]

TORQUES AT CRANKSHAFT | 4 cylinder 3017 min-1 | 34.5 kW (brake) | 90 % 450.0000 400.0000 350.0000 300.0000 250.0000 200.0000 150.0000 100.0000 50.0000 0.0000 -50.0000 -100.0000 -200

2.0000 1.0000 0.0000 -1.0000 -2.0000 -3.0000

-100

0

100

200

300

400

CRANK ANGLE [ deg. ]


109

-4.0000 Total 1 cyl. -5.0000 Total incl. valve gear -6.0000 Total w/t friction engine 500Total 600 Friction Torque Valve gear torque at cra

Mechanické ztráty - rozdělení v klikovém mechanismu ŠKODA 3EA111 Indicator dia. (VW EA111 - PROMO calculation) Friction acc. to CR_SK4M*.XLS, valve gear SK_R_VAS,V (1.92 N.m for the whole 3 cyl. eng. at crank - with hydraulic followers) and accessories 2.5 N.m at 5000 min-1 Rings Piston Skirt Pist.Pin

Main Bearing CR Bearing 9% 9%

CR Bearing Main Bearing

Rings 50%

Pist.Pin 17% Piston Skirt 15%


110

Mechanické ztráty - rozložení ztrát v ústrojích motoru

OVERALL MECHANICAL LOSSES 3000 min-1 | 25.6 kW (brake) | 90 %

Crank Mechanism Engine Cycle Valve Gear Injection Pump Luboil Pump

0%

7%

8%

0% 1%

Water Pump Cooling Fan Vent. Loss Losses=const.

21% 63% 0%


111

Využití modelu motoru pro stanovení okrajových podmínek; mechanické ztráty motoru Mechanické ztráty Shares of Loss Sources at Maximum Torque Characteristics (1 dm3 engine) 100%

Valve Gear [N.m] 1 cyl. Luboil Pump [N.m] 1 cyl.

Share of Torque Source

80%

Water Pump [N.m] 1 cyl. Piston Rings [N.m] 1 cyl.

60%

Piston Skirt [N.m] 1 cyl. Piston Pin [N.m] 1 cyl.

40%

Conrod Crank Bearing [N.m cyl. Main Bearings [N.m] 1 cyl.

20%

0% 901

999

2000

3000

Matematické modelování oběhu motorů[min-1] Engine Speed

4004

5001

112

Mechanika motoru - moment a setrvačné síly INERTIAL EFFECTS AND THEIR BALANCING 4 cylinder 3017 min-1 | 34.5 kW (brake) |

350

2500

300

2000

Torques [N.m]

1000

200

500

150 0

100

-500

50

-1000

0 -200.00

-100.00

0.00

100.00

200.00

300.00

400.00

-50

500.00

-1500

600.00

-2000

Crank Angle [deg]


113

Inertial Force [N]

1500

250

Využití modelu motoru pro stanovení okrajových podmínek. Nerovnoměrnost otáčení. INSTANTANEOUS ANGULAR SPEED 4 cylinder 3017 min-1 | 34.5 kW (brake) |

317.00

Angular Speed [rad.s-1

316.50 316.00 315.50 315.00 314.50 314.00 -200

-100

0

100

200

300

400

Crank Angle [deg]


114

500

Mechanické ztráty motoru - model dle OBEH • Zjednodušená závislost na časově středním tlaku ve válci a na otáčkách; parametry možno zjistit regresí na výsledky VYVAZ.XLS

(

p z = p z ,0 + K1 pVA − pVA,ref

)

   n nM M + K2  − 1  + K 3   n M ,mech ,ref   n M ,mech ,ref   

2    − 1    

• časově střední tlak respektuje jednoduše tření na pístu a v ložiskách, které není jednoznačně závislé ani na pe, ani na pmax • pVA,ref se s výhodou využije, pokud je znám z protáčení při nM,ref ztrátový tlak (moment) pz,0; ten je pro přepočet pro jiné zatížení potřebné korigovat právě o vliv středního tlaku ve válci; při běžném protáčení se střední tlak spočte pro cyklus komprese-expanze-výfuk-sání (např. pomocí OBEH pro velmi malou (ale nenulovou) dávku paliva), při protáčení se zavřenými ventily se vypočte cyklus komprese-expanze s takovým počátečním tlakem, aby pVA,ref byl zhruba atmosférický. Matematické modelování oběhu motorů

115

Výstupy v souvislosti s turbodmychadlem - příklad interpretace výsledků: tlaky v závislosti na úhlu kliky pro válec a potrubí, „rovnotlaký“ výfuk osmiválce •


116

Výstupy: indikátorový diagram pro válec a potrubí, „rovnotlaký“ výfuk 1*8 válců

•


117

Výstupy: průtoky ventily a turbinou, výkon turbiny; „rovnotlaký“ výfuk 1*8 válců •


118

Výstupy - průběh parametrů turbiny (součinitel hltnosti, účinnost, tlakový poměr, rychlostní poměr); - „rovnotlaký“ výfuk osmiválce 1*8 válců •


119

Výstupy - parametry turbiny v rychlostní charakteristice turbiny (součinitel hltnosti, účinnost, tlakový poměr); „rovnotlaký“ výfuk osmiválce 1*8 válců •


120

Výstupy: tlaky v závislosti na úhlu kliky pro válec a potrubí, „rovnotlaký“ výfuk šestiválce PRESSURES ENGINE:

• 0.5 0.45

p EXH

PRESSURE [MPa]

0.4

p CY p IN

0.35 0.3 0.25 0.2 0.15 0.1 100

150

200

250

300

350

400

450

Crank Angle [deg from CTDC]


121

500

550

600

Výstupy: indikátorový diagram pro válec a potrubí, „rovnotlaký“ výfuk šestiválce CYLINDER CHARGE EXCHANGE p-V ENGINE:

• 0.8

0.7

PRESSURE [MPa]

0.6

0.5

0.4

0.3

0.2

0.1 0

0.005

0.01

0.015

0.02

CYL. VOLUME [m3] p CYL


p IN

p EXH

122

0.025

Výstupy: průtoky ventily a turbinou, výkon turbiny; „rovnotlaký“ výfuk šestiválce Flow Rates and Turbine Power ENGINE:

•

0 600.0

4.000 m.T 3.500

m.S

PT

400.0

P T ise

2.500 2.000

300.0

1.500

200.0

1.000 100.0 0.500 0.0

0.000 -0.500 100

150

200

250

300

350

400

450



123

500

550

-100.0 600

Turbine Power [kW]

3.000

Flow Rate [kg.s-1]

500.0

m.V

Výstupy - průběh parametrů turbiny (součinitel hltnosti, účinnost, tlakový poměr, rychlostní poměr); - „rovnotlaký“ výfuk šestiválce TURBINE PARAMETERS = f(CA) ENGINE:

• u/c s

mu T

pi T

3

1.1

2.5

1 2

0.9

1.5

0.8 0.7

1

0.6 0.5

0.5 0.4 100

0 150

200

250

300



124

350

400

Pressure Ratio pi T [1]

Isentropic Efficiency eta T; Discharge Coefficient mu T; Velocity Ratio x=u/c s [1]

eta T1.2

Výstupy - parametry turbiny v rychlostní charakteristice turbiny (součinitel hltnosti, účinnost, tlakový poměr); „rovnotlaký“ výfuk šestiválce 1*6 válců Turbine Characteristics

• 3

1.2

Discharge Coefficient Isentropic Efficiency; Discharge Coefficient [1]

1

2.5

Isentropic Efficiency

0.8

2

0.6

1.5

0.4

1

0.2

0.5

0

0 -0.1

0.1

0.3

0.5

0.7

0.9

1.1

Velocity Ratio [1]


125

1.3

Pressure ratio [1]

Pressure Ratio

Výstupy: tlaky v závislosti na úhlu kliky pro válec a potrubí, pulsační výfuk 2*3 válce

•


126

Výstupy: indikátorový diagram pro válec a potrubí, pulsační výfuk 2*3 válce

•


127

Výstupy: průtoky ventily a turbinou, výkon turbiny; pulsační výfuk 2*3 válce

•


128

Výstupy - parametry turbiny v závislosti na úhlu kliky (součinitel hltnosti, účinnost, tlakový poměr); pulsační výfuk 2*3 válce •


129

Výstupy - parametry turbiny v rychlostní charakteristice turbiny (součinitel hltnosti, účinnost, tlakový poměr); pulsační výfuk 2*3 válce •


130

Dodatek

Přechodové režimy - viz Transient.ppt

Dvouzónové a vícezónové modely pro spalování

1-D nestacionární proudění v potrubích - MKO, charakteristiky, počáteční a okrajové podmínky, CFL kriterium

.


131


1.




Vícezónové (Q-D) modely Zákony zachování

Úvod


6. Závěry

d {s m}j

Zachování látek: • přítok, • odtok, • chemická kinetika.

N

=

∑ i

j

= ∑ V&P i

 {s m}i  .γ j , i . V  i

j, i +

{s m}j

d {s m}CH d t

Zákon zachování hybnosti ve vektorové formě platné pro kartézské souřadnice Zákon zachování energie pro entalpii otevřené soustavy

d t

Nj

Vj

j

=

 d {s m}CH .δ j, i  +  d t 

s*r C .

j

d {r m}CH

j

d t

r Nj    r mi  dw m r j j r r = ∑ V& j , i . wi . .γ j, i + w j . .δ j, i  − p i .Ai , j .ni , j − τ i, j . Aτ i, j  − mj.   d t Vi Vj  i     r r r d mj − wj. + ∑ Fi , k , kapka → plyn + ∑ Fi , p , prekázka → plyn + ... d t ki pi dH dt

j

=

  H Ki .γ  V& j , i .  Vi  

∂H ∂T j ,i

+

. H

dT

j

dt Kj

.δ

Vj

+ j ,i

∂H ∂p   +  


.

∑ i

dp dt

j

T  ∂ H  d {s m }j + s =  . m dt ∂  

 α Q .A i , j .(T i - T j )  + 

132

∑ k

Q&

k, j

+ V j.

dp dt

j


Úprava základních rovnic Q-D modelu šíření plamene s homogenním tlakovým polem

Teplo z chemické reakce

1.

Úvod





6. Závěry

Teplo přenesené do plynu

{H u }

• d {r mi } •  •  − ∑ α i , j S i , j Ti − TCH ,i , j + m i −1,i  {i i −1 } − m i ,i +1  {i i } = E i = r*x C dt j    

= {i i }

d {mi }  +  c pi dt 

T

T

(

{ }

T

T

)

Ti  dp {mi } − Vi  + c p ,i p  dt

{ }

Derivace společného tlaku

T

T

({ } {m })

Ti dVi {mi } − c pi Vi dt

T

i

Ti {r}

d {mi } {r}T {mi } dt T

Korekce změny entalpie (Kirchhoff)

•     T T   Ei {ii } {r}  d {mi } − dV  dp 1   = V − −  ∑ i  T T T  dt  dt dt 1  c { m } T c { m } T { r } { m } i   Vi i i i i i   p ,i   p ,i    ∑ Vi −  T p c p ,i {mi }Ti  i   Celková změna objemu

{ }

{ }

{ }

• 1 dTi dp  T d {mi } = + Vi   Ei − {ii } T dt dt dt  c pi {mi } 

{ }

Teplotní změny v zónách Matematické modelování oběhu motorů

Změna složení - DR neobsahují další derivace 133


1.





6. Závěry

Vazebné vztahy pro 2 zónový model plamene: neshořelá-1/shořelá-2 zóna

Úvod

toky mezi vloženými zónami

transformační matice stechiometrických koeficientů r složek a x reakcí

d {m i }  •   •  d {r m i } = m i −1,i  − m i ,i +1  + r * x C ; dt dt     2 zones 1 ≡ unburned 2 ≡ burning & burned

hlavní složka(y) reakce(í)

hmotnostní konc. {m1 }  •   •  A f w f ( relative ) čerstvé směsi m 0 ,1  = {0} ; m1, 2  = − V1     d { r m1 } d {r m 2 } dx = 0 ; for r = 1 =− 1 m 2 .... ROHR dt dt dt d {r m 2 }  •  =  r m1, 2  .... immediate burning of newly coming fuel dt   w d {r m 2 } {r m 2 } = = { r m 2 } b .... delayed burning (e.g ., due to flame dt τb wf rychlost postupu “zadní fronty” plamene další stupeň volnosti


134


1.

Úvod





6. Závěry

Geometrické údaje (plochy, objemy) zón v závislosti na poloze geometricky definovaného (např. hemisférického) plamene a poloze pístu Automaticky generovaná geometrická data pomocí CAD programu Úhel kliky Plocha fronty plamene

(poloha pístu

Vzdálenost od svíčky

Spalovací prostor 4 ventilové hlavy


135

Simulace a analýza hoření pomocí zónového modelu válce zážehového motoru 1.

Úvod





6. Závěry

1.0E+00 9.0E-01 8.0E-01

4 ventilová hlava

Normovaný vývin tepla [1]

relativní množství shořelé směsi [-]

Průběhy hoření (2 zónový model) pro stálou rychlost plamene

7.0E-01 6.0E-01 5.0E-01

2V hlava válců 4V hlava válců

4.0E-01 3.0E-01

2 ventilová hlava

2.0E-01 1.0E-01

úhel natočení klikového hřídele a [ o]

Úhel kliky [°]


136

50

45

40

35

30

25

20

15

10

5

0

-5

-10

-15

-20

-25

-30

0.0E+00


1.




Rychlost vývinu tepla

Úvod

1,1E-05


6. Závěry

1,0E-05

Normovaná rychlost hmotnost hořící směsi [kg] vývinu tepla [1/°]

9,0E-06

2V hlava válců

4 ventilová hlava

8,0E-06

4V hlava válců

7,0E-06 6,0E-06

2 ventilová hlava

5,0E-06 4,0E-06 3,0E-06 2,0E-06 1,0E-06

úhel natočení klikového hřídele a [o]

Úhel kliky [°]


137

50

45

40

35

30

25

20

15

10

5

0

-5

-10

-15

-20

-25

-30

0,0E+00

1-D bez pohyblivých hranic = 1-D proudění v potrubích: Základní toky do bilancí w+d w x

p +dxp ρ +dx ρ i +dxi A +dA

hmotnostní tok •

•

•

∂m m = ρ wF A = ρ w A ; dx m = dx ; wG = 0 ∂x •

•

dt (ρ) Adx

m+ dx m = (ρ + dx ρ )( w + dx w )( A + dA ) = = ρ w A + ρ w dA + wA dx ρ + ρ A dx w

• hybnostní tok •

w m = wρ A w • d x  w m  = ρ w 2 dA + w 2 A d x ρ + 2 w ρ A d x w =   •

dt (ρ w) Adx ρ w p τ i A

w +dxw

w*dm

dt (ρ i) Adx

αn

τ

n; |n|=1 dx

•

= m d x w + wd x m

• působící síly = tlak a vazké napětí na stěnách

x

   dA w  w r r F = ∫∫ ( pnx + τt x ) dS + = ∫∫  p cosα n, x − τ sinα n, x dS = ∫  p + τO dx w dx w  L δV δV   d x p cosα n, x dS w O w O dT ρw2 d x F = d x ( pA) = − p dA − A d x p + p cosα n, x dS + = − A d x p ; d x pz = = − τ dx = − λ dx 2 A w A w 4A 2


138

1-D proudění v potrubích: Parciální diferenciální rovnice pro zachování hmotnosti a hybnosti • zachování hmotnosti ∂ (ρ Adx ) ∂ρ w A   dx  = ρ wA−ρ wA+ ∂t ∂x   ∂ρ ∂w ∂ρ ∂ (ρ w A ) ∂ρ d ln A A + =0; +w = −ρ − ρw dx ∂t ∂x ∂t ∂x ∂x

• zachování hybnosti • tok hybnosti (přítok zvyšuje akumulaci v kontrolním objemu, odtok snižuje) • tlakové síly • síly tření (stac. 1-D) • konzervativní tvar • po dosazení kontinuity

 • d x  w m  = ρ w 2 dA + w 2 A d x ρ + 2 w ρ A d x w ;  

∂ (ρ ∂ ( Aρw ) = −A ∂t ∂x

∂ (ρ w ) ∂w ∂ ( Aρw ) ∂p O ρ ww + A ρw +w = −A − Aλ ∂t ∂x ∂x ∂x 4A 2 ∂w ∂w 1 ∂p O ww +w =− −t ; t = λ ρ ∂x ∂t ∂x 4A 2 A


139

)

1-D proudění v potrubích: Parciální diferenciální rovnice zachování energie •

∂T  Odx ∂T  d x Q = q Odxdt + d x  − λQ A α Q (Tst − T )dt − λQ d x  A  dt =  dt ; ∂ x sin ∂ x α     dW = pd tV = 0 pro w G = 0

• toky energie

• sdílení tepla ze stěn a vedením; práce nulová • tok klidové entalpie (přítok zvyšuje akumulaci, ...)  w2  dI = m i0 dt = mc p (T − Tref ) +  dt 2  •

•

• • •  w 2  w2  • d x I = d x m c p (T − Tref ) +  = c p m d xT + m wd x w + c p (T − Tref ) +  d x m 2  2    •

• zachování energie pro u0=u+w2/2 ∂ (ρu0 Adx ) ∂ (ρ w A i0 ) dt = − dxdt + ∂t

O ∂  ∂T  dxα Q (T − Tst )dt − λQ A  sin α ∂x  ∂x 

∂x

 w2   ∂  ρ c vT − rTref + ρ •  2 • ∂T • 2  ∂w w ∂m O ∂ 2T ∂T dA   αQ (T − Tst ) − λQ A 2 − λQ A = −c p m − mw − c p T − Tref + +  ∂t ∂x ∂x  2  ∂x sin α ∂x dx ∂x

(

)

(

)

∂T ∂ρ ∂w ∂T ∂ρ O ∂ 2T ∂T d ln A 2 ∂w cv ρ + c vT − rTref + ρw = − c p ρw + c p T − Tref − ρw + αQ (T − Tst ) − λQ 2 − λQ ∂t ∂t ∂t ∂x ∂t ∂x A sin α ∂x dx ∂x

(

)

[ (

)]


140

1-D proudění v potrubích: Parciální diferenciální rovnice zachování energie • zachování energie pro cp=const.; cp-cv=r ; Tref = 0 K (lze jen v tomto případě konstantních měrných tepelných kapacit) ∂T ∂T p ∂ρ ∂w 2 ∂w cv + cpw − +w +w =q ∂t ∂x ρ 2 ∂t ∂t ∂x λQ ∂ 2T λQ ∂T d ln A O z hybnosti α Q (T − Tst ) − q= − 2 ρ ∂x ρ ∂x dx Aρ sin α r ∂T p ∂ρ w ∂p ∂T κr w + − 2 − = wt + q κ − 1 ∂t κ − 1 ∂x ρ ∂t ρ ∂x • dosazení ze stavové rovnice a její derivace pro vyloučení teploty dp dρ dT p ∂ρ ∂p ∂T − = ; = − rρ výkon vazkých p T ρ ρ ∂t ∂t ∂t napětí („třecí r ∂T ∂T ∂T 1 ∂p w ∂p κr w − = q + wt teplo“ - ve 3D +r + − ∂t κ − 1 ∂x ρ ∂t ρ ∂x κ − 1 ∂t jen část tohoto členu!)  ∂p + w ∂p  − κp  ∂ρ + w ∂ρ  = ρ (κ − 1)(q + wt )     ∂ t ∂ x ∂ t ∂ x ρ  141 Matematické modelováníoběhumotorů

1-D proudění v potrubích: Parciální diferenciální rovnice – rekapitulace Použitelné v této formě, pokud není sdílení tepla vedením v plynu ve směru x ∂ρ + ∂t

−

κp ∂ρ + ρ ∂t

w ∂w + ∂t

∂ρ + ∂x

κp ∂ρ ∂p + − w + ρ ∂x ∂t

Pozn.: Tedy pro vektor derivací s koeficienty, které jsou jen funkcí (nederivovaných) neznámých:

∂w ∂x ∂w w + ∂x

ρ

= 1 ∂p ρ ∂x ∂p w ∂x

=

− ρw

d ln A dx

−t

= ρ (κ − 1)(q + tw)

λQ ∂ 2T λQ ∂T d ln A O − q= α Q (T − Tst ) − Aρ sin α ρ ∂x 2 ρ ∂x dx O ww t =λ 4A 2

A {yt }+

B {y x } =

{P}

Předpokládejme, že známe hledané funkce podél jisté křivky, na níž leží bod x0, t0. Křivka určuje vazbu dt a dx v okolí x0, t0. Hledáme podmínky pro zjištění přírůstku hledaných funkcí v libovolném směru v okolí bodu x0, t0 (Cauchyho úloha). Matematické modelování oběhu motorů

142

1-D proudění v potrubích: Parciální diferenciální rovnice - pojem a rovnice charakteristiky Soustavu pro neznámé derivace dt ∂ρ + ∂t pak doplníme o podmínku ∂w přírůstku na křivce, určující dx a dt + ∂ t dt:

dx

= dρ dx

dt

∂p + ∂t

∂w ∂x

= dw dx

A {yt }

tedy (E je jednotková matice):

∂ρ ∂x

B {y x }

∂p ∂x

=

=

{P}

dt E {yt }+ dx E {y x } =

{dy}

+

dp

Neznámé derivace v okolí křivky lze určit jen tehdy, když má lineární soustava nenulový determinant soustavy (Cramerovo pravidlo). Na mezi řešitelnosti leží případ, kdy je křivka tzv. charakteristikou soustavy PDR. Pak musí být současně

det

A

B

dt E

dx E A

hodnost dt E Matematické modelování oběhu motorů

=0 B dx E

143

{P}

{dy} = hodnost

A dt E

B dx E

1-D proudění v potrubích: Parciální diferenciální rovnice – podmínky na charakteristikách A det dt E

B =0 dx E A

hodnost dt E

B dx E

{P} {dy} = hodnost

A

B

dt E

dx E

Rozvojem determinantu soustavy podle části s jednotkovými maticemi se dostane kubická rovnice charakteristik Rozvojem determinantů rozšířené matice (podmínka stejné hodnosti obou matic) se dostanou podmínky na charakteristice, pro případ bez vedení tepla plynem

− w 3 dt 3 + 3w 2 dt 2 dx − 3 wdtdx 2 + dx 3 − 2

  (dx − wdt )3 −  κp  dt 2 ( dx − wdt ) = 0  ρ  dx =w 1. dt κp dx 2.,3. =w± = w ± κrT dt ρ


144

κp 2 dt ( dx − wdt ) = 0 ρ

1-D proudění v potrubích: Parciální diferenciální rovnice – Riemannovy invarianty Rozvojem determinantů rozšířené matice (podmínka stejné hodnosti obou matic) se dostanou podmínky na charakteristice, pro isoentropický případ po integraci Riemannovy invarianty. Další viz např. Jenny. Pro intensivní vedení tepla se soustava PDR musí rozšířit o další rovnice, které vzniknou rozpisem derivací druhého řádu. Vedení tepla rychlost zvuku ovlivňuje, tření a přestup v 1-D případě ne (ovšem mění rychlost zvuku po částech v jednotlivých elementech).

det

A

B

dt E

dx E A

hodnost dt E

=0 B dx E

{P}

A {dy} = hodnost dt E


B dx E

2 dw = ± da κ −1 2 (a − a0 ) w − w0 = ± κ −1 145

Dodatek: 1-D proudění v potrubích: Parciální diferenciální rovnice zachování energie ve formě pro entropii (uplatní se pro podmínku na 3. charakteristice w v obecném případě) • dosazení z Gibbsovy rovnice do zákona zachování energie dT p 1 r  dp dρ  rρ r dp κr dρ   − 2 dρ = + d = − − T T ρ κ −1 p ρ  ρ κ −1 p κ −1 ρ ∂p  κ r  ∂ρ ∂ρ  rρ (q + wt ) r  ∂p +w − +w   = (κ − 1)p  ∂t ∂x  (κ − 1)ρ  ∂t ∂x  p  ∂s + w ∂s  = rρ (q + wt )   ∂ t ∂ x p  

ds = c v

• entropie se mění v důsledku (vratného) přívodu a (nevratného) vedení tepla i (nevratné) třecí ztráty • neznámé: ρ, w, p (pokud se počítá s vnitřním vedením tepla, pak nutno upravit ještě i člen s derivací teploty) Matematické modelování oběhu motorů

146

Numerické metody: Cauchyho úloha pro hyperbolické PDE V rovnicích MKO jsou prostorové derivace ukryty v integrálech konvektivních a difusních členů. Pro PDE jde o s elegantním fyzikálním významem pro hyperbolické (zobecněné vlnové) PDE. Počáteční úlohu v KO lze řešit, pokud není jeho obsah zasažen za časový krok interferencí tlakových vln, vycházejících v čase počáteční podmínky z jeho nejvzdálenějších míst. Pro 1-D izoentropický případ je tedy (CFL vyjadřuje „bezpečnost“ proti interferenci): max ( w ± a ) ∆x CFL = < 1 ; ∆t < .CFL t ∆x max ( w ± a ) ∆t ∆t

t0

∆x

x


147

Nestacionární teplotní pole při periodickém průběhu teploty na povrchu polomasivu Teplotní vlny ve stěně - rovnice Nestacionární 1-D vedení v tuhé fázi: ZZE bez konvekce uvnitř a prakticky bez konání práce v důsledku malých objemových změn.    Q p − Qo dt = dH − Vdp ; dp = 0   • ∂T  • •  ρ A q − A q + d q dt = c Adx dt   x p  dt    •

•

!

•

∂q ∂T = cpρ ∂x dt • ∂T ∂ 2T ∂T q = −λ Q ⇒ λQ 2 = c p ρ dx dt ∂x ! λ ∂ 2T ∂T Q a 2 = ; a= dt cpρ ∂x

−

Rovnice je lineární, tedy libovolná lineární kombinace partikulárních řešení bude také řešením. Rozdělme hledané part. řešení na dvě nezávislé funkce t, x. Pak jejich poměr musí být konstantní, konstantu zvolíme tak, aby byla splněna periodická okrajová podmínka. ! ∂θ ∂T θ = T − T∞ =τ (t )ξ ( x ) ; =


∂...

∂...

d 2ξ ( x ) dτ (t ) a ( t ) = τ ξ (x ) 2 dt dx d 2ξ ( x ) dτ (t ) a dx 2 = dt ⇒ nutne = const. ξ (x ) τ (t ) 148

Nestacionární sdílení tepla - odhad teploty povrchu součástí Nestacionární jednorozměrné vedení tepla, Fourierova rovnice : Teplotní vlny v povrchové vrstvě součástí spalovacího prostoru. Obecně až pro okrajovou podmínku 3. druhu, tj. daný přestup tepla do povrchu s kolísající teplotou 2 ∂Θ ∂ Θ = a. 2 ; Θ = T - T∞ ; ∂t ∂x ∞

•

j.i.ω .t = A ; q okoli = α Q (Θokolí − Θ´ x =0 ) Θokolí ∑ i e i =0 Matematické modelování oběhu motorů

149

Nestacionární teplotní pole při periodickém průběhu teploty na povrchu polomasivu Teplotní vlny ve stěně – partikulární a úpné řešení ! dτ (t ) e jωt + e − jωt e jωt − e − jωt kt jωt − Kτ (t ) = 0 ; τ = e ; k = K ; if K = jω then τ = e = + j = cos ωt + j sin ωt dt 2 2j

ODE 1. řádu pro časovou funkci musí mít řešení s goniometrickými funkcemi, tedy konstanta musí být imaginární číslo. ! d 2ξ ( x ) K − ξ (x ) = 0 ; ξ = e px ; 2 a dx

p1, 2 = ±

K =± a

(

)

jω ω jπ 2 ω ω (1 + j ) =± e =± cos π + j sin π = ± 4 4 a a a 2a

Z obou kořenů nepřipadá v úvahu kladný, neboť by amplituda teploty s rostoucí hloubkou pod povrchem rostla, což neodpovídá zkušenosti.

j ωt

θ part = e e

−

ω 2a

(1+ j ) x

Pak zřejmě teplota na povrchu kmitá s frekvencí ω a směrem do hloubky se šíří v harmonických vlnách s exponenciálně klesající amplitudou. Libovolné okrajové podmínce na povrchu lze vyhovět ve formě Fourierovy řady, jejíž koeficienty se určí z okrajové podmínky. ω − x   ω  ω  2a    x  + B cos ωt − x  θ =e  A sin  ωt − 2a  2a    


150

Nestacionární teplotní pole při periodickém průběhu teploty na povrchu polomasivu Teplotní vlny ve stěně – druhy okrajových podmínek θ (x = 0) = f per (ωt ) = ∑ ( Ai′ cos i ωt + Bi′ sin i ωt ) i

 ∂θ  q i ( x ) = −λQ   =  ∂x  •

= λQ

iω − e 2a

iω x 2a

   iω  iω   iω − x  + B cos ωt − x  + λQ e  A sin ωt − 2 a 2 a 2 a      

iω x 2a

   iω  iω   x  − B sin  i ωt − x   A cos i ωt − 2 a 2 a     

iω  ∂θ  [( A − B )sin (i ωt ) + ( A + B ) cos(i ωt )] q i ( x = 0 ) = −λQ   = λQ 2a  ∂x  x =0 •

Pro danou teplotu povrchu v čase – OP 1. druhu, pro daný tepelný tok jde o okrajovou podmínku OP 2. druhu. OP 1. druhu se výhodně použije při měření filmovým teploměrem – z průběhu teploty lze určit fluktuující tepelný tok. Pro analytické řešení motoru se často používá přibližně OP 2. druhu, protože se ukazuje, že teplotní kmity na povrchu tepelně vodivé součásti jsou malé (desítky K), zatímco teplota plynu kolísá o 1500 – 2000K. Tedy  ∂θ  q( x = 0 ) = −λQ   = α g (Tg − T∞ − θ ) ≈ α g (Tg − T x =0 ) = ∑ (C i cos i ωt + Di sin i ωt )  ∂x  x =0 i •


151

Nestacionární teplotní pole při periodickém průběhu teploty na povrchu polomasivu Teplotní vlny ve stěně – druhy okrajových podmínek •

n

q ( x = 0 ) = λQ ∑ i =0

n iω  [( A − B )sin (i ωt ) + ( A + B ) cos(i ωt )] = α g  Tg − T∞ − ∑ [A sin (i ωt ) + B cos(i ωt )] 2a i =0  

OP 3. druhu s daným součinitelem přestupu a proměnnou teplotou povrchu při dané periodicky proměnné teplotě plynu Tg je však jedině přesná pro motor. Tuto podmínku lze řešit analyticky srovnáním s Fourierovým rozvojem teploty plynu jen při konstantním součiniteli přestupu tepla, což však není případ motoru. Pak nezbývá než iterační výpočet – viz dále. Povrchová teplota by činila při αg= 500 W.m-2.K-1=const. pro 500 min-1 a dvoudobý motor (cca 1000 min-1 pro čtyřdobý motor) jen 0.45% amplitudy Tg (např. 4,5 K z 1000 K) pro stěnu ze šedé litiny. Pro litinu je a=60/7280/540=1,5e-5 m2.s-1. Pro hloubku 1 mm je při otáčkách dvoudobého motoru 500 min-1 útlum fluktuace povrchové teploty na 27%, pro 3000 min-1 na 4%, tedy 1-D předpoklad bude zřejmě použitelný.


152

Nestacionární sdílení tepla - odhad teploty povrchu součástí s okrajovou podmínkou 3. druhu při konstantním součiniteli přestupu tepla Θ = Θ A . e j.i.ω .t .

cosh( µ .x) - sinh( µ .x) 1+ µ.

µ = (1 + j).

Θ x=0,max = Θ A .

λ steny αQ

;

i.ω ; 2.a 1

1+ µ.

λ αQ

Dosavadní analytická řešení neuspokojivá s ohledem na současnou proměnlivost teploty a součinitel přestupu tepla. Jen numerické řešení nebo krajní odhady jsou pak možné. Vliv adiabatizace vrstvou úsad by neměl být podceňován. Matematické modelování oběhu motorů

153

Nestacionární teplotní pole při periodickém průběhu teploty na povrchu polomasivu Teplotní vlny ve stěně – numerický výpočet pro okrajovou podmínku 3. druhu s proměnlivým součinitelem přestupu OP 3. druhu s

Tepelný tok stěnou

daným součinitelem

2250.00

3 500 000

přestupu a 3 000 000

proměnnou teplotou

Zadaný průběh q°

povrchu při dané proměnné teplotě plynu Tg řeší iterační výpočet pomocí podmínky 2. druhu s opravou toku podle

Hustota tepelného toku [W.m-2]

periodicky

Fouriérův polynom q° Teplota plynu

2 500 000

1750.00

Teplota stěny

2 000 000

Součinitel přestupu

1250.00

1 500 000

750.00

1 000 000

500 000 250.00 0

předpokládaného průběhu povrchové

-500 000 -180

-250.00 -90

0

90

180

270

Úhel kliky [st.]

teploty. Konverguje i pro špatně vodivé materiály.


154

360

450

540

Nestacionární teplotní pole při periodickém průběhu teploty na povrchu polomasivu Teplotní vlny ve stěně – výpočet pro okrajovou podmínku 3. druhu Průběh teplotu na

Teplota stěny v různé hloubce

povrchu litinové

480.00

součásti pro otáčky

478.00

skutečného 476.00

zážehového čtyřdobého motoru hloubce 1 mm pod povrchem.

Teplota stěny [K]

1000 min-1 a v

474.00

Teplota stěny x=0 mm Teplota stěny x

472.00

470.00

T stěny v hloubce 1.00 mm

468.00 Harmonická analýza: Tepelný tok stěnou 500000

466.00

464.00

toku včetně průměrné hodnoty

462.00 -180

-90

0

(0. harmonická) po ½ harmonických vůči otáčkám.

90

Harmonická složka: Hustota tepelného toku [W.m-2]

400000

Spektrum tepelného

300000

200000

180

270

360

450

100000

0

-100000 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

6.5

7

7.5

8

8.5

Řád harmonické [cos,sin] amplituda


540

Úhel kliky [st.]

155

9

9.5

10

10.5

11

11.5

12

12.5

13

13.5

14

14.5

Sdílení tepla do stěn motoru Energetická bilance – vliv adiabatizace a snížení odvodu tepla (MSOT) zvýšením tepelného odporu stěn 837 828 773

780

900

856

706

856

799 603

800 700

633

600 624

500

722

400

200

217 599

267

300

600

100 200

100

0 0

100

100

0

1

0 0

AM

1 1

LI A Z M 3

[1] rL

[1] rP

[W /m 2/K] ag

[oC ] tg

0 [oC ] tT1

[oC ] tL 1

[oC ] tP1

M SO T 2


156

M SO T 1

Energetická bilance – adiabatizace: redukce spotřeby měrné paliva 2 1 0 -1 lambda [1] -2 del m pe [%] -3 -4 -5

l [1]

-6 -7 LIAZ M3

del mpe [%] MSOT 1

MSOT 2


AM

157

Náčrt postupu kalibrace 0-D modelu přeplňovaného motoru • znalost parametrů v úplné nebo alespoň na vnější charakteristice – pe, mpe, pK2 , m°k, TWCH motoru i mezichladiče, pa, Ta, α1v nebo 1z – pozor na motor management, tlaky pK1 a pT2, další regulované parametry - λ, průřez wastegate nebo VTG, Ts, TT1(špatně přímo využitelná); znalost geometrie motoru včetně přesného ε a typu turbodmychadla i mezichladiče. znalost indikátorového diagramu (alespoň pmax, ideálně průběh hoření – pozor na HU a korekci teplotního ovlivnění snímače talku i plavání nuly tlaku), mechanická účinnost (protáčení vyžaduje důkladné vyhodnocení termodynamických ztrát, jednoduchá a vyhovující bývá při určité zkušenosti extrapolace spotřeby, indikace při korekci HU je nejlepší); průběh vstřiku u vznětového motoru nebo způsob regulace λ u zážehového motoru; charakteristiky klapky, pokud je použita a má se s ní počítat; Matematické modelování oběhu motorů

158

Náčrt postupu kalibrace 0-D modelu přeplňovaného motoru znalost průběhu zdvihu a průtokových vlastností ventilů (jinak odhad podle podobného provedení); znalost středních nebo alespoň místních teplot charakteristických míst stěn spalovacího prostoru, výfukového potrubí před turbinou a účinnosti mezichladiče; znalost charakteristiky kompresoru a turbiny (alespoň při u/cs=opt., pak použití univerzální charakteristiky) a její mechanické účinnosti; údaje o statických vlastnostech vestavěných regulátorů – obvod λ-sondy, karburátor nebo směšovač včetně regulace tlaku paliva, look-up tables pro elektronickou regulaci, charakteristika waste-gate nebo VTG, pokud nejsou zahrnuty v základních údajích z charakteristik motoru atp.

další speciální údaje např. pro přechodové režimy nebo výpočty nerovnoměrnosti chodu (hmotnosti, momenty setrvačnosti, tepelné kapacity stěn), údaje o připojené zátěži (statika, dynamika), chemii (složení paliva, vlastnosti katalyzátoru), údaje o dynamických vlastnostech vestavěných regulátorů atp. Matematické modelování oběhu motorů

159

Náčrt postupu kalibrace 0-D modelu přeplňovaného motoru Postup plnění vstupního souboru a postupných aproximací vstupů • geometrie a počítaný režim (parametry na vnější charakteristice – pe, nM ,pa, Ta, TWCH ,α1v nebo 1z, resp. počátek hoření, λ pro zážeh motor) • pK1 a pT2 podle podobných motorů, pozor na otáčkovou závislost (obecně s průtokem kvadratická, u katalyzátoru a některých mezichladičů lineární)

• • • •

další regulované parametry na základní hodnotu – wastegate zavřená průběh vstřiku z dávky paliva a rozumného úhlu vstřiku; pmax, ∆αs - z podobných motorů; mechanická účinnost – z VYVAZ.XLS pro podobný motor; průběh zdvihu a průtokových vlastností ventilů - z podobného motoru; střední teploty dílů – tepelné odpory z podobného motoru; charakteristika kompresoru – pro max. účinnost, a turbiny (alespoň při u/cs=opt., pak použití univerzální charakteristiky); Matematické modelování oběhu motorů

160

Náčrt postupu kalibrace 0-D modelu přeplňovaného motoru m´p = peVz1 (ηeHu )

• pe, nM, pa, Ta, TWCH

( )

′ lneni mpLt • dodržet λ pomocí dávky paliva a pK2 λ =( pK2 (rTs ))Vz1ηnap • m°k ovlivnit pK2, tedy průřezem turbiny a účinností turbodmychadla; u proplachovaných motorů překrytím ventilů, jinak rozvodem obvykle málo, pokud rozumné otáčky (SZ, překrytí, průřez ventilů); pK1 a pT2 mohou být důležité; při vstřiku do sání vnitřní chlazení - Ts, zejména u bohaté směsi; • optimalizovat α1v nebo 1z, resp. počátek hoření a • nastavit ∆αs, při bohaté směsi ηchem; na tvaru přívodu tepla tolik nezáleží; • nastavit pz pro mechanické ztráty podle nM a středního tlaku ve válci; • mpe (např. velká)– pokud pe v pořádku a spotřeba paliva ne, pak • kombinovaný vliv • naplnění válce (např. velké – pak se musí projevit na m°k) a • mechanická účinnost (např. malá) a/nebo • indikovaná účinnost vysokotlaké části (např. malá – hoření rozvleklé nebo špatně umístěné – pak p max malý, teploty stěn malé, přechlazené); • velká spotřeba práce na výměnu náplně (např. velká - nevhodné ventily, malá účinnost turbodmychadla, velké odpory v potrubích, wastegate otevřená). 161 Matematické modelování oběhu motorů

Co dál Inverzní algoritmy (především pro hoření). Empirická doplnění modelů nižší hloubky • • • • • •

zákony hoření na algebraické i na diferenciální - zónové úrovni, modely zařízení pro tvorbu směsi, kinetika tvorby škodlivin, modely průtoku ventily, modely turbodmychadla, napojení na konstrukční výpočty motoru...

Bilanční rovnice pro 3-D modely – RANS a další průměrování NS. Modely turbulence, podrobnější přehled konstitutivních vztahů pro transportní veličiny. Modely chemické kinetiky a dvoufázového proudění. Více o numerických metodách pro hyperbolické PDE. Matematické modelování oběhu motorů

162

Struktura kombinovaného modelu s využitím MKO. Zdrojové členy z lagrangeovského modelu. Lagrangeovský model skupiny kapek Individuální popis skupiny kapek, vazba na eulerovský model plynu pomocí zdrojových členů hybnost=odporová síla; látka=odpařování; energie=tepelný tok konvekcí


163

Struktura kombinovaného modelu s využitím MKO. Zobecněný lagrangeovský model plamene s porézní stěnou a „dvoufázovým“ modelem KO (metoda Volume Of Fluid). Určení části KO zasaženého plamenem.


164

Numerické metody: eulerovský model a numerický transport. Odhad numerického transportu v důsledku velikosti KO. r r m& ∆c ∆x j = = D ∇ c = D. ; max( jnum ) = .∆c Př. Eulerovský popis A ∆x ∆t difuse látky v KO s homogenním obsahem. Dosáhne-li látka hranici KO, je pro další časový krok rovnoměrně rozptýlena v celém KO. Obdobně je tomu s tepelnou vodivostí. Viskozita - lze použít jen pro velmi hrubý odhad

r Q& ∆T ∆x r q& = λ ∇ T = = λ. ; max(q& num ) = .ρ .c p .∆T A ∆x ∆t 2   λ ν ∆ x ⇒ min  Dturb , turb = aturb  >> ; turb = Prturb ≈ 1 ρ .c p ∆t aturb   ∆w y ∝ ρ .w′x .w′y τ yx = ρ .ν turb . ∆x w′x2 + w′y2 + w′z 2

ν turb ∝ L.w′ = L.

2 = L. I turb .w ; 2 2 w + w y + wz 2 x

2

ν num ∝ ∆x.I turb .w ; ν turb >> ∆x.I turb .w Matematické modelování oběhu motorů

165

Numerické metody: Volba kroků řešení • Numerickou vazkost lze snížit použitím metod vyššího řádu pro náhradu prostorových derivací. • Nevýhody: • náročnější na formulaci okrajových podmínek, potíže v rozích oblasti; • částečná ztráta přesnosti při přechodu na sítě s obecnými tvary prvků (neregulární, protáhlé, trojúhelníkové / tetrahedrické...); • sklon k oscilacím řešení. • Odstranění oscilací řízenými tlumivými členy (ENO,...). • Prostorový krok shora omezen požadovanou přesností, zdola zaokrouhlovacími chybami. • Časový krok shora omezen CFL kriteriem, zdola numerickým transportem. Matematické modelování oběhu motorů

166

The End.


167

Tříválcový motor 1,2 l 2V 40 kW Cíle vývoje

Vysoký točivý moment v nízkých otáčkách Nízká spotřeba Splnění emisních norem EU IV Možná zástavba FSI Nízké náklady na údržbu Kultivovaný běh motoru Dlouhodobá kvalita Nízká hmotnost Efektivní a výhodná sériová výroba Minimalizace doby vývoje použitím CAD/CAE Matematické modelování oběhu motorů

168

Tříválcový motor 1,2 l 2V 40 kW Pracovní oběh motoru Předpověď kolísání tlaku v sacím potrubí benzinového motoru pomocí programu GT Power - srovnání s experimentem

Model motoru a průběh tlaku v sacím potrubí


169

Tříválcový motor 1,2 l 2V 40 kW Pracovní oběh motoru Předpověď změny hmotnostního toku v motoru pomocí programu GT Power


170

fyzikální principy modelování - zákony zachování a konstitutivní vztahy, empirické uzávěry; druhy modelů podle hloubky (počet souřadnic) a šířky

Recommend Documents