Simulace teplotn´ıho pole za proudov´ym motorem lehk´eho letounu

(1)

Simulace teplotn´ıho pole za proudov´ ym motorem lehk´ eho letounu

Bakal´ aˇrsk´ a pr´ ace

Studijn´ı program: B3901 – Aplikované vˇedy v inˇzenýrstv´ı Studijn´ı obor: 3901R055 – Aplikované vˇedy v inˇzenýrstv´ı Autor práce: Gleb Pokatilov

Vedouc´ı pr´ace: doc. Ing. ˇSidlof Petr, CSc.

(2)

(3)

(4)

Prohl´ aˇsen´ı

Prohlaˇsuji, ˇze svou bakaláˇrskou práci jsem vypracoval samostatnˇe jako p˚uvodn´ı d´ılo s pouˇzit´ım uvedené literatury a na základˇe konzultac´ı s vedouc´ım mé bakaláˇrské práce a konzultantem.

Jsem si vˇedom toho, ˇze na mou bakaláˇrskou práci se plnˇe vztahuje zákon ˇc. 121/2000 Sb., o právu autorském, zejména § 60 – ˇskoln´ı d´ılo.

Beru na vˇedom´ı, ˇze Technická univerzita v Liberci nezasahuje do mých autorských práv uˇzit´ım mé bakaláˇrské práce pro vnitˇrn´ı potˇrebu Tech- nické univerzity v Liberci.

Uˇziji-li bakaláˇrskou práci nebo poskytnu-li licenci k jej´ımu vyuˇzit´ı, jsem si vˇedom povinnosti informovat o této skuteˇcnosti Technickou univerzitu v Liberci; v tomto pˇr´ıpadˇe má Technická univerzita v Liberci právo ode mne poˇzadovat úhradu náklad˚u, které vynaloˇzila na vytvoˇren´ı d´ıla, aˇz do jejich skuteˇcné výˇse.

Souˇcasnˇe ˇcestnˇe prohlaˇsuji, ˇze text elektronické podoby práce vloˇzený do IS STAG se shoduje s textem tiˇstˇené podoby práce.

Beru na vˇedom´ı, ˇze má bakaláˇrská práce bude zveˇrejnˇena Technickou univerzitou v Liberci v souladu s § 47b zákona ˇc. 111/1998 Sb., o vy- sokých ˇskolách a o zmˇenˇe a doplnˇen´ı dalˇs´ıch zákon˚u (zákon o vysokých ˇskolách), ve znˇen´ı pozdˇejˇs´ıch pˇredpis˚u.

Jsem si vˇedom následk˚u, které podle zákona o vysokých ˇskolách mohou vyplývat z poruˇsen´ı tohoto prohláˇsen´ı.

31. 5. 2020 Gleb Pokatilov

(5)

Simulace teplotn´ıho pole za proudov´ ym motorem lehk´ eho letounu

Abstrakt

Tato bakaláˇrská práce se zabývá CFD simulac´ı a analýzou teplotn´ıho pole za proudovým motorem lehkého letounu a porovnán´ım teplotn´ıch pol´ı za r˚uzných provozn´ıch podm´ınek jako rychlost letu a úhel nábˇehu.

V práci jsou popsány základy mechaniky tekutin, modelován´ı turbulentn´ıho proudˇen´ı a popis numerických metod.

K ˇreˇsen´ı byl pouˇzit výpoˇcetn´ı program ANSYS Fluent. Simulace pouˇz´ıvá výpoˇcetn´ı s´ıt’ o 360 000 elementech. ˇReˇsen´ı bylo nalezeno pomoc´ı turbulentn´ıho k–ω modelu a Density-based ˇreˇsiˇce. V této simulaci se pˇredpokládá stacionárn´ı charakter proudˇen´ı.

Kl´ıˇcová slova: CFD, ANSYS, Fluent, turbulentn´ı proudˇen´ı, poˇc´ıtaˇcová simulace, proudový motor, teplotn´ı pole.

Simulation of the temperature field from a light aircraft jet engine

Abstract

This bachelor thesis operates with CFD simulation and analysis of the temperature field from a light aircraft jet engine and comparsion of these temperature fields under different operating conditions such as flight speed and angle of attack.

The work describes the basics of fluid mechanics, modeling of turbulent flow and a description of numerical methods.

The computer program ANSYS Fluent was used for the solution. The simulation uses a mesh of 360 000 elements. The solution was found using the turbulent k–ω model and the Density-based solver. Stationary flow is assumed in this simulation.

Keywords: CFD, ANSYS, Fluent, turbulent flow, computer simulation, jet engine, temperature field.

(6)

Podˇ ekov´ an´ı

T´ımto dˇekuji vˇsem, kteˇr´ı mi na mé práci pomohli. Pˇredevˇs´ım dˇekuji vedouc´ımu práce doc. Ing. Petru ˇSidlofovi, Ph.D za jeho neocenitelné rady, trpˇelivost, maximáln´ı vstˇr´ıcnost a perfektn´ı veden´ı pˇri tvorbˇe této práce. V neposledn´ıˇradˇe chci podˇekovat za vˇecné rady a kritický pohled na tuto práci Ing. Jiˇr´ımu Pokatilovi a Ing. Alexandrovi Filippovi.

(7)

Obsah

Seznam obr´azk˚u. . . 11

Seznam tabulek . . . 12

Seznam zkratek . . . 12

1 Uvod´ 13 2 Matematick´y popis 14 2.1 Dynamika kontinua . . . 14

2.1.1 Kontinuum . . . 14

2.1.2 Viskozita. . . 14

2.1.3 Euler˚uv popis . . . 14

2.1.4 Definice tekutiny . . . 15

2.1.5 Z´akon zachov´an´ı hmotnosti . . . 16

2.1.6 Eulerova pohybov´a rovnice . . . 17

2.1.7 Navier-Stokesovy rovnice . . . 17

2.2 Poˇc´ateˇcn´ı a okrajov´e podm´ınky . . . 18

2.2.1 Poˇc´ateˇcn´ı podm´ınky . . . 18

2.2.2 Okrajov´e podm´ınky . . . 18

2.3 Matematick´e modelov´an´ı turbulentn´ıho proudˇen´ı. . . 19

2.3.1 DNS - Direct Numerical Simulation . . . 19

2.3.2 LES - Large Eddy Simulation . . . 20

2.3.3 RANS - Reynolds Average Navier-Stokes . . . 20

3 Numerick´e ˇreˇsen´ı 24 3.1 V´ypoˇcetn´ı s´ıt’ . . . 24

3.1.1 Strukturovan´a s´ıt’ . . . 25

3.1.2 Nestrukturovan´a s´ıt’ . . . 25

3.1.3 Hybridn´ı s´ıt’ . . . 26

3.2 Metoda s´ıt´ı . . . 26

3.3 Metoda koneˇcn´ych prvk˚u . . . 27

3.4 Metoda koneˇcn´ych objem˚u . . . 27

3.5 Pressure-based solver . . . 28

3.6 Density-based solver . . . 28

4 Geometrie a okrajov´e podm´ınky 30 4.1 Geometrie 2D ˇrezu . . . 30

4.2 V´ypoˇcetn´ı s´ıt’ 2D modelu . . . 30

4.3 Parametry proudov´eho motoru PBS TJ100 . . . 31

(8)

4.4 Okrajov´e podm´ınky 2D modelu . . . 32

4.5 Geometrie 3D modelu . . . 33

4.6 Okrajov´e podm´ınky 3D modelu . . . 34

5 V´ysledky 37 5.1 2D simulace . . . 37

5.1.1 Rychlostn´ı pole . . . 37

5.1.2 Teplotn´ı pole . . . 37

5.2 3D simulace . . . 39

5.3 Vliv v´ypoˇcetn´ı s´ıtˇe . . . 40

5.4 Porovnán´ı r˚uzných nastaven´ı systému . . . 43

5.4.1 Simulace C, maximáln´ı rychlost, maximáln´ı úhel . . . 43

5.4.2 Simulace B, cestovn´ı rychlost, maxim´aln´ı ´uhel . . . 47

5.4.3 Simulace F, minimáln´ı rychlost, nulový úhel . . . 49

5.4.4 Simulace G, cestovn´ı rychlost, nulov´y ´uhel . . . 51

6 Z´avˇer 53 7 Pˇr´ılohy 55 7.1 Simulace A . . . 55

7.2 Simulace B . . . 60

7.3 Simulace C . . . 65

7.4 Simulace D . . . 70

7.5 Simulace E . . . 75

7.6 Simulace F. . . 80

7.7 Simulace G . . . 86

7.8 Simulace H . . . 92

(9)

Seznam obr´ azk˚ u

2.1 Spojitá deformace v tekutinˇe, zp˚usobená teˇcným napˇet´ım τ [3]. . . 14

2.2 Tok kontroln´ı plochou plochou S [3] . . . 16

2.3 Porovn´an´ı (zleva) DNS, LES a uRANS modelu turbulence proudˇen´ı [6]. . . 20

3.1 Typy bunˇek v´ypoˇcetn´ıch s´ıt´ı [7] . . . 24

3.2 Pˇr´ıklady strukturovan´ych s´ıt´ı [7] . . . 25

3.3 Pˇr´ıklady nestrukturovan´ych s´ıt´ı . . . 25

3.4 Hybridn´ı s´ıt’ [7] . . . 26

3.5 Diskretizace prostoru na mˇr´ıˇzky indexovan´ych bod˚u (i,j) v os´ach x,y [10]. . 27

3.6 N´akres kontroln´ıho objemu k ilustraci diskretizace transportn´ıch rovnic [7] 28 3.7 Algoritmus Density-based solveru [7] . . . 29

4.1 V´ypoˇcetn´ı s´ıt’ 2D modelu . . . 30

4.2 Výpoˇcetn´ı prostor s proudovým motorem a proudˇen´ım okoln´ıho vzduchu [11] 31 4.3 Názvy okrajových oblast´ı. . . 32

4.4 Geometrie letadla importovan´a do prostˇred´ı Fluent . . . 33

4.5 Kl´ıˇcov´e vstupn´ı veliˇciny . . . 34

4.6 Grafické znázornˇen´ı hranic výpoˇcetn´ı oblasti . . . 35

5.1 Rychlostn´ı pole 2D modelu . . . 37

5.2 Grafy pr˚ubˇehu teplot na povrchu ˇrezu v´yˇskov´eho kormidla . . . 38

5.3 Teplotn´ı pole 2D modelu . . . 38

5.4 Proudnice v rovinˇe R1 viz obr. 5.5. Úhel nábˇehu 25^◦, letová rychlost 77ms⁻¹ 39 5.5 Charakteristické roviny (ˇrezy) . . . 40

