Numerické modelování proudění ve vodních turbínách

www.KMA.zcu.cz

ˇ ve vodních Numerické modelování proudení turbínách Problémy a otázky Bohumír Bastl, Marek Brandner, Jiˇrí Egermaier, Hana Horníková, Kristýna Michálková, Jan Šourek, Eva Turnerová Katedra matematiky, Evropské centrum excelence NTIS – Nové technologie pro informaˇcní ˇ Západoˇceská univerzita v Plzni spoleˇcnost, Fakulta aplikovaných ved, cˇ erven 2016

ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

1 / 78

Obsah

www.KMA.zcu.cz

1

Úvod

2

Geometrické modelování

3

ˇ Numerické modelování proudení Navierovy-Stokesovy rovnice Galerkinova metoda ˇ Casová diskretizace ˇ Rešení nelineární soustavy (stacionární úlohy) Nekonformní síteˇ – DGFEM

4

Isogeometrická analýza B-spline, NURBS Použití v IgA

5

Numerické simulace

6

ˇ Záver ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

2 / 78

Úvod

www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

3 / 78

Úvod

www.KMA.zcu.cz

Úvod

Problematika ˇrešena v rámci dvou projektu: ˚ I

TA03011157 Inovativní postupy pro zvyšování užitných vlastností vodních turbín s využitím tvarové optimalizace založené na moderních metodách geometrického modelování

I

H2020 678727 MOTOR - Multi-ObjecTive design Optimization of fluid eneRgy machines


cˇ erven 2016

4 / 78

Úvod

www.KMA.zcu.cz

Úvod

ˇ Rešený problém má tyto cˇ ásti: I

Geometrické modelování vodní turbíny.

I

Numerické modelování dynamiky vazké nestlaˇcitelné tekutiny s využitím isogeometrické analýzy (s modelováním turbulence).

I

ˇ Tvarová optimalizace nekterých cˇ ástí turbíny. Cílem je optimalizace – pˇredchozí úlohy jsou její souˇcástí.


cˇ erven 2016

5 / 78

Úvod

www.KMA.zcu.cz

Kaplanova turbína


cˇ erven 2016

6 / 78

Úvod

www.KMA.zcu.cz

Kaplanova turbína


cˇ erven 2016

7 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

8 / 78


www.KMA.zcu.cz

Lopatkový kanál I I

I

I

ˇ plášt’ turbíny – vymezuje tak prostor, reprezentuje vnitˇrní a vnejší ve kterém proudí voda uvnitˇr turbíny vnitˇrní plášt’ ukrývá zaˇrízení související s pˇrevodem odebrané ˇ plášt’ je samotným pláštem ˇ energie na elektrickou energii, vnejší celého stroje ˇ plášt’ je rotaˇcní plochou – lze jej jak vnitˇrní, tak vnejší reprezentovat rovinnou kˇrivkou v prostoru složenou z úseˇcek, cˇ ástí kružnic a parabol a osou rotace pˇri tvorbeˇ modelu je tato kˇrivka reprezentována jako NURBS kˇrivka, výsledná rotaˇcní plocha je NURBS plochou


cˇ erven 2016

9 / 78


www.KMA.zcu.cz

ˇ lopatka Rozvádecí

I

ˇ jedná se o zobecnenou kuželovou plochu

I

rovinný profil lopatky v rovinné mˇríži se nejprve ztransformuje do kruhové mˇríže

I

ˇ „stoˇcena“ do kuželové plochy a následneˇ je kruhová mˇríž opet profil v kruhové mˇríži je tak nanesen na kužel

I

dále vytvoˇríme kuželovou plochu urˇcenou profilem naneseným na ˇ lopatky (osa kuželové ploše a pruseˇ ˚ cíkem osy rotace rozvádecí ˇ lopatku) a osou turbíny výseˇce oblasti pro rozvádecí

I

posledním krokem je oˇríznutí získané kuželové plochy kulovými plochami


cˇ erven 2016

10 / 78


www.KMA.zcu.cz

ˇ lopatka Rozvádecí


cˇ erven 2016

11 / 78


www.KMA.zcu.cz

ˇ Lopatka obežného kola

I

typicky urˇcena 3 − 7 rovinnými profily a odtokovou hranu (nebo osou lopatky)

I

každý rovinný profil je v rovineˇ pˇríslušneˇ natoˇcen a nanesen na ˇ (osou je vždy osa turbíny) válcovou plochu urˇcitého polomeru

I

pokud je zadána osa lopatky, je profil na válcové ploše posunut tak, aby osa lopatky procházela pˇríslušným bodem daného profilu

I

následneˇ jsou všechny profily propojeny do výsledné plochy lopatky s cílem dosáhnout co nejvyšší spojitosti

I

ˇ je získaná plocha opet ˇ oˇríznuta z obou stran kulovými na záver ˇ ˇ turbíny plochami, urˇcenými vnitˇrním a vnejším pláštem


cˇ erven 2016

12 / 78


www.KMA.zcu.cz

ˇ Lopatka obežného kola


cˇ erven 2016

13 / 78


www.KMA.zcu.cz

Savka I

I

