Numerické modelování problémů mechaniky kontinua s pomocí výpočetního balíku FEniCS

Numerické modelování problémů mechaniky kontinua s pomocí

výpočetního balíku FEniCS

Bakalářská práce

Studijní program: B2612 Elektrotechnika a informatika Studijní obor: Elektronické informační a řídicí systémy

Autor práce: Martin Votýpka

Vedoucí práce: doc. Ing. Petr Šidlof, Ph.D.

Ústav nových technologií a aplikované informatiky

Zadání bakalářské práce

Numerické modelování problémů mechaniky kontinua s pomocí výpočetního balíku FEniCS

Jméno a příjmení: Martin Votýpka

Osobní číslo: M17000065

Studijní program: B2612 Elektrotechnika a informatika Studijní obor: Elektronické informační a řídicí systémy

Zadávající katedra: Ústav nových technologií a aplikované informatiky Akademický rok: 2019/2020

Zásady pro vypracování:

1. Seznamte se se základy mechaniky elastických těles.

2. Nastudujte základy dynamiky tekutin, zejména z oblasti externího proudění (obtékání těles).

3. Seznamte se se základními principy metody konečných prvků, osvojte si jednoduché koncepty programování v C++ nebo Python.

4. Na základě benchmarkového problému specifikovaného v publikaci [4] připravte geometrii a výpočetní síť

5. S pomocí volně dostupné knihovny FEniCS realizujte zvlášť numerickou simulaci deformace elastického nosníku a simulaci obtékání tuhého tělesa ve 2D.

6. Navrhněte možné přístupy, jak do budoucna ve FEniCS sestavit numerickou simulaci sdruženého problému interakce proudění a pružného tělesa.

Rozsah grafických prací: dle potřeby Rozsah pracovní zprávy: 30-40 stran

Forma zpracování práce: tištěná/elektronická

Jazyk práce: Čeština

Seznam odborné literatury:

[1] White F. M. (2006), Fluid Mechanics, McGraw-Hill.

[2] Brennen C. E. (2006), Internet Book on Fluid Dynamics, Caltech, http://brennen.caltech.edu/fluidbook/content.htm (Online).

[3] Schäfer M. (2006), Computational Engineering Introduction to Numerical Methods, Springer- Verlag.

[4] Turek S., Hron J. (2006), Proposal for Numerical Benchmarking of Fluid-Structure Interaction between an Elastic Object and Laminar Incompressible Flow, in: Fluid-Structure Interaction, pp.

371-385, DOI 10.1007/3-540-34596-5_15.

[5] FEniCS Project, https://fenicsproject.org/ (Online).

Vedoucí práce: doc. Ing. Petr Šidlof, Ph.D.

Ústav nových technologií a aplikované informatiky

Datum zadání práce: 9.října 2019 Předpokládaný termín odevzdání: 18.května 2020

L.S.

prof. Ing. Zdeněk Plíva, Ph.D. Ing. Josef Novák, Ph.D.

děkan vedoucí ústavu

Prohlášení

Prohlašuji, že svou bakalářskou práci jsem vypracoval samostatně jako původní dílo s použitím uvedené literatury a na základě konzultací s vedoucím mé bakalářské práce a konzultantem.

Jsem si vědom toho, že na mou bakalářskou práci se plně vztahuje zákon č. 121/2000 Sb., o právu autorském, zejména § 60 – školní dílo.

Beru na vědomí, že Technická univerzita v Liberci nezasahuje do mých autorských práv užitím mé bakalářské práce pro vnitřní potřebu Technické univerzity v Liberci.

Užiji-li bakalářskou práci nebo poskytnu-li licenci k jejímu využití, jsem si vědom povinnosti informovat o této skutečnosti Technickou univerzitu v Liberci; v tomto případě má Technická univerzita v Liberci právo ode mne požadovat úhradu nákladů, které vynaložila na vytvoření díla, až do jejich skutečné výše.

Současně čestně prohlašuji, že text elektronické podoby práce vložený do IS/STAG se shoduje s textem tištěné podoby práce.

Beru na vědomí, že má bakalářská práce bude zveřejněna Technickou univerzitou v Liberci v souladu s § 47b zákona č. 111/1998 Sb., o vysokých školách a o změně a doplnění dalších zákonů (zákon o vysokých školách), ve znění pozdějších předpisů.

Jsem si vědom následků, které podle zákona o vysokých školách mohou vyplývat z porušení tohoto prohlášení.

1. června 2020 Martin Votýpka

Abstrakt

Hlavním cílem této bakalářské práce je numerická simulace obtékání tělesa tekutinou a simulace deformace elastických těles v prostředí FEniCS. Dlouhodobým cílem je řešení problému interakce proudění s pružnými tělesy. V první části této práce je uveden velmi zjednodušený popis metody konečných prvků a také představení výpočetního balíku FEniCS, který je použit pro numerické řešení parciálních diferenciálních rovnic. Práce se dále zabývá numerickou simulací deformací spojitě zatíženého nosníku a numerickou simulací obtékání tělesa tekutinou. V případě ohybu spojitě zatíženého nosníku jsou výsledky simulací ověřeny pomocí analyticky získaného řešení, v případě simulace obtékaní tělesa tekutinou jsou výsledky ověřeny pomocí benchmarkových dat. V poslední části je uveden možný postup při řešení problému interakce proudění s pružnými tělesy.

Klíčová slova: metoda konečných prvků, FEniCS, numerické simulace, deformace pružných těles, Navierovy-Stokesovy rovnice, interakce proudění a pružného tělesa

Abstract

The main goal of this bachelor's thesis is numerical simulation of fluid flow around a body and a simulation of deformation of elastic bodies in the FEniCS environment. The long-term goal is to solve the problem of fluid-structure interaction. The first part of this work presents a very simplified description of the finite element method and also the introduction of the FEniCS libraries, which are used for the numerical solution of partial differential equations. This thesis further deals with the numerical simulation of deformations of a continuously loaded beam and numerical simulation of the fluid flow around a body. In the case of bending of a continuously loaded beam, the results of simulations are verified using an analytically obtained solution, in the case of simulating the fluid flow around a body, the results of simulations are verified using benchmark data. The last part presents a possible method for solving the problem of fluid-structure interaction.

Keywords: finite element method, FEniCS, numerical simulations, deformation of elastic bodies, Navier-Stokes equations, fluid-structure interaction

Obsah

Seznam použitých zkratek a symbolů...9

Úvod...10

1 Metoda konečných prvků...11

1.1 FEniCS Project...13

2 Numerická simulace ohybu spojitě zatíženého nosníku...14

2.1 Rovnice popisující deformace pružného tělesa...14

2.2 Parametry výpočtu a simulace...16

2.2.1. Geometrie nosníku a zvolené konstanty...16

2.2.2. Výpočetní sítě...16

2.3 Analytický výpočet deformace nosníku...17

2.4 Implementace ve FEniCSu...18

2.5 Výsledky numerických simulací...22

3 Numerická simulace obtékání tělesa tekutinou...23

3.1 Navierovy-Stokesovy rovnice...23

3.2 Parametry simulace...26

3.2.1. Geometrie oblasti...26

3.2.2. Parametry proudící tekutiny...27

3.2.3. Okrajové podmínky...27

3.2.4. Výpočetní sítě...28

3.3 Implementace ve FEniCSu...29

3.4 Výsledky numerických simulací...37

3.4.1. Výsledky simulace při použití hrubší výpočetní sítě...37

3.4.2. Výsledky simulace při použití jemnější výpočetní sítě...40

4 Interakce proudění a pružného tělesa...45

Seznam obrázků

Obrázek 1: nestrukturovaná síť (858 prvků)...15

Obrázek 2: strukturovaná síť (712 prvků)...15

Obrázek 3: průhyb spojitě zatíženého nosníku...21