5.6 Porovn´an´ı proudnic model˚u o (shora) 50, 360 a 900 tis´ıc´ıch prvk˚u . . . 41

5.7 Teplotn´ı pole pro r˚uzn´e okrajov´e podm´ınky, viz tabulka 4.2 . . . 42

5.8 Simulace C. Teplota na povrchu letounu, pohled shora . . . 43

5.9 Simulace C. Teplota na povrchu letounu, pohled zdola. . . 44

5.10 Simulace C. Pr˚ubˇeh teploty v rovinˇe R3 . . . 44

5.11 Simulace C. Rychlostn´ı pole v rovinˇe R1 . . . 45

5.12 Simulace C. Proudnice veden´e rovinou R1 . . . 45

5.13 Simulace C. Proudnice veden´e rovinou boˇcn´ı roviny (outlet bottom) . . . . 46

5.14 Simulace C. Proudnice rychlosti v horizont´aln´ı rovinˇe R4 . . . 46

5.15 Simulace B. Teplotn´ı pole v prostoru . . . 47

5.16 Simulace B. Rychlostn´ı pole v rovinˇe R1 . . . 48

5.17 Simulace B. Proudnice v rovinˇe R1 . . . 48

5.18 Teplotn´ı pole v prostoru pro simulaci F . . . 49

5.19 Simulace F. Teplotn´ı pole v rovinˇe R1. . . 50

(10)

5.20 Simulace F. Rychlostn´ı pole v rovinˇe R1 . . . 50

5.21 Teplotn´ı pole v prostoru pro simulaci G . . . 51

5.22 Simulace G. Teplotn´ı pole v rovinˇe R1 . . . 52

5.23 Simulace G. Rychlostn´ı pole v rovinˇe R1 . . . 52

7.1 Simulace A. Teplotn´ı pole v prostoru . . . 55

7.2 Simulace A. Teplotn´ı pole v rovinˇe R1 . . . 55

7.3 Simulace A. Teplotn´ı pole na povrchu letounu, pohled shora . . . 56

7.4 Simulace A. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 56

7.5 Simulace A. Teplotn´ı pole v rovinˇe R3 . . . 57

7.6 Simulace A. Rychlostn´ı pole v rovinˇe R1 . . . 57

7.7 Simulace A. Proudnice v rovinˇe R1 . . . 58

7.8 Simulace A. Proudnice v rovinˇe R4 . . . 58

7.9 Simulace A. Tlakov´e pole v rovinˇe R1 . . . 59

7.10 Simulace A. Pole hustot v rovinˇe R1 . . . 59

7.11 Simulace B. Teplotn´ı pole v prostoru . . . 60

7.12 Simulace B. Teplotn´ı pole v rovinˇe R1 . . . 60

7.13 Simulace B. Teplotn´ı pole na povrchu letounu, pohled shora . . . 61

7.14 Simulace B. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 61

7.15 Simulace B. Teplotn´ı pole v rovinˇe R3 . . . 62

7.16 Simulace B. Rychlostn´ı pole v rovinˇe R1 . . . 62

7.19 Simulace B. Tlakov´e pole v rovinˇe R1 . . . 64

7.20 Simulace B. Pole hustot v rovinˇe R1. . . 64

7.21 Simulace C. Teplotn´ı pole v prostoru . . . 65

7.22 Simulace C. Teplotn´ı pole v rovinˇe R1 . . . 65

7.23 Simulace C. Teplotn´ı pole na povrchu letounu, pohled shora . . . 66

7.24 Simulace C. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 66

7.25 Simulace C. Teplotn´ı pole v rovinˇe R3 . . . 67

7.26 Simulace C. Rychlostn´ı pole v rovinˇe R1 . . . 67

7.27 Simulace C. Proudnice v rovinˇe R1 . . . 68

7.28 Simulace C. Proudnice v rovinˇe R4 . . . 68

7.29 Simulace C. Tlakov´e pole v rovinˇe R1 . . . 69

7.30 Simulace C. Pole hustoty v rovinˇe R1 . . . 69

7.31 Simulace D. Teplotn´ı pole v prostoru . . . 70

7.32 Simulace D. Teplotn´ı pole v rovinˇe R1 . . . 70

7.33 Simulace D. Teplotn´ı pole na povrchu letounu, pohled shora . . . 71

7.34 Simulace D. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 71

7.35 Simulace D. Teplotn´ı pole v rovinˇe R3 . . . 72

7.36 Simulace D. Rychlostn´ı pole v rovinˇe R1 . . . 72

7.37 Simulace D. Proudnice v rovinˇe R1 . . . 73

7.38 Simulace D. Proudnice v rovinˇe R4 . . . 73

7.39 Simulace D. Tlakov´e pole v rovinˇe R1 . . . 74

7.40 Simulace D. Pole hustoty v rovinˇe R1 . . . 74

7.41 Simulace E. Teplotn´ı pole v prostoru . . . 75

7.42 Simulace E. Teplotn´ı pole v rovinˇe R1. . . 75

(11)

7.43 Simulace E. Teplotn´ı pole na povrchu letounu, pohled shora . . . 76

7.44 Simulace E. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 76

7.45 Simulace E. Teplotn´ı pole v rovinˇe R3. . . 77

7.46 Simulace E. Rychlostn´ı pole v rovinˇe R1 . . . 77

7.47 Simulace E. Proudnice v rovinˇe R1 . . . 78

7.48 Simulace E. Proudnice v rovinˇe R4 . . . 78

7.49 Simulace E. Tlakov´e pole v rovinˇe R1 . . . 79

7.50 Simulace E. Pole hustoty v rovinˇe R1 . . . 79

7.51 Simulace F. Teplotn´ı pole v prostoru . . . 80

7.53 Simulace F. Teplotn´ı pole v rovinˇe R1. Kontrastn´ı mˇeˇr´ıtko . . . 81

7.54 Simulace F. Teplotn´ı pole v rovinˇe R2. Kontrastn´ı mˇeˇr´ıtko . . . 81

7.55 Simulace F. Teplotn´ı pole na povrchu letounu, pohled shora . . . 82

7.56 Simulace F. Teplotn´ı pole na povrchu letounu, pohled zdola. . . 82

7.58 Simulace F. Rychlostn´ı pole v rovinˇe R1 . . . 83

7.59 Simulace F. Proudnice v rovinˇe R1 . . . 84

7.60 Simulace F. Proudnice v rovinˇe R4 . . . 84

7.61 Simulace F. Tlakov´e pole v rovinˇe R1 . . . 85

7.62 Simulace F. Pole hustoty v rovinˇe R1 . . . 85

7.63 Simulace G. Teplotn´ı pole v prostoru . . . 86

7.65 Simulace G. Teplotn´ı pole v rovinˇe R1. Kontrastn´ı mˇeˇr´ıtko . . . 87

7.66 Simulace G. Teplotn´ı pole v rovinˇe R2. Kontrastn´ı mˇeˇr´ıtko . . . 87

7.67 Simulace G. Teplotn´ı pole na povrchu letounu, pohled shora . . . 88

7.68 Simulace G. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 88

7.70 Simulace G. Rychlostn´ı pole v rovinˇe R1 . . . 89

7.71 Simulace G. Proudnice v rovinˇe R1 . . . 90

7.72 Simulace G. Proudnice v rovinˇe R4 . . . 90

7.73 Simulace G. Tlakov´e pole v rovinˇe R1 . . . 91

7.74 Simulace G. Pole hustoty v rovinˇe R1 . . . 91

7.75 Simulace H. Teplotn´ı pole v prostoru . . . 92

7.76 Simulace H. Teplotn´ı pole v rovinˇe R1 . . . 92

7.77 Simulace H. Teplotn´ı pole na povrchu letounu, pohled shora . . . 93

7.78 Simulace H. Teplotn´ı pole na povrchu letounu, pohled zdola . . . 93

7.79 Simulace H. Teplotn´ı pole v rovinˇe R3 . . . 94

7.80 Simulace H. Rychlostn´ı pole v rovinˇe R1 . . . 94

7.81 Simulace H. Proudnice v rovinˇe R1 . . . 95

7.82 Simulace H. Proudnice v rovinˇe R4 . . . 95

7.83 Simulace H. Tlakov´e pole v rovinˇe R1 . . . 96

7.84 Simulace H. Pole hustoty v rovinˇe R1 . . . 96

(12)

Seznam tabulek

2.1 Turbulentn´ı konstanty . . . 23

4.1 Z´akladn´ı parametry nastaven´ı okrajov´ych podm´ınek 2D modelu . . . 33

4.2 Vstupn´ı parametry pro simulace A-H . . . 34

4.3 N´azvy jednotliv´ych simulac´ı . . . 34

4.4 Okrajov´e podm´ınky simulace C . . . 36

4.5 Okrajov´e podm´ınky simulace G . . . 36

5.1 Z´akladn´ı vstupn´ı parametry pro simulaci C . . . 43

5.2 Z´akladn´ı vstupn´ı parametry pro simulaci B . . . 47

5.3 Z´akladn´ı vstupn´ı parametry pro simulaci F . . . 49

5.4 Z´akladn´ı vstupn´ı parametry simulace G . . . 51

Seznam zkratek

CFD Computation Fluid Dynamics DNS Direct Numerical Simulation LES Large Eddy Simulation

RANS Reynolds Average Navier-Stokes SST Shear Stress Transport

FDM Finite-Difference Methods FEM Finite-Element Methods FVM Finite-Volume Methods

(13)

1 Uvod ´

Pˇri návrhu stroj˚u z r˚uzných sfér (automobilový pr˚umysl, stroj´ırenstv´ı, tˇeˇzaˇrský pr˚umysl, letectv´ı, kosmonautika, atd.) se dnes nelze obej´ıt bez poˇc´ıtaˇcových simulac´ı. Napˇr´ıklad statické zat´ıˇzen´ı nosn´ıku je moˇzné vypoˇc´ıtat na pap´ıˇre, ale ˇcasto se mus´ı ˇreˇsit daleko kom- plikovanˇejˇs´ı problémy jako napˇr´ıklad deformace trupu, nebo turbulentn´ı proudˇen´ı. V dobˇe, kdy byl výpoˇcetn´ı výkon poˇc´ıtaˇc˚u niˇzˇs´ı neˇz dnes, byly tyto úlohy jen stˇeˇz´ı ˇreˇsitelné, nebo se povaˇzovaly dokonce za numericky neˇreˇsitelné vzhledem k poˇctu operac´ı, které bylo tˇreba provádˇet. Výkon dneˇsn´ıch poˇc´ıtaˇc˚u je tak vysoký, ˇze je moˇzné ˇreˇsit i kompliko- vaná turbulentn´ı proudˇen´ı. D´ıky simulac´ım jsme schopni lépe optimalizovat stroj, nebo analyzovat jeho r˚uzné provozn´ı stavy. T´ım zaniká potˇreba vytváˇret mnoho zkuˇsebn´ıch prototyp˚u. Výroba a zkouˇsen´ı prototyp˚u je ˇcasovˇe i finanˇcnˇe velmi nároˇcné oproti tvorbˇe poˇc´ıtaˇcových simulac´ı.