ˇ poslední cˇ ást turbíny za obežným kolem, hlavním smyslem je odebrání co nejvíce ze zbývající kinetické energie vody, také pˇrispívá k tomu, že turbínu není nutné umístit tak hluboko do zemeˇ ˇ složena ze 3 cˇ ástí – za obežným kolem následuje kuželová cˇ ást, která pˇrechází v „kvádr se zaoblenými rohy“


cˇ erven 2016

14 / 78


www.KMA.zcu.cz

ˇ Objem vnitˇrní cˇ ásti turbíny pro simulaci proudení

I

B-spline model turbíny je následneˇ využit pro generování B-spline objemové parametrizace vnitˇrní cˇ ásti turbíny

I

tato objemová parametrizace dále slouží jako výpoˇcetní sít’ pro ˇ výpoˇcty provádené metodou isogeometrické analýzy typicky používáme dva typy objemových parametrizací:

I

I

I

ˇ objemová parametrizace vnitˇrní cˇ ásti turbíny pouze s rozvádecími ˇ lopatkami, ale bez obežných lopatek – slouží k získání rychlostního pole za ˇ rozvádecími lopatkami, využívá se periodicity oblasti ˇ ˇ objemová parametrizace obežného kola, resp. prostoru okolo jedné obežné ˇ ˇ na nekolik ˇ lopatky – rozdelen opet cˇ ástí, využívá se periodicity oblasti


cˇ erven 2016

15 / 78


www.KMA.zcu.cz

ˇ Objem vnitˇrní cˇ ásti turbíny pro simulaci proudení


cˇ erven 2016

16 / 78

ˇ Numerické modelování proudení

www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

17 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

18 / 78


www.KMA.zcu.cz

ˇ Úloha vazkého nestlaˇcitelného proudení I

ˇ Navierovy-Stokesovy rovnice pro vazké nestlaˇcitelné proudení ∂u + ∇p + u · ∇u − ν∆u = 0 ∂t ∇·u = 0

I

in Ω × (0, T ) in Ω

okrajové a poˇcáteˇcní podmínky u(x, t) = w(x, t) u(x, 0) = u0 (x)

on ∂Ω × [0, T ] (okrajové podmínky) in Ω (poˇcáteˇcní podmínky)

ˇ kde u = u(x, t) rychlost proudení, p = p(x, t) tlak, ν kinematická viskozita. I

Reynoldsovo cˇ íslo Re =

Ud ν


cˇ erven 2016

19 / 78


www.KMA.zcu.cz

ˇ Úloha vazkého nestlaˇcitelného proudení

I

Vazba mezi konvektivním cˇ lenem a vazkým cˇ lenem (Reynoldsovo cˇ íslo).

I

ˇ ˇ dat (DNS – výpoˇcetní Turbulentní proudení, prudké zmeny nároˇcnost, modely turbulence, stabilizace).

I

Tlak vystupuje pouze v cˇ lenu s první derivací (formálneˇ konvektivní cˇ len).

I

V rovnici kontinuity nevystupuje cˇ len s cˇ asovou derivací.


cˇ erven 2016

20 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

21 / 78


www.KMA.zcu.cz

Galerkinova metoda I

Prostorová diskretizace realizována pomocí isogeometrického pˇrístupu (jiné bázové funkce než ve FEM)

I

ˇ Definujeme koneˇcne-dimenzionální podprostory V h ⊂ V , V0h ⊂ V0 , W h ⊂ L2 (Ω)

I

Hledáme uh ∈ V h a ph ∈ W h takové, že pro všechny vh ∈ V0h a qh ∈ W h platí R R R ∂uh ∂t · vh + ν ∇uh · ∇vh + (uh · ∇uh )vh − Ω R Ω Ω (1) − ph ∇ · vh = 0, Ω

Z qh ∇ · uh = 0,

(2)

Ω


cˇ erven 2016

22 / 78


www.KMA.zcu.cz

Galerkinova metoda

I

Pˇribližné ˇrešení uh a ph je lineární kombinací bázových funkcí ϕui , ϕpi , tyto lineární kombinace dosadíme do (1) a (2) u

n X uh = (u1i , u2i )> ϕui , i=1 I

p

ph =

n X

pi ϕpi ,

i=1

Bázové funkce – B-splajny, NURBS


cˇ erven 2016

23 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

24 / 78


www.KMA.zcu.cz

ˇ Casová diskretizace

I

Aplikace cˇ asové diskretizace na celou soustavu (ˇrešení nelineární soustavy PDR se složitou vnitˇrní strukturou v každém cˇ asovém kroce).

I

Segregaˇcní metody (metoda korekce tlaku, metoda korekce ˇ rychlosti, metody štepení – spojitý a algebraický pˇrístup). V každém cˇ asovém kroku dekompozice na jednodušší úlohy.


cˇ erven 2016

25 / 78


www.KMA.zcu.cz

ˇ Casová diskretizace

Jednoduchá verze metody projekce, un ≈ u(x, tn ) : n+1 −un I u + ∇pn+1 + un · ∇un ∆t I un+1 = −∆t∇pn+1 + gn

I

∇ · un+1 = 0

I

∆t∆pn+1 = ∇ · gn

− ν∆un = 0

Problém s okrajovou podmínkou pro Poissonovu rovnici.


cˇ erven 2016

26 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

27 / 78


www.KMA.zcu.cz