Obrázek 4: geometrie oblasti...25

Obrázek 5: hrubá výpočetní síť (3794 prvků)...27

Obrázek 6: jemná výpočetní síť (30850 prvků)...27

Obrázek 7: jemná síť v okolí obtékaného tělesa...27

Obrázek 8: rychlost tekutiny obtékající těleso v čase 9 sekund...36

Obrázek 9: tlak tekutiny obtékající těleso v čase 9 sekund...36

Obrázek 10: graf aerodynamické vztlakové síly v čase 0 až 10 sekund...37

Obrázek 11: graf aerodynamické vztlakové síly v čase 9 až 9,6 sekund...37

Obrázek 12: graf aerodynamické odporové síly v čase 0 až 10 sekund...38

Obrázek 13: graf aerodynamické odporové síly v čase 9 až 9,6 sekund...38

Obrázek 14: rychlost tekutiny obtékající těleso v čase 9 sekund...39

Obrázek 15: tlak tekutiny obtékající těleso v čase 9 sekund...39

Obrázek 16: graf aerodynamické vztlakové síly v čase 0 až 10 sekund...40

Obrázek 17: graf aerodynamické vztlakové síly v čase 9 až 9,6 sekund...40

Obrázek 18: graf aerodynamické odporové síly v čase 0 až 10 sekund...41

Obrázek 19: graf aerodynamické odporové síly v čase 9 až 9,6 sekund...41 Obrázek 20: srovnání výsledků aerodynamické odporové síly z benchmarku [1] a ze simulací.42 Obrázek 21: srovnání výsledků aerodynamické vztlakové síly z benchmarku [1] a ze simulací. 42

Seznam použitých zkratek a symbolů

ALE Arbitrary Lagrangian-Eulerian

Δ

Laplaceův operátor

E

Youngův modul

ϵ

symetrický gradient posunutí, tenzor rychlosti deformace

I

jednotková matice

J

moment setrvačnosti

λ

Lamého konstanta

Mo

ohybový moment

μ

smykový modul, dynamická viskozita

n

vnější jednotková normála

∇

nabla operátor

ν

Poissonovo číslo, kinematická viskozita

p

tlak

R e

Reynoldsovo číslo

ρ

^hustota

σ

tensor napětí

t

čas

tr

stopa matice

u vektorové pole posunutí, vektor rychlosti

v

testovací funkce

Úvod

Proudění tekutin i deformace pružných těles jsou témata s širokým spektrem uplatnění. S aplikací proudění tekutin se setkáváme například v meteorologii, při testech aerodynamických vlastností těles, například letadel, nebo při proudění tekutin v potrubích. S pružnými deformacemi se setkáváme u každého stroje. Často dochází také ke sdruženému problému, kdy v důsledku proudění tekutiny dochází k deformaci těles a tyto vzniklé deformace následně ovlivňují proudící tekutinu.

Cílem této práce je numerická simulace pružné deformace spojitě zatíženého nosníku a simulace obtékání tělesa tekutinou. Dlouhodobým cílem, tedy možným navázáním na tuto práci, je simulace sdruženého problému. K těmto simulacím byl využit výpočetní balík FEniCS, který disponuje všemi potřebnými funkcemi pro řešení parciálních diferenciálních rovnic metodou konečných prvků.

V první části se věnuji metodě konečných prvků. Je zde uveden velmi zjednodušený popis metody a následně je představena open-source knihovna FEniCS, obsahující funkce založené na metodě konečných prvků.

Druhá část je věnována pružným deformacím spojitě zatíženého vetknutého nosníku. Nejprve jsou uvedeny rovnice popisující deformace pružného tělesa a zároveň úpravy rovnic pro případ rovinné napjatosti. Následně je ukázán analytický výpočet a kód implementující simulaci ohybu spojitě zatíženého nosníku. Nakonec jsou srovnány výsledky získané ze simulací při využití různých výpočetních sítí s analyticky získaným výsledkem. Geometrie i vlastnosti materiálu nosníku odpovídají hodnotám uvedeným v benchmarku [1].

Třetí část je věnována simulaci obtékání tělesa nestlačitelnou vazkou tekutinou. Nejprve jsou představeny rovnice popisující proudění nestlačitelné vazké tekutiny a jejich následné úpravy.

Poté jsou uvedeny parametry simulací, které se shodují s parametry v benchmarku, a také výpočetní sítě použité při numerických simulacích. Následně je popsán kód implementující výpočet pomocí knihoven FEniCS. Nakonec je uvedeno srovnání dosažených výsledků s výsledky v benchmarku.

Ve čtvrté a poslední části je uveden možný postup při řešení sdruženého problému, tedy interakce proudění s pružnými tělesy.

1 Metoda konečných prvků

Metoda konečných prvků je numerická metoda pro řešení parciálních diferenciálních rovnic. V současné době nachází využití v širokém spektru aplikací jako je například meteorologie či simulace deformace namáhaných částí strojů. Její princip spočívá v diskretizaci spojité oblasti a převedení diferenciálních rovnic na soustavu lineárních rovnic. Informace použité v této kapitole byly získány z [4] a [8].

Princip metody konečných prvků lze demonstrovat například na Laplaceově rovnici ve 2D

−Δ u=f

^v

Ω

^{, (1)}

s okrajovými podmínkami

u=0

na

∂Ω

. (2)

Nejprve je nutné získat slabou formulaci. Rovnici (1) násobíme libovolnou funkcí

v

s hodnotou rovnou nule na hranicích oblasti. Výsledná rovnice po úpravě je

−Δ uv=f v

^{. (3)}

Dále bude potřeba rovnici (3) zintegrovat

− ∫

Ω

Δ uv dx= ∫

Ω

f v dx

. (4)

Pomocí Greenovy věty

∫

Ω

f ∂ g

∂ x

i

= ∫

∂Ω

f g n

_i

− ∫

Ω

g ∂ f

∂ x

i (5)

je možné dále rovnici upravit do tvaru

∑

i=1

2