C´ılem této bakaláˇrské práce je analyzovat teplotn´ı pole za proudovým motorem a popsat termické ovlivnˇen´ı výˇskového kormidla lehkého letounu. Tato práce je souˇcást´ı pro- jektu, jehoˇz podstatou je upravit lehký vrtulový letoun na proudový letoun. Ten by slouˇzil pˇredevˇs´ım pro cviˇcné úˇcely d´ıky n´ızkým náklad˚um na výrobu, provoz i jednoduˇs´ı pilo- tován´ı oproti bˇeˇzným proudovým letadl˚um. Letoun je sestrojen z kompozitn´ıch materiál˚u a kl´ıˇcové je znát teplotn´ı a proudové pole z trysky motoru, aby nedoˇslo k poˇskozen´ı výˇskového kormidla.

V term´ınu práce se nepodaˇrilo z´ıskat reálnou geometrii letounu a proto byla k analýze pouˇzita geometrie ˇceskoslovenského akrobatického a sportovn´ıho motorového letounu Zl´ın Z-50. Geometrii jsem importoval z bezplatného serveru GRABCAD [1]. Model letounu jsem zmˇenil pˇridán´ım proudového motoru PBS TJ100 [2], tento konkrétn´ı typ motoru byl zadán konzultantem V ÚTS. Výstupn´ı rychlost plyn˚u z trysky motoru bylo tˇreba vypoˇc´ıtat z údaj˚u zjiˇstˇených v podrobném popisu provozn´ıch parametr˚u motoru od výrobce. Pro potˇreby simulace byla geometrie letounu zjednoduˇsena odebrán´ım ménˇe podstatných prvk˚u pro konkrétn´ı simulaci jako je anténa, podvozek. Výsledná geometrie je ˇreˇsena pomoc´ı turbulentn´ıho k–ω modelu a Density-based ˇreˇsiˇcem.

(14)

2 Matematick´ y popis

2.1 Dynamika kontinua

2.1.1 Kontinuum

V mechanice tekutin pokl´ad´ame tekutiny za kontinuum. Je tˇreba si tento pojem definovat.

Jedná se o takový pohled na látku, ˇze jej´ı jednotlivé ˇcástice, molekuly a vztahy mezi nimi interpolujeme do nekoneˇcnˇe malých rozmˇer˚u. D´ıky tomu m˚uˇzeme k poˇc´ıtán´ı v mechanice tekutin vyuˇz´ıt infinitesimáln´ıho poˇctu.

2.1.2 Viskozita

Obrázek 2.1: Spojitá deformace v tekutinˇe, zp˚usobená teˇcným napˇet´ım τ [3].

Teˇcné napˇet´ı τ vzniká v tekutinˇe pˇri pohybu kontinua (obrázek2.1). Dynamická viskozita µ je konstantou úmˇernosti mezi teˇcným napˇet´ım τ a rychlostn´ım gradientem, v pˇr´ıpadˇe proudˇen´ı ve smˇeru x (vyjádˇreno sloˇzkou rychlosti u)

τ = µdu

dy. (2.1)

Dynamická viskozita µ je materiálová vlastnost tekutiny a závis´ı pˇredevˇs´ım na teplotˇe.

Obecn´e vyj´adˇren´ı tenzoru napˇet´ı bude uveden v dalˇs´ı kapitole.

2.1.3 Euler˚ uv popis

Jedná se o takový popis systému, kdy jsou pevnˇe zvolené souˇradnice v prostoru. Tekutinu (systém) lze pozorovat, jak proud´ı okolo. Jedná se o tzv. popis systému ”z venku”. Vyuˇz´ıvá

(15)

se Eulerových promˇenných x, y, z, t. Napˇr´ıklad tlak, nebo vektor rychlosti je vyjádˇren jako funkce prostorových promˇenných a ˇcasu p(x, y, z, t), ~v(x, y, z, t).

2.1.4 Definice tekutiny

Tekutina se definuje pomoc´ı pol´ı termodynamických veliˇcin a to tlakové pole, teplotn´ı pole, pole hustoty a rychlostn´ı pole. Rovnováˇzný stav tekutiny popisuje rovnice rovnováhy tekutin.

Pole stavov´ych veliˇcin

Mechanika tekutin operuje s rychlostn´ım, tlakovým, hustotn´ım a teplotn´ım polem, pˇri- ˇcemˇz plat´ı stavová rovnice ideáln´ıho plynu p = ρRT , kde R je plynová konstanta a lze ji vyjádˇrit jako rozd´ıl mˇerné tepelné kapacity plynu pˇri konstantn´ım tlaku a mˇerné tep. kap.

pˇri konst. objemu R = Cp− CV. Stavová rovnice dává do souvislosti tlakové, hustotn´ı a teplotn´ı pole.

p(x, y, z, t) ρ(x, y, z, t) T(x, y, z, t)







T = p

R ρ (2.2)

Rychlostn´ı pole

Jedn´a se o z´asadn´ı promˇennou v mechanice tekutin. Rychlostn´ı pole se definuje

~v(~r, t) = ~iu(x, y, z, t) + ~jv(x, y, z, t) + ~kw(x, y, z, t), (2.3) kde u, v a w jsou skalárn´ı funkce rychlosti závislé na poloze a ˇcase. Sloˇzky ~i, ~j, ~k jsou vektory standardn´ı báze prostoru R³. Vektorové pole zrychlen´ı ~a lze vyjádˇrit následuj´ıc´ım zp˚usobem

~a= d~v

dt = ~idu

dt + ~jdv

dt + ~kdw

dt . (2.4)

Pole zrychlen´ı pro x-ovou sloˇzku ax lze popsat tot´aln´ım diferenci´alem a_x = du(x, y, z, t)

dt = ∂u

∂t + u∂u

∂x + v∂u

∂y + w∂u

∂z = ∂u

∂t + (~∇ · ~v)u. (2.5) Nyn´ı se na celé pole zrychlen´ı (2.4) aplikuje totáln´ı diferenciál (2.5)

~a= d~v dt = ∂~v

∂t + (u∂~v

∂x + v∂~v

∂y + w∂~v

∂z) = ∂~v

∂t + (~∇ · ~v)~v. (2.6) Sloˇzka ∂~v/∂t se nazývá lokáln´ı zrychlen´ı. Je nulová, pokud je tok stabiln´ı. Sloˇzka (~∇ · ~v)~v se nazývá konvektivn´ı zrychlen´ı. [3]

Rovnice rovnov´ahy tekutin

K odvozen´ı se uvaˇzuje ideáln´ı tekutina. Taková tekutina má vˇsechna smyková napˇet´ı rovna nule a je dokonale tekutá. K odvozen´ı se vyjde z rovnice rovnováhy elastického kontinua (vyuˇz´ıvám Einsteinovu sumaci)

∂τij

∂x + Fi = 0. (2.7)

(16)

F_i jsou sloˇzky vektoru objemových hustot vnˇejˇs´ıch objemových sil. τij jsou sloˇzky tenzoru napˇet´ı a ten je pro ideáln´ı tekutiny definován takto:

τ =





−p 0 0

0 −p 0

0 0 −p



 (2.8)

Rovnice (2.7) a (2.8) lze vyj´adˇrit vektorovˇe v n´asleduj´ıc´ım tvaru

−~∇p + ~F = ~0. (2.9)

respektive po vydˇelen´ı hustotou ρ

−1 ρ

∇p + ~~ G= ~0 , (2.10)

kde G je hmotnostn´ı hustota objemových sil. Toto je výsledný tvar rovnice rovnováhy tekutin.

2.1.5 Z´ akon zachov´ an´ı hmotnosti

Ubytek hmotnosti v kontroln´ım objemu V je roven toku hmotnosti pˇres kontroln´ı povrch´ S. Jinými slovy to je hmotnost tekutiny, která pˇreteˇce kontroln´ı plochou S za jednotku ˇcasu smˇerem ven z kontroln´ıho objemu V , viz obrázek 2.2. Vektor rychlosti se prom´ıtne

Obr´azek 2.2: Tok kontroln´ı plochou plochou S [3]

do smˇeru normály, skalárnˇe se vynásob´ı vektor ~v a ~n. Pr˚umˇet vektoru rychlosti toku do normály k ploˇse dS se násob´ı hustotou tekutiny ρ a tou samou plochou dS, tento souˇcin je infinitezimáln´ı hmotnostn´ı tok. Celkový hmotnostn´ı tok kontroln´ı plochou se z´ıská seˇcten´ım vˇsech infinitesimáln´ıch hmotnostn´ı tok˚u pˇres celou kontroln´ı plochu S, tedy integrac´ı souˇcinu pˇres uzavˇrenou plochu S. To se mus´ı rovnat úbytku hmotnosti kontroln´ıho objemu v ˇcase

‹

ρ(~v· ~n) dS = −d dt

˚

ρ dV. (2.11)

(17)

Gauss-Ostrogradského vˇetou se uprav´ı levá strana rovnice. Na pravé stranˇe rovnice se integruje pˇres kontroln´ı objem, který je v ˇcase nemˇenný. V Eulerovˇe popisu lze zamˇenit ˇcasovou derivaci s prostorovou integrac´ı

˚

V

∇ · (ρ~v) dV = −~

˚

V

∂ρ

∂t dV. (2.12)

Obˇe strany rovnice (2.12) se integruj´ı pˇres ten sam´y objem. Rovnice se pˇrep´ıˇse do tvaru s nulovou pravou stranou a oba integr´aly spoj´ı do jednoho

˚

V

∇ · (ρ~v) +~ ∂ρ

∂t

dV = 0. (2.13)

Pokud se tyto veliˇciny povaˇzuj´ı za spojité, pak je moˇzné prohlásit, ˇze integrand rovnice (2.13) je roven nule. Rovnice v diferenciáln´ı formˇe má následuj´ıc´ı tvar

∂ρ

∂t + ~∇ · (ρ~v) = 0 . (2.14)

Tato rovnice vyjadˇruje zákon zachován´ı hmotnosti stlaˇcitelného a nestacionárn´ıho prou- dˇen´ı. ˇCasto se uvád´ı jako rovnice kontinuity. Veliˇcina ρ~v se nazývá hustota toku tekutiny, jednotka [kg m⁻²s⁻¹].

2.1.6 Eulerova pohybov´ a rovnice