ˇ Rešení nelineární soustavy (stacionární úlohy)

I

ˇ stlaˇcitelnosti (hledáme ustálený stav evoluˇcního Metoda umelé ˇ problému, využití pˇrístupu˚ pro stlaˇcitelné vazké proudení).

I

Picardova linearizace, Newtonova metoda (s vhodným ˇrešiˇcem získané soustavy lineárních algebraických rovnic s nesymetrickou maticí)

I

ˇ metody typu SIMPLE (Picardova linearizace a štepení na ˇ Poissonovu rovnici pro tlak a advekˇcne-difuzní rovnice pro složky rychlostí).

I

Segregaˇcní metody (hledáme ustálený stav, napˇr. metoda projekce)


cˇ erven 2016

28 / 78


www.KMA.zcu.cz

Stabilita I I

I I

I I

I

Standardní FEM, isogeometrická analýza, centrální diference ˇ v prostorové promenné v rámci FDM, centrální numerické toky v rámci FVM – lineárneˇ nestabilní metody (pro nulovou viskozitu). Upwind (FDM, FVM, FEM), Taylor-Galerkin FEM, SUPG, LW, BW – L2 stabilní metody. ˇ ˇrešení. Pokud pˇrevládá konvekce, (ne vždy) prudce se menící ˇ ˇ Nekteré z výše uvedených metod použitelné i pro nekteré úlohy ˇ s prudce se menícími daty. Nelze však použít obecneˇ (metody nezachovávají monotonii ˇrešení). ˇ V nekterých speciálních pˇrípadech lze dosáhnout L1 a L∞ stability (LW3, QUICK, schémata lichého ˇrádu). Metody, které zachovávají monotonii, jsou TVD, TVB, monotonní, ˇ jsou nelineární (stabilizace i „vuˇ ˚ ci prudkým zmenám dat, nikoliv pouze stabilizaˇcní metody pro konvekˇcní rovnici s nulovou viskozitou“). ˇ Stabilizace zaˇclenením modelu˚ turbulence. ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

29 / 78


www.KMA.zcu.cz

Stabilita II – LBB I

Standardní implementace FEM, FDM a FVM (stejné bázové funkce pro složky rychlosti a tlaku) vede v pˇrípadeˇ Navierových-Stokesových rovnic k nestabilnímu algoritmu. Nezávisle na Reynoldsovu cˇ íslu.

I

Problém nastává i u Stokesovy úlohy (problém nesouvisí s konvekcí, tlakový cˇ len má formálneˇ konvektivní tvar).

I

ˇ LBB podmínky (volba ruzného Splnení ˚ stupneˇ bázových funkcí).

I

Použití stabilizaˇcních pˇrístupu˚ (podobneˇ jako u advekce).

I

ˇ Segregaˇcní metody lze dát do souvislosti se stabilizací, umelou ˇ vazkostí a stlaˇcitelností (pˇresto je také nekdy nutná také ˇ LBB; obvykle se u metod projekce stabilizace nebo splnení ˇ nezminuje).


cˇ erven 2016

30 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

31 / 78


www.KMA.zcu.cz

Nekonformní síteˇ – DGFEM

I

Jednotlivé „patche“ na sebe nelze navázat konformneˇ (konformní sítí nelze diskretizovat celý objem, kterým proudí tekutina).

I

ˇ Potˇrebujeme použít metodu, která pˇripouští nekonformní síte.

I

Volíme metodu DGFEM – nespojitá Galerkinova metoda.

I

V pˇrípadeˇ DGFEM se používá nespojitá aproximace ˇrešení.

I

V našem pˇrípadeˇ nespojitost pouze na rozhraních patchu. ˚


cˇ erven 2016

32 / 78


www.KMA.zcu.cz

Nekonformní síteˇ – DGFEM Zavedeme oznaˇcení Z (p, q) =

pqdΩ,

pro skalární funkce,

ZΩ u · vdΩ,

(u, v) =

pro vektorové funkce.

Ω

1 1 }i + }j on Γ, 2 }i + 2 }j on Γ, {}} = } on ∂Ω. } on ∂Ω, pi ni + pj nj = ni (pi − pj ) on Γ, JpnK = pn on ∂Ω

J}K =


cˇ erven 2016

33 / 78


www.KMA.zcu.cz

Nekonformní síteˇ – DGFEM Hledáme uh ∈ V h and ph ∈ W h takovou, že ∀v ∈ V h a ∀q ∈ W h ∂uh , v + aIP (uh , v) + c(uh ; uh , v)− ∂t −b(v, ph ) + ({ph }, Jn · vK)Γ = lIP (v),

(3)

b(uh , q) − ({q}, Jn · uh K)Γ = 0,

where γ (Jn ⊗ uK, Jn ⊗ vK)Γ − h −(ν{∇u}, Jn ⊗ vK)Γ − (Jn ⊗ uK, ν{∇v})Γ , 1 1 c(u; w, v) = ((u · ∇)w, v) + ( {u · n}JwK − |{u · n}|JwK, JvK)Γ , 2 2 Z aIP (u, v) = (ν∇u, ∇v) +

q∇ · vdΩ, lIP (v) = (f , v).

b(v, p) = Ω


cˇ erven 2016

34 / 78


www.KMA.zcu.cz