( ^∫

^∂Ω

^{v ( − ∂ ∂ x u}

ⁱ

^{) n}

ⁱ

^{− ∫}

^Ω

^{( ∂ x ∂ v}

ⁱ

^{)( − ∂ ∂ x u}

ⁱ

⁾ ) ^{= ∫}

^Ω

^{f v}

^{. (6)}

Jelikož má v tomto případě funkce

v

nulovou hodnotu na hranicích oblasti, je výsledek křivkového integrálu vždy rovný nule a je tedy možné rovnici dále upravit do tvaru

∫

Ω

∇ v⋅∇ u dx= ∫

Ω

f v dx

. (7)

Následně provedeme aproximaci slabé formulace. Budeme hledat po částech lineární funkci

u

_hs hodnotou rovnou nule na hranicích oblasti tak, aby

∫

Ω

∇ u

h

⋅∇ v

h

dx = ∫

Ω

f v

_h

dx

, (8) kde

v

_hjsou po částech lineární funkce splňující okrajovou podmínku.

Funkce

v

_ha

u

_hjsou elementy prostoru

V

_h

={v

_h

: Ω →ℝ }

. (9)

V

_hje prostor s konečnou dimenzí, jehož bázi lze zapsat jako

{φ

1

, ... , φ

n

}

, (10)

kde

n

je dimenze prostoru

V

_h(pro lineární prvky rovný počtu vnitřních uzlů výpočetní sítě).

Funkci

u

_hje možné vyjádřit jako

∑

j=1 n

α

j

φ

j , (11)

kde

α

1

,... , α

njsou neznámé.

Výsledkem této aproximace je soustava lineárních rovnic

∑

j=1 n

α

j

a

_ij

=b

_i ,

i=1 ,..., n

, kde

a

_ij

= ∫

Ω

∇ φ

j

⋅∇ φ

i

dx

,

b

_i

= ∫

Ω

f φ

i

dx

.

(12)

Soustavu je možné také zapsat maticově jako

( ^{a a ⋮ ...}

^{n 111}

^{... a ... a ⋮}

^{1 nnn}

) ^{( α α ⋮}

¹ⁿ

^{) = ( b b ⋮}

¹ⁿ

⁾

^{. (13)}

1.1 FEniCS Project

FEniCS je výpočetní balík pro řešení parciálních diferenciálních rovnic, který uživatelům umožňuje rychlé převedení rovnic do kódu konečných prvků. Je možné jej použít jak na běžných pracovních počítačích, tak na výkoných strojích využívajících paralelního výpočtu.

FEniCS byl původně vytvořen v roce 2003 a je vyvíjen za účasti univerzit a výzkumných ústavů po celém světě. Bližší informace ohledně knihoven FEniCS, jejich instalaci a použití lze nalézt v [4].

Přímo lze FEniCS instalovat pouze v operačním systému Linux. Pokud má být FEniCS nainstalován na zařízení s jiným operačním systémem, je nutné využít dalších programů, jako je například Docker, který instalaci a následné používání knihovny FEniCS umožní.

Po úspěšné instalaci balíku FEniCS je možné přistoupit k tvorbě zdrojového kódu. Ten je možné psát buď v programovacím jazyku C++ a nebo v jazyku Python. Osobně jsem pro tvorbu zdrojového kódu zvolil Python. Učinil jsem tak hlavně z důvodu většího množství materiálů, ze kterých bylo možné čerpat potřebné informace.

2 Numerická simulace ohybu spojitě zatíženého nosníku

2.1 Rovnice popisující deformace pružného tělesa

V této kapitole jsou uvedeny rovnice popisující deformace pružného tělesa a úprava rovnic pro případ rovinné napjatosti. Informace byly čerpány z [6] a [4].

Pružné těleso je takové těleso, které se po odstranění vnějších působících sil navrátí do původního stavu. Jeho deformace jsou popsány pomocí rovnice

−∇⋅ σ =f

, (14)

kde

f

je síla působící na jednotku objemu deformovaného tělesa a

σ

je tenzor napětí. Jelikož se zde jedná o případ přímého tenkého nosníku, je možné zanedbat kolmou složku napětí a uvažovat případ rovinné napjatosti, čímž dojde k zjednodušení výsledných rovnic. Pro případ rovinné napjatosti lze tenzor napětí

σ

zapsat jako

σ = [ ^{σ σ 0}

¹¹²¹

^{σ σ 0}

¹²²²

^{0 0 0} ]

^{. (15)}

Tento tenzor je dán vztahem

σ =λ⋅tr (ϵ )⋅I+2μ⋅ϵ

^, (16)

kde symbol

λ

a

μ

jsou Lamého konstanty definované jako

λ = E ⋅ ν

(1+ν )(1−2ν )

^a

μ= E

2 (1+ν )

^{, (17)}

symbol

I

ve vztahu (16) je jednotková matice a

ϵ

je symetrický gradient vektorového pole posunutí definovaný pro případ rovinné napjatosti jako

ϵ= [ ^{ϵ ϵ 0}

¹¹²¹

^{ϵ ϵ 0}

¹²²²

^{ϵ 0 0}

³³

]

^{. (18)}

Symetrický gradient posunutí je dán vztahem

ϵ (u)= 1

2 (∇ u+(∇ u)

^T

)

, (19)

kdeuje vektorové pole posunutí.

Celý problém může být následně zredukován pouze na dva rozměry. Nelze však ignorovat vliv parametru

ϵ

₃₃. Aby byla redukce možná, je nutné rovnici (16) upravit.

[ ^{σ σ 0}

¹¹²¹

^{σ σ 0}

¹²²²

^{0 0 0} ] ^=λ⋅(ϵ

¹¹

^+ϵ

²²

^+ϵ

³³

^{)⋅ [ 1 0 0 0 1 0 0 0 1 ] +2 μ⋅} [ ^{ϵ ϵ 0}

¹¹²¹

^{ϵ ϵ 0}

¹²²²

^{ϵ 0 0}

³³

]

λ⋅(ϵ

11

+ϵ

22

+ϵ

33

)+2 μ⋅ϵ

33

=0 ϵ

33

=− λ λ +2μ (ϵ

¹¹

+ϵ

22

)

(20)

Dosazením výsledku úpravy (20) do rovnice (16) lze následně získat

σ =λ⋅(ϵ

11

+ϵ

22

− λ λ +2μ (ϵ

¹¹

+ϵ

22

))⋅I+2 μ⋅ϵ

σ =λ⋅((1− λ λ +2 μ )(ϵ

¹¹

+ϵ

22

))⋅I+2 μ⋅ϵ

σ = ( ^{λ +2μ 2 μ⋅λ} ) ^(ϵ

¹¹

^+ϵ

²²

^{)⋅I+2 μ⋅ϵ}

^.

(21)

Při označení

λ

^∗

= ( ^{λ +2 μ 2 μ⋅λ} )

⁽²²⁾

je výsledná rovnice pro dvourozměrný problém

σ =λ

^∗

⋅tr (ϵ )⋅I+2μ⋅ϵ

^{. (23)}

2.2 Parametry výpočtu a simulace

2.2.1.Geometrie nosníku a zvolené konstanty

Parametry simulace ohybu spojitě zatíženého nosníku vycházejí z parametrů daných v benchmarku [1]. Jedná se o nosník s délkou 35 centimetrů, výškou 2 centimetry, vyrobený z polypropylenu.

Tabulka 1: hodnoty parametrů výpočtu a simulace ohybu spojitě zatíženého nosníku

Délka

L=0 ,35

m

Výška

H =0 ,02

m

Youngův modul

E=9⋅10

⁸

kg

ms

²

Hustota

ρ =1100 kg

m

³

Poissonovo číslo

ν =0,42

2.2.2.Výpočetní sítě

Při vlastních simulacích pružných deformací nosníku jsem využil dvě výpočetní sítě, každou s jiným uspořádáním prvků.

První síť použitá pro simulace byla nestrukturovaná síť, tedy síť s nepravidelným rozložením jednotlivých prvků. Výslednou nestrukturovanou síť je možné vidět na obrázku 1.

Druhá použitá síť byla strukturovaná, tedy síť s pravidelným uspořádáním prvků, kterou je možné vidět na obrázku 2.

Různé výpočetní sítě pro simulaci ohybu nosníku byly použity pro následné porovnání obdržených výsledků.

Obrázek 2: strukturovaná síť (712 prvků)

Obrázek 1: nestrukturovaná síť (858 prvků)

2.3 Analytický výpočet deformace nosníku

Při výpočtu průhybu spojitě zatíženého nosníku budu vycházet ze vzorce

w

^{' '}

(x)= Mo (x)

E ⋅Jy

^{, (24)}

kde

w (x)

je průhybová křivka,

Mo (x)

je ohybový moment,

E

je Youngův modul pružnosti a

Jy

je moment setrvačnosti vzhledem k ose y definovaný jako

Jy = ∫

₋^b

2 b 2

∫

₋^h

2 h

2

y

²

dy dz= ∫

₋^b

2 b 2

h

³

12 dz= b ⋅h

³

12

^{, (25)}

kde

b

je šířka nosníku a

h

je výška nosníku.

Ohybový moment

Mo(x)

je definovaný pro spojitě zatížený nosník

Mo (x)=m⋅g(x− x

²

2 L − L

2 )

^{, (26)}

kde

m ⋅g

je působící zatížení a

L

je délka nosníku.

Rovnici (24) je následně možné dvakrát zintegrovat, čímž se získá rovnice

w (x)= m ⋅g

E ⋅Jy ( ^{−x 24 L}

⁴

^{+ x 6}

³

^{− L ⋅x 4}

²

)

^{. (27)}

Při výpočtu maximálního průhybu, tedy průhybu na volném konci nosníku, je možné rovnici (27) dále upravit do tvaru

L ⋅h⋅b⋅ρ⋅g E⋅ 1

12 b ⋅h

³

⋅ ( ^{− 24 3 L}

³

) ^{=− 3 2 ρ⋅g⋅L E} ⋅h

^{24 . (28)}

Po dosazení zadaných parametrů je výsledný průhyb na volném konci nosníku

w (0 ,35)≈−0 ,000674719

^m. (29)

Vzorec (24) a jeho následné úpravy byly získány z [7].

2.4 Implementace ve FEniCSu

V této kapitole bude popsán kód implementující výpočet pod vlastní vahou deformovaného pružného nosníku. Velká část kódu byla převzata z [6]. Informace ohledně použitých funkcí byly čerpány z [3] a [6]. Parametry použité při simulaci lze nalézt v kapitole 2.2.

Aby bylo možné pro výpočet diferenciálních rovnic použít funkce FEniCSu, je nejprve nutné knihovnu FEniCS importovat.

from fenics import *

Knihovna fenics obsahuje všechny potřebné funkce pro výpočet parciálních diferenciálních rovnic, založených na metodě konečných prvků. Po jejím importování již bude možné potřebné funkce využít.

Dále bude potřeba nadefinovat konstanty použité při výpočtu.

L = 0.35 H = 0.02

E = Constant(9e8) nu = Constant(0.42) rho = 1100

g = 9.81 rho_g = rho*g

f = Constant((0,-rho_g)) mu = E/2/(1+nu)

lmbda = E*nu/(1+nu)/(1-2*nu) lmbda = 2*mu*lmbda/(lmbda+2*mu)

Proměnná L je délka nosníku, H je výška nosníku, konstanta E je Youngův modul v pružnosti a nu je Poissonovo číslo. Proměnná rho je hustota, g je gravitační zrychlení a konstanta f je zatížení nosníku. mu a lmbda jsou Lamého konstanty odpovídající předpisu (17). Z odvození v kapitole 2.1 je zřejmé, že aby mohla být provedena redukce problému čistě na dva rozměry, je nutné rovnici upravit do tvaru (23). Výsledná hodnota proměnné lmbda odpovídá hodnotě koeficientu

λ

^∗.

Následně je třeba načíst, případně vytvořit výpočetní síť. Při načítání již vytvořené sítě je nutné, aby síť byla uložena ve formátu xml. Pokud je síť uložena v jiném formátu, je nutné nejprve provést konverzi. K tomu je možné využít nástroj dolfin-convertr, který je implementován v knihovnách FEniCS. Pokud síť není předem vytvořena, je možné ji vytvořit přímo ve FEniCS.

Tímto způsobem lze vytvořit jednodušší sítě. Aby mohla být síť vytvořena pomocí knihoven FEniCS, je nejprve nutné importovat komponentu mshr, která potřebné funkce pro vytvoření sítě obsahuje.

mesh = Mesh("beam.xml")

Ve chvíli, kdy je načtena, případně vytvořena výpočetní síť, je možné přistoupit k definici funkčního prostoru a následně bázové a testovací funkce.

V = VectorFunctionSpace(mesh, 'Lagrange', degree=2) u_ = TrialFunction(V)

v_ = TestFunction(V)

Funkce VectorFunctionSpace přebírá tři vstupní parametry. První parametr je výpočetní síť, druhý parametr značí typ prvku a třetí parametr je stupeň prvku. Proměnná V je tedy vektorový funkční prostor s kvadratickými prvky typu Lagrange. Funkce TrialFunction vytvoří bázovou funkci daného funkčního prostoru. Funkce TestFunction vytvoří testovací funkci daného funkčního prostoru.

Aby byl program snáze čitelný je dobré předem nadefinovat dlouhé a často používané výrazy. Z tohoto důvodu je výpočet symetrického gradientu posunutí (19) a tenzoru napětí (16) realizován jako funkce.

def epsilon(u):

return 0.5*(nabla_grad(u)+nabla_grad(u).T)

Funkce nabla_grad je jednou z funkcí obsažených v knihovnách FEniCS. Vstupním parametrem je vektor, výstupním pak jeho gradient. Parametr T značí transponovanou matici gradientu vstupního vektoru.

def sigma(u):

return lmbda*tr(epsilon(u))*Identity(2) + 2.0*mu*epsilon(u)

lmbda je již definovaná proměnná odpovídající hodnotě koeficientu (22). Funkce tr vrací stopu vstupní čtvercové matice, tedy součet prvků na hlavní diagonále. Funkce Identity vytvoří jednotkovou matici o velikosti závislé na vstupním parametru. V tomto případě se jedná o matici 2x2. Funkce epsilon je funkce pro výpočet tenzoru rychlosti deformace (16) popsaná výše.

Dále bude potřeba definovat hranice oblasti po pozdější definici okrajových podmínek. V případě deformace vetknutého nosníku bude definován jeden okraj v místě vetknutí.

def left(x, on_boundary):

return near(x[0],0.)

Funkce near vrací booleovskou hodnotu udávající, zda se hodnota prvního parametru dostatečně blíží hodnotě druhého. Parametr x[0] značí souřadnice na ose x.

Po definování okraje je možné definovat okrajovou podmínku.

bc = DirichletBC(V, Constant((0.,0.)), left)

Funkce DirichletBC je funkce pro definici Dirichletovy okrajové podmínky. Prvním argumentem je funkční prostor, druhý argument je popis chování funkce na okraji a třetí argument je okraj.

Po definování všech potřebných konstant, vztahů a okrajových podmínek je možné přistoupit k výpočtu. Nejprve je nutné získat slabou formulaci. Tu lze získat vynásobením rovnice (14) testovací funkcí, integrací a následnou aplikací Greenovy věty a okrajové podmínky. Rovnice po této úpravě je ve tvaru

∫

Ω

σ (u)⋅∇ v dx= ∫

Ω

f v dx

, (30) kde

σ (u)

je tenzor napětí,

v

je testovací funkce, a

f

je zatížení nosníku. Jelikož je tenzor

σ (u)

symetrický, tak po skalárním součinu s tenzorem

∇ v

budou mít vliv na výsledek pouze symetrické části tenzoru

∇ v

. Rovnici (30) je tedy dále možné upravit do tvaru

∫

Ω

σ (u)⋅ϵ (v)dx= ∫

Ω

f v dx

, (31) kde

ϵ (v)

je symetrická část gradientu funkce

v

. Implementace této rovnice ve FEniCSu vypadá následovně:

a = inner(sigma(u_), epsilon(v_))*dx l = dot(f, v_)*dx

Aby mohla být rovnice vyřešena, musí nejdříve být rozdělena na pravou a levou stranu.

Proměnná a značí levou stranu rovnice (31), proměnná l pravou stranu rovnice (31). Funkce inner je funkce pro skalární součin tenzorů, funkce sigma a epsilon jsou již definované funkce popsané výše. Proměnná u_ je bázová funkce, proměnná v_ je testovací funkce. Parametr dx je parametr značící integraci. Funkce dot je funkcí pro skalární součin vektorů a proměnná f je zatížení nosníku.

Následně je již možné pomocí funkcí knihovny fenics rovnici vyřešit.

u = Function(V, name="Displacement") solve(a == l, u, bc)

Proměnná u je funkce, do které bude uložen výsledek simulace. function vytvoří funkci v závislosti na vstupním parametru. Vstupní parametr může být funkční prostor, nebo jiná funkce. Pokud je funkce tvořena z jiné funkce, je možné přidat další parametr určující číslo dílčí funkce, která má být extrahována. Parametr name je název funkce, který umožní snazší práci při zpracování výsledků. K řešení rovnice je využita funkce solve. Tato funkce má celkem čtyři různé způsoby použití. V tomto případě bude řešena rovnice

a =l

, kde a je levá strana rovnice (31) a l pravá strana. Funkce u bude výsledkem výpočtu a bc jsou okrajové podmínky.

Po dokončení výpočtu bude možné výsledky uložit a dále zpracovat. Aby mohly být uloženy hodnoty napětí je nejprve nutné promítnout tenzor napětí do vhodného funkčního prostoru.

Vsig = TensorFunctionSpace(mesh, "DG", degree=0) sig = Function(Vsig, name="Stress")

sig.assign(project(sigma(u), Vsig))

Funkce TensorFunctionSpace je funkce pro vytvoření tenzorového funkčního prostoru. První vstupní parametr funkce je výpočetní síť, druhý parametr je typ prvků, třetí parametr je stupeň prvků. Funkce project vrací projekci funkce či předpisu v prvním parametru do prostoru konečných prvků v parametru druhém. Výsledek projekce je nakonec přiřazen do proměnné sig pomocí funkce assign.

Výsledky jsou nakonec pro další zpracování uloženy.

file_results = XDMFFile("elasticity_results.xdmf") file_results.parameters["flush_output"] = True

file_results.parameters["functions_share_mesh"] = True file_results.write(u, 0.)

file_results.write(sig, 0.)

2.5 Výsledky numerických simulací

Při simulacích průhybu nosníku byly použity dvě sítě s odlišným uspořádáním prvků. Výsledky získané při použití jednotlivých sítí se od sebe mírně liší. Výsledné hodnoty průhybu nosníku získané při simulacích a z analytického výpočtu je možné najít v tabulce 2. Grafické zpracování výsledků simulací je možné vidět na obrázku 3.

Tabulka 2: Výsledné hodnoty průhybu na volném konci nosníku

Nestrukturovaná síť: 0,0006752933 m

Strukturovaná síť: 0,0006750631 m

Analytický výpočet: 0,0006747185 m

Výsledky numerických simulací souhlasí s analyticky získaným výsledkem. Hodnota vychýlení nosníku na jeho konci se při využití nestrukturované sítě liší od analyticky získaného výsledku o necelých 0,09%. U strukturované sítě se výsledná hodnota simulace od analytického výsledku liší o přibližně 0,05%. Ze získaných hodnot je patrné, že byť strukturovaná síť obsahuje nižší počet prvků a šlo by ji tedy považovat za méně detailní, tak díky vhodném uspořádání prvků vzhledem k řešenému problému podává v tomto případě přesnější výsledky než síť nestrukturovaná.

Obrázek 3: průhyb spojitě zatíženého nosníku

3 Numerická simulace obtékání tělesa tekutinou

3.1 Navierovy-Stokesovy rovnice

S jistou formou proudění se každý z nás setkává denně a není tedy divu, že se o popsání tohoto jevu vědci dlouhá léta snažili. První, kdo rovnice popisující proudění vazké nestlačitelné tekutiny odvodil, byl Claude Louis Navier a o několik let později George Gabriel Stokes.

Navierovy-Stokesovy rovnice jsou nejznámější model proudění nestlačitelné viskózní tekutiny a mají velmi široké využití od zkoumání počasí po proudění tekutin v potrubí [5].

Navierovy-Stokesovy rovnice vyplývají z užití Newtonova druhého zákona. Ten lze pro proudící tekutinu zapsat ve tvaru ( viz [2], [4] a [5] )

ρ ( ∂

u

∂t +u⋅∇ u)=f +∇⋅ σ (u , p)

^, (32)

Kde symbol

ρ

je hustota tekutiny, symboluje vektor rychlosti,

t

je čas,

f

představuje objemovou sílu, která je často volena rovna nule a

σ (u , p)

je tensor napětí, který je pro Newtonovské kapaliny dán vztahem

σ (u , p)=2μϵ (u)− pI

^{. (33)}

Symbol

I

značí jednotkovou matici, symbol

μ

je dynamická viskozita proudící tekutiny a symbolem

ϵ (u)

je označen tensor rychlosti deformace, který je dán vztahem

ϵ (u)= 1

2 (∇ u+(∇ u)

^T

)

(34)

neboli

ϵ (u)= 1

2 ( [ ^{∂u ∂u ∂ u ∂ x ∂ x ∂ x}

^xyz

^{∂u ∂u ∂ u ∂ y ∂ y ∂ y}

^xyz

^{∂ u ∂u ∂u ∂ z ∂ z ∂ z}

^xyz

] ⁺ [ ^{∂ u ∂ u ∂ u ∂ y ∂ x ∂ z}

^xxx

^{∂u ∂u ∂u ∂ y ∂ x ∂ z}

^yyy

^{∂u ∂u ∂u ∂ y ∂ x ∂ z}

^zzz

] ⁾

^{. (35)}

Po dosazení vztahu (34) do rovnice (33), lze tenzor napětí zapsat

Výsledný tvar tenzoru napětí (36) je poté možné dosadit do rovnice (32), čímž se získá vektorový tvar Navierových-Stokesových rovnic.

Vektorový tvar Navierovy-Stokesovy rovnice je

∂u

∂ t +u⋅∇ u=ν Δ u− 1

ρ ∇ p+f

^{, (37)}

kde symbol

ν

je kinematická viskozita

ν = μρ

^{. (38)}

Neznámými v rovnici (37) jsou rychlostua tlak

p

. Soustavu je nutné doplnit o rovnici kontinuity, která vyjadřuje zákon zachování hmoty. V diferenciálním tvaru lze rovnici kontinuity zapsat

∇⋅u=0

. (39)

Výsledná soustava rovnic popisuje proudění nestlačitelné newtonovské tekutiny, což znamená, že libovolná část tekutiny zachovává svůj objem [5].

Soustavu rovnic (37) a (39) je nutné ještě doplnit o počáteční a okrajové podmínky.

Jako počáteční podmínku volíme hodnotu rychlosti a tlaku v počátečním čase. Tato hodnota je často volena rovna nule.

Okrajové podmínky lze charakterizovat jako popis hledané funkce v mezních bodech.

Nejčastější předepisovanou okrajovou podmínkou je Dirichletova okrajová podmínka (okrajová podmínka prvního typu)

v =g

, (40)

kde

v

je neznámá funkce a

g

je známá skalární funkce udávající chování funkce neznámé na okrajích (často

g =0

).

Dalším příkladem okrajové podmínky je Neumannova okrajová podmínka (okrajová podmínka druhého typu)

∂v

∂n =g

, (41)

kde

n

je vnější jednotková normála k hranici oblasti.

Další tvar Navierových-Stokesových rovnic, který bude v této práci uveden, je tzv. bezrozměrný tvar. Ten je možné získat zavedením parametru charakteristického rozměru

L

a parametru charakteristické rychlosti

U

a následnou úpravou jednotlivých členů rovnice na bezrozměrné proměnné [2]

u '= u

U ∇'=L ∇ x ' = x

L y ' = y

L z '= z L

^.

t '= t U

L p' = p

ρ⋅U

²

(42)

Výsledné bezrozměrné rovnice jsou

∇'⋅u'=0

∂u '

∂t ' =−∇ ' p' + 1

R e ∇ '

²

u '

, (43)

kde parametr

R e

je Reynoldsovo číslo dané vztahem

R e= ρ U L

μ

^{. (44)}

Reynoldsovo číslo je bezrozměrný parametr, který charakterizuje chování proudící tekutiny.

Pomocí znalosti Reynoldsova čísla lze usuzovat, zda bude proudění laminární či turbulentní.

Hodnota Reynoldsova čísla, při které dochází k přechodu z laminárního proudění do turbulentního, se nazývá kritická hodnota. Tato hodnota se pro různé druhy uspořádání liší a zjišťuje se experimentálně.

Reynoldsovo číslo má velký význam při studiu proudění. U některých těles, jako jsou například letadla, není možné kvůli rozměru tělesa měřit působení odporových sil a je nutné tyto síly získat z modelu. Získané hodnoty budou použitelné pro originál za předpokladu, že Reynoldsovo číslo zůstane stejné [9].

3.2 Parametry simulace

Geometrie výpočetní oblasti a okrajové podmínky jsou nastaveny podle specifikace benchmarkové úlohy [1].

3.2.1.Geometrie oblasti

Geometrie oblasti, ve které tekutina proudí, je dvourozměrný kanál šířky 0,41 metrů, délky 2,5 metrů. V kanálu se na souřadnicích 0,2 metrů a 0,2 metrů nachází překážka válcového tvaru s průměrem 0,05 metrů s uprostřed vetknutým nosníkem délky 0,35 metrů a šířky 0,02 metrů.

Výslednou oblast je možné vidět na obrázku 4. Hodnoty všech parametrů jsou v tabulce 3.

Tabulka 3: hodnoty parametrů simulace obtékání tělesa tekutinou

Šířka kanálu

H =0 ,41

m

Délka kanálu

L=2 ,5

m

Poloměr válce

R =0 ,05

m

Poloha středu válce na vodorovné ose

C 1=0 ,2

m Poloha středu válce na svislé ose

C 2 =0,2

m

Šířka nosníku

h=0 ,02

m

Délka nosníku

l =0 ,35

m

Obrázek 4: geometrie oblasti

3.2.2.Parametry proudící tekutiny

Jako proudící tekutina byly zvolena tekutina s hustotou

ρ =1000 kg

m

³ a kinematickou viskozitou

ν =0,001 m s

² . Jednotlivé parametry lze vidět v tabulce 4.

Tabulka 4: parametry proudící tekutiny

Hustota proudící kapaliny

ρ =1000 kg

m

³ Kinematická viskozita

ν =0,001 m s

²

3.2.3.Okrajové podmínky

Jako rychlostní profil na vstupu byl zvolen kvadratický profil s maximem ve středu kanálu popsaný vztahem

u(0 , y)=3 y (H − y)

( ^{H 2} )

^{2 . (45)}

Pro zadanou šířku kanálu

H

je výsledná rychlost na vstupu

u (0 , y)= 12

0 ,1681 y (0 ,41− y)

. (46)

Jako okrajová podmínka na výstupu se volí referenční hodnota tlaku. V případě simulace proudění nestlačitelné kapaliny může tato hodnota být volena libovolně. V případě proudění stlačitelné tekutiny má tato hodnota vliv na napětí [1]. V tomto případě má zvolená referenční hodnota tlaku na výstupu nulovou průměrnou hodnotu.

Na zbylých okrajích, tedy na horní a spodní stěně a na okrajích obtékaného tělesa, je volena nulová rychlost tekutiny vzhledem k hranicím oblasti.

3.2.4.Výpočetní sítě

Výpočetní sítě, na kterých vlastní simulace proudění nestlačitelné viskózní tekutiny probíhaly, byly dvě. První síť, kterou jsem pro simulace použil, byla příliš hrubá a tak bylo nutné síť výrazně zjemnit. Druhá použitá síť obsahovala přibližně osminásobné množství prvků sítě první, čímž bylo dosaženo výrazně vyšší přesnosti simulace, avšak za cenu výrazného zvýšení doby trvání celé simulace. Hrubou výpočetní síť je možné vidět na obrázku 5, jemnou výpočetní síť na obrázku 6.

Obrázek 5: hrubá výpočetní síť (3794 prvků)

Obrázek 6: jemná výpočetní síť (30850 prvků)

Obrázek 7: jemná síť v okolí obtékaného tělesa

3.3 Implementace ve FEniCSu

V následujícím textu bude popsán kód implementující výpočet obtékání válce s deskou nestlačitelnou newtonovskou kapalinou. Velká část kódu byla převzata z [4]. Informace ohledně výpočtu aerodynamických sil a implementaci tohoto výpočtu ve FEniCSu byly získány z [11]. Informace ohledně použitých funkcí byly čerpány z [3]. Parametry simulace je možné nalézt v kapitole 3.2.

Nejprve je nutné importovat potřebné knihovny:

from fenics import * from mshr import *

Knihovna fenics obsahuje všechny potřebné funkce pro výpočet parciálních diferenciálních rovnic založených na metodě konečných prvků. mshr je komponenta FEniCSu pro generování sítí.

Dále je potřeba nadefinovat konstanty využité při výpočtu:

T = 10

num_steps = 50000 dt = T / num_steps mu = 1

rho = 1000

inflow_profile = ('0.2*1.5*4*x[1]*(0.41 - x[1]) / 0.1681', '0')

Konstanta T je celkový čas simulace, num_steps je počet kroků, dt je délka jednoho kroku, mu je dynamická viskozita, rho je hustota a inflow_profile je rychlost na vstupu do oblasti.

Následně je třeba načíst, případně vytvořit výpočetní síť. Jak již bylo v kapitole 2.4 zmíněno, při načítání již vytvořené sítě je nutné, aby síť byla uložena ve formátu xml. Pokud je síť uložena v jiném formátu, je nutné nejprve provést konverzi pomocí nástroje dolfin-convertr, který je implementován v knihovnách FEniCS. Pokud síť není předem vytvořena, je možné ji vytvořit přímo ve FEniCS.

mesh = Mesh("meshLD.xml")

bndry = MeshFunction("size_t", mesh) s = Measure("ds", subdomain_data=bndry) I = Identity(mesh.geometry().dim())

Pomocí funkce MeshFunction bude možné přistupovat k jednotlivým součástem sítě. To bude třeba při definování hranic pro budoucí výpočet aerodynamických sil. K výpočtu aerodynamických sil bude dále potřeba funkce Measure, která je použita pro integraci přes povrch. Proměnnou s tak bude možné použít pro integraci přes povrch, kde vstupní parametr bude značit index plochy. Funkce Identity vytvoří jednotkovou matici. Rozměr matice je dán dimenzí úlohy. V tomto případě se jedná a matici 2x2.

Po načtení, popřípadě vytvoření sítě je třeba vytvořit funkční prostory. Proměnná V je vektorový funkční prostor pro rychlost, proměnná Q je skalární funkční prostor pro tlak. Pro rychlost jsou zvoleny po částech kvadratické prvky, pro tlak po částech lineární prvky.

V = VectorFunctionSpace(mesh, 'P', 2) Q = FunctionSpace(mesh, 'P', 1)

Funkce FunctionSpace přebírá tři vstupní parametry. První parametr je výpočetní síť, druhý parametr značí typ prvku a třetí parametr je stupeň prvku.

Aby bylo možné simulovat proudění nestlačitelné viskózní tekutiny, je nutné správně definovat hranice oblasti, ve které bude kapalina proudit. Celkově se jedná o čtyři hranice, z nichž první je přítok, druhá je odtok, třetí jsou stěny kanálu a poslední je těleso obtékané kapalinou.

Hranice obtékaného tělesa může být zadána přesně, ale lze použít funkce FEniCSu, u kterých stačí zadat část výpočetní sítě a funkce hranice rozpozná automaticky.

inflow = 'near(x[0], 0)' outflow = 'near(x[0], 2.5)'

walls = 'near(x[1], 0) || near(x[1], 0.41)'

cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.61 && x[1]>0.1 && x[1]<0.3'

Funkce near vrací booleovskou hodnotu udávající, zda se hodnota prvního parametru dostatečně blíží hodnotě druhého. Parametr x[0] značí souřadnice na ose x, parametr x[1] je souřadnice na ose y. on_boundary při definici hranice obtékaného tělesa dává funkci najevo, že hranice má hledat v síti v definované oblasti.

Další možný způsob definování hranic je pomocí přímého přístupu k součástem sítě. Zde je možné k jednotlivým částem oblasti přiřadit indexy. Toho je využito při výpočtu aerodynamických sil působících na obtékané těleso.

eps = 1e-10

for f in facets(mesh):

mp = f.midpoint() if near(mp[0], 0.0):

bndry[f] = 1 elif near(mp[0], 2.5):

bndry[f] = 2

elif near(mp[1], 0.0) or near(mp[1], 0.41):

bndry[f] = 3

elif mp.distance(Point(0.2, 0.2)) <= (0.05 + eps) or (0.20 <= mp[0] <= (0.6 + eps) and (0.19 - eps) <= mp[1] <= (0.21 + eps)):

bndry[f] = 5

Proměnná eps slouží k eliminaci chyb způsobených v důsledku zaokrouhlování čísel v počítačové aritmetice. Další možný způsob eliminace takto vzniklých chyb je užití funkce near.

Pokud by případné chyby nebyly eliminovány, mohlo by dojít k vynechání některých bodů a tím ke znehodnocení výsledků simulace.

Po definování jednotlivých okrajů a počáteční podmínky pro rychlost je možné přistoupit k definování okrajových podmínek.

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow) bcu_walls = DirichletBC(V, Constant((0, 0)), walls)

bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder) bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]