Eulerovu rovnici popisuj´ıc´ı dynamiku ideáln´ı (nevazké) tekutiny lze z´ıskat pomoc´ı rovnice rovnováhy (2.10), d’Alembertova zákona a zrychlen´ı tekutiny (2.6) v této formˇe:

∂~v

∂t + (~∇ · ~v) ~v + 1 ρ

∇p − ~~ G= ~0 (2.15)

2.1.7 Navier-Stokesovy rovnice

Navier-Stokesovy rovnice popisuj´ı proudˇen´ı newtonovské tekutiny. Pˇri odvozován´ı se mus´ı vyj´ıt z obecného vyjádˇren´ı tenzoru napˇet´ı tekutin

τ =−p I + λ~∇ · ~v I + 2µ ˙e, (2.16)

kde µ je dynamick´a viskozita. ˙e je tenzor rychlosti deformace. I je jednotkov´a matice.

Parametru λ se ˇr´ıká tzv. druhá viskozita, projev´ı se v matici na diagonále spolu s tlakem.

Navier-Stokesovy rovnice (2.18) lze z´ıskat dosazen´ım tenzoru napˇet´ı (2.16) do obecn´e pohybov´e rovnice (2.17)

∂~v

∂t + (~∇ · ~v)~v = ~G+1 ρ

∇τ~ (2.17)

∂~v

∂t + (~∇ · ~v)~v = ~G1 ρ

∇(λ~∇ · ~v) +~ 1 ρ

∇ · (2µ ˙e)~ (2.18)

(18)

2.2 Poˇ c´ ateˇ cn´ı a okrajov´ e podm´ınky

Poˇcáteˇcn´ı a okrajové podm´ınky urˇcuj´ı jednoznaˇcné ˇreˇsen´ı soustavy parciáln´ıch diferenci-

´aln´ıch rovnic.

2.2.1 Poˇ c´ ateˇ cn´ı podm´ınky

Jedná se o inicializaci poˇcáteˇcn´ıch hodnot systému, nastaven´ı jednotlivých hodnot v ˇcase 0.

Takto se nastav´ı rychlostn´ı pole, tlakov´e pole, teplotn´ı pole a pole hustot:

~

v(x, y, z, 0) = ~v₀(x, y, z) (2.19)

p(x, y, z, 0) = p0(x, y, z) (2.20)

T(x, y, z, 0) = T0(x, y, z) (2.21)

ρ(x, y, z, 0) = ρ0(x, y, z) (2.22)

Pokud se poˇc´ıtá nestacionárn´ı proudˇen´ı, tedy ˇcasovˇe závislé, pak tyto rovnice urˇc´ı star- tovn´ı hodnoty, kterými výpoˇcet zaˇcne. Kdyˇz se uvaˇzuje proudˇen´ı stacionárn´ı, tedy ne- závislé na ˇcase, pak se tyto podm´ınky zadávat nemus´ı. V praxi je vˇsak dobré poˇcáteˇcn´ı podm´ınky zadat i u stacionárn´ıho proudˇen´ı proto, ˇze výpoˇcetn´ı software zaˇcne iterovat bl´ızko ˇzádouc´ıch hodnot a nemus´ı zbyteˇcnˇe ztrácet ˇcas nacházen´ım správných podm´ınek.

Nav´ıc by se zvyˇsovalo riziko, ˇze nˇekteré metody bud’ v˚ubec nezkonverguj´ı, nebo jejich výsledek se bude pˇr´ıliˇs liˇsit od skuteˇcné hodnoty. Obzvláˇstˇe d˚uleˇzité je to pˇri poˇc´ıtán´ı turbulentn´ıho proudˇen´ı.

2.2.2 Okrajov´ e podm´ınky

Zadávaj´ı se na hranici vymezené oblasti ∂Ω. Nastaven´ı okrajových podm´ınek pˇr´ımo ovli- vˇnuje ˇreˇsen´ı systému. Mezi okrajové podm´ınky napˇr´ıklad patˇr´ı nastaven´ı vektoru rychlosti proudˇen´ı na hranici vstupu (inlet), nebo nastaven´ı podm´ınek nepropustné stˇeny. Na hranic´ıch lze definovat r˚uzné typy okrajových podm´ınek. V dalˇs´ıch kapitolách je popsáno matematické vyjádˇren´ı základn´ıch dvou podm´ınek

Neumanova podm´ınka

Tato podm´ınka urˇcuje, jakých hodnot bude nabývat derivace neznámé veliˇciny φ na dané hranici. Neumanova podm´ınka vzhledem k Laplaceovˇe rovnici (speciáln´ımu pˇr´ıpadu Po- issonovy rovnice)

∆φ(x) = 0 ∀x ∈ Ω (2.23)

je psan´a

∂φ(x)

∂~n = f (x) ∀x ∈ ∂Ω. (2.24)

kde φ je neznámá funkce, x je nezávisle promˇenná, definiˇcn´ı obor funkce se znaˇc´ı Ω.

(19)

Dirichletova podm´ınka

Tato podm´ınka definuje, jakých hodnot má neznámá veliˇcina x nabývat na dané hranici.

Dirichletova podm´ınka vzhledem k Laplaceovˇe rovnici je psan´a

φ(x) = f (x) ∀x ∈ ∂Ω, (2.25)

Skalárn´ı funkce f je definována na okrajové oblasti ∂Ω. V numerických simulac´ıch se podm´ınky zakomponovávaj´ı pˇr´ımo do soustavy algebraických rovnic. Algebraický systém odvozený z numerického algoritmu je uveden v následuj´ıc´ı maticové rovnici







k_1,1 k_1,2 · k_1,m−1 k_1,m k_2,1 k_2,2 · k_2,m−1 k_2,m

· · · · ·

k_m−1,1 k_m−1,2 · k^m−1,m−1 k_m−1,m k_m,1 k_m,2 · k_m,m−1 k_m,m







·





 x₁ x₂

· x_m−1

x_m







=





 a₁ a₂

· a_m−1

a_m







, (2.26)

kde ki,j jsou prvky algebraického operátoru, xijsou neznámé v závislosti na stupni volnosti n, ai jsou známé podm´ınky. Nejjednoduˇsˇs´ı cestou k zanesen´ı Dirichletovy podm´ınky (pro n-tou sloˇzku) do systému je následuj´ıc´ı úprava soustavy rovnic







k1,1 · · · k^1,m

· · · · ·

0 0 1 0 0

· · · · ·

k_m,1 · · · k^m,m







·





 x1

· x_n

· x_m







=





 a1

· f

· a_m







, (2.27)

kde f je hodnota, kterou mus´ı nabýt n-tá sloˇzka vektoru neznámých ~x.

2.3 Matematick´ e modelov´ an´ı turbulentn´ıho proudˇ en´ı

Matematický model proudˇen´ı vycház´ı ze základn´ıch zákon˚u zachován´ı. Pomoc´ı nˇej lze popsat laminárn´ı i turbulentn´ı proudˇen´ı. Pohybové rovnice laminárn´ıho proudˇen´ı tvoˇr´ı uzavˇrenou soustavu rovnic, ˇreˇsen´ı je jednoznaˇcné a jeho nalezen´ı mnohem jednoduˇsˇs´ı oproti turbulentn´ımu proudˇen´ı. Pˇri hledán´ı ˇreˇsen´ı turbulentn´ıho proudˇen´ı vystupuj´ı jed- notlivé promˇenné tlak p, hustota ρ, vektor rychlosti ~v, termodynamická teplota T jako okamˇzité hodnoty náhodných nestacionárn´ıch veliˇcin. Zat´ım neum´ıme naj´ıt pˇresné ˇreˇsen´ı turbulentn´ıho proudˇen´ı a proto vˇzdy hledáme pouze jeho aproximaci [4].

Dnes se pouˇz´ıvá nˇekolik pˇr´ıstup˚u k výpoˇctu turbulentn´ıch proudˇen´ı. Budou popsány tˇri základn´ı metody, jedna z nich se pak vyuˇzije k ˇreˇsen´ı této bakaláˇrské práce.

2.3.1 DNS - Direct Numerical Simulation

DNS znamená ˇcesky pˇr´ımé numerické simulace Navier-Stokesových rovnic. Tato simulace vyˇzaduje velmi hustou výpoˇcetn´ı s´ıt’, aby bylo pokryto celé spektrum v´ırových struktur.

Diskretizaˇcn´ı s´ıt’ mus´ı být schopna zachytit v´ıry o rozmˇeru Kolmogorova mˇeˇr´ıtka. Pomˇer velikosti nejvˇetˇs´ıch l0 a nejmenˇs´ıch η struktur v proudovém poli udává poˇcet stupˇn˚u volnosti n. Tento pomˇer také souvis´ı s Reynoldsovým ˇc´ıslem ^l_η⁰ ≈ Re^3/4 (Kolmogorova

(20)

teorie pro izotropn´ı turbulenci). Pro poˇcet stupˇn˚u volnosti n dan´e simulace v prostoru plat´ı:

n ≈ Re⁹⁴ (2.28)

Z tohoto je patrné, ˇze výpoˇcetn´ı nároˇcnost této simulace bude velice záviset na velikosti Re a bude ji moˇzné vyuˇz´ıvat pouze pro n´ızká Reynoldsova ˇc´ısla [5].

2.3.2 LES - Large Eddy Simulation

LES, neboli metoda simulace velkých v´ır˚u je zaloˇzena na principu filtrován´ı Navier- Stokesových rovnic. Turbulentn´ı proudˇen´ı se skládá z v´ır˚u r˚uzných mˇeˇr´ıtek, od malých v´ır˚u, ve kterých prob´ıhá vazká disipace aˇz po velké v´ıry velikostnˇe srovnatelné s mˇeˇr´ıtkem celého systému. Tyto velké v´ıry zprostˇredkovávaj´ı pˇreváˇznou ˇcást turbulentn´ıho trans- portu hmotnosti, hybnosti a energie. Filtrován´ım Navier-Stokesových rovnic se oddˇel´ı malé v´ıry, které jsou menˇs´ı neˇz ˇs´ıˇrka zvoleného filtru. Výsledné rovnice pak popisuj´ı dynamiku pohybu velkých v´ır˚u [5].

2.3.3 RANS - Reynolds Average Navier-Stokes

Reˇsen´ı stˇredovan´ˇ ych Navier-Stokesových rovnic ke svému numerickému ˇreˇsen´ı vyˇzaduje z uvedených pˇr´ıstup˚u nejmenˇs´ı výpoˇcetn´ı výkon. V praxi se jedná o nejpouˇz´ıvanˇejˇs´ı metodu. Soustava tˇechto rovnic nen´ı uzavˇrená a mus´ı být doplnˇena tzv. modelem turbulence.