Nekonformní síteˇ – DGFEM I

ˇ per partes, integrand cˇ lenu RNa konvektivní cˇ len uplatneno (u · n)w · vdΓ nahrazen vhodným numerickým tokem (pˇribližný Γ

Riemannuv ˚ ˇrešiˇc). I

ˇ Na vazký cˇ len uplatneno per partes (získáme cˇ leny na rozhraní), dále je pˇridána penaltová bilineární forma (spojitost) a cˇ len, který zajistí symetrii.

I

Na tuto soustavu lze uplatnit projekˇcní metodu. Druhá varianta: nejdˇríve uplatníme metodu projekce, poté na Poissonovu rovnici a ˇ advekˇcne-difuzní rovnice aplikujeme DGFEM.

I

Pro ˇrešení jednotlivých soustav lineárních algebraických rovnic použity v souˇcasnosti pˇrímé metody.


cˇ erven 2016

35 / 78

Isogeometrická analýza

www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

36 / 78


www.KMA.zcu.cz

Isogeometrická analýza – motivace I

inspirováno CAD (Computer Aided Design) – praxi je design typicky reprezentován pomocí CAD metod, výpoˇcetní sít’ je generována z tohoto CAD popisu

I

ˇ což muže sít’ obecneˇ nepopisuje oblast pˇresne, ˚ vnášet nepˇresnosti do výpoˇctu

I

ˇ generování síteˇ je vetšinou poloautomatický proces, který vyžaduje ˇradu ruˇcních zásahu˚ a je velmi cˇ asoveˇ nároˇcný

I

ˇ vyžadují nové generování síteˇ jakékoli designové zmeny

I

isogeometrická analýza - hlavním cílem je geometrická pˇresnost, nezávisle na hrubosti diskretizace

I

ˇ ˇ dalším cílem je zjednodušení procesu zjemnování síteˇ a zmeny tvaru oblasti


cˇ erven 2016

37 / 78


www.KMA.zcu.cz

Isogeometrická analýza I I I I I I I I

metoda pro analýzu problému˚ popsaných pomocí parciálních diferenciálních rovnic ˇ mnoho prvku˚ spoleˇcných s metodou koneˇcných prvku, ˚ nekteré prvky spoleˇcné s „bezmeshovými“ metodami je založena na NURBS objektech, které jsou standardem v souˇcasných CAD systémech název isogeometrická analýza proto, že bázové funkce prostoru ˇrešení a funkce použité pro reprezentaci geometrie jsou stejné výpoˇcetní sít’ je urˇcena cˇ ástmi NURBS objektu, ˚ máme tedy pˇresnou reprezentaci oblasti ˇ zjemnování ˇ lze snadno provádet síteˇ (h-refinement) nebo zvyšování stupneˇ bázových funkcí (p-refinement) ˇ ˇ novou možností je tzv. k-refinement (vysvetlíme pozdeji) ˇ zustává podstatnou vlastností je, že v každé úrovni zjemnení ˚ ˇ objekt reprezentován pˇresne ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

38 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

39 / 78


www.KMA.zcu.cz

B-spline báze I I

B-spline a NURBS objekty jsou založeny na tzv. B-spline bázi k jejímu zkonstruování potˇrebujeme vektor parametrizace T = (t1 , t2 , . . . , tn+p+1 ),

I

I I

ti ≤ ti+1 ,

ti nazýváme uzly, p je stupenˇ B-spline báze a n je poˇcet bázových funkcí ˇ eˇ – uniformní vektor parametrizace, uzly rozloženy rovnomern jinak neuniformní obecneˇ mohou být uzly ve vektoru parametrizace násobné definice B-spline báze 1 ti ≤ t < ti+1 Ni,0 (t) = 0 jinde ti+p+1 − t t − ti Ni,p−1 (t) + Ni+1,p−1 (t) Ni,p (t) = ti+p − ti ti+p+1 − ti+1 ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

40 / 78


www.KMA.zcu.cz

B-spline báze I

I

I

B-spline bázové funkce stupneˇ p mají obecneˇ p − 1 spojitých derivací, tzn. kvadratické bázové funkce obecneˇ spadají do tˇrídy C 1 apod. násobné uzly ve vektoru parametrizace snižují spojitost bázových funkcí (a následneˇ i odpovídající B-spline/NURBS kˇrivky) – jestliže vektor parametrizace obsahuje k-násobný uzel (k > 1 ∧ k <= p), potom odpovídající bázové funkce jsou tˇrídy C p−k pro neperiodický vektor parametrizace bázové funkce interpolují kraje intervalu parametru (podobneˇ pro uzel násobnosti p)

Duležité ˚ vlastnosti B-spline bázových funkcí: Pn I ∀t : i=1 Ni,p (t) = 1 I

support každé bázové funkce Ni,p je kompaktní (lokální) – jde o interval [ti , ti+p+1 ]

I

bázové funkce jsou nezáporné (Ni,p (t) ≥ 0, ∀t) ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

41 / 78


www.KMA.zcu.cz

B-spline báze – periodický vektor parametrizace I

bázové funkce pro uniformní periodický vektor T = (0, 1, 2, . . . , 7) stupneˇ 0 až 3 – vždy stejné, pouze posunuté bázové funkce


cˇ erven 2016