bcp = [bcp_outflow]

Funkce DirichletBC je funkce pro definici Dirichletovy okrajové podmínky. Prvním argumentem je funkční prostor, druhý argument je popis chování funkce na okraji a třetí argument je hranice. Všechny okrajové podmínky jsou následně uloženy do listu pro snazší přístup a úspornější zápis.

Jelikož jsou definovány dva různé funkční prostory, je nutné také vytvořit dvě sady bázových a testovacích funkcí

u = TrialFunction(V) v = TestFunction(V) p = TrialFunction(Q) q = TestFunction(Q)

Funkce TrialFunction vytvoří bázovou funkci daného funkčního prostoru. Funkce TestFunction vytvoří testovací funkci daného funkčního prostoru.

Dále je nutné vytvořit funkce reprezentující výsledky jednotlivých výpočtů. Jak pro rychlost tak pro tlak bude potřeba vytvořit dvě funkce. První funkce v sobě ponese informaci o předchozím výsledku, druhá funkce bude výsledkem v novém kroku.

u_i = Function(V) u_ = Function(V) p_i = Function(Q) p_ = Function(Q)

function vytvoří funkci v závislosti na vstupním parametru. Vstupní parametr může být funkční prostor nebo jiná funkce. Pokud je funkce tvořena z jiné funkce, je možné přidat druhý parametr určující číslo dílčí funkce, která má být extrahována.