Pomoc´ı RANS vˇsak nelze modelovat v´yvoj mal´ych poruch v ˇcase a prostoru, tˇemito meto- dami tedy nelze spolehlivˇe pˇredpovˇedˇet pˇrechod do turbulence ani odtrˇzen´ı mezn´ı vrstvy.

Tato skuteˇcnost je patrn´a na obr´azku 2.3.

RANS vyuˇz´ıvá ke svému výpoˇctu stˇredn´ıch hodnot turbulentn´ıch veliˇcin. Okamˇzité hodnoty turbulentn´ıch veliˇcin jsou souˇctem stˇredn´ıch a fluktuaˇcn´ıch sloˇzek. V pˇr´ıpadˇe RANS je fluktuaˇcn´ı sloˇzka a⁰⁰ rovna nule [5].

A(xi, t) = ¯A(xi) + a⁰⁰(xi, t). (2.29)

Obr´azek 2.3: Porovn´an´ı (zleva) DNS, LES a uRANS modelu turbulence proudˇen´ı [6].

(21)

Model k–

Standardn´ı model k– je pouˇzitelný v dostateˇcné vzdálenosti od obtékané stˇeny, kde je vysoké turbulentn´ı Reynoldsovo ˇc´ıslo Ret. Tato vlastnost vyplývá z pr˚ubˇehu rychlosti disipace turbulentn´ı kinetické energie v bl´ızkosti stˇeny. Tento model spolu s modelem k–ω patˇr´ı k takzvaným dvourovnicovým model˚um a oba se pouˇz´ıvaj´ı u metod RANS.

Témˇeˇr vˇsechny turbulentn´ı modely pracuj´ı s veliˇcinou turbulentn´ı kinetické energie k, která je definována vztahem

k = ρu⁰_iu⁰_i

2¯ρ . (2.30)

Clen ρuˇ ⁰_iu⁰_j vyjadˇruje vliv turbulentn´ıch fluktuac´ı na pˇrenos hybnosti v tekutinˇe. Veliˇcina

¯

ρ je stˇredovan´a hustota [4].

Dvourovnicové modely vyuˇz´ıvaj´ı pˇri ˇreˇsen´ı dalˇs´ı veliˇcinu, která se dostane ˇreˇsen´ım transportn´ı rovnice, rychlost disipace . Turbulentn´ı viskozita daného modelu je pak dána vztahem

µ_t= Cµρ¯k²

, (2.31)

kde Cµ je konstanta. Pro úplnost je tˇreba uvést transportn´ı rovnice modelu k– imple- mentované v programu ANSYS bez podrobné analýzy jednotlivých rovnic a jejich ˇclen˚u.

∂

∂t(ρk) + ∂

∂x_i(ρkui) = ∂

∂x_j

µ+ µ_t σ_k

∂k

∂x_j

+ Gk+ Gb− ρ − Y^M + Sk (2.32)

∂

∂t(ρ) + ∂

∂x_i(ρui) = ∂

∂x_j

µ+µ_t σ

∂

∂x_j

+ C1

k(Gk+ C3G_b)− C2ρ²

k + S (2.33) V tˇechto rovnic´ıch Gk pˇredstavuje generován´ı turbulence kinetické energie v d˚usledku pr˚umˇerných gradient˚u rychlosti. Gb je generován´ı turbulence kinetické energie v d˚usledku vztlaku. YM pˇredstavuje pˇr´ıspˇevek kol´ısavé dilatace ve stlaˇcitelné turbulenci k celkové rychlosti rozptylu. C1, C2 a C3 jsou konstanty. σka σ jsou turbulentn´ı Prandtlova ˇc´ısla pro k a . Sk a S jsou zdroje definované uˇzivatelem [7].

Model k–ω

Tato sekce je vˇenov´ana dvˇema z´akladn´ım model˚um. Standardn´ı k–ω model a SST (Shear Stress Transport) k–ω model.

Standardn´ı model k–ω

Tento model kvalitnˇe popisuje chován´ı proudˇen´ı v bl´ızkosti stˇen a nejˇcastˇeji se vyuˇz´ıvá k výpoˇctu smykového proudˇen´ı v bl´ızkosti obtékaných stˇen. Specifická rychlost disipace se znaˇc´ı ω (ω = _k). Turbulentn´ı kinetická energie k se z´ıská z následuj´ıc´ıch transportn´ıch rovnic implementovaných v programu ANSYS.

∂

∂t(ρk) + ∂

∂xi

(ρkui) = ∂

∂xj

(Γk

∂k

∂xj

) + Gk− Y^k+ Sk, (2.34)

∂

∂t(ρω) + ∂

∂x (ρωui) = ∂

∂x (Γω

∂ω

∂x ) + Gω− Yω+ Sω, (2.35)

(22)

kde veliˇcina Gk (Gω) znamen´a generov´an´ı k (ω) v d˚usledku stˇredn´ıho gradientu rychlosti.

Y_k (Yω) reprezentuje disipaci k (ω) v d˚usledku turbulence. Sk a Sω jsou pˇreddefinované zdrojové ˇcleny. Γk (Γω) je tzv. efektivn´ı difuzivita veliˇciny k (ω). Efektivn´ı difuzivitu modelu k–ω nám definuj´ı tyto vztahy:

Γk = µ + µ_t

σ_k, (2.36)

Γω = µ + µ_t

σ_ω, (2.37)

kde σk a σω jsou turbulentn´ı Prandtlova ˇc´ısla pro k a ω. Turbulentn´ı viskozita µt pro k–ω model se z´ısk´a pomoc´ı n´asleduj´ıc´ıho vztahu

µ_t= α^∗ρ k

ω . (2.38)

Pomoc´ı koeficientu α^∗ se potlaˇcuje turbulentn´ı viskozita vzniklá pˇri n´ızkých hodnotách Reynoldsova ˇc´ısla. Pˇri vysokých Reynoldsových ˇc´ıslech plat´ı α^∗ = α^∗_∞= 1. Dále plat´ı:

α^∗ = α^∗_∞ α^∗₀ +^Re_R^t

k

1 + ^Re_R^t

k

!

(2.39)

Ret= ρ k

µ ω Rk = 6

α₀^∗ = β_i

3 β_i = 0, 072

SST model k–ω

SST kombinuje oba dva typy model˚u (k– a k–ω). Ve vnitˇrn´ı oblasti mezn´ı vrstvy (bl´ızko stˇeny) poˇc´ıtá se standardn´ım k–ω modelem a ve vnˇejˇs´ı oblasti mezn´ı vrstvy (dál od stˇen) poˇc´ıtá s modelem k–.

Pohybové rovnice SST jsou stejné jako u standardn´ıho k–ω modelu. U následuj´ıc´ıch promˇenných, nastala zmˇena. Turbulentn´ı viskozita tohoto modelu µt je dána vztahem

µ_t= ρ k ω

1 max[_α¹∗,^{S F}_α ²

1ω] (2.40)

σk = F₁

σ_k,1 +1− F¹ σ_k,2

−1

(2.41)

σ_ω = F₁

σ_ω,1 + 1− F¹ σ_ω,2

−1

(2.42) veliˇcina S je velikost deformace. F1 a F2 pˇredstavuj´ı smˇeˇsovac´ı funkce.

α = α∞

α^∗

α₀+^Re_R^t

ω

1 + ^Re_R^t

ω

!

(2.43) Re_ω = 2,95

α₀ = 1

(23)

Konkr´etn´ı nastaven´ı pro simulace

V tabulce 2.1 jsou uvedeny hodnoty konstant turbulentn´ıho modelu. Tyto hodnoty jsou nastaveny pro SST k–ω model turbulentn´ıho proudˇen´ı pro danou simulaci.

Konstanta Hodnota

α^∗_∞ 1

α∞ 0,52

β_∞^∗ 0,09

a1 0,31

β_i,1 0,075

β_i,2 0,0828

σ_k,1 1,176

σ_k,2 1

σ_ω,1 2

σ_ω,1 1,168

P rt 0,85

P r_t,w 0,85

PLCF 10

Tabulka 2.1: Turbulentn´ı konstanty

(24)

3 Numerick´ e ˇreˇsen´ı

Tato kapitola se zabývá postupy, jak numericky vypoˇc´ıtat rovnice zákonu zachován´ı hmotnosti, hybnosti a energie, jak správnˇe diskretizovat spojitý problém, jak ho co neje- fektivnˇeji a s co nejvˇetˇs´ı pˇresnost´ı vyˇreˇsit. Zat´ım nejde nalézt ˇreˇsen´ı diferenciáln´ı rovnice turbulentn´ıho proudˇen´ı jinak, neˇz rozdˇelen´ım na diskrétn´ı podmnoˇziny a numericky vyˇc´ıslovat jednotlivé promˇenné. Otázkou d˚ukazu existence spojitého ˇreˇsen´ı Navier- Stokesových rovnic se zabývá mnoho matematik˚u a patˇr´ı to také k jednomu ze sedmi matematických problému tis´ıcilet´ı.

3.1 V´ ypoˇ cetn´ı s´ıt’

Diskretizac´ı výpoˇcetn´ıho prostoru se vytvoˇr´ı výpoˇcetn´ı s´ıt’. V praxi se tomuto procesu ˇr´ıká jeho anglickým ekvivalentem

”meshing“ (s´ıt’ován´ı). Dynamika proudˇen´ı tekutin se zat´ım neum´ıˇreˇsit spojitˇe, proto se vyuˇz´ıvá procesu diskretizace k rozdˇelen´ı výpoˇcetn´ıho prostoru na malé ˇcásti a pro kaˇzdou z nich se vypoˇc´ıtávaj´ı urˇcité diferenciáln´ı rovnice. Poˇcáteˇcn´ı podm´ınky se pˇreb´ıraj´ı ze sousedn´ıch ˇcást´ı.

Obr´azek 3.1: Typy bunˇek v´ypoˇcetn´ıch s´ıt´ı [7]

Dvourozmˇerné výpoˇcetn´ı s´ıtˇe se mohou standardnˇe skládat z trojúhel- n´ıkových a ˇctyˇrúheln´ıkových bunˇek (obrázek 3.1). Prostorové výpoˇcetn´ı s´ıtˇe se mohou skládat ze ˇctyˇrstˇen˚u aˇz obrazc˚u s v´ıce stˇenami (mno- hostˇeny). Tyto obrazce nemus´ı být pravidelné, mohou se pˇrizp˚usobovat geometrii výpoˇcetn´ıho prostoru v závislosti na pouˇzitém zp˚usobu diskretizace. Rozdˇelujeme nˇekolik druh˚u s´ıt´ı. Nestrukturovanou, hybridn´ı a strukturovanou. V této práci jsem k výpoˇctu pouˇzil nestrukturovanou výpoˇcetn´ı s´ıt’. [8]