42 / 78


www.KMA.zcu.cz

B-spline báze – otevˇrený (neperiodický) vektor parametrizace I

vektor parametrizace T = (0, 0, 0, 0, 1, 2, 3, 4, 4, 4, 4), tj. uniformní vektor parametrizace, p = 3, n = 7, výsledný spline je tˇrídy C 2


I

vektor parametrizace T = (0, 0, 0, 1, 2, 2, 3, 3, 3), tj. neuniformní vektor parametrizace, p = 2, n = 6, výsledný spline je obecneˇ tˇrídy C0

cˇ erven 2016

43 / 78


www.KMA.zcu.cz

B-spline kˇrivky I I I

získáme jako lineární kombinaci B-spline bázových funkcí, koeficienty bázových funkcí nazýváme ˇrídící body lomená cˇ ára spojující ˇrídící body se nazývá ˇrídící polygon pro ˇrídící body Pi , i = 1, . . . , n, a bázové funkce Ni,p , i = 1, . . . , n, je B-spline kˇrivka stupneˇ p dána vztahem n X C(t) = Pi Ni,p (t) i=1

Duležité ˚ vlastnosti B-spline kˇrivek: I

B-spline kˇrivka stupneˇ p je obecneˇ (p − 1)-krát spojiteˇ diferencovatelná, pokud vektor parametrizace neobsahuje ˇ ˇrídící body násobné uzly nebo nesplývají nekteré

I

uzel násobnosti k ve vektoru parametrizace snižuje spojitost B-spline kˇrivky o k − 1

I

afinní invariantnost – je jedno, jestli transformujeme afinní cˇ erven 2016 transformací kˇrivku nebo její ˇrídící polygon, výsledek je stejný


44 / 78


www.KMA.zcu.cz

B-spline kˇrivka – pˇríklad I

B-spline kˇrivka urˇcená ˇrídícím polygonem P1 = [0, 0], P2 = [1, 4], P3 = [3, 5], P4 = [5, −1], P5 = [7, −1], P6 = [9, 2] a vektorem parametrizace T = (0, 0, 0, 1, 2, 2, 3, 3, 3):

I

kˇrivka je složena ze 3 oblouku˚ a je tˇrídy C 0 , jak naznaˇcoval vektor parametrizace ˇ první cˇ ást kˇrivky je ovlivnena pouze prvními tˇremi bázovými funkcemi (žlutá, zelená, tyrkysová), ostatní jsou na [0, 1] nulové, a ˇ tuto cˇ ást kˇrivky, podobneˇ tedy pouze první tˇri ˇrídící body ovlivnují pro ostatní cˇ ásti ˇrídícího první poslední cují teˇcny ke kˇr2016 ivce na ˇ hrana Numerickéa modelování proudení ve vodních turbínách polygonu urˇ cˇ erven 45 / 78

I

I


www.KMA.zcu.cz

B-spline plochy

I

pro danou ˇrídící sít’ {Pi,j }, i = 1, . . . , n, j = 1, . . . , m, a vektory parametrizace U = (u1 , . . . , un+p+1 ), V = (v1 , . . . , vm+q+1 ) je B-spline plocha stupneˇ p, q definována vztahem S(u, v) =

n X m X

Pi,j Ni,p (u)Mj,q (v),

i=1 j=1

kde Ni,p (u), Mj,q (v) jsou B-spline bázové funkce odpovídající vektorum ˚ parametrizace U , V I

tímto zpusobem ˚ mužeme ˚ popsat i rovinné oblasti

I

ˇ analogicky mužeme ˚ popsat i telesa


cˇ erven 2016

46 / 78


www.KMA.zcu.cz

NURBS – racionální B-spline kˇrivky I

I I

ˇ nekteré standardní objekty (napˇr. kružnice, elipsa), nelze popsat jako B-spline kˇrivky, proto se zavádí racionální B-spline kˇrivky ˇ (NURBS kˇrivky), které tento nedostatek odstranují ˇ ˇrídících bodu˚ NURBS objekt je urˇcen ˇrídícími body, váhami techto a potˇrebným poˇctem vektoru˚ parametrizace NURBS objekt v Rd popsaný parametrizací CR (u) je možné získat projektivní transformací B-spline objektu v Rd+1


cˇ erven 2016

47 / 78


www.KMA.zcu.cz

Racionální B-spline kˇrivky/plochy I

pokud Pi a wi , i = 1, . . . , n jsou rˇídící body a váhy NURBS kˇrivky, T je vektor parametrizace, potom parametrizace NURBS kˇrivky stupneˇ p je Pn Ni,p (u)wi Pi C(u) = Pi=1 (4) n i=1 Ni,p (u)wi

I

analogicky pro plochy, NURBS plocha stupneˇ p, q pro ˇrídící sít’ {Pi,j }, s váhami wi,j , i = 1, . . . , n, j = 1, . . . , m, a vektory parametrizace U = (u1 , . . . , un+p+1 ), V = (v1 , . . . , vm+q+1 ) je dána vztahem Pn Pm i=1 j=1 wi,j Pi,j Ni,p (u)Mj,q (v) S(u, v) = Pn Pm i=1 j=1 wi,j Ni,p (u)Mj,q (v)


cˇ erven 2016

48 / 78


www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