Dále bude nutné definovat výrazy použité při výpočtu.

U = 0.5*(u_i + u) n = FacetNormal(mesh) f = Constant((0,0)) k = Constant(dt) mu = Constant(mu) rho = Constant(rho)

U je proměnná určující rychlost v polovině intervalu dána vzorcem

u

^{i+ 12}

≈ 1

2 (u

ⁱ

+u

^∗

)

, (47)

kde

i

je předchozí časový krok a

u

^∗je hodnota předběžné rychlosti, tedy hodnota získaná výpočtem z hodnot tlaku v předchozím časovém kroku.

Proměnná n je vnější jednotková normála, f je konstanta udávající objemovou sílu, k je délka jednoho časového kroku, mu je dynamická viskozita a rho je hustota.

Pro snazší práci s programem je dobré nadefinovat často používané výrazy. Výsledný program je poté kratší, snáze čitelný a při případné chybě snáze opravitelný. Z těchto důvodů je výpočet tenzoru rychlosti deformace (34) realizován jako funkce, která může být v případě potřeby využita a není tak nutné psát celý výraz znovu.

def epsilon(u):

return 0.5*(nabla_grad(u)+nabla_grad(u).T)

Funkce nabla_grad je jednou z funkcí obsažených v knihovnách FEniCS. Vstupním parametrem je vektor, výstupním pak jeho gradient. Parametr T značí transponovanou matici gradientu vstupního vektoru.

