ÚVOD DO MODELOVÁNÍ V MECHANICE Přednáška č. 10
LAMINÁRNÍ PROUDĚNÍ NESTLAČITELNÝCH VAZKÝCH KAPALIN – aplikace v biomechanice Ing. Jan Vimmr, Ph.D.
Úvod do modelování v mechanice (UMM)
Obsah přednášky: 1. Základní pojmy 2. Proudění tekutin a jeho rozdělení 3. Laminární izotermické proudění nestlačitelné Newtonovy kapaliny 4. Speciální případy analytického řešení laminárního proudění 5. Aplikace laminárního proudění nestlačitelné kapaliny v biomechanice - proudění krve 2D a 3D modely femorálního a koronárního bypassu - ukázky vybraných numerických výsledků
Úvod do modelování v mechanice (UMM) Základní pojmy Všechny látky se skládají z atomů, které se sdružují v molekuly. Tuhé látky – veliké mezimolekulární síly → pravidelné uspořádání atomů do krystalické mřížky (krystalická struktura látky) – mezi molekulami, popř. atomy v krystalické mřížce působí síly přitažlivé (kohézní) a odpudivé (adhézní) → částice kmitají kolem rovnovážné polohy Tekutiny – látky, které nemají vlastní tvar a přijímají tvar nádoby, v níž se nacházejí – kapaliny – vytvářejí kapky (voda, olej, …), nemění samovolně svůj objem (molekuly netvoří stálou mřížku, ale působí mezi nimi ještě přitažlivé síly, které způsobují soudržnost kapaliny), jsou obecně málo stlačitelné, při pohybu (proudění) kladou odpor proti pohybu, tj. jsou vazké – plyny ( i páry) – soudržnost mezi molekulami téměř nulová → molekuly plynu se snaží vyplnit prostor, v němž se nacházejí → jsou rozpínavé – vzdálenosti mezi molekulami plynů jsou velké oproti kapalinám → jsou stlačitelné, málo vazké.
Úvod do modelování v mechanice (UMM) Vedle reálné (skutečné) tekutiny, která je stlačitelná a vazká, zavádíme pojem ideální (dokonalá) kapalina, která je nestlačitelná a nevazká, tj. bez vnitřního tření. Ideální tekutinu chápeme jako aproximaci reálné tekutiny. Tekutinu považujeme za spojité prostředí – izotropické kontinuum (stejné vlastnosti ve všech směrech) → parametry tekutiny (tlak, hustota, rychlost) se mění spojitě → vyjadřujeme je spojitými funkcemi. Z hlediska matematického můžeme vyjádřit stav elementárního objemu tekutiny dV a zákonitosti integrálně rozšířit na celý objem tekutiny V. Síly působící na tekutinu – vnitřní (dány vzájemným působením atomů a molekul) – vnější (vyvolány vnějším silovým polem) – objemové (setrvačné síly, např. odstředivá síla, gravitační síly) – plošné (tlakové síly, tečné (smykové) síly, kapilární síly) Viskozita tekutin – projevuje se při proudění reálných tekutin odporem proti pohybu – první formulace viskozity: Newton (1687) – potvrzena experimentálně
Úvod do modelování v mechanice (UMM) Představme si proudění ve vodorovném směru x podél desky jako pohyb tenkých vrstev tekutiny o tloušťce dy, rovnoběžných s deskou. Takové proudění ve vrstvách se nazývá laminární. Na desce je rychlost tekutiny nulová (ulpívá na ní). Rychlost ostatních vrstev se zvětšuje se vzdáleností od desky (brzdící účinek desky se zmenšuje). Jednotlivé vrstvy tekutiny vzájemně po sobě klouzají → dochází k jejich vzájemnému posuvu. Mezi vrstvami působí smykové (třecí) síly vyvolané viskozitou tekutiny.
tgα =
dx du ≡ dy dy
Tečné (smykové) napětí τ [Pa ] od viskozity je podle Newtona určeno vztahem
τ =η
du dy
Úvod do modelování v mechanice (UMM)
[ ]
du … gradient rychlosti s −1 v kolmém směru na pohyb tekutiny dy −1 −1 η … dynamická viskozita tekutiny Pa ⋅ s ≡ kg ⋅ m s
[
Zavádí se pojem: kinematická viskozita ν =
]
η ρ
[m
2
s −1
]
Viskozita tekutin je definována Newtonovým vztahem za předpokladu laminárního proudění. Dynamická a kinematická viskozita závisí na teplotě tekutiny. U plynů roste viskozita s teplotou. U kapalin s rostoucí teplotou viskozita klesá. Newtonské kapaliny – vyhovují Newtonovu zákonu viskozity Nenewtonské kapaliny – závislost smykového napětí τ na gradientu rychlosti du / dy nelze vyjádřit Newtonovým vztahem (např. krev při proudění nízkými rychlostmi v menších arteriích) .
Úvod do modelování v mechanice (UMM) Proudění tekutin a jeho rozdělení Proudění – pohyb tekutiny Hydrodynamika – nauka o proudění kapalin Aerodynamika (vnitřní, vnější) – nauka o proudění plynů Rozdělení proudění – podle fyzikální vlastnosti tekutin: 1. proudění ideální kapaliny a) potenciální proudění (nevířivé) – částice tekutiny se pohybují po křivočarých trajektoriích tak, že se vůči pozorovateli neotáčejí kolem vlastní osy
b) vířivé proudění – částice tekutiny se vůči pozorovateli natáčejí kolem vlastních os
Úvod do modelování v mechanice (UMM) 2. proudění reálných (vazkých) tekutin a) laminární – částice tekutiny se pohybují ve vrstvách (lamina – vrstva)
b) turbulentní – částice tekutiny mají kromě postupné rychlosti turbulentní (fluktuační) rychlost, jíž se přemisťují po průřezu – promíchávají se. – podle kinematických hledisek: 1. uspořádání proudění v prostoru a) proudění třírozměrné (prostorové) – rychlost u = u (x, y, z) b) proudění dvourozměrné (rovinné) – rychlost u = u (x, y) c) proudění jednorozměrné – u = u (x) 2. rozložení rychlosti v prostoru a) rovnoměrné proudění – vyvinuté proudění v trubici b) nerovnoměrné proudění – rychlost proudění v prostoru se mění, např. obtékání profilu v jeho blízkosti
Úvod do modelování v mechanice (UMM) 3. Závislost proudění na čase a) ustálené (stacionární) proudění – veličiny proudového pole (rychlost, tlak, hustota, teplota) se nemění s časem b) neustálené (nestacionární) proudění – veličiny proudového pole se mění s časem Částice tekutiny – elementární objem tekutiny vymezený uzavřenou kontrolní plochou Při popisu pohybu tekutiny můžeme užít dva přístupy: Lagrangeův popis – sledujeme pohyb určité částice tekutiny (analogie k vyšetřování pohybu hmotného bodu v mechanice tuhých těles) Eulerův popis – sledujeme proudění tekutiny v určitém místě (např. změnu rychlosti a tlaku). Tímto místem protékají různé částice tekutiny, což vede ke složitějšímu vyjádření zrychlení částice tekutiny ve sledovaném místě. Tento přístup se v mechanice tekutin užívá častěji. Výchozí systém rovnic popisujících proudění reálných tekutin ve 3D je vyvozen ze základních fyzikálních zákonů zachování:
Úvod do modelování v mechanice (UMM) a) zákon zachování hmotnosti → rovnice kontinuity (1 rce) b) zákon zachování hybnosti → pohybové Navierovy – Stokesovy rovnice (3 rce) c) zákon zachování celkové energie → energetická rovnice (1 rce)
Systém pěti nelineárních PDR
Rozdíl v kinematice laminárního a turbulentního proudění → plyne z časových průběhů rychlostí Laminární proudění – nedochází k promíchávání sousedních vrstev tekutiny Turbulentní proudění – časově střední hodnota rychlosti u s – turbulentní (fluktuační) složka rychlosti u′ (je malá, časově proměnná velikostí i směrem) → turbulence je nahodilý jev, který se vyhodnocuje statickými metodami
Úvod do modelování v mechanice (UMM) Řešení laminárního proudění – jednodušší ve srovnání s turbulentním – uplatňuje se Newtonův vztah pro smykové napětí – obecně pomocí numerických metod (metoda konečných objemů nebo metoda konečných diferencí) – speciální případy lze řešit exaktně (analyticky) – viz dále Výskyt laminárního proudění – proudění v úzkých plochých kanálech (malé průtokové rychlosti), např. zařízení hydraulických mechanismů a strojů – těsnící mezery, ložiska s hydrodynamickým mazacím filmem, … – proudění krve v arteriích
Úvod do modelování v mechanice (UMM) Laminární izotermické proudění nestlačitelné Newtonovy kapaliny Nestlačitelná kapalina → hustota kapaliny ρ = konst Izotermické proudění Newtonovy kapaliny → dynamická viskozita kapaliny η = konst Matematický model proudění ve 3D je tvořen soustavou rovnic (1) – (4): - rovnice kontinuity
- pohybové Navierovy – Stokesovy rovnice
∂u ∂v ∂w + + =0 ∂x ∂y ∂z
(1)
( )
∂u ∂ u 2 1 ∂p ∂ (uv ) ∂(uw) η ∂ 2u ∂ 2u ∂ 2u + + + + = 2 + 2 + 2 ∂x ∂y ∂z ∂z ρ ∂x ρ ∂x ∂y ∂t ∂v ∂ (uv ) ∂ v 2 1 ∂p ∂ (vw) η ∂ 2 v ∂ 2v ∂ 2v + + + + = 2 + 2 + 2 ∂x ∂y ∂z ρ ∂y ρ ∂x ∂y ∂z ∂t ∂w ∂ (uw) ∂(vw) ∂ w2 1 ∂p η ∂ 2 w ∂ 2 w ∂ 2 w + + + = 2 + 2 + 2 + ρ ∂z ρ ∂x ∂x ∂y ∂z ∂y ∂z ∂t
( )
( )
(2) (3) (4)
kde t je čas, p je tlak, v = (u, v, w) je vektor rychlosti kapaliny a y = (x, y, z ) je vektor prostorových souřadnic. T
T
Nelineární systém PDR (1) – (4) obecně nazýváme systém Navierových-Stokesových (NS) rovnic pro izotermické proudění nestlačitelné Newtonovy kapaliny.
Úvod do modelování v mechanice (UMM) Omezíme se dále na ustálené (plně vyvinuté) proudění mezi dvěma rovnoběžnými nekonečně širokými a nekonečně dlouhými deskami ve vodorovném směru. Vertikální vzdálenost desek je 2H. Proudění ve vodorovném směru → u ≠ 0, v = w → 0. Dosadíme do rovnice
∂u = 0 ⇒ u = u ( y, z ) , neboť ve směru x je složka rychlosti u ∂x konstantní. ∂u =0 Proudění je ustálené (stacionární) → ∂t
kontinuity (1) →
Dosazením uvedených předpokladů do NS rovnic (2) – (4) dostáváme:
∂ 2u ∂ 2u ∂p ∂2 p ∂p = η 2 + 2 ⇒ = 0 ⇒ = konst 2 ∂x ∂ y ∂ z ∂ x ∂ x ∂p ⇒ p není funkcí y =0 ∂y ⇒ p = p (x ) ∂p ⇒ p není funkcí z =0 ∂z
Úvod do modelování v mechanice (UMM) Jedná se o rovinné proudění (v rovinách rovnoběžných s rovinou xy jsou stejné rychlostní poměry) ⇒ u = u ( y ) … rychlost u není funkcí z Konečně dostáváme
dp d 2u = η 2 = konst , dx dy
p = p(x )
(5)
Rovnice (5) představuje matematický model nejjednoduššího laminárního proudění nestlačitelné Newtonovy kapaliny.
Úvod do modelování v mechanice (UMM) Příklad: Couetteovo proudění – proudění způsobené pouze pohybem horní desky dp rychlostí U = konst ⇒ =0 dx d 2u =0 dy 2
Rovnice (5) přejde na tvar: Řešení: u ( y ) = ± U 1 + y
2
u (− H ) = 0 , okrajové podmínky: u (H ) = ±U (6)
H
Couetteovo proudění je způsobené pohybem jedné stěny a rychlostní profil je lineární. Dále určíme hodnotu smykového napětí na stěně (WSS – wall shear stress):
τw =η
du dy
y =± H
= ±η
U 2H
Úvod do modelování v mechanice (UMM) Příklad: Proudění způsobené pouze tlakovým gradientem dp ≠ 0 , kdy obě desky dx jsou fixovány. Řešíme rovnici (5): Protože
dp d 2u = η 2 = konst dx dy
dp = konst , musí být rozložení dx
tlaku ve směru osy x přímkové. Kapalina proudí ve směru tlakového spádu ⇒ p1 > p2 Tedy: p(x ) = −
p1 − p2 x + p1 l
dp p − p2 =− 1 <0 dx l
, okrajové podmínky:
u (− H ) = 0 u (H ) = 0
Úvod do modelování v mechanice (UMM) Řešení:
u( y ) = −
(
)
(
)
1 dp 1 p − p2 H 2 − y2 ≡ H 2 − y2 1 2η dx 2η l
(7)
Rychlostní profil je v tomto případě, kdy proudění je způsobeno pouze tlakovým gradientem, parabola. Rychlost u je nezávislá na poloze x ⇒ ve všech průřezech je stejné parabolické rozložení rychlosti ⇒ plně vyvinuté proudění.
du 1 dp =0= y ⇒ y=0 , dy η dx
u ( y = 0) ≡
[
u max
H 2 dp H 2 p1 − p2 =− ≡ 2η dx 2η l
]
Průtočné množství (průtočný objem) Q m 3 s −1 kapaliny v mezeře definujeme vztahem:
Q = ∫ u ( y ) dA
(8)
( A)
Pro střední rychlost u avg kapaliny v mezeře platí: H
uavg
(
)
1 1 2 2 dp = − H − y ⋅ b dy = ∫ 2 Hb − H 2η dx
u avg =
Q 1 = ∫ u ( y ) dA A A ( A)
2 H 2 dp 2 − ⇒ uavg = umax 3 2η dx 3 u max
(9)
Úvod do modelování v mechanice (UMM) Smykové napětí na stěně mezery:
τw =η
du dy
y =± H
=y
dp dx
y =± H
= ±H
2η dp H ≡ ∓ ( p1 − p2 ) ⇒ τ w = ∓ u max dx l H
(10)
Úvod do modelování v mechanice (UMM) Příklad: Zobecněné Couetteovo proudění - proudění v mezeře je způsobeno kromě tlakového gradientu (dp / dx ≠ 0 ) také pohybem horní desky rychlostí U = konst.
dp d 2u u (− H ) = 0 =η = konst Řešíme rovnici (5): , okrajové podmínky: dx dx 2 u (H ) = ±U Řešení:
u( y ) = −
(
)
1 dp U y H 2 − y2 ± 1 + 2η dx 2 H
(11)
Rychlostní profil je v případě zobecněného Couetteova proudění sečtením rychlostních profilů (6) a (7) z předchozích příkladů
Úvod do modelování v mechanice (UMM) Pro střední rychlost uavg platí:
uavg
(
)
H 1 1 1 U y 2 2 dp = ∫ u ( y ) dA = − H − y ± 1 + bdy ∫ A ( A) 2 Hb − H 2η dx 2 H
2 U ⇒ uavg = umax ± 3 2
(12)
Smykové napětí na stěně mezery:
τw =η
du dy
y=± H
y dp U = η ± η dx 2 H
y =± H
⇒ τ w = ±H
dp η ± U dx 2 H
(13)
Úvod do modelování v mechanice (UMM) Příklad: Ustálené laminární proudění Newtonovy kapaliny ve vodorovné válcové trubce s vnitřním poloměrem R a s nedeformovatelnými stěnami V tomto případě platí: u = u ( y, z ), v = 0, w = 0
∂ 2u ∂ 2u ∂p ⇒ η 2 + 2 = = konst ∂ y ∂ z ∂ x
p1 > p2
dp p − p2 =− 1 <0 dx 1
Řešení tohoto problému provedeme ve válcových souřadnicích x, r,ϕ, kde r 2 = y 2 + z 2 ⇒ u = u (r ) Rovnováha sil působících na vytknutý objemový element proudící tekutiny:
dp 2 r dp πr ⇒ τ = dx 2 dx dp = konst Tečné napětí τ je podle (14) přímo úměrné poloměru r pro dx du τ = η Dosadíme do (14) Newtonův vztah a dostaneme: dr du r dp = , okrajové podmínky: u (R ) = 0 dr 2η dx
τ ⋅ 2π r ⋅1 − p2 ⋅ π r 2 + p1 ⋅ π r 2 = 0 ⇒ τ ⋅ 2π r =
(14)
Úvod do modelování v mechanice (UMM) Řešení: u (r ) = −
(
)
(
)
1 2 dp 1 2 p − p2 R − r2 ≡ R − r2 1 4η dx 4η l
(15)
Rychlostní profil je v tomto případě rotační paraboloid. Maximální rychlost umax v trubce je: u (r = 0) ≡ umax
R 2 dp R 2 p1 − p2 =− ≡ 4η dx 4η l
Průtočný objem (průtok kapaliny trubicí):
π R 4 dp Q = ∫ u (r ) dA = ∫ u (r ) 2π r dr ⇒ Q = − 8η dx ( A) 0 R
(16)
Vztah (16) odvodil francouzský lékař Poiseuille (1840), který studoval proudění krve v žilách. Nezávisle na něm odvodil tento vztah německý fyzik Hagen (1839) → vztah (16) pro průtočný objem označujeme jako Hagenova-Poiseuilleova formule. Z této formule plyne, že největší vliv na změnu proudění kapaliny má poloměr trubice. Má-li být průtočný objem kapaliny trubicí konstantní, tak 1%-ní zmenšení poloměru trubice vyžaduje 4 %-ní přírůstek tlakového spádu. Př.: Hypertenze (vysoký krevní tlak) je vyvolán zúžením krevních cév.
Úvod do modelování v mechanice (UMM) Pro střední rychlost u avg kapaliny v trubici platí:
Q R 2 dp uavg = =− ⋅ π R2 8η dx du Smykové napětí na stěně trubice: τ w = η dr
1 uavg = umax 2 R dp = 2 dx
⇒ r =R
(17) (18)
Smyková rychlost Newtonovy kapaliny na stěně trubice:
γɺ =
du dr
r =R
[ ]
4 uavg −1 R dp = =− s R 2η dx
Rozběhová dráha laminárního proudu kapaliny v trubici Vstup do trubice – kapalina má rychlostní profil odpovídající dokonalé kapalině. Vzdálenost xr , na níž se vyvine rychlostní profil ve tvaru paraboloidu, se nazývá rozběhová dráha laminárního proudu a platí pro ní Boussinesqův vztah
xr ≥ 0,065 Re , 2R
Re =
2 Rρ uavg
η
,
(19)
kde Re je Reynoldsovo číslo. Ze vztahu (19) je zřejmé, že k ustálení rychlostního profilu dojde daleko od vstupního průřezu → v krátkých trubkách se laminární rychlostní profil nevyvine, a proto u nich Hagenův – Poiseuilleův vztah neplatí.
Úvod do modelování v mechanice (UMM) Aplikace laminárního proudění nestlačitelné kapaliny v biomechanice – proudění krve v cévách přemostěných bypassovým štěpem Kardiovaskulární choroby (infarkt, mrtvice) – jsou příčinou 50% předčasných úmrtí v ČR Ateroskleróza – kornatění tepen vlivem usazování cholesterolu ve vnitřní vrstvě cévy → zúžení průsvitu cévy → snížení průtoku krve → nedostatečné prokrvení tkáně (ischemie). V případě srdce (ischemická choroba srdeční) může tento stav vést k infarktu myokardu (lokálnímu odumření srdečního svalu). – její výskyt ovlivněn lokálním charakterem proudění, nejčastěji v místech větvení cévy (bifurkace)
Úvod do modelování v mechanice (UMM) – při 50 – 70% zúžení průsvitu cévy (stenóza) – léčba medikamenty a úpravou životosprávy – při větší stenóze – balónková angioplastika (nechirurgický zákrok) – aplikace bypassových štěpů (chirurgický zákrok) – syntetické (polymery) – autologní (žilní, arteriální)
stehenní (femorální) bypass
kyčelní bypass
sekvenční aortokoronární bypass
Sekvenční bypass – vícenásobné přemostění nativní artérie
dvojitý aorto-koronární bypass
Úvod do modelování v mechanice (UMM) Základní typy anastomóz (spojení mezi bypassovým štěpem a arterií)
end-to-end
end-to-side
side-to-side
Životnost bypassových štěpů je omezena. Příčiny selhání: ztráta průchodnosti implantovaného bypassového štěpu v oblasti anostomózy při hojení narušené cévní stěny chirurgickým zákrokem
Ztrátu průchodnosti bypassového štěpu, resp. proces narušení cévní stěny ovlivňuje lokální charakter proudění (výskyt recirkulačních zón, nízké hodnoty smykového napětí na stěně – WSS, …)
Úvod do modelování v mechanice (UMM) → Snaha o lepší pochopení souvislostí mezi ztrátou průchodnosti bypassových štěpů a lokálním charakterem proudění vede k dlouhodobému výzkumu této problematiky. Proto je žádoucí zabývat se matematickým modelováním a numerickými simulacemi proudění krve v modelech bypassu a ověřovat numerické výsledky experimentálně s cílem optimalizace tvaru bypassových štěpů vedoucí k prodloužení jejich životnosti.
Úvod do modelování v mechanice (UMM) Ukázka vybraných numerických výsledků proudění krve idealizovanými modely bypassu Výpočetní sítě – 2D modelů bypassu (se stenózou, s okluzí, dvoucestný):
Výpočetní sítě – 3D modelů bypassu (s okluzí, se stenózou):
Úvod do modelování v mechanice (UMM)
2D koronární bypass s okluzí (Re=230, D=6,8mm, α=45°) – profily rychlostí
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=45°) – profily rychlostí
Úvod do modelování v mechanice (UMM)
2D koronární bypass s okluzí (Re=230, D=6,8mm, α=45°) – izo čáry rychlosti
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=45°) – izo čáry rychlosti
Úvod do modelování v mechanice (UMM)
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=20°) – profily rychlostí
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=45°) – profily rychlostí
2D koronární bypass se stenózou (Re=230, D=6,8mm, α =70°)
– profily rychlostí
Úvod do modelování v mechanice (UMM)
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=20°) – izo čáry rychlosti
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=45°) – izo čáry rychlosti
2D koronární bypass se stenózou (Re=230, D=6,8mm, α =70°)
– izočáry rychlosti
Úvod do modelování v mechanice (UMM)
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=20°, d(št ěp)=0,5D) – profily rychlostí
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=20°, d(št ěp)=D) – profily rychlostí
2D koronární bypass se stenózou (Re=230, D=6,8mm, α =20°, d(št ěp)=1,5D) – profily rychlostí
Úvod do modelování v mechanice (UMM)
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=20°, d(št ěp)=0,5D) – izočáry rychlosti
2D koronární bypass se stenózou (Re=230, D=6,8mm, α=20°, d(št ěp)=D) – izočáry rychlosti
2D koronární bypass se stenózou (Re=230, D=6,8mm, α =20°, d(št ěp)=1,5D) – izočáry rychlosti
Úvod do modelování v mechanice (UMM) 2D koronární bypass s dvoucestným štěpem (Re=230, D=6,8mm, α=45°) – profily rychlostí
Úvod do modelování v mechanice (UMM) 2D koronární bypass s dvoucestným štěpem (Re=230, D=6,8mm, α=45°) – izo čáry rychlosti
Úvod do modelování v mechanice (UMM) 3D koronární bypass s okluzí (Re=230, D=6,8mm, α=45°), newtonská kapalina
3D koronární bypass s okluzí (Re=230, D=6,8mm, α =45°), nenewtonská kapalina
3D femorální bypass s okluzí (Re=125, D=3,3mm, α =45°), newtonská kapalina
3D femorální bypass s okluzí (Re=125, D=3,3mm, α =45°), nenewtonská kapalina