49 / 78


www.KMA.zcu.cz

„NURBS elementy“ I

I

I

analogií ke „koneˇcným prvkum“ ˚ v MKP jsou v isogeometrické analýze tzv. NURBS elementy (prvky), které vyplývají z vektoru˚ parametrizace v 1D jsou NURBS elementy cˇ ásti NURBS kˇrivky pro intervaly [ui , ui+1 ] nenulové délky ve vektoru parametrizace U = (u1 , . . . , un+p+1 ), napˇr. pro U = (0, 0, 0, 1, 2, 3, 3, 3) máme 3 NURBS elementy pro intervaly [0, 1], [1, 2], [2, 3] ve 2D, resp. 3D jsou NURBS elementy cˇ ásti NURBS plochy, resp. objemu, pˇríslušné oblastem [ui , ui+1 ] × [vj , vj+1 ], resp. [ui , ui+1 ] × [vj , vj+1 ] × [wk , wk+1 ] v parametrické oblasti

→ ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

50 / 78


www.KMA.zcu.cz

h-refinement I I I

I

ˇ analogem zjemnování síteˇ v MKP je v isogeometrické analýze vkládání uzlu˚ do vektoru parametrizace ˇ vložení uzlu je možné provést bez toho, aby se zmenila geometrie kˇrivky nebo dokonce její parametrizace ˇ pˇríklad: mejme kvadratický B-spline urˇcený ˇrídícími body P1 = [0, 0], P2 = [1/2, 1], P3 = [1, 0] a vektorem parametrizace T = (0, 0, 0, 1, 1, 1) a chceme vložit uzel t¯ = 1/2

nové ˇrídící body po vložení uzlu jsou P¯1 = [0, 0], P¯2 = [1/4, 1/2], P¯3 = [3/4, 1/2], P¯4 = [1, 0], parametricky a geometricky je ale kˇrivka modelování stejná proudení ˇ ve vodních turbínách Numerické cˇ erven 2016 51 / 78


www.KMA.zcu.cz

p-refinement I

I

I

I

ˇ možné provést beze zvýšení stupneˇ bázových funkcí je opet ˇ geometrie nebo parametrizace kˇrivky zmeny vektor parametrizace se upraví tak, že násobnost každého uzlu se zvýší o 1, abychom zachovali tˇrídu spojitosti; je potˇreba zvýšit i poˇcet ˇrídících bodu˚ ˇ pˇríklad: mejme stejný kvadratický B-spline jako v pˇredchozím a chceme zvýšit stupenˇ o 1

nový vektor parametrizace je T¯ = (0, 0, 0, 0, 1, 1, 1, 1) a nové ˇrídící body jsou P¯1 = [0, 0], P¯2 = [1/3, 2/3], P¯3 = [2/3, 2/3], P¯4 = [1, 0] ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

52 / 78


www.KMA.zcu.cz

k-refinement I I

I

využívá faktu, že procesy vložení uzlu a zvýšení stupneˇ nekomutují pokud nejprve vkládáme nové uzly a až poté zvyšujeme stupenˇ bázových funkcí, dostáváme pˇríliš mnoho bázových funkcí s nízkou tˇrídou spojitosti ˇ proto je nejprve zvýšit stupenˇ bázových funkcí na co výhodnejší ˇ nejhrubší síti a až poté vkládat uzly a zjemnovat sít’ – dostáváme mnohem méneˇ bázových funkcí, které navíc spadají do vyšší tˇrídy pojitosti T=(0,0,1,1)

→

T=(0,0,1/3,2/3,1,1)

→

T=(0,0,0,1/3,1/3,2/3,2/3,1,1,1)

T=(0,0,1,1)

→

T=(0,0,0,1,1,1)

→

T=(0,0,0,1/3,2/3,1,1,1)


cˇ erven 2016

53 / 78


www.KMA.zcu.cz

Rovinné a prostorové NURBS oblasti I jedním z hlavních geometrických problému˚ v isogeometrické analýze je nalezení

I

I I I

I I

popisu rovinných a prostorových NURBS oblastí z jejich hranice (která je obvykle dána CAD modelem) základním požadavkem je, aby zobrazení popsané parametrizací NURBS plochy/objemu, která zobrazuje parametrickou oblast do reálného prostoru, bylo prosté obvykle se navíc vyžaduje, aby nad celou parametrizací platilo det(J) > 0, kde J je Jakobián daného zobrazení (parametrizace) kvalita parametrizace oblasti muže ˚ mít znaˇcný vliv na výsledky výpoˇctu˚ ˇ pro generování rovinných NURBS oblastí existuje nekolik ruzných ˚ pˇrístupu, ˚ jak vygenerovat vnitˇrní body – jednodušší jsou založeny na ˇrešení lineární soustavy ˇ obvykle poˇcáteˇcní rozmístení ˇ vnitˇrních ˇrídících bodu˚ (napˇr. rovnic, složitejší pomocí diskrétních Coonsových plátu) ˚ modifikují na základeˇ daného optimalizaˇcního procesu (typicky nelineární úlohy) hlavní problémem (a v souˇcasnosti stále otevˇreným problémem) je generování objemových NURBS parametrizací pro prostorové NURBS oblasti urˇcené hranicí ˇ prozatím existují metody pouze pro nekteré speciální tˇrídy ploch, jako jsou napˇr. rotaˇcní plochy a sweep plochy ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