Ze stejných důvodů jako u tenzoru rychlosti deformace (34) je i výpočet tenzoru napětí (33) realizován jako funkce.

def sigma(u, p):

return 2*mu*epsilon(u) – p*I

Funkce epsilon je funkce pro výpočet tenzoru rychlosti deformace (34) popsaná výše.

Proměnná I je již definovaná jednotková matice.

Stejně jako výpočet tenzoru rychlosti deformace (34) a tenzoru napětí v kapalině (33) je i výpočet aerodynamických sil a jejich následné uložení realizováno pomocí funkce. Vstupním parametrem do této funkce je soubor, do kterého mají být hodnoty uloženy a čas, ve kterém síly na obtékané těleso působí.

Síla, která působí na povrch obtékaného tělesa, je dána jako integrál

F= ∫

S

T⋅n

, (48)

kde

S

je povrch obtékaného tělesa,

T

je tensor napětí a

n

je jednotkový normálový vektor k ploše

S

.

def Lift_Drag(LDFile,t):

T2 = sigma(u_,p_) force = -dot(T2, n) D = (force[0])*s(5) L = (force[1])*s(5) drag = assemble(D) lift = assemble(L)

LDFile.write("%e;%e;%e\n"%(t,lift,drag)) return

Proměnná T2 je tenzor napětí. K jeho výpočtu je využita již definovaná funkce popsaná výše.