Kvalitu v´ypoˇcetn´ı s´ıtˇe lze definovat

tvarem bunˇek. Tedy jak moc se objem buˇnky liˇs´ı od ide´aln´ıho objemu buˇnky (pravideln´e).

Nebo vzájemnou ortogonalitou sousedn´ıch bunˇek. Kvalitu lze posoudit i velikost´ı jednot- livých bunˇek, tedy hustotou s´ıtˇe. Obecnˇe plat´ı, ˇze ˇc´ım hustˇs´ı s´ıt’ je, t´ım pˇresnˇejˇs´ı výsledek lze oˇcekávat. To vˇsak m˚uˇze znamenat pˇr´ıliˇs velkou výpoˇcetn´ı nároˇcnost. Proto se mus´ı naj´ıt výpoˇcetn´ı s´ıt’ s co nejmenˇs´ım poˇctem element˚u, ale aby byl pˇri tom výsledek správný.

To se ovˇeˇr´ı tak, ˇze pˇri vˇetˇs´ım zhuˇst’ov´an´ı s´ıtˇe se uˇz ˇreˇsen´ı kvalitativnˇe nezmˇen´ı.

(25)

3.1.1 Strukturovan´ a s´ıt’

Kaˇzdý element této s´ıtˇe lze oznaˇcit pomoc´ı index˚u i, j, k v závislosti na prostorové orien- taci. D´ıky tomu nen´ı tˇreba ukládat souˇradnice element˚u a t´ım sn´ıˇzit výpoˇcetn´ı nároˇcnost pro nˇekteré numerické metody. Tento typ s´ıt´ı je vhodný pro jednoduˇsˇs´ı geometrie. Vyuˇz´ıvá se u mezn´ıch vrstev. Na obrázc´ıch 3.2a a 3.2b jsou pˇr´ıklady strukturovaných s´ıt´ı.

(a) Strukturovaná ˇctyˇrúheln´ıková s´ıt’ (b) Strukturovaná multibloková 3D s´ıt’

Obr´azek 3.2: Pˇr´ıklady strukturovan´ych s´ıt´ı [7]

3.1.2 Nestrukturovan´ a s´ıt’

Nestrukturovanou výpoˇcetn´ı s´ıt’ lze pouˇz´ıt i na sloˇzitou geometrii. Jednotlivé elementy jsou uspoˇrádány v prostoru nepravidelnˇe a libovolnˇe. Je potˇreba ukládat konektivitu jed- notlivých bunˇek, t´ım se zvyˇsuje výpoˇcetn´ı nároˇcnost procesu. V praxi se tento typ s´ıt´ı pouˇz´ıvá nejˇcastˇeji. Na obrázc´ıch 3.3a a 3.3b jsou pˇr´ıklady nestrukturovaných s´ıt´ı.

(a) Trojúheln´ıková s´ıt’ [7] (b) 3D s´ıt’, kterou jsem sestrojil pomoc´ı programu ANSYS Obrázek 3.3: Pˇr´ıklady nestrukturovaných s´ıt´ı

(26)

3.1.3 Hybridn´ı s´ıt’

Obsahuje oblasti strukturované i nestrukturované s´ıtˇe (obrázek 3.4).

Obr´azek 3.4: Hybridn´ı s´ıt’ [7]

3.2 Metoda s´ıt´ı

Anglicky FDM (Finite-Difference Methods, metoda koneˇcných diferenc´ı) je metoda jej´ımˇz principem je diskretizace prostoru viz obrázek3.5, diskretizace parciáln´ıch diferenciáln´ıch rovnic popisuj´ıc´ıch daný systém a následný popis hlavn´ı stavové promˇenné v jednotlivých uzlech. Tato metoda je nekonzervativn´ı.

Derivace se tedy vyjádˇr´ı pomoc´ı diferenc´ı. Rozliˇsuje se tzv. zpˇetná diference (3.1), dopˇredná diference (3.2) a centráln´ı diference (3.3)

f⁰(x0) = f(x0)− f(x⁰− h)

h −1

2h· f⁰⁰(ξ1), (3.1)

f⁰(x0) = f(x0+ h)− f(x⁰)

h − 1

2h· f⁰⁰(ξ1), (3.2) f⁰(x0) = f(x0+ h)− f(x⁰− h)

h − h²

6 [f⁰⁰⁰(ξ1) + f⁰⁰(ξ1)], (3.3) délka kroku se znaˇc´ı h. Tato metoda má jednoduchou implementaci, ale funguje pouze na velmi jemných s´ıt´ıch [9]. U FDM se rozliˇsuj´ı 2 základn´ı pˇr´ıstupy, explicitn´ı a implicitn´ı.

Explicitn´ı

• K výpoˇctu je pouˇzita dopˇredná a centráln´ı diference.

• Oblast stability, kter´y zaruˇcuje, ˇze chyba nebude divergovat, je velice mal´y.

• Tato metoda je výpoˇcetnˇe ménˇe nároˇcná.

(27)

Obrázek 3.5: Diskretizace prostoru na mˇr´ıˇzky indexovaných bod˚u (i,j) v osách x,y [10].

Implicitn´ı

• K výpoˇctu je pouˇzita zpˇetná a centráln´ı diference.

• Oblast stability je mnohem vˇetˇs´ı neˇz u explicitn´ı metody.

• V kaˇzdém kroku se mus´ı ˇreˇsit soustava lineárn´ıch rovnic. Jedná se o výpoˇcetnˇe nároˇcnˇejˇs´ı metodu.

3.3 Metoda koneˇ cn´ ych prvk˚ u

S´ıt’ metody FEM (Finite-Element Methods) nemus´ı být pravidelná, lze ji aplikovat na r˚uzné tvary. Výpoˇcetn´ı prostor (spojité kontinuum) se rozdˇel´ı do urˇcitého (koneˇcného) poˇctu prvk˚u. Tyto malé elementy se nazývaj´ı koneˇcné prvky. FEM se primárnˇe pouˇz´ıvá k analýze mechanického napˇet´ı a deformace pevných tˇeles. Dokáˇze dosáhnout vysoké pˇresnosti i na hrubých s´ıt´ıch a lze ji pouˇz´ıt i k simulaci laminárn´ıch proudˇen´ı. Nevýhodou je vysoká výpoˇcetn´ı nároˇcnost u velkých oblast´ı.

Metoda koneˇcných prvk˚u spoˇc´ıvá v ˇreˇsen´ı neznámých funkc´ı, popisuj´ıc´ı spojité ˇreˇsen´ı problému. Metoda hledá lineárn´ı kombinace pˇredem zvolených bázových funkc´ı a nezná- mých parametr˚u ˇreˇsen´ı. Tyto parametry je pak schopna vypoˇc´ıtat z posun˚u uzl˚u s´ıtˇe.

Celý proces se tedy transformuje z hledán´ı spojitých funkc´ı na hledán´ı koneˇcného poˇctu parametr˚u (posun˚u v uzlech s´ıtˇe). [9]

3.4 Metoda koneˇ cn´ ych objem˚ u

Výpoˇcetn´ı software ANSYS, který se vyuˇz´ıvá v této práci, pouˇz´ıvá k výpoˇctu metodu koneˇcných objem˚u FVM (Finite-Volume Methods). Pˇrevád´ı obecnou transportn´ı rovnici na algebraickou rovnici, která jiˇz m˚uˇze být numericky vyˇc´ıslena. Transportn´ı rovnice se integruje kolem kaˇzdého kontroln´ıho objemu. T´ım se z´ıská diskrétn´ı rovnice, která vyjadˇruje zákon zachován´ı v kontroln´ım objemu. [10]

Pro nestacionárn´ı (ˇcasovˇe závislé) simulace mus´ı být ˇr´ıd´ıc´ı rovnice diskretizovány v prostoru i ˇcase. Prostorová diskretizace pro ˇcasovˇe závislé rovnice je identická s pˇr´ıpadem ustáleného stavu. ˇCasová diskretizace zahrnuje integraci kaˇzdého ˇclenu do diferenciáln´ıch

(28)

Obr´azek 3.6: N´akres kontroln´ıho objemu k ilustraci diskretizace transportn´ıch rovnic [7]

rovnic v ˇcasovém kroku ∆t. Jak je ukázáno n´ıˇze, integrace pˇrechodných ˇclen˚u je pˇr´ımá.

Obecný výraz pro ˇcasový vývoj promˇenné φ je dán:

∂φ

∂t = F (φ). (3.4)

Funkce F pˇredstavuje jakoukoli prostorovou diskretizaci. Pokud je ˇcasová derivace diskre- tizována pomoc´ı zpˇetné diference, je pˇresná ˇcasová diskretizace prvn´ıho ˇrádu následuj´ıc´ı [7]:

φ_n+1− φⁿ

∆t = F (φ). (3.5)

3.5 Pressure-based solver

Tento solver byl historicky vyvinut pro n´ızkorychlostn´ı nestlaˇcitelné toky. Dnes jsou vˇsak oba typy solver˚u (tlakový i hustotn´ı) rozˇs´ıˇreny a pˇreformulovány tak, aby ˇreˇsily a fungovaly pro ˇsirokou ˇskálu tokových podm´ınek nad jejich p˚uvodn´ı zámˇer. [7]

3.6 Density-based solver

Tento solver byl historicky vyvinut naopak pro vysokorychlostn´ı stlaˇcitelné toky. Rovnice kontinuity se pouˇz´ıvá k z´ıskán´ı pole hustoty, zat´ımco tlakové pole se urˇcuje ze stavové rovnice. Kaˇzdá iterace výpoˇcetn´ıho algoritmu se skládá z krok˚u znázornˇených na obrázku (3.7) a nast´ınˇených n´ıˇze [7].

1. Aktualizace vlastnost´ı tekutiny na základˇe aktuáln´ıho ˇreˇsen´ı. Pokud výpoˇcet právˇe zaˇcal, vlastnosti tekutiny budou aktualizovány na základˇe inicializace poˇcáteˇcn´ıch podm´ınek.

2. ˇReˇsen´ı spojen´ych rovnic kontinuity, hybnosti a (kde je to vhodn´e) rovnice energie.

3. ˇReˇsen´ı rovnic pro skaláry (turbulence) pomoc´ı jiˇz dˇr´ıve aktualizovaných hodnot ostatn´ıch promˇenných.

4. Kontrola konvergence sady rovnic.

(29)

Obr´azek 3.7: Algoritmus Density-based solveru [7]

(30)

4 Geometrie a okrajov´ e podm´ınky

Pˇresná geometrie letounu Zl´ın Z-50 byla z´ıskána z bezplatné databáze GRABCAD [1].

3D model byl vytvoˇren dle v´ykres˚u, ale nebyl vhodn´y pro pouˇzit´ı v programu Fluent.