54 / 78


www.KMA.zcu.cz

Oblast urˇcená 4 okrajovými NURBS kˇrivkami I mejme ˇ jednoduše souvislou oblast

danou 4 okrajovými NURBS kˇrivkami a chceme najít ˇrídící body a pˇríslušné uzlové vektory pro rovinnou NURBS plochu, popisující tuto oblast: 1. pro každou dvojici „protilehlých“ kˇrivek potˇrebujeme, aby tyto kˇrivky ˇ stejný stupenˇ (zvýšení stupne) ˇ mely 2. ... a stejné uzly (vkládání uzlu) ˚ 3. potom najdeme vnitˇrní ˇrídící body – existují ruzné ˚ pˇrístupy 4. získáme tak NURBS plochu popisující danou oblast, která zachovává zadané hraniˇcní kˇrivky


cˇ erven 2016

55 / 78


www.KMA.zcu.cz




cˇ erven 2016

55 / 78


www.KMA.zcu.cz




cˇ erven 2016

55 / 78


www.KMA.zcu.cz




cˇ erven 2016

55 / 78


www.KMA.zcu.cz




cˇ erven 2016

55 / 78


www.KMA.zcu.cz




cˇ erven 2016

55 / 78

Numerické simulace

www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

56 / 78

Numerické simulace

www.KMA.zcu.cz

Numerické simulace - nestacionární úloha

I I I

I I I

Re = 200, nestacionární úloha nulová poˇcáteˇcní podmínka ˇ metoda projekce, zpetná Eulerova metoda pro cˇ asovou diskretizaci báze stupneˇ 2 a 3 nekonformní sít’ - DGFEM na rozhraních patchu˚ 18844 stupnˇ u˚ volnosti, ∆t = 0.1s, CF L = 0.5 ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

57 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I

Re = 200, nestacionární úloha, t = 1s nulová poˇcáteˇcní podmínka ˇ metoda projekce, zpetná Eulerova metoda pro cˇ asovou diskretizaci báze stupneˇ 2 a 3 nekonformní sít’ - DGFEM na rozhraních patchu˚ 18844 stupnˇ u˚ volnosti, ∆t = 0.1s, CF L = 0.5 ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

58 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

59 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

60 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

61 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

62 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

63 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

64 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

65 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

66 / 78

Numerické simulace

www.KMA.zcu.cz


I I I

I I I


cˇ erven 2016

67 / 78

Numerické simulace

www.KMA.zcu.cz

Numerické simulace - stacionární úloha

I

Re = 200, hledáme stacionární ˇrešení,

I

zastavovací podmínka:

I

metoda projekce 18844 stupnˇ u˚ volnosti, ∆t = 1.5s , 80 iterací (ˇcas. kroku) ˚

I

||un+1 −un || ||un+1 ||


< 10−5

cˇ erven 2016

68 / 78

Numerické simulace

www.KMA.zcu.cz

Numerické simulace I

18844 stupnˇ u˚ volnosti

I

nestacionární úloha I

I

I

segregaˇcní metodou (projekce): 157s (z toho 80% cˇ asu ˇrešení soustav - sparseLU), 900 MB RAM ˇrešení nesegregované soustavy: ≈ 8000s (z toho >99% cˇ asu ˇrešení soustav), 2, 2 GB RAM

stacionární úloha I

I

segregaˇcní metodou (projekce): 80 iterací, 119s ˇrešení nesegregované soustavy pomocí Picardových iterací: 16 iterací, 1621s


cˇ erven 2016

69 / 78

ˇ Záver

www.KMA.zcu.cz

1

Úvod

2


3


4


5

Numerické simulace

6


cˇ erven 2016

70 / 78

ˇ Záver

www.KMA.zcu.cz

ˇ Záver

I

Efektivní ˇrešení optimalizaˇcní úlohy

I

Využití knihovny G+Smo

I

Úˇcinná stabilizace, model turbulence

I

ˇ Efektivní iteraˇcní ˇrešiˇce (s pˇredpodmínením).


cˇ erven 2016

71 / 78

ˇ Záver

www.KMA.zcu.cz

ˇ Dekujeme za pozornost.


cˇ erven 2016

72 / 78

ˇ Záver

www.KMA.zcu.cz


I

K popisu geometrie a k popisu ˇrešení použijeme stejnou skupinu funkcí.

I

K popisu ˇrešení volíme funkce, které „jej dobˇre vystihují“. Tento ˇ po pˇrístup je obvyklý i v oblasti FVM (polynomy vyššího stupne, cˇ ástech konstantní funkce v okolí nespojitosti, T-splajny ve více dimenzích apod.).

I

Souˇcasneˇ však chceme konstruovat metodu použitelnou i pro ˇ nonkonformní síte.


cˇ erven 2016

73 / 78

ˇ Záver

www.KMA.zcu.cz