Proměnná force značí vnitřek integrálu (48). Proměnná n je opačně orientovaná vnější normála. Proměnná D odpovídá aerodynamické odporové síle, proměnná L aerodynamické vztlakové síle. Parametr s značí integraci přes povrch, přičemž číslo 5 je identifikátor plochy.

Daný předpis pro výpočet sil musí být nejdřív sestaven. K tomu slouží funkce assemble, která předpis sestaví a vrátí odpovídající tenzor. Následně už je jen třeba získané výsledky uložit. K tomu slouží funkce write, která do předem otevřeného souboru uloží dané hodnoty.

Dále bude potřeba definovat rovnice popisující daný problém. Celý výpočet rychlosti a tlaku je rozdělen do tří částí. V první části je z předchozí hodnoty rychlosti a tlaku vypočtena předběžná rychlost. Výsledná rovnice pro první krok je

∫

_Ω

^{ρ ( u}

^∗

_{Δ t} ^−u

ⁱ

^{) ⋅v dx+} ∫

_Ω

^{σ (u}

^{i+ 12}

^{, p}

ⁱ

^{)⋅ ϵ (v)dx=} ∫

_Ω

^f

ⁱ⁺¹

^{⋅v dx}

− ∫

Ω

ρ⋅u

ⁱ

⋅∇ u

ⁱ

⋅v dx− ∫

∂Ω

p

ⁱ

⋅n⋅v ds+ ∫

∂Ω

μ⋅∇ u

^{i+ 12}

⋅n⋅vds

,

(49)

kde symbol

u

^∗je neznámá předběžná hodnota rychlosti,

u

ⁱa

p

ⁱ jsou hodnoty rychlosti a tlaku v předchozím časovém kroku,

Δ t

je velikost jednoho časového kroku,

v

je testovací funkce a

u

^{i+ 12} je hodnota rychlosti v polovině intervalu aproximovaná aritmetickou hodnotou (47). Celý výraz lze získat vynásobením rovnice (32) testovací funkcí

v

a následnou integrací. Výsledná implementace prvního kroku ve FEniCSu vypadá následovně.

F1 = rho*dot((u - u_i) / k, v)*dx \