M´ısto sestavy objemových prvk˚u byl tvoˇren jediným prvkem tvoˇreným trojrozmˇernými plochami. Musela se nejprve optimalizovat geometrie. V programu Solid Works se troj- rozmˇerné plochy pˇremodelovaly na objemové prvky, se kterými bylo moˇzné pracovat v progamu DesignModeler (modul od ANSYS). Posléze byl pˇridán proudový motor PBS TJ100. Upravený 3D model byl exportován ve standardn´ım formátu IGS a dále se s n´ım pracovalo pˇr´ımo v prostˇred´ı ANSYS Fluent, Design Modeler.

Byl zvolen model stacionárn´ıho proudˇen´ı. Tato práce se zaj´ımá pˇredevˇs´ım o ustálené proudˇen´ı (v ˇcase nemˇenné). Dále model poˇc´ıtá s vazkým, stlaˇcitelným proudˇen´ım ideáln´ıho plynu a jeho turbulenc´ı v prostoru. Ovlivnˇen´ı gravitac´ı bylo zanedbáno. Na výstupn´ıch hranic´ıch bylo zamezeno zpˇetnému (vratnému) proudˇen´ı zpátky do výpoˇcetn´ıho prostoru.

4.1 Geometrie 2D ˇrezu

2D model má ve srovnán´ı s 3D modelem mnohem ménˇe prvk˚u. V této práci slouˇz´ı pˇredevˇs´ım k odladˇen´ı simulace a okrajových podm´ınek systému. Odladˇen´ı na 2D modelu zabere daleko ménˇe výpoˇcetn´ıho ˇcasu.

3D geometrie byla importována do modulu Design Modeler a vytvoˇren podélný ˇrez vedený stˇredem proudového motoru. Výˇskové kormidlo bylo pro úˇcely ladˇen´ı modelu ve 2D geometrii umˇele posunuto smˇerem dol˚u do výstupn´ıho proudu z motoru.

4.2 V´ ypoˇ cetn´ı s´ıt’ 2D modelu

Obr´azek 4.1: V´ypoˇcetn´ı s´ıt’ 2D modelu

(31)

2D geometrie byla zas´ıt’ována pomoc´ı nestrukturované trojúheln´ıkové s´ıtˇe. Hrubˇs´ı s´ıt’

je rozdˇelena na trojúheln´ıkové buˇnky o stranˇe 50 mm. Podél pevných stˇen modelu, kˇr´ıdla, motoru a výˇskového kormidla se mˇr´ıˇzka lineárnˇe zjemˇnuje aˇz na buˇnky o stranˇe 15mm, viz Obrázek 4.1. Celkem tato 2D výpoˇcetn´ı s´ıt’ obsahuje 15 900 element˚u. Pˇri výpoˇctu ve 2D byl uvaˇzován pouze nulový úhel nábˇehu, proto je geometrie zjednoduˇsená a je u n´ı odstranˇena pˇredn´ı ˇcást kˇr´ıdla.

4.3 Parametry proudov´ eho motoru PBS TJ100

Parametry proudového motoru PBS TJ100 jsou pouˇzity z informac´ı z dostupných na inter- netu [2]. Tato práce se zaob´ırá pˇredevˇs´ım nejv´ıce nepˇr´ıznivou situac´ı, která m˚uˇze nastat.

Byly proto nastaveny maximáln´ı moˇzné a udrˇzitelné hodnoty výkonu a teplot. K nastaven´ı okrajových podm´ınek bylo nutné vypoˇc´ıtat hodnotu rychlosti výstupn´ıch plyn˚u ue prou- dového motoru, která nebyla v dokumentaci uvedená. Na obrázku4.2je vidˇet nákres situace za provozu. Rychlost okoln´ıho vzduchu je vyjádˇrena ~V0, pa je tlak okoln´ıho prostˇred´ı.

Obrázek 4.2: Výpoˇcetn´ı prostor s proudovým motorem a proudˇen´ım okoln´ıho vzduchu [11]

Hmotnostn´ı tok vzduchu na vstupu proudového motoru je ˙m_a, výstupn´ı hmotnostn´ı tok plyn˚u z trysky je ˙m_e, ˙m_f hmotnostn´ı tok paliva pˇricházej´ıc´ıho do spalovac´ıho procesu.

Výtoková rychlost plyn˚u z trysky (exhaust velocity) je ~ue. Tah je oznaˇcován jako FR. Z údaj˚u dokumentace motoru [2] jsou známy maximáln´ı pˇribliˇzná hodnota ˙m_a vstup- n´ıho toku plyn˚u do trysky. Tento údaj je promˇenlivý a nen´ı zcela pˇresný. Dále je známa maximáln´ı teplota výstupn´ıch plyn˚u z trysky tmax. Byl zvolen operaˇcn´ı tah FR, který lze udrˇzovat po dlouhou dobu. Tyto hodnoty jsou povaˇzovány za konstanty a pouˇzij´ı se k dalˇs´ım výpoˇct˚um i simulac´ım.

hmotnostn´ı tok na vstupu m˙a= 1,7 kg·s⁻¹ teplota na v´ystupu t_max = 780^◦C

tah F_R= 680 N

(32)

V dokumentaci nen´ı pˇr´ımo uveden hmotnostn´ı tok paliva, ale tzv. SF C (Specific Fuel Consumption, ˇcesky mˇern´a spotˇreba paliva). Tuto hodnotu je tˇreba pˇrepoˇc´ıtat do poˇzadovan´e formy:

SF C ≤ 1, 07 kg · daN⁻¹· h⁻¹ (4.1)

˙

m_f = SF C · F^R 3600· 10

= 0, 02 kg. · s⁻¹ (4.2)

Je vidˇet, ˇze hmotnostn´ı tok paliva je v tomto pˇr´ıpadˇe témˇeˇr zanedbatelný a oproti hmotnostn´ımu toku vstupn´ıch plyn˚u ˇcin´ı 1-2 %. Pˇresto je správné v tomto kroku spoˇc´ıtat cel- kový výstupn´ı hmotnostn´ı tok plyn˚u z trysky, který je ve zjednoduˇseném tvaru souˇctem hmotnostn´ıch tok˚u na vstupu proudového motoru [11]

˙

m_e = ˙m_a+ ˙m_f (4.3)

Nyn´ı jsou vˇsechny potˇrebné údaje k urˇcen´ı výtokové rychlosti plyn˚u z trysky ~u_e. Dále se pouˇzije vztah pro výpoˇcet tahu proudového motoru [12]:

F~_r= ˙m_e~u_e− ˙m_aV~₀ (4.4) Proto, ˇze se v tomto pˇr´ıpadˇe jedná o pohyb v jednom smˇeru, lze pˇrevést vektory rychlosti na skaláry a dále poˇc´ıtat:

u_e = F_R+ ˙m_aV₀

˙

m_a+ ˙m_f = 396, 2 ms⁻¹ .

= 400 ms⁻¹ (4.5)

4.4 Okrajov´ e podm´ınky 2D modelu

Obrázek 4.3: Názvy okrajových oblast´ı

Okrajové podm´ınky jednotlivých ˇcást´ı hranice jsou znázornˇeny na obrázku 4.3 a jsou shrnuté v tabulce 4.1.

(33)

n´azev hranice rychlost (ve smˇeru x) [ms⁻¹] tlak [kPa] teplota [^◦C]

inlet 70 100 26,85

inlet turbojetOUT 400 100 780

outlet voln´y vektor 100 ∂T /∂~n= 0

outlet borders voln´y vektor 100 ∂T /∂~n= 0

wall elevator 0 ∂p/∂~n = 0 ∂T /∂~n= 0

wall wing 0 ∂p/∂~n = 0 ∂T /∂~n= 0

wall turbojet 0 ∂p/∂~n = 0 ∂T /∂~n= 0

Tabulka 4.1: Z´akladn´ı parametry nastaven´ı okrajov´ych podm´ınek 2D modelu

4.5 Geometrie 3D modelu

Obr´azek 4.4: Geometrie letadla importovan´a do prostˇred´ı Fluent

Typ letadla, který byl zvolen pro tuto práci prezentuje moˇzný stav. Skuteˇcná situace mus´ı m´ıt optimalizovanou polohu proudového motoru. Optimalizace geometrie letadla a um´ıstˇen´ı motor˚u nen´ı souˇcást´ı této práce. Typ letadla také nen´ı naprosto vhodný k podobnému um´ıstˇen´ı pˇr´ıdavných proudových motor˚u. Hlavn´ı c´ıl spoˇc´ıvá ve správném nastaven´ı programu a podm´ınek simulace.

Pˇri optimalizaci geometrie byly odstranˇeny ˇcásti a detaily letadla, které témˇeˇr neo- vlivˇnuj´ı poˇzadované ˇreˇsen´ı. Byla odstranˇena anténa, podvozek letadla, vztlakové klapky.

Um´ıstˇen´ı proudového motoru bylo odhadnuto tak, aby mˇelo význam pro analýzu daného problému. Zjednoduˇsen´ı geometrie má jeden z hlavn´ıch c´ıl˚u, sn´ıˇzit výpoˇcetn´ı nároˇcnost a zlepˇsit konvergenci výpoˇcetn´ıch metod.

(34)

4.6 Okrajov´ e podm´ınky 3D modelu

V tabulce 4.2 jsou popsány simulace A-H, které se liˇs´ı v nastaven´ı kl´ıˇcových vstupn´ıch veliˇcin a to vektoru letové rychlosti, respektive jeho velikosti |~v0|, úhlu nábˇehu α, teploty výstupn´ıch plyn˚u z trysky te a vektoru rychlosti výstupn´ıch plyn˚u z trysky, respektive jeho velikosti |~uê|, viz obrázek 4.5. Tyto veliˇciny definuj´ı o jaký typ pohybu letounu v okoln´ım prostˇred´ı se jedná. Napˇr´ıklad model B definuje let cestovn´ı rychlost´ı s maximáln´ım

´

uhlem nábˇehu. Simulace C má maximáln´ı velikost vstupn´ıho vektoru okoln´ıho prostˇred´ı a pˇredstavuje maximáln´ı rychlost letounu s maximáln´ım úhlem nábˇehu. Pádovou rychlost pˇri nulovém nábˇehovém úhlu definuje simulace F. Simulace G definuje nekritickou situaci, tedy cestovn´ı rychlost a nulový úhel nábˇehu. Kompletn´ı okrajové podm´ınky pro simulaci C jsou v tabulce 4.4 a G v tabulce 4.5.

Obr´azek 4.5: Kl´ıˇcov´e vstupn´ı veliˇciny

oznaˇcen´ı simulace A B C D E F G H