NURBS plochy I NURBS plocha stupneˇ p, q je dána ˇrídící sítí P (ˇrídícími body Pi,j , i = 0, . . . , n,

ˇ j = 0, . . . , m), váhami wi,j rˇídících bodu˚ a dvema uzlovými vektory U = (u0 , . . . , un+p+1 ), V = (v0 , . . . , vm+q+1 ) Pn Pm n X m X i=0 j=0 wi,j Pi,j Ni,p (u)Mj,q (v) S(u, v) = Pn Pm = Pi,j Ri,j (u, v) i=0 j=0 wi,j Ni,p (u)Mj,q (v) i=0 j=0

I B-splajnové bázové funkce Ni,p (t) jsou dány uzlovými vektory T a stupnem ˇ p

Ni,0 (t)

=

Ni,p (t)

=

1 ti ≤ t < ti+1 0 jinak t − ti ti+p+1 − t Ni,p−1 (t) + Ni+1,p−1 (t) ti+p − ti ti+p+1 − ti+1

I Volíme kombinace bázových funkcí splnujících ˇ LBB podmínku (pˇrípadneˇ

stabilizovanou verzi metody).


cˇ erven 2016

74 / 78

ˇ Záver

www.KMA.zcu.cz

ˇ Modelování turbulentního proudení I

I

I

I

I

Pˇrímá numerická simulace (volíme tak jemnou sít’, abychom postihli dynamiku všech víru) ˚ je extrémneˇ výpoˇcetneˇ nároˇcná (konec 21. století?). Existuje celá ˇrada pˇrístupu, ˚ pomocí kterých turbulenci ˇ škály než je modelujeme (modelujeme jevy, které mají jemnejší diskretizaˇcní sít’). ˇ Mezi cˇ asto zminované patˇrí RANS (Reynolds Averaged Navier-Stokes), LES (Large Eddy Simulation), ILES (Implicit LES), VMS (Variational MultiScale) . . . ˇ V rámci techto postupu˚ jsou obvykle pˇríslušné veliˇciny (rychlost, ˇ ˇ tlak) rozdeleny do více složek (dvou, tˇrí). Nekteré jsou stanovovány pomocí numerické metody pro ˇrešení NS rovnic, další jsou odpovídajícím zpusobem ˚ modelovány. ˇ ˇ V pˇrípadeˇ RANS pˇrístupu delíme veliˇciny na v cˇ ase prum ˚ erovanou složku a složku fluktuaˇcní ¯ (x) + u0 (x, t) . u (x, t) = u ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

75 / 78

ˇ Záver

www.KMA.zcu.cz

RANS method I

apply the decomposition and time - averaged Navier – Stokes equations ¯ ∂u ∂t

I

¯ · ∇¯ +u u = −∇¯ p + ν∆¯ u − u0 · ∇u0 , ¯ = 0. ∇·u

Extra new term appears in the system u0 · ∇u0 which is usually approximated using Boussinesq assumption 2 T 0 0 −u · ∇u = ∇ · νT ∇¯ u + (∇¯ u) − kI 3 where νT in turbulent viscosity and k is kinetic energy. ¯ ∂u ∂t

I

¯ · ∇¯ +u u = −∇¯ p + ∇ [(ν + νT ) ∇¯ u] + ∇(νT (∇¯ u)T ) − 23 ∇k, ¯ = 0. ∇·u

How to determine νT and k? ˇ ve vodních turbínách Numerické modelování proudení

cˇ erven 2016

76 / 78

ˇ Záver

www.KMA.zcu.cz

SST k − ω model I

we implement k − ω SST model which include two more equations for kinetic energy and turbulent specific dissipation ∂k ¯ · ∇k = Pk + ∇ · [(σk νT + ν) ∇k] − β ∗ kω, +u ∂t ∂ω 1 +¯ u·∇ω = αf +∇·[(σω νT +ν)∇ω]−βω 2 +2 (1 − F1 ) σω2 ∇k·∇ω, ∂t ω where " " ! ##4  √ k 500ν 4σ k ω2  F1 = tanh  min max , , β ∗ ωy y 2 ω CDkω y 2 1 ∇¯ u + (∇¯ u)T , 2 - constants dependent on wall distance y Then k νT = max (ω, SF2 )

Pk = min (νT f, 10β ∗ kω), β ∗ , β, α, σk , σω , σω2


f=

cˇ erven 2016

77 / 78

ˇ Záver

www.KMA.zcu.cz

Turbulent model I

Turbulent initial boundary value problem is described ¯ · ∇¯ +u u = −∇¯ p + ∇ [(ν + νT ) ∇¯ u] + ∇(νT (∇¯ u)T ) − 23 ∇k, ¯ = 0, ∇·u ∂k ¯ + u · ∇k = Pk + ∇ · [(σk νT + ν) ∇k] − β ∗ kω, ∂t ∂ω 1 +¯ u·∇ω = αf +∇·[(σω νT +ν)∇ω]−βω 2 +2 (1 − F1 ) σω2 ∇k · ∇ω. ∂t ω initial and boundary condition Discretization: similar to NS problem - FEM using quadratic basis ¯ and linear basis functions for pressure p, functions for velocity u kinetic energy k and specific dissipation ω Linearization: Picard’s method for both RANS and k − ω equations Disadvantage: k − ω models do not provide satisfactory results through near-wall regions - very fine grid is necessary near walls or the so called wall functions are applied ¯ ∂u ∂t

I I

I I


cˇ erven 2016

78 / 78

Numerické modelování proudění ve vodních turbínách

Recommend Documents