+ rho*dot(dot(u_i, nabla_grad(u_i)), v)*dx \ + inner(sigma(U, p_i), epsilon(v))*dx \

+ dot(p_i*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \ - dot(f, v)*dx

a1 = lhs(F1) L1 = rhs(F1)

Proměnná u je neznámá rychlost, proměnná v je testovací funkce. Proměnná u_i je rychlost v minulém časovém kroku, proměnná p_i je tlak v předchozím čase, U je hodnota rychlosti v polovině intervalu (47), k je definovaná velikost jednoho časového kroku, n je vnější jednotková normála a f je objemová síla. Parametr dx značí integraci. Pro integraci přes povrch je použit parametr ds. Jelikož jsou napětí a rychlost deformace tensory, nemůže být pro výpočet integrálu

∫

_Ω

^{σ (u}

^{i+ 12}

^{, p}

ⁱ

^{)⋅ ϵ (v)dx}

použita funkce dot. Místo ní je použita funkce inner. Pro vektory vrací funkce inner stejnou hodnotu jako funkce dot. Nakonec je nutné separovat levou a pravou stranu rovnice. K tomu jsou užity funkce lhs a rhs, které rozdělení rovnice provedou automaticky.

V dalším kroku bude z výsledné rychlosti z prvního kroku vypočten tlak v novém čase. Výsledná rovnice pro tlak v čase i+1 je

∫

Ω

∇ p

ⁱ⁺¹

⋅∇ q dx= ∫

Ω

∇ p

ⁱ

⋅∇ q dx− 1

Δ t ∫

Ω

∇⋅u

^∗

⋅q dx

. (50) Symbol

q

v rovnici (50) je skalární testovací funkce z prostoru tlaku.

a2 = dot(nabla_grad(p), nabla_grad(q))*dx

L2 = dot(nabla_grad(p_i), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

V předchozím kroku byly jednotlivé strany rovnice získány funkcí lhs a rhs. V tomto případě je levá a pravá strana rovnice separována rovnou.

V posledním kroku je spočtena výsledná rychlost v čase i+1. Vzorec pro výpočet rychlosti je

∫

Ω

u

^{i+ 1}

⋅v dx= ∫

Ω

u

^∗

⋅v dx−Δ t ∫

Ω

∇( p

ⁱ⁺¹

− p

ⁱ

)⋅v dx

. (51)

a3 = dot(u, v)*dx

L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_i), v)*dx

Jelikož jsou levé strany všech tří rovnic časově nezávislé, není nutné je sestavovat v každém časovém kroku a stačí tak učinit pouze jednou.

A1 = assemble(a1) A2 = assemble(a2) A3 = assemble(a3)

Nakonec je nutné na výsledné matice použít definované okrajové podmínky. Pro matici získanou v první kroku jsou aplikovány okrajové podmínky určující rychlost tekutiny na stěnách a okraji obtékaného předmětu. Na matici získanou v druhém kroku je aplikována okrajová podmínka určující referenční tlak na výstupu. Třetí matice je počítána z hodnot získaných v prvním a druhém kroku a není tedy již nutné okrajové podmínky aplikovat.

[bc.apply(A1) for bc in bcu]

[bc.apply(A2) for bc in bcp]

Aby bylo možné s výsledky simulace dále pracovat, bude nutné je uložit. Výsledky budou ukládány do vybraných složek a souborů.

ufile = File("results/velocity.pvd") pfile = File("results/pressure.pvd")

File('navier_stokes_cylinder/cylinder.xml.gz') << mesh LDFile = open("results/LiftDrag.txt","a")

Před začátkem cyklu, ve kterém bude počítán průběh proudící kapaliny, je vytvořena proměnná udávající aktuální čas simulace.

t = 0

for x in range(num_steps):

Nejprve je k hodnotě aktuálního času přičtena hodnota časového kroku.

t += dt

Následně je zjištěn aktuální pokrok v simulaci a při dosažení další mezní hodnoty je tato skutečnost signalizována. Tato část kódu slouží čistě pro signalizace pokroku v simulaci.

if (x % (num_steps/100)) == 0 :

print("Done:",100*x/num_steps,"%")

V prvním kroku je nejprve sestavena pravá strana rovnice (49). Poté jsou na výsledný vektor použity okrajové podmínky. Nakonec je možné rovnici vyřešit.

b1 = assemble(L1)

[bc.apply(b1) for bc in bcu]

solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

K řešení rovnice je využita funkce solve. Tato funkce má celkem čtyři různé způsoby použití. V tomto případě je využita pro řešení lineárního systému, který je pro tento určitý případ

A 1 ⋅u

−

. vector ()=b 1

^, (52)

kde

A 1

je matice a

u

₋

. vector ()

^a

b 1

jsou vektory. Ostatní parametry jsou volitelné a určují způsob a podmínky řešení, což má vliv na dobu výpočtu a využití paměti.

Stejným způsobem jako probíhal první krok výpočtu probíhají i kroky další.

b2 = assemble(L2)

[bc.apply(b2) for bc in bcp]

solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg') b3 = assemble(L3)

solve(A3, u_.vector(), b3, 'cg', 'sor')

Po vyřešení všech rovnic je možné výsledky uložit. Jelikož je však jednotlivých kroků velké množství, nejsou výsledky ukládány v každém kroku. V tomto případě jsou výsledky uloženy v každém dvacátém průchodu cyklem.

if (x % 20) == 0 :

Lift_Drag(LDFile,t) ufile << u_

pfile << p_

Funkce Lift_Drag je funkce pro výpočet a uložení hodnot aerodynamických odporových sil popsaná výše. ufile je soubor pro ukládání hodnot rychlosti, pfile je soubor pro ukládání hodnot tlaku.

Nakonec je ještě nutné aktualizovat hodnoty rychlosti a tlaku v minulém kroku.

u_i.assign(u_) p_i.assign(p_)

Funkce assign přiřadí funkci, která je vstupním parametrem funkci druhé.

3.4 Výsledky numerických simulací

3.4.1. Výsledky simulace při použití hrubší výpočetní sítě

Na obrázcích 8 a 9 je možné vidět proudová pole získaná při použití hrubé sítě .

Z obrázku je patrné, že téměř nedochází k oscilacím úplavu za obtékaným tělesem. Tato skutečnost je zapříčiněna příliš nízkým počtem prvků výpočetní sítě. Při užití sítě jemnější, tedy sítě s více jak osminásobným počtem prvků než má síť použitá v tomto případě, již k oscilacím dochází. Výsledky numerických simulací s využitím jemnější sítě lze nalézt v kapitole 3.4.2.

Obrázek 8: rychlost tekutiny obtékající těleso v čase 9 sekund

Obrázek 9: tlak tekutiny obtékající těleso v čase 9 sekund

Na obrázcích 10 až 13 jsou vykresleny průběhy aerodynamické odporové a vztlakové síly působící na obtékané těleso.

Z grafu 10 je patrné, že i když z proudového pole by se dalo soudit, že k žádným oscilacím úplavu za obtékanou překážkou nedochází, tak úplav lehce osciluje. Tyto kmity mají mírně vzestupnou tendenci. Z detailnějšího zobrazení průběhu aerodynamické vztlakové síly v grafu 11 je patrné, že střední hodnota průběhu síly v čase 9 až 9,6 sekund se pohybuje přibližně kolem 13 N.

0 1 2 3 4 5 6 7 8 9 10

-15 -10 -5 0 5 10 15 20 25

t(s)

F_L(N)

Obrázek 10: graf aerodynamické vztlakové síly v čase 0 až 10 sekund

9 9,1 9,2 9,3 9,4 9,5 9,6

10 11 12 13 14 15 16

t(s)

F_L(N)

Obrázek 11: graf aerodynamické vztlakové síly v čase 9 až 9,6 sekund