letov´a rychlost |~v⁰| [ms⁻¹] 35 70 83 70 70 35 70 70

´

uhel n´abˇehu α 20^◦ 20^◦ 20^◦ 20^◦ 10^◦ 0^◦ 0^◦ 20^◦

v´ystupn´ı r. toku z trysky |~u⁰| [ms⁻¹] 400 400 400 400 400 400 400 70 v´ystupn´ı teplota z trysky te [^◦C] 780 780 780 570 780 780 780 26,85

Tabulka 4.2: Vstupn´ı parametry pro simulace A-H

B Cestovn´ı rychlost, maximáln´ı úhel C Maximáln´ı rychlost, maximáln´ı úhel F Pádová rychlost, nulový úhel

G Cestovn´ı rychlost, nulový úhel Tabulka 4.3: Názvy jednotlivých simulac´ı

(35)

Na obrázku 4.6 jsou znázornˇeny oblasti okrajových podm´ınek. Povrch letadla je zde symbolizován oblastmi typu ”wall”, plat´ı pro nˇe Dirichletova podm´ınka (2.25) pro rychlost, teplotu a nulová Neumanova podm´ınka pro tlak (^∂p_∂~_n = 0) (2.24). Výstup proudˇen´ı z prostoru je oˇsetˇren oblastmi typu ”outlet”na kterých plat´ı Dirichletova podm´ınka pro tlak. Zvláˇstn´ı oblast je ”outlet symetry”která vytváˇr´ı symetrické proudˇen´ı. Pro tlak na této oblasti plat´ı nulová Neumanova podm´ınka ^∂p_∂n = 0. Pro vektor rychlosti zde plat´ı

~

v· ~n = ~0.

Obrázek 4.6: Grafické znázornˇen´ı hranic výpoˇcetn´ı oblasti

(36)

oblast ~v [ms⁻¹] p [kPa] t [^◦C]

inlet front (0 28,4 -78) ∂p/∂~n= 0 26,85

inlet bottom (0 28,4 -78) ∂p/∂~n= 0 26,85

inlet turbojetOUT (0 0 -400) ∂p/∂~n= 0 780

outlet turbojet hm. tok ˙m_a=1,7 kg s⁻¹ ∂p/∂~n= 0 ∂T /∂~n= 0

outlet symetry voln´y vektor 100 ∂T /∂~n= 0

outlet side voln´y vektor 100 ∂T /∂~n= 0

outlet upper voln´y vektor 100 ∂T /∂~n= 0

outlet back voln´y vektor 100 ∂T /∂~n= 0

wall airplane (0 0 0) ∂p/∂~n= 0 ∂T /∂~n= 0

wall elevator (0 0 0) ∂p/∂~n= 0 ∂T /∂~n= 0

Tabulka 4.4: Okrajov´e podm´ınky simulace C

oblast ~v [ms⁻¹] p[kPa] t [^◦C]

inlet front (0 0 -70) ∂p/∂~n= 0 26,85

inlet bottom (outlet) voln´y vektor 100 ∂T /∂~n= 0 inlet turbojetOUT (0 0 -400) ∂p/∂~n= 0 780

outlet turbojet hm. tok ˙m_a=1,7 kg s⁻¹ ∂p/∂~n= 0 ∂T /∂~n= 0

outlet symetry voln´y vektor 100 ∂T /∂~n= 0

outlet side voln´y vektor 100 ∂T /∂~n= 0

outlet upper voln´y vektor 100 ∂T /∂~n= 0

outlet back voln´y vektor 100 ∂T /∂~n= 0

wall airplane (0 0 0) ∂p/∂~n= 0 ∂T /∂~n= 0

wall elevator (0 0 0) ∂p/∂~n= 0 ∂T /∂~n= 0

Tabulka 4.5: Okrajov´e podm´ınky simulace G

(37)

5 V´ ysledky

Pro 2D model probˇehla simulace na 200 iterac´ı. Vstupn´ı podm´ınky okoln´ıho prostˇred´ı jsou nastaveny na 26,85^◦C a |~v₀| =70 ms⁻¹. Vstupn´ı podm´ınky trysky jsou nastaveny na te =780^◦C a |~ue| =400 ms⁻¹.

5.1 2D simulace

5.1.1 Rychlostn´ı pole

Obr´azek 5.1: Rychlostn´ı pole 2D modelu

Z obrázk˚u 5.1 je vidˇet, ˇze proud z motoru se v této konfiguraci pˇripoj´ı ke spodn´ımu povrchu výˇskového kormidla.

5.1.2 Teplotn´ı pole

V grafu 5.2 nahoˇre je barevnˇe znázornˇena teplota v jednotlivých bodech na povrchu výˇskového kormidla. Graf5.2dole je tvoˇren ”rozvinem”plochy výˇskového kormidla a uka- zuje pr˚ubˇeh teploty na povrchu. Poˇcátek je zvolen na konci klapky a smˇer rozvinu je v obou grafech barevnˇe znázornˇen. Maxima teploty se dosahuje u spodn´ıho okraje výˇskového kormidla. Maximáln´ı teplota je tmax = 281^◦C.

(38)

Obrázek 5.2: Grafy pr˚ubˇehu teplot na povrchu ˇrezu výˇskového kormidla

Obr´azek 5.3: Teplotn´ı pole 2D modelu

(39)

Tato simulace operuje s dvourozmˇerným modelem, který nedokáˇze spolehlivˇe vypoˇc´ıtat chován´ı turbulentn´ıho proudˇen´ı. Hlavn´ım problémem je absence tˇret´ı dimenze. Výsledky 2D simulace mohou ukazovat správné chován´ı pouze jednoduché geometrie, která je ve smˇeru z nemˇenná. Napˇr´ıklad kruh ve dvourozmˇerné geometrii znamená ve skuteˇcnosti ˇreˇsen´ı obtékán´ı dostateˇcnˇe vysokého válce. Z toho plyne, ˇze v této simulaci napˇr´ıklad nen´ı zahrnut rotaˇcnˇe symetrický tvar proudového motoru, nebo rozd´ılná geometrie kˇr´ıdel v r˚uzných ˇrezech i pˇr´ıtomnost trupu letadla, který bude hrát nemalou roli v charak- teru proudového pole. Tento výsledek je proto tˇreba brát jako zkuˇsebn´ı, nedemonstruj´ıc´ı reálnou situaci.

5.2 3D simulace

Pro správný popis daného problému je nutné poˇc´ıtat se 3D geometri´ı. V praxi se pˇri výpo- ˇctech ˇcasto vyuˇz´ıvá symetrie, pokud to jde. V simulaci se poˇc´ıtá s vstupn´ım rychlostn´ım polem rovnobˇeˇzným s podélnou osou letadla a tak lze rozdˇelit celou geometrii podle roviny symetrie. Na polovinˇe geometrie se sn´ıˇz´ı poˇcet stupˇn˚u volnosti 2x, ˇc´ımˇz se výpoˇcetn´ı nároˇcnost implicitn´ıho iteraˇcn´ıho ˇreˇsiˇce sn´ıˇz´ı (2^p)krát, kde p bude v intervalu 1-2.

Letoun Zl´ın Z-50 je projektován tak, ˇze pˇri letu s n´ızkým úhlem nábˇehu bude výˇskové kormidlo proudem z trysky motoru velmi málo ovlivnˇeno. K významnˇejˇs´ımu ovlivnˇen´ı dojde pˇri vyˇsˇs´ıch hodnotách úhlu nábˇehu. Je tˇreba dbát na to, ˇze kritický úhel nábˇehu je pˇredevˇs´ım v závislosti na geometrii pˇribliˇznˇe 20^◦. Tuto skuteˇcnost jsem ovˇeˇril simulac´ı, nastaven´ım vstupn´ıch podm´ınek na nábˇehový úhel 25^◦, který by mˇel být za hodnotou kritického úhlu nábˇehu.

Obrázek 5.4: Proudnice v rovinˇe R1 viz obr. 5.5. Úhel nábˇehu 25^◦, letová rychlost 77ms⁻¹

(40)

Z obrázku5.4vid´ıme odtrˇzen´ı proudˇen´ı v bl´ızkosti nábˇeˇzné hrany s masivn´ı separaˇcn´ı zónou nad kˇr´ıdlem. V tomto okamˇziku kˇr´ıdlo ztrác´ı vztlak. Jedná se o kritický stav, kdy letadlo zaˇcne padat.

Byly vytvoˇreny grafické výstupy v nˇekolika charakteristických rovinách (ˇrezech), viz obrázky 5.5. V tˇechto rovinách jsou vykreslovány teplotn´ı pole, rychlostn´ı pole, pole tlak˚u a pole hustot a také proudnice rychlosti.

(a) Rovina R1 (b) Rovina R2

(c) Rovina R3 (d) Rovina R4

Obr´azek 5.5: Charakteristick´e roviny (ˇrezy)

5.3 Vliv v´ ypoˇ cetn´ı s´ıtˇ e

Pˇri volbˇe r˚uznˇe husté výpoˇcetn´ı s´ıtˇe lze dosáhnout odliˇsných výsledk˚u. Obecnˇe plat´ı, ˇze ˇc´ım vˇetˇs´ı poˇcet prvk˚u ve výpoˇcetn´ı s´ıti je, t´ım je simulace pˇresnˇejˇs´ı. V praxi se snaˇz´ıme vytváˇret prvk˚u co nejménˇe, ale aby výsledek byl stále správný. Hustota výpoˇcetn´ı s´ıtˇe pˇr´ımo ovlivˇnuje poˇcet operac´ı pˇri simulaci a celkový ˇcas.

Obrázek5.6 zahrnuje simulaci B s výpoˇcetn´ımi s´ıtˇemi o r˚uzných poˇctech prvk˚u. Jsou tam vykresleny proudnice pro s´ıtˇe o hustotˇe 50, 360 a 900 tis´ıc element˚u. Je zˇretelnˇe vidˇet, ˇze v prvn´ım modelu s 50 000 prvky simulace vypoˇc´ıtala odtrˇzen´ı proudu za kˇr´ıdlem.

Simulace s výpoˇcetn´ı s´ıt´ı o 360 000 prvc´ıch jiˇz odtrˇzen´ı neurˇcila. Výsledek simulace s nejvˇetˇs´ı výpoˇcetn´ı s´ıt´ı 900 000 element˚u témˇeˇr nepˇrinesl zmˇenu, pouze vˇetˇs´ı výpoˇcetn´ı nároˇcnost. Pro úˇcely této práce postaˇc´ı 3D model se s´ıt´ı o poˇctu 360 000 prvk˚u.

(41)

Obr´azek 5.6: Porovn´an´ı proudnic model˚u o (shora) 50, 360 a 900 tis´ıc´ıch prvk˚u