Mechanika nenewtonských tekutin. Josef Málek

Mechanika nenewtonských tekutin

Josef Málek

1. přednáška, 5. října 2011 Otázky: Q1) Co se rozumí mechanikou? Q2) Co je tekutina? Q3) Co je newtonská tekutina? Q4) Co je nenewtonská tekutina? Q5) Proč je NDIR057 důležitá přednáška? Q6) Jaké jsou jevy, které nelze klasickými (newtonskými) tekutinami popsat? Q7) V čem je tento kurz výjimečný (oproti jiným kurzům z mechaniky tekutin možným absolvovat na jiných technických školách)? Q8) Neměli bychom používat jiné přístupy, např. teorie směsí, statistické přístupy, diskrétní přístupy? Odpovědi: Ad Q1) Mechanika nenewtonských tekutin je součást mechaniky kontinua (termodynamiky kontinua). Základní koncept – kontinuum

Q: Jak popsat odezvu materiálu na externí zatížení (mechanické, tepelné, elektrické, magnetické)? A: Dva typy prostředků k popisu odezev: 1) Obecnějšího charakteru – bilanční rovnice hmoty, hybnosti, momentu hybnosti, energie a formulace druhého zákona termodynamiky. 2) Konstitutivní vztahy (rovnice) – charakterizují odezvy „idealizovaných materiálůÿ. Příklady konstitutivních vztahů: (a) Elastický materiál (po vypnutí zatížení se materiál okamžitě vrátí do původního stavu) Pokud je vztah mezi napětím (silou) a deformací (změny úhlů, délek) lineární, pak mluvíme o lineárně pružné pevné látce. (b) Viskózní tekutina Pokud je napětí přímo úměrné rychlosti změny délek a úhlů, pak mluvíme o lineárně viskózní tekutině. Ad Q2) Tekutina je materiál, který neudrží smykové napětí v časovém měřítku pozorovatele. Ad Q3) Konstituvní vztah pro nestlačitelnou newtonskou tekutinu: T = −pI + 2µD(v),

div v = 0,

kde p je tlak, µ je dynamická viskozita, D je symetrická část gradientu rychlosti, v je rychlost tekutiny. Stlačitelná newtonská tekutina má tenzor napětí ve tvaru T = −p(ρ, θ)I + 2µ(ρ, θ)D(v) + λ(ρ, θ)(div v)I, kde ρ je hustota tekutiny, θ je teplota. Ad Q4) Nenewtonská tekutina je taková tekutina, kterou nelze popsat výše uvedenými vztahy. Ad Q5) Aplikace v biomateriálech (efekty v medicíně), geomateriálech (silnice, ranveje), chemickém průmyslu, potravinovém průmyslu (kečup ;-)). Ad Q6) Budeme mluvit o tekutinách, které se nemusí na první pohled jevit jako tekutiny, jak je známe – třeba voda. Ukážeme tekutiny, které vykazují chování, které newtonské tekutiny popsané Navierovými-Stokesovými rovnicemi nevykazují. Ukážeme například viskoelastické tekutiny, které se chovají částečně viskozně (jako newtonská tekutina – voda) a částěčně elasticky (jako elastická látka – guma). Nenewtonské jevy (ukázána videa):

1

• Weissenberg effect (rod climbing – šplhání po tyči) • Fano flow (tekutina se brání přerušení proudu ze stříkačky) • Open channel extensional flow (tekutina teče z kádinky dál, i když by už neměla) • Barus effect (Die swell) (rozšíření proudu při výtoku z trubice) • Kaye effect (tekutina volně padající do misky s tou tekutinou tvoří „tlustá stříkající vláknaÿ) • Electro-rheological fluid (tekutina mění ohromně svou viskozitu v elektrickém poli) (křeslo pro vozíčkáře) • Viscoealstic solid-like fluid (viskoelastická tekutina, která se chová více jako pevná látka) • Silly Putty (hračka s různým chováním, na krátké časové škále křehká, postupně s delším časem se chová více jako tekutina) • Shear thickening non-newtonian fluid (tekutina se při rychlém pohybu chová jako pevná látka, při pomalém pohybu jako tekutina) Příkladem nenewtonské tekutiny je třeba krev, je to viskoelastická shear-thinning tekutina. V tepnách proudí rychle a má nízkou viskozitu, ve vlásečnicích pomalu a má vysokou viskozitu (krev proudí okolo zranění pomalu s vysokou viskozitou, více ulpívá, zranění se snadněji zahojí). Ad Q7) Tento kurz je výjimečný v tom, že si zde ukážeme nové přístupy (Rajagopal), které se jinde neučí. Ad Q8) Budeme se soustředit na konstitutivní rovnice, zůstaneme ve standardním kontinuu.

2

2. přednáška, 12. října 2011

Rámec Newtonských tekutin Nechť B je abstraktní těleso, K0 (B) konfigurace v počátečním stavu a Kt (B) konfigurace v čase t (viz obrázek).

Pak definujme pohyb χ zobrazující z K0 (B) do Kt (B) : x = χ(X, t). Dále definujme další kinematických veličiny: vektor posunutí rychlostní pole deformační gradient

u = x − X = χ(X, t) − X v = ∂χ ∂t ∂χ FK0 = ∂X

Bilanční rovnice ρ˙ = −ρ div v ρv˙ = div T + ρb ρE˙ = div(Tv − q) + ρb · v + ρr,

bilance hmoty bilance hybnosti bilance energie

(I)

kde E = e + (ρ|v|2 )/2 a ∂z + ∇z · v, ∂t ∂z z˙ = + (∇z)v. ∂t z˙ =

Dále platí (v poslední rovnosti použijeme symetrii T = TT ) (I)3 − (I)2 · v ⇔ ρe˙ − div q = T · ∇v + ρr = T · D + ρr. Domácí úkol č. 1

Přepište rovnice do tvaru ∂ρ + div(ρv) = 0, ∂t ∂(ρv) + div(ρv ⊗ v) − div T = ρb, ∂t

(III)

∂(ρE) + div(ρEv) + div q − div(Tv) = ρb · v + ρr. ∂t Domácí úkol č. 2 Definujme operátor časoprostorové divergence Divt,x ve třech prostorových proměnných takto: ∂ ∂ ∂ ∂ ∂ Divt,x = + + + = + div . ∂t ∂x1 ∂x2 ∂x3 ∂t Najděte vektory Hi , i = 1, 2, 3, 4, 5, a přepište soustavu z domácího úkolu č. 1 do tvaru Divt,x Hi = gi ,

i = 1, 2, 3, 4, 5.

(Tij )3i,j=1

Naše neznámé jsou ρ, v, e a T = + předpoklad platnosti bilance momentu hybnosti T = TT ⇒ 6 veličin tenzoru napětí a znalosti q. Celkem tedy máme i pro tento jednoduchý případ 9 konstitutivních vztahů pro T, q. 3

Navierovy-Stokesovy-Fourierovy rovnice Stlačitelná tekutina T = −p(ρ, θ)I + 2µ(ρ, θ)D(v) + λ(ρ, θ)(div v)I, q = K(ρ, θ)∇θ, kde µ je smyková viskozita, λ je objemová viskozita a K(ρ, θ) je pozitivně definitní matice. Nestlačitelná tekutina Objem tekutiny je roven Z V (Pt ) =

dx. Pt

Materiál je nestlačitelný, když d V (Pt ) = 0 dt a z toho plyne det FK0 = 1 ⇒ div v = 0. Nestlačitelná tekutina splňuje div v = 0, T = −pI + 2µ(ρ, θ)D(v), q = K(ρ, θ)∇θ. Zde tlak p má jiný význam než v předchozí rovnici. U nestlačitelné tekutiny se jedná o neznámou (není způsob jak konstitutivně určit sférickou část tenzoru napětí), u stlačitelné musí být popsán konstitutivním vztahem. Z podmínky div v = 0 neplyne, že by hustota byla konstantní, tj. ρ(x, t) = ρ∗ , kde ρ∗ ∈ (0, ∞). Neboť z (I)1 plyne ∂ρ + ∇ρ · v = 0 ⇔ ρ je konstantní podél charakteristik, ∂t kde charakteristika je χ. Když v K0 (B)∃X1 , X2 : ρ(X1 ) 6= ρ(X2 ), pak takový materiál nazveme nehomogenní nestlačitelná tekutina. Je-li ovšem na počátku hustota konstantní, pak zůstává konstantní po celou dobu a mluvíme o homogenní nestlačitelné tekutině. Dále v případě izotermálního procesu, kde θ(x, t) = θ∗ ∈ (0, ∞) a ∇θ = 0, máme bilanční rovnice pro nehomogenní nestlačitelné Navierovy-Stokesovy rovnice div v = 0, ∂ρ + ∇ρ · v = 0, ∂t ∂(ρv) + div (ρv ⊗ v) = −∇p + div(2µ(ρ)D(v)) + ρb ∂t a bilanční rovnice pro homogenní nestlačitelné Navierovy-Stokesovy rovnice

ρ∗



div v = 0,  ∂v + div (v ⊗ v) = −∇p + µ∗ ∆v + ρ∗ b. ∂t

Newton (1687): ”The resistence arising from the want of lubricity in parts of the fluid, other things being equal, is proportional to the velocity with which the parts of the fluid are separated from one another.”

4

Síla je úměrná rozdílu rychlostí a nepřímoúměrná výšce h v(h) − v(0) F ∼ . A h |{z} Txy

Máme-li tedy jednoduché smykové proudění s rychlostí ve tvaru v = (v(y), 0, 0), pak Txy = µv 0 (y) = 2µDxy (v). Někdo může chápat slovíčko proportional obecněji ve smyslu libovolné závislosti, pak můžeme chápat Newtonův výrok ve smyslu implicitní vazby G(Txy , Dxy , . . . ) = 0. Všimněme si nyní fyzikálních rozměrů dynamické viskozity [µ] = a kinematické viskozity

kg m s−2 N m = s = kg m−1 s−1 m2 m s−1 m2 [µ∗ ] = m2 s−1 . [ρ∗ ]

Příklady velikostí viskozit Tekutina vzduch (18 ◦ C) voda olivový olej motorový olej SAE50 med kečup burákové máslo asfalt zemská kůra

Viskozita [cP] 0,02638 1 84 540 2000 – 3000 50000 – 70000 150000 – 250000 3 × 1010 3 × 1025

Domácí úkol č. 3 Jaká jednotka je označována cP? Má blíže k dynamické (µ) nebo kinematické viskozitě (µ∗ /ρ)? Jméno jakého vědce je „schovánoÿ v této jednotce? Jaká byla jeho oblast zájmu a co studoval?

5

3. přednáška, 19. října 2011 Převod do bezrozměrného tvaru Máme charakteristické veličiny T ∗ , L∗ , V∗ , p∗ a charakteristické proměnné x ˜=

t v p x ˜ ˜ = ∗ , p˜ = ∗ , , t = ∗, v ∗ L T V p

pak ∂ 1 ∂ 1 ∂ ∂ = ∗ , = ∗ . ∂t T ∂ t˜ ∂xi L ∂x ˜i Pokud dosadíme do Navierových-Stokesových rovnic, dostaneme ˜ = 0, div v ∗

∗

∗

˜ µ V V ∂v p∗ V ˜ ˜ div (˜ v ⊗ v ) − ∆ ∇x˜ p˜, v = − + x ˜ x ˜ T ∗ ∂ t˜ L∗ ρ∗ (L∗ )2 L∗ ∗

pokud V∗ = L∗ /T ∗ , pak přepíšeme na ˜ = 0, div v ∗

∗

˜ L ∂v µ 1 p∗ L ˜ ) − ∗ ∗ ∗ ∆x˜ v ˜ = − ∗ ∇x˜ p˜. v⊗v + ∗ 2 divx˜ (˜ ∗ 2 (T ) ∂ t˜ (T ) ρ T L L ∗

Nyní vynásobíme (T ∗ )2 /L∗ a dostáváme ˜ = 0, div v ∗

∗

˜ ∂v µ T p∗ (T ∗ )2 ˜ ) − ∗ ∗ 2 ∆x˜ v ˜=− + divx˜ (˜ v⊗v ∇x˜ p˜. ρ (L ) (L∗ )2 ∂ t˜ Nazveme nyní Re :=

ρ∗ (L∗ )2 µ∗ T ∗

Reynoldsovým číslem. Pokud (L∗ )2 , (T ∗ )2 pak dostáváme Navierovy-Stokesovy rovnice v bezrozměrném tvaru p∗ =

div v = 0, ∂v 1 + div(v ⊗ v) − ∆v = −∇p. ∂t Re Domácí úkol č. 4 Ukažte, že platí: Je-li (v, p) řešením bezroměrných Navierových-Stokesových rovnic, pak ∀λ > 0 vλ (t, x) = λv(λ2 t, λx), pλ = λ2 p(λ2 t, λx) řeší také bezrozměrné Navierovy-Stokesovy rovnice. Okrajové podmínky, speciální oblasti Mějme mechanicky izolované těleso, tj. nic neprotéká hranicí oblasti Ω, v · n = 0 na ∂Ω. Pokud platí vτ = 0 na ∂Ω, kde zτ = z − (z · n)n, mluvíme o no-slip podmínce – dokonalém ulpívání na hranici. Dále, pokud (Tn)τ = 0, jedná se o slip podmínku – dokonalý skluz na hranici. Kombinací těchto podmínek je Navierova1 podmínka skluzu λ(Tn)τ + (1 − λ)αvτ = 0, kde α je koeficient tření a λ ∈ [0, 1]. Mezní hodnoty dávají no-slip a slip podmínku. Cílem kurzu bude ukázat, že i tyto podmínky jsou konstitutivní rovnice a že podmínky na hranici lze také odvodit. 1 Navier

má ve Francii špatnou pověst, protože most, který stavěl, spadnul.

6

Jednoduché typy proudění Známe dva jednoduché typy proudění. Poiseuilleovo proudění je řízeno gradientem tlaku při proudění ve válci. Couetteovo proudění mezi dvěma soustřednými válci je řízeno rychlostí otáčení jednoho z válců. Rovinné Poiseuilleovo proudění je proudění mezi dvěma deskami řízené gradientem tlaku. Rovinné Couetteovo proudění v mezikruží je řízeno rychlostí otáčení jedné z kružnic. Domácí úkol č. 5 Hledejte řešení pro ustálené proudění ve tvaru v = (u(y), 0, 0) pro bezrozměrné NavierovyStokesovy rovnice v geometrii rovinného Poiseuilleova proudění.

Použijte tyto okrajové podmínky. Na hranici Γ1 u(1) = 0 a na Γ−1 uvažujte Navierovu okrajovou podmínku s λ ∈ [0, 1]. Podmínky na vstupu a výstupu Na Γin obvykle předepisujeme (i) Dirichletova podmínka pro rychlost v = α (ii) (Tn)n = p1 , vτ = 0 Na Γout obvykle předepisujeme (i) (Tn)n = 0, vτ = 0 (ii) (Tn)n = p2 , vτ = 0 Speciální okrajová podmínka je ”do-nothing” podmínka −pn + ν

∂v = 0. ∂n

Poznámka. Pro následující příklad s no-slip okrajovou podmínkou všude kromě modrého výstupu, kde je ”donothing” podmínka

by člověk očekával, že existuje jen jediné řešení této úlohy, a to nulové. Pokud to tak není, pak by to znamenalo, že ”do nothing” podmínka nemusí být dobrá.

Entropie a newtonské tekutiny Definice 1. Existuje specifická hustota entropie η, která je funkcí stavových veličin η = η˜(y0 , y1 , . . . ), přičemž y0 = e je vnitřní energie, a platí: (i)

∂ η˜ ∂˜ e > 0 ⇒ e = e˜(η, y1 , y2 , . . . ) označme teplotu θ := ∂e ∂η 7

(ii) η → 0+ , pokud θ → 0+ Z (iii) S(t) = (ρη)(t, x) dx → Smax pro t → +∞ pokud Ω je mechanicky a tepelně (energeticky) izolovaná Ω

Pro Newtonské tekutiny platí η = η˜(e, ρ) ⇔ e = e˜(η, ρ). Bilanční rovnice: ρ˙ = −ρ div v, ρv˙ = div T + ρb, ρE˙ = div(Tv + q) + ρb · v. Zderivujeme e ρe˙ = ρ

∂˜ e ∂˜ e η˙ + ρ ρ˙ ∂η ∂ρ

a dosadíme z bilančních rovnic ρb · v + div(Tv + q) − (div T) · v − ρb · v = ρ

∂˜ e ∂˜ e η˙ − ρ2 div v. ∂η ∂ρ

Nyní označme θ=

∂˜ e , ∂η

p = ρ2

∂˜ e ∂ρ

a dostáváme ρθη˙ = T · ∇v + div q + p div v symetrie T

=

T · D + div q + p div v.

Dále rozepíšeme pomocí deviatorických částí 1 Ad = A − (tr A)I 3 a získáváme 1 1 ρθη˙ = Td · Dd + (tr T) (div v) |{z} I · I + div q + p div v 3 3 3   1 = Td · Dd + (tr T) + p div v + div q. 3 Podělíme teplotou θ ρη˙ − div

q θ

    1 1 q · ∇θ d d = T ·D + (tr T) + p div v + . θ 3 θ | {z } ξ

Druhý zákon termodynamiky říká, že rychlost produkce entropie ξ ≥ 0, tedy ζ :=

q ξ = ρη˙ − div ≥ 0. θ θ

Označme m = (tr T)/3 průměr normálových napětí, pak platí      d 1 q · ∇θ ∇θ d d d ξ =T ·D + (tr T) + p div v + = T , m + p, q · D , div v, 3 θ θ | {z } | {z } termodynamické toky

termodynamické afinity

8

Jak vypadá ξ pro Navierovy-Stokesovy rovnice? Stlačitelné Navierovy-Stokesovy-Fourierovy rovnice T = −p(ρ, η)I + 2µ(ρ, η)D + λ(ρ, η)(div v)I, q = K(ρ, η)∇θ. Platí m+p=

2µ + 3λ div v 3

a Td = 2µDd . Dosaďme nyní do ξ. Varianta 1

2µ + 3λ |∇θ|2 (div v)2 + K . 3 θ Pokud µ ≥ 0, 2µ + 3λ ≥ 0 a K ≥ 0, pak ξ ≥ 0. ξ = 2µ|Dd |2 +

Varianta 2 ξ=

1 d2 3 1 |q|2 |T | + (m + p)2 + . 2µ 2µ + 3λ K θ

Pokud µ > 0, 2µ + 3λ > 0 a K > 0, pak ξ > 0. Ptáme se teď: Kdy ξ = 0? Varianta 1 Nechceme omezení na rychlost nebo teplotu, proto µ = 2µ + 3λ = K = 0. (To je ale divné, protože kdybych měřil viskozitu, tak ji nenaměřím nulovou.) Varianta 2 Lze jen, když Td = 0, m = −p, q = 0. 4. přednáška, 26. října 2011 Nestlačitelné Navierovy-Stokesovy-Fourierovy rovnice Varianta 1 ξ = 2µ|Dd |2 + K

|∇θ|2 . θ

Pokud µ ≥ 0 a K ≥ 0, pak ξ ≥ 0. Varianta 2 ξ=

1 d2 1 |q|2 |T | + . 2µ K θ

Pokud µ > 0 a K > 0, pak ξ > 0. Poznámka. θ= Gibbsova rovnice

∂e , ∂η

p = ρ2

∂e . ∂ρ

  1 , θ dη = de + p d ρ

lze z ní odvodit předchozí vztah pro tlak. Tato rovnice platí pro stlačitelné plyny, není důvod, proč by měla obecně platit pro komplikovanější materiály. Otázka: T a q?

Existuje způsob, jak ze znalosti konstitutivních vztahů pro η a pro ξ určit konstitutivní vztahy pro

9

Princip maximalizace rychlosti produkce entropie Našim cílem je nyní odvodit Navierovy-Stokesovy rovnice, nechť tedy platí, že ξ = Td · Dd + (m + p) div v +

q · ∇θ , θ

(1)

a vybereme konstitutivní vztah pro ξ˜ 2 ˜ d , div v, ∇θ) = 2µ|Dd |2 + 2µ + 3λ (div v)2 + K |∇θ| . ξ = ξ(D 3 θ

˜ vybereme tu maximální, tj. maximalizujeme ξ˜ přes všechny Dd , div v, ∇θ tak, Je mnoho možností, kdy ξ = ξ, aby bylo splněno (1) ˜ d , div v, ∇θ). max ξ(D (1)+Dd ,div v,∇θ

Dalo by se říci, že materiál dodržuje „princip lenostiÿ, snaží se deformovat co nejrychleji, aby pak mohl odpočívat. Tento princip nazýváme princip maximalizace rychlosti produkce entropie. Maximalizaci provedeme pomocí Lagrangeových multiplikátorů. Definujeme tedy Lagrangeovu funkci    ˜ d , div v, ∇θ) + λ ξ(D ˜ d , div v, ∇θ) − Td · Dd + (m + p) div v + q · ∇θ L = ξ(D , θ pro kterou nyní hledáme extrém, ∂L ∂L ∂L = 0, = 0, = 0. ∂Dd ∂ div v ∂∇θ Dostáváme ∂ ξ˜ = λTd , ∂Dd ∂ ξ˜ (1 + λ) = λ(m + p), ∂ div v q ∂ ξ˜ =λ , (1 + λ) ∂∇θ θ (1 + λ)

(2)

nyní vypočítejme λ tak, že násobme první rovnici Dd , druhou div v a třetí ∇θ 1+λ = λ

ξ˜ ∂ ξ˜ ∂Dd

·

Dd

+

∂ ξ˜ ∂ div v

div v +

∂ ξ˜ ∂∇θ

= · ∇θ

1 . 2

Dosadíme-li nyní do (2), provedeme derivace a dostáváme Navierovy-Stokesovy rovnice Td = 2µDd , 2µ + 3λ div v, m+p= 3 q ∇θ =K . θ θ

10

Domácí úkol č. 6

Pomocí principu maximalizace rychlosti produkce entropie max

(1)+Td ,m,q

˜ d , m, q), ξ(T

kde ξ˜ je dáno vztahem 1 d2 3 1 |q|2 ξ˜ = |T | + (m + p)2 + , 2µ 2µ + 3λ K θ odvoďte Navierovy-Stokesovy rovnice. Předpoklad maximalizace ξ˜ pro nestlačitelné materiály Pro izotermální nestlačitelné materiály platí tr D = div v = 0, θ = θ∗ = konst. Vazba pro ξ je ξ = T · D. Nejprve, nechť ξ˜ je dáno pomocí

(3)

ξ˜ = 2µ|D|2 .

Použijme opět princip maximalizace rychlosti produkce entropie max (3)+D+tr D=0

˜ ξ(D),

definujeme Lagrangeovu funkci   ˜ ˜ L = ξ(D) + λ1 ξ(D) − (T · D) + λ2 tr D a tu maximalizujeme vzhledem k D, (1 + λ1 )

∂ ξ˜ = λ1 T + λ2 I. ∂D

Vynásobme tuto rovnici D a získáváme 1 + λ1 T·D 1 = ˜ = , ∂ξ λ1 2 ∂D · D dále proveďme stopu na tutéž rovnici a získáváme −3λ2 = λ1 tr T ⇒

λ2 1 = − tr T. λ1 3

Máme tedy 1 tr T. 3 Odvoďme ještě nyní Navierovy-Stokesovy rovnice dalším způsobem. Nechť nyní T = 2µD −

1 d2 ξ˜ = |T | 2µ

(q = 0)

a maximalizujeme ˜ d ). max ξ(T

(4)+Td

Dostaneme

1 1 d T = D =D= 2µ 2µ d



 1 T − (tr T)I , 3

což je méně obvyklý, ale ekvivalentní zápis Navierových-Stokesových rovnic.

11

(4)

Cauchyův model konečné pružnosti Kinematické veličiny

∂χ , ∂X = LFK0 ⇒ BK0 = FK0 FT K0 ,

FK0 = ˙K F 0

kde BK0 je

levý Cauchy-Greenův tenzor, lze vypočítat jeho derivaci ˙ K = LBK + BK LT a tr B ˙ K = 2D · BK . B 0 0 0 0 0 Mějme nyní entropii η = η˜(e, BK0 ) ⇔ e = e˜(η, BK0 ), například neo-Hookeův materiál

µ (tr BK0 − 3). 2

e˜ = e0 (η) + Dosadíme za e do bilance energie

ρθη˙ = div(Tv + q) − div T · v − µBK0 · D = T · D − µBK0 · D + div q = (T − µBK0 ) · D + div q. Dostáváme ρη˙ + div

q θ

=

  1 q · ∇θ (T − µBK0 ) · D + . θ θ | {z } produkce entropie

Elastický materiál nedisipuje, proto produkce entropie musí být nulová, a tedy T = µBK0 , pro nestlačitelný materiál navíc platí T = −ΦI + µBK0 . 5. přednáška, 2. listopadu 2011 Příklad (varující): Máme proudění mezi dvěma rovnoběžnými deskami, hledáme řešení pro ustálené proudění NavierovýchStokesových rovnic, tekutina ulpívá na obou stěnách, tj. u(1) = u(−1) = 0, předpokládáme rychlostní pole ve tvaru v = (u(y), 0, 0). Mechanika V případě čistě mechanickém máme tyto rovnice div v = 0, div(v ⊗ v) = div T, T = −pI + 2µD. Definujme průtok Q=

1 2

Z

1

u(y) dy. −1

Řešení úlohy je pak u(y) =

3Q (1 − y 2 ). 2

12

Termodynamika V termodynamickém případě pak máme tyto rovnice div v = 0, div(v ⊗ v) = div T, T = −pI + 2µD, div(Ev) = div(Tv + q)

a teplo neprostupuje hranicí, q · n = 0 na ∂Ω ⇒ q2 (±1) = 0. Dále předpokládáme, že p = p(x, y, z), D = D(y), e = e(y), q = q(y), S = 2µD. Dostáváme

  u(−p + S11 )  S12 u Tv =  S13 u

a platí ∂p ∂S12 = ⇒ S12 = C1 y + C2 , p = C1 x, ∂y ∂x ∂S22 ∂p = ⇒ p nezávisí na y, ∂y ∂y ∂q2 ∂ ∂ − = ((−p + S11 )u) + (S12 u). ∂y ∂x ∂y

(5) (6) (7)

Dosadíme výsledek z (5) a (6) do (7) −

∂ ∂q2 = −C1 u + (S12 u) ∂y ∂y

a zintegrujeme od -1 do 1 Z

1

q2 (−1) − q2 (1) = −C1

u(y) dy + S12 (1)u(1) − S12 (−1)u(−1), −1

Použijeme okrajové podmínky a dostaneme Z

1

C1

u(y) dy = 0. −1

To znamená, že buď Z

1

C1 = 0, nebo

u(y) dy = 0. −1

První část, když C1 = 0, znamená, že u = 0. Druhá část znamená, že u0 musí změnit na [-1,1] alespoň dvakrát znaménko a musí platit 0 ≤ ξ = T · D = S12 u0 , přičemž S12 mění znaménko nejvýše jednou, což v součinu s u0 nedává nezáporné číslo a porušuje druhý zákon termodynamiky. Tedy u musí být identicky rovno nule. Závěr: Ustálené Poisseuillovo proudění není možné řešení mezi dvěma rovnoběžnými deskami, pokud desky jsou dokonale tepelně izolované a tekutina ulpívá na hranici.

Nenewtonské jevy a modely pro nenewtonské tekutiny Definice 2. Tekutina je nenewtonská, pokud její chování nelze popsat modelem Navierovy-Stokesovy tekutiny stlačitelné

T = −p(ρ, θ)I + 2µ(ρ, θ)D + λ(ρ, θ)(div v)I

nestlačitelné

T = −pI + 2µ(ρ, θ)D.

Tato definice není moc užitečná. Popíšeme nenewtonské jevy. 13

Nenewtonské jevy (1) Schopnost tekutiny zesílit/zeslabit (rychlost) smyk(u) (shear thinning, shear thickening) (2) Schopnost tekutiny zesílit tlak (pressure thickening) (3) Přítomnost aktivačního/deaktivačního kritéria (spojeného s napětím nebo s rychlostí smyku) (4) Přítomnost nenulových rozdílů normálových napětí v jednoduchém smykovém poli (normal stress differences) (5) Napěťová relaxace (stress relaxation) (6) Nelineární tečení (non-linear creep2 ) Pro první tři body platí, že vztah mezi smykovým napětím a rychlostí smyku není lineární. Uvažujme tekutinu proudící rychlostí v = (u(y), 0, 0) s tenzorem napětí ve tvaru T = −pI + S. Ad (1) Schopnost tekutiny zesílit/zeslabit (rychlost) smyk(u) (shear thinning, shear thickening) Newtonská tekutina Pro newtonskou tekutinu pak platí S12 = µD12 a D12 = u0 (y). V literatuře se používají různé symboly pro rychlost smyku |D12 |, např. κ, nebo γ. ˙ Pro smykové napětí |S12 | se používá jako značení σ, nebo τ . Graf závislosti |S12 | na |D12 | je v případě newtonské tekutiny lineární s koeficientem úměrnosti µ (viz obr. 1). Zadefinujme nyní zobecněnou viskozitu µg (|D12 |) =

|S12 (|D12 |)| . |D12 |

V případě newtonské tekutiny je µg = µ = konst., graf závislosti µg na γ˙ je na obr. 2.

µg

σ

µ

σ = µγ˙

0

γ˙

Obrázek 1: Newtonská tekutina, (smykové napětí / rychlost smyku)

0

γ˙

Obrázek 2: Newtonská tekutina, (zobecněná viskozita / rychlost smyku)

Nenewtonská tekutina Schopnost zesílit smyk (shear thickening, (dilatant fluids)) znamená, že funkce je superlineární (obr. 3). Zoběcněná viskozita µg je rostoucí funkce, na obr. 4 je příklad, kde zobecněná viskozita degeneruje, obvykle je zobecněná viskozita v počátku kladná. Schopnost zeslabit smyk (shear thinning, (pseudoplastic fluids)) znamená, že funkce je sublineární (obr. 5). Zoběcněná viskozita µg je klesající funkce, na obr. 6 je příklad, kde je zobecněná viskozita singulární v počátku, obvykle je zobecněná viskozita v počátku konečná. Obecněji není důvod, aby to byla pěkná rostoucí funkce3 . Nejznámější modely jsou power-law tekutiny T = −pI + 2µ∗ |D(v)|r−2 D(v). {z } | µg (|D|)=˜ µg (|D|2 )

Pro r = 2 se jedná o newtonskou tekutinu, pro r > 2 jde o shear thickening tekutinu a pro 1 < r < 2 o shear thinning tekutinu. 2 V literatuře se někdy mluví o creeping flow. Tím se myslí něco úplně jiného, jde o pomalé stokesovo proudění, kdy lze zanedbat konvektivní člen. 3 Nemusí jít ani o funkci.

14

µg

σ

zobecnˇen´ a viskozita degeneruje 0

γ˙

0

γ˙

Obrázek 3: Shear thickening, (smykové napětí / rychlost smyku)

Obrázek 4: Shear thickening, (zobecněná viskozita / rychlost smyku)

σ

µg zobecnˇen´ a viskozita je singul´ arn´ı

0

γ˙

Obrázek 5: Shear thinning, (smykové napětí / rychlost smyku) Domácí úkol č. 7

0

γ˙

Obrázek 6: Shear thinning, (zobecněná viskozita / rychlost smyku)

Uvažujte dva modely tekutin, kde vztah pro tenzor napětí je následující T = −pI + 2µ0 A + |D(v)|2


r−2 2

D,
r−2 T = −pI + 2µ0 + 2µ1 |D(v)|2 D,

kde A, µ0 , µ1 > 0. Určete, jaké nenewtonské jevy tekutiny vykazují v závislosti na různých r ∈ R. Domácí úkol č. 8

Uvažujte model tekutiny T = −pI + ν∞ +

ν0 − ν∞ (1 + Γ2 |D|2 )

1−n 2

D,

kde ν0 , ν∞ , Γ > 0. Určete, jaké nenewtonské jevy tekutina vykazuje pro různá n ∈ R a r ∈ R. Příklad tekutiny vykazující zeslabení smyku je krev, zesílení smyku vykazují např. barvy, různé omáčky nebo kečupy. U průmyslových produktů se uvádí viskozita jako kvalita výrobku. Ad (3) Přítomnost aktivačního/deaktivačního kritéria (spojeného s napětím nebo s rychlostí smyku). Tekutina začne téct až po dosažení kritické hodnoty napětí τ , anglicky yield stress. V případě následné lineární závislosti nazýváme tyto tekutiny Binghamovy tekutiny4 . Pokud je závislost nelineární, mluvíme o Herschel-Bulkleyho tekutinách. Standardní zápis je následující D = 0 ⇔ |S| ≤ τ, D T = −pI + 2τ +µ ˜g (|D|2 )D ⇔ |S| ≤ τ. |D| Pokud je µ ˜g konstantní, jde o Binghamovu tekutinu, pokud nekonstatní, jde o Herschel-Bulkleyho tekutinu (viz obr. 7). Dalším příkladem je tekutina, kde je odezva spojená s chemickými procesy. Při určité hodnotě |D12 | se tekutina uzamkne (angl. locking), viz obr. 8. 4 Nazýváme

tekutinou, ačkoli v definici tekutiny máme, že neumí udržet smykové napětí.

15

σ

σ

uzamknut´ı τ

τ 0

0

γ˙ Obrázek 7: Herschel-Bulkley

γ˙ Obrázek 8: Uzamkutí

Kromě těchto explicitních vazeb mezi smykovým napětím |S12 | a rychlostí smyku |D12 | lze uvažovat i zcela obecnou implicitní vazbu ve tvaru g(|S12 |, |D12 |) = 0. 6. přednáška, 9. listopadu 2011 Domácí úkol č. 9

Ukažte, že odezva zobecněné Binghamovy tekutiny |S| ≤ τ ∗ ⇔ D = 0, τ ∗D + 2µ(|D|2 )D |S| > τ ∗ ⇔ S = |D|

je ekvivalentní konstitutivní rovnici
2µ(|D|2 ) τ ∗ + (|S| − τ ∗ )+ D = (|S| − τ ∗ )+ S, + + kde τ ∗ > 0, µ(·) : R+ 0 → R0 a x = max{0, x}.

Kromě toho, že budeme tenzor S dosazovat do bilance hybnosti a rozhodovat o jeho hodnotě na základě velikosti tenzoru D, je možné definovat tenzorovou funkci G
G(D, S) = 2µ(|D|2 ) τ ∗ + (|S| − τ ∗ )+ D − (|S| − τ ∗ )+ S zobrazující d×d d×d d×d Rsym × Rsym → Rsym

a tu přidat k systému bilančních rovnic. Výhodou je, že G je spojitá a že není třeba používat žádné speciální nástroje. Nevýhodou, kterou za to zaplatíme, je, že se nám velmi zvětší systém rovnic. Domácí úkol č. 10

Projděte si pořádně všechny PDF soubory na stránkách kurzu

http://www.karlin.mff.cuni.cz/~malek/new/index.php/NDIR057_Mechanika_nenewtonovských_tekutin.

Ad (2) Schopnost tekutiny zesílit tlak (pressure thickening). Viskozita není konstantní, jako to je u newtonských tekutin5 (viz obr. 9), ale je funkcí tlaku p (viz obr. 10). Pozoruje se, že viskozita je jen rostoucí funkcí tlaku, ne klesající. Kolem roku 1890 navrhl Barus model µ(p) = µ0 exp(αp), α > 0. Za experimentální potvrzení viskozity závislé na tlaku byla udělena Nobelova cena. Tyto modely se používají například v ložiskách, kde jsou různě velké tlaky a viskozita se s nimi ohromně mění. 5 Už

Stokes ve své stati psal o závislosti viskozity na tlaku. Pro jednoduchost ovšem dále uvažoval pouze konstantní viskozitu.

16

µg

µg

µ0 µ0

0

p

0

Obrázek 9: Newtonská tekutina

p Obrázek 10: Viskozita závislá na tlaku

Stlačitelná vs. nestlačitelná tekutina Pro stlačitelnou tekutinu je p = p˜(ρ), pokud si budeme myslet, že tento vztah lze obrátit, ρ = ρ˜(p), pak zápis T = −˜ p(ρ)I + 2µ(ρ)D + λ(ρ)(div v)I je ekvivalentní se zápisem ˜ T = −pI + 2˜ µ(p)D + λ(p)(div v)I, což je model s viskozitou závislou na tlaku. I když je materiál stlačitelný (ovšem ne nějak extrémně), používá se nestlačitelný model s viskozitou závislou na tlaku. Modely, kde viskozita je funkcí tlaku, jsou obtížně jak numericky, tak analyticky. Ad (4) Přítomnost nenulového rozdílu normálových napětí v jednoduchém smykovém poli (normal stress differences). Definujme tři (jsme v trojrozměrném prostoru) viskometrické funkce: µ N1 := T11 − T22 N2 := T22 − T33

viskozita první rozdíl normálových napětí druhý rozdíl normálových napětí

V jednoduchém smykovém poli předpokládejme rychlostní pole v = (u(y), 0, 0). Newtonská tekutina Vypočítejme, jak vypadá tenzor napětí,     0 u0 (y) 0 −p µu0 (y) 0 1 0 u (y) 0 0 , T = −pI + 2µD = µu0 (y) −p 0 . D= 2 0 0 0 0 0 −p Vidíme, že N1 = N2 = 0, a tedy newtonská tekutina nevykazuje přítomnost nenulového rozdílu normálových napětí v jednoduchém smykovém poli. Nenewtonská tekutina Nenewtonské tekutiny, jako je třeba škrob rozpuštěný ve vodě, vykazují přítomnost nenulového rozdílu normálových napětí v jednoduchém smykovém poli. Působíme-li na tekutinu silou v jednom směru, tekutina začne reagovat i v jiném (obvykle kolmém) směru. S tímto jevem jsou spojované efekty: • Rod climbing (šplhání po tyči) • Die swell (rozšíření proudu při výtoku z trubice) Označme D průměr trubice a DE průměr tekutiny v místě největšího „nabobtnáníÿ tekutiny při výtoku z trubice. Pro newtonské tekutiny bylo experimentálně zjištěno, že nejvíce je DE = 1,13. D Pro nenewtonské tekutiny dochází až ke čtyřnásobnému rozšíření proudu a platí "    2 #1/6 DE 1 N1 = 0,13 + 1 + , D 2 2S w kde S značí smykové napětí a w označuje, že se jedná hodnotu na stěně. 17

• Proudění nakloněným kanálem V případě newtonské tekutiny je tekutina na povrchu hladká (volná hranice), v případě nenewtonské tekutiny se tekutina vybouluje nad povrchem. • Obrácené sekundární proudění Mějme kádinku s tekutinou, z vrchu ji těsně přikryjeme otáčivým víčkem. Tekutina přilne k víčku a otáčením víčka ji uvedeme do pohybu. V tekutině se vytvoři sekundární proudění, které lze vidět na řezu kádinky. V případě nenewtonské tekutiny má toto sekundární proudění opačný směr než v případě newtonské tekutiny. Tento směr je daný velikostí N2 6 . Ad (5) Napěťová relaxace (stress relaxation). Test napěťové relaxace je následující (viz obr. 11): Deformujeme v čase t = 0 materiál konstantní deformací ε a pak sledujeme, jak se chová smykové napětí σ.

σ

ε

ε0

???

⇒

0

0

t

t

Obrázek 11: Test napěťové relaxace Uvažujme dva základní materiály, jeden se chová jako lineární pružina a jeden jako lineární tlumič (dva soustředné válce o téměř stejném poloměru, kde meziválči je vyplněno newtonskou tekutinou o viskozitě µ). Lineární pružina je popsaná Hookeovým zákonem σ = Eε, kde E je modul pružnosti. Lineární tlumič je popsán vztahem σ = µε. ˙ Výsledek testu napěťové relaxace je pro lineární pružinu na obr. 12, pro lineární tlumič na obr. 13. V případě tlumiče je napětí singulární v nule.

σ

σ

Eε0

0

0

t Obrázek 12: Lineární pružina

t Obrázek 13: Lineární tlumič

Většina materiálů je kombinací dvou základních materiálů, na obr. 14 je výsledek testu napěťové relaxace pro materiál blízký viskoelastické pevné látce, na obr. 15 je výsledek testu napěťové relaxace pro materiál blízký viskoelastické tekuté látce. Ad (6) Nelineární tečení (non-linear creep). Test nelineárního tečení je následující (viz obr. 16): V čase t = 0 zatížíme materiál konstantním napětím σ, v čase t = t∗ napětí vypneme a pak sledujeme, jak se chová deformace ε. Výsledek testu nelineárního tečení je v případě lineární pružiny na obr. 17, v případě lineárního tlumiče na obr. 18. Opět, většina materiálů je kombinací dvou základních materiálů, na obr. 19 je výsledek testu nelineárního tečení pro materiál blízký viskoelastické pevné látce, na obr. 20 je výsledek testu nelineárního tečení pro materiál blízký viskoelastické tekuté látce. 6V

literatuře se uvádí, že N2  N1 , řádově je pak N2 ∼ 0,1N1 .

18

OBRÁZEK CHYBÍ

OBRÁZEK CHYBÍ

BUDE DOPLNĚN

BUDE DOPLNĚN

Obrázek 14: Viskoelastická pevná látka

Obrázek 15: Viskoelastická tekutá látka

ε

σ

σ0

???

⇒

0

t∗

0

t

t

Obrázek 16: Test nelineárního tečení

ε

ε

σ0 E

σ0 µ

0

t∗

0

t

t∗

Obrázek 17: Lineární pružina

Obrázek 18: Lineární tlumič

OBRÁZEK CHYBÍ

OBRÁZEK CHYBÍ

BUDE DOPLNĚN

BUDE DOPLNĚN

Obrázek 19: Viskoelastická pevná látka

Obrázek 20: Viskoelastická tekutá látka

t

Kromě modelů algebraických existují také modely rychlostního typu, integrálního typu a stochastického typu (více viz PDF soubor ’Od vodíku po asfalt’). 7. přednáška, 16. listopadu 2011

Princip objektivity a materiálová symetrie Cílem je určit, na čem závisí a jak konkrétně vypadá vztah pro tenzor napětí.

Princip objektivity Chceme vědět, funkcí jakých proměnných je např. tenzor napětí T = T(ρ, v, ∇v, ∇ρ, . . . ). Lze závislost na některých proměnných vyeliminovat? 19

Začneme mechanikou hmotného bodu a druhým Newtonovým pohybovým zákonem m

d2 x = F, dt2

provedeme-li transformaci polohy (pouze translaci) a času x∗ = x + vt,

t∗ = t − t0 ,

pak se nezmění tvar Newtonova zákona m

d2 x∗ = F∗ . dt2

Přidáme-li rotaci, pak x∗ = Q(x − x0 ) + vt,

t∗ = t − t0 ,

kde Q je ortogonální matice, QQT = QT Q = I. V této transformaci se opět nemění tvar Newtonových zákonů. Princip invariance Newtonových zákonů vůči rovnoměrnému pohybu byla popsána Galileem. V mechanice kontinua nás bude zajímat tenzor napětí. Tenzor napětí se z ryze geometrického důvodu musí transformovat jako T = QT T∗ Q, z Galileiho principu relativity platí Q = konst. Nyní si představme, že jsme v klidu a chceme znát sílu, která je potřeba k protažení pružiny o délku ∆l. Síla je pak rovna F = k(∆l). Nyní pokus zopakujeme s rotující pružinou (i nerovnoměrně Q = Q(t), úhlová rychlost ω = ω(t)). Na pravou stranu musíme přidat zdánlivé síly – odstředivou, Eulerovu a Corriolisovu. Změna polohy, tj. rychlost je rovna dx dx∗ = +ω×x dt dt a zrychlení d2 x∗ d2 x dω dx dx = + ×x+ω× +ω× + ω × (ω × x) = 2 2 dt dt dt dt dt d2 x dx dω = × x + 2ω × + ω × (ω × x). + dt2 dt dt Dosadíme do Newtonova zákona m

d2 x∗ = F∗ − m dt2



 dω dx × x + 2ω × + ω × (ω × x) . dt dt

Vidíme, že koeficient pružiny k je transformací nedotčený. Budeme tedy dále požadovat, aby se tvar nezměnil ani pro Q = Q(t) závislé na čase. Pro tenzor napětí tedy požadujeme T(x, t) = Q(t)T T∗ (x∗ , t∗ )Q(t), chceme tedy více než Galileův princip relativity, povolujeme navíc závislost na čase, nazveme princip objektivity. Co nám z toho tedy plyne? Zkusme nejprve, že tenzor napětí je funkcí T = T(ρ, v, ∇v). Hned na začátek můžeme říct, že tenzor napětí nemůže být funkcí rychlosti v. Uvažujme např. konstitutivní vztah T = −pI + 2µ(|v|)D. První pozorovatel je v soustavě pohybující se rychlostí v = v, druhý pozorovatel v soustavě pohybující se v = v + w, kde w = konst. Počítejme D∗ =


1
1 ∇x∗ v∗ + (∇x∗ v∗ )T = ∇x v + (∇x v)T = D, 2 2

neboť ∇x∗ v∗ = ∇x v a symetrická část tenzoru napětí je stejná v obou soustavách. Pak tedy máme dva tenzory napětí T = −pI + 2µ(|v|)D T∗ = −pI + 2µ(|v + w|)D,

20

a tedy jeden pozorovatel naměří jinou viskozitu než druhý a to porušuje Galileiho princip relativity (použili jsme zatím bez závislosti na čase). Nyní použijeme princip objektivity (Q = Q(t)). Gradient rychlosti rozdělíme na symetrickou a antisymetrickou část ∇v = D + W =


1
1 ∇v + (∇v)T + ∇v − (∇v)T , 2 2

každá část má význam, D má význam deformace (vzájemný posun částic), W je zbytek, který rotuje. x = Q(t)X,

QQT = I,

počítejme rychlost a její gradient dx ˙ = QX, dt

v= deformační gradient je pak roven F=

∂x =Q ∂X

a platí pro něj ˙ = LF = (∇v)F ⇒ ∇v = QQ ˙ T. F Dosaďme nyní do ∇v ∇v = D + W =

   1˙ T ˙ T )T + 1 QQ ˙ T − (QQ ˙ T )T = QQ ˙ T=W QQ + (QQ | {z } 2| 2 {z } ˙ T =−QQ

=0

a rotace kontinua je obsažena jen v antisymetrické části W. Ukážeme nyní, že T je funkcí D, ne celého gradientu T(ρ, ∇v) = QT T∗ (ρ∗ , (∇v)∗ )Q QT(ρ, ∇v)QT = T∗ (ρ∗ , (∇v)∗ ). Použijeme x∗ = Q(t)(x − x0 ) a pomocí derivace a algebraickými úpravami dostaneme ˙ T + Q∇vQT , ∇x∗ v∗ = QQ my ale chceme získat ∇x∗ v∗ = Q∇vQT . Rozdělením na symetrickou a antisymetrickou část ∇x ∗ v ∗ = D ∗ + W ∗ =

QDQT | {z }

symetrická část

˙ T. + QWQT + QQ {z } | antisymetrická část

zjistíme, že potřebujeme jen závislost na D T(ρ, v, ∇v) = T(ρ, D).

Materiálová symetrie Z materiálové symetrie vyplývá7 , že platí (tenzor napětí T = G(D)) G(D) = QT G(QDQT )Q. Nechť G je polynom G(A) = α0 I + α1 A + α2 A2 + α3 A3 + . . . , pak z materiálové symetrie platí G(A) = QT G(QAQT )Q = α0 I + α1 A + α2 A2 + α3 A3 + . . . , polynomy jsou tedy dobré funkce. Speciálně v R3 se nám výsledek ještě více zjednoduší. Připomeňme, že pro charakteristický polynom platí det(B − λI) = 0 3

−λ + i1 λ2 − i2 λ + i3 = 0, 7 Voda

ve sklenici se musí chovat stejně jako voda ve sklenici natočená o určitý úhel.

21

kde


1 (tr B)2 − tr(B2 ) , i3 = det B. 2 Navíc dle Cayley-Hamiltonovy věty platí, že matice B je kořenem svého polynomu i1 = tr B, i2 =

B3 = i1 B2 − i2 B + i3 I. Víme tedy, že B3 = f (i1 , i2 , i3 , B, B2 ), pak B4 = Bf (i1 , i2 , i3 , B, B2 ) = g(i1 , i2 , i3 , B, B2 ). Můžeme tedy každou mocninu vyšší než dva napsat pomocí i1 , i2 , i3 , B a B2 a platí G(A) = β0 (i1 , i2 , i3 )I + β1 (i1 , i2 , i3 )A + β2 (i1 , i2 , i3 )A2 . V dimenzi tři můžeme přesně nahradit nekonečný rozvoj polynomem druhého stupně. Zatím tedy víme T(ρ, D) = β0 (i1 , i2 , i3 , ρ)I + β1 (i1 , i2 , i3 , ρ)D + β2 (i1 , i2 , i3 , ρ)D2 . Jaké jsou důsledky tohoto tvaru tenzoru napětí? Co když připustíme jen lineární členy? Pak ponecháme jen i1 a i2 , i3 vypustíme a ponecháme jen lineární funkce D. T(ρ, D) = γ0 (ρ)I + γ1 (ρ)(tr D)I + γ2 (ρ)D máme stlačitelnou tekutinu, pokud budeme mít nestlačitelnou tekutinu, dostaneme T = γ0 I + γ2 D a přeznačíme-li koeficienty T = −pI + 2µD, dostaneme Navierovy-Stokesovy rovnice. Vraťme se zpátky k obecnému vztahu pro nestlačitelnou tekutinu T(D) = −pI + 2µ(D)D + 2µ(D)D + 2˜ µ(D)D2 . Nemůžeme ale takhle získat T = −pI + 2µ(p, D)D. Vraťme se zpátky k explicitnímu konstitutivnímu vztahu pro tenzor napětí T = G(D), lze napsat také jako implicitní vztah mezi tenzorem napětí T a symetrickou částí gradientu rychlosti D H(T, D) = 0. Provedeme-li celý postup jako doteď dostáváme α0 I + α1 T + α2 D + α3 (TD + DT) + α5 T2 + α6 D2 + α7 (T2 D + DD2 ) + α9 (D2 T + TD2 ) + α11 (D2 T2 + T2 D2 ) = 0, kde α0 , . . . , α11 je funkcí tr T, tr T2 , tr T3 , tr D, tr D2 , tr D3 , tr(TD), tr(T2 D), tr(TD2 ), tr(T2 D2 ). Označíme-li p = −(tr T)/3, dostaneme po linearizaci α0 I + α1 T + α2 D = 0. Domácí úkol č. 11

Ukažte, že v R3 platí
1 (tr B)2 − tr(B2 ) = (det B) tr(B−1 ). 2

22

Domácí úkol č. 12 Ukažte, že v R2 platí. Provedeme-li polární dekompozici deformačního gradientu na rotaci R a prodloužení (right stretch tensor) U F = RU, pak pro U platí

 √   1 C+ det C I , U= p √ tr C + 2 det C

kde C = FT F je pravý Cauchy-Greenův tenzor napětí. 8. přednáška, 23. listopadu 2011 Bude doplněno. Domácí úkol č. 13 přes Td ∈ A, kde

˜ d , Dd ) Použijte princip maximalizace rychlosti produkce entropie a maximalizujte ξ = ξ(T  d d A := Td ∈ R3×3 , sym , ξ = T · D

kde ξ˜ je dáno vztahy: 1. ξ˜ = |Td |r/(r−1) , 2. ξ˜ = |Dd |2−r |Td |2 . Domácí úkol č. 14

Uvažujte Stokesovu tekutinu s konstitutivním vztahem T = −pI + 2µD + αD2 ,

kde α a µ jsou konstanty. Určete, jaké nenewtonské jevy model je schopen zachytit a jaké není. Vysvětlete. 9. přednáška, 30. listopadu 2011 Obecněji. Pro nestlačitelný materiál platí tr D = 0, uvažujme dále isotermální proces θ = θ0 = konst. Rychlost produkce entropie v tom případě, jak již víme, splňuje ξ = Td · Dd a nechť máme obecný konstitutivní vztah pro rychlost produkce entropie ˜ ξ = ξ(m, Td , Dd ) ≥ 0 konvexní v Td a Dd . Z objektivity platí ˜ ˜ ξ(m, Td , Dd ) = ξ(m, QTd QT , QDd QT ) ∀Q ∈ SO3. Pak z Cayley-Hamiltonovy věty vyplývá ˆ ξ = ξ(m, tr D, tr D2 , tr D3 , tr Td , tr(Td )2 , tr(Td )3 , tr(Td D), tr((Td )2 D), tr(Td D2 ), tr((Td )2 D2 )) a použijeme princip maximalizace produkce entropie max

Td +ξ=Td ·Dd

ˆ . . ). ξ(.

Lagrangeova funkce má tvar ˆ . . ) + λ(ξ(. ˆ . . ) − Td · Dd ) L = ξ(. Zderivujeme ji a položíme rovnu nule 1 + λ ∂ ξˆ = Dd λ ∂Td a vypočítáme 1+λ ξˆ = λ ∂ ξˆ · Td ∂Td a dosadíme zpět do (8) Dd =

ξˆ

∂ ξˆ . ∂Td ∂ ξˆ d ·T ∂Td 23

(8)

Vypočítáme derivace a získáváme Dd = a0 Td + a1 Dd + a2 D2 + a3 Td Dd + a4 Td D2 + a5 (Td )2 + a6 I, kde ai je funkcí m a všech invariantů8 .

Tekutiny Kortewegova typu Ukážeme si, jak se dostat k modelu, který odvodil v roce 1901 holandský matematik Korteweg. Hranice mezi plynem a tekutinou není ostrá. Tento model je používán k zachycení popisu kapilárních efektů, vícefázových materiálů, granulovaných materiálů. Model je popsán následujícími rovnicemi ρ˙ = −ρ div v, ρv˙ = div T, T = −p(ρ)I + 2µ(ρ)D(v) + λ(ρ)(div v)I + K, κ K = (∆(ρ2 ) − |∇ρ|2 )I − κ(∇ρ ⊗ ∇ρ). 2 Ve článku Dunn, Serrin (1985) je ukázáno, že tento model je termodynamicky konzistentní. Víme, že pro stlačitelné Navierovy-Stokesovy-Fourierovy rovnice, pokud předpokládáme konstitutivní vztah η = η˜(e, ρ) dostaneme ρη˙ + div

q θ

  1 q · ∇θ d d T · D + (m + p) div v + = θ θ   1 q · ∇θ dis d d = Tdis · D + tdis div v + , θ θ

kde zavedeme disipativní členy (s indexem dis) a určíme konstitutivní vztah pro rychlost produkce entropie 2 1 ˜ d , tdis , qdis ) = 1 |Td |2 + t2 + |qdis |2 ξ = ξ(T dis 2µ dis 2µ + 3λ dis K a pomocí principu maximalizace rychlosti produkce entropie dostaneme Navier-Stokes-Fourierův systém. Domácí úkol č. 15 Uvažujte Kortewegův model v jednoduchém smykovém poli v = (u(y), 0, 0). Ukažte, že rozdíly normálových napětí jsou nenulové T11 − T22 6= 0, T22 − T33 6= 0. Pro Kortewegův model na začátku uvažujeme entropii navíc závislou na gradientu hustoty η = η˜(e, ρ, ∇ρ) ↔ e = e˜(η, ρ, ∇ρ). Derivujme ρe˙ = ρ

∂˜ e ∂˜ e ˙ ∂˜ e η˙ + ρ ρ˙ + ρ ∇ρ, ∂η ∂ρ ∂∇ρ

˙ víme, že ρ˙ = −ρ div v, potřebujeme vypočítat ∇ρ ˙ = −∇(ρ div v) − (∇v)∇ρ. ∇ρ Označme θ= 8 Porovnejte

∂˜ e ∂˜ e , p = ρ2 ∂η ∂ρ

výsledek s obecným vztahem G(Td , Dd , m) = 0.

24

máme tedy ρ˙ ∂˜ e ˙ ρθη˙ = ρE˙ − ρv˙ · v − p − ρ ∇ρ = ρ ∂∇ρ   ∂˜ e ∂˜ e = div(Tv + q) − div T · v + p div v + ρ · ∇(ρ div v) + ρ ⊗ ∇ρ · ∇v = ∂∇ρ ∂∇ρ       ∂˜ e ∂˜ e ∂˜ e = T · ∇v + ρ ⊗ ∇ρ · ∇v + div q + p div v + div ρ2 div v − ρ(div v) div ρ . ∂∇ρ ∂∇ρ ∂∇ρ Princip objektivity říká, že e je funkcí velikosti ∇ρ e = e˜(η, ρ, |∇ρ|), z čehož vyplývá, že ∂˜ e ⊗ ∇ρ ∂∇ρ je symetrický tenzor. Dostáváme tedy      ∂˜ e ∂˜ e ρθη˙ = T + ρ ⊗ ∇ρ · D + p − ρ div ρ div v+ ∂∇ρ ∂∇ρ    d ! ∂˜ e ∂˜ e div q + ρ2 (div v) = Td + ρ ⊗ ∇ρ ·Dd + ∂∇ρ ∂∇ρ {z } | Td dis

      ∂˜ e ∂˜ e e ρ ∂˜ 2 . p − ρ div ρ · ∇ρ div v + div q + ρ (div v) +m+ ∂∇ρ 3 ∂∇ρ ∂∇ρ {z } {z } | | tdis

Dostáváme rovnici ρη˙ − div

q

dis

θ



=

qdis

  qdis · ∇θ 1 Tddis · Dd + tdis div v + . θ θ

ξ = 0 pokud Tddis = 0, tdis = 0, qdis = 0, pak d     ∂˜ e ρ ∂˜ e ∂˜ e ⊗ ∇ρ − p(. . . )I − · ∇ρ I + ρ div ρ I ∂∇ρ 3 ∂∇ρ ∂∇ρ     ∂˜ e ∂˜ e = −p(e, ρ, ∇ρ)I + ρ div ρ I−ρ ⊗ ∇ρ . ∂∇ρ ∂∇ρ

T = Td + mI = −ρ



Speciálně volbou e = e˜(η, ρ, ∇ρ) = e0 (η, ρ) +

β ∂˜ e |∇ρ|2 ⇒ ρ = β∇ρ 2ρ ∂∇ρ

dostaneme Kortewegův model T = −p(. . . )I + ρβ div ∆ρI − β(∇ρ ⊗ ∇ρ)
β ∆(ρ2 ) − |∇ρ|2 I − β(∇ρ ⊗ ∇ρ). = −pI + 2 Volbou konstitutivní rovnice pro rychlost produkce entropie 3 1 ˜ d , tdis , qdis ) = 1 |Td |2 + ξ = ξ(T t2dis + |qdis |2 dis dis 2µ 2µ + 3λ K a volbou9 Tddis = 2µDd 2µ + 3λ tdis = div v 3 qdis = K∇θ, 9 Je možné použít princip maximalizace produkce entropie, pro jednoduchost snadno volíme díky tomu, že rychlost produkce entropie je kvadratická.

25

vidíme, že ξ ≥ 0 a je splněn druhý zákon termodynamiky. Pak dostáváme ∂˜ e = K∇θ − βρ(div v)∇ρ ∂∇ρ
β T = −pI + 2µD + λ(div v)I + ∆(ρ2 ) − |∇ρ|2 I − β(∇ρ ⊗ ∇ρ) 2 q = K∇θ − ρ2 div v

Kortewegův model. 10. přednáška, 7. prosince 2011

Modelování vazkopružných (viskoelastických) materiálů Chceme modelovat materiály schopné zachytit napěťovou relaxaci a nelineární creep (tečení), příp. rozdíl normálových napětí v jednoduchém smykovém poli. O jaké materály jde? Jde o geofyzikální materiály, biologické materiály, polymery a jiné chemické produkty, potraviny. Budou se kombinovat vlastní tekutiny a pevné látky, např. v plazmě se nacházejí pevné částečky a ty dodávají elasticko-plastický charakter. Dalším příkladem je asfalt, jeho vlastnosti závisí velmi na teplotě. Jde o směs materiálu s tekutou amorfní částí a makromolekul, které dodávají pružnost. Zopakujme si, co to je test napěťové relaxace a nelineárního tečení. Test napěťové relaxace Test napěťové relaxace je následující (viz obr. 21): Deformujeme po čás t materiál konstantní deformací ε, a pak sledujeme, jak se chová smykové napětí σ. Pro lineární pružinu platí

σ

ε

ε0

???

⇒

0

0

t

t

Obrázek 21: Test napěťové relaxace σ = Eε, pro lineární tlumič σ = µε. ˙ Pro lineární pružinu a tlumič pak dostáváme Definujme funkci napěťové relaxace

σ

σ

Eε0

0

0

t Obrázek 22: Lineární pružina

t Obrázek 23: Lineární tlumič

G(t) = reálné materiály vykazují následující chování 26

σ(t) , ε0

OBRÁZEK CHYBÍ

OBRÁZEK CHYBÍ

BUDE DOPLNĚN

BUDE DOPLNĚN

Obrázek 24: Viskoelastická pevná látka

Obrázek 25: Viskoelastická tekutá látka

Creep test (pro tečení) Test nelineárního tečení je následující (viz obr. 26): V čase t = 0 zatížíme materiál konstantním napětím σ a v čase t = t∗ napětí vypneme, a pak sledujeme, jak se chová deformace ε.

ε

σ

σ0

???

⇒

0

t∗

0

t

t

Obrázek 26: Test nelineárního tečení Výsledek testu napěťové relaxace je v případě lineární pružiny na obr. 27. V případě lineárního tlumiče na obr. 28.

ε

ε

σ0 E

σ0 µ

0

t∗

0

t

Obrázek 27: Lineární pružina

t∗ Obrázek 28: Lineární tlumič

Definujme funkci nelineárního creepu J(t) =

ε(t) . σ0

Pro reálné materiály je výsledek následující OBRÁZEK CHYBÍ

OBRÁZEK CHYBÍ

BUDE DOPLNĚN

BUDE DOPLNĚN

Obrázek 29: Viskoelastická pevná látka

Obrázek 30: Viskoelastická tekutá látka

Původní přístupy k modelování Modely začali zkoumat Maxwell, Burgers, Boltzmann, Kelvin, Voigt, Oldroyd, . . .

27

t

Maxwell

Maxwell dospěl k modelu

dσ dε σ =E − , dt dt λ kde parametr λ má rozměr času, pošleme-li λ → ∞ dostaneme vzorec pro lineární pružinu, pošleme-li λ → 0 a Eλ → µ dostáváme Newtonskou tekutinu. Meyer (Kelvin-Voigt) dε , dt v každém bodě tělesa máme přítomnou elastickou i vazkou část. σ = Gε + η

Boltzmann Z

t

σ(t) =

G(t − t0 )

−∞

dε 0 (t ) dt0 , dt

kde materiál má paměť skrz G.

Tvorba jednoduchých modelů kombinací lineárních pružinek a lineárních tlumičů Maxwellův prvek Maxwellův prvek vznikne sériovým zapojením lineární pružiny (linear spring) a lineárního tlumiče (linear dashpot) Pro lineární pružinu platí

Spring

Dashpot

l Obrázek 31: Maxwellův prvek

FS = E∆S , tedy σS = EεS , pro lineární tlumič ˙ D , tedy σD = µε˙D . FD = µ∆ Pro Maxwellův prvek provedeme výpočet podrobně. Nejprve odvodíme konstitutivní vztah. V sériovém zapojení je síla v pružině i tlumiči stejná, tj. σD = σS , celkové prodloužení je součtem prodloužení pružiny a tlumiče, tj. ∆ = ∆S + ∆D . Dosadíme do konstitutivních vztahl pro lineární pružinu a lineární tlumič ˙ ˙ =∆ ˙S+∆ ˙ D = FS + FD , ∆ E µ a získáváme konstitutivní vztah pro Maxwellův prvek E ε˙ = σ˙ +

E σ ⇔ p1 σ˙ + p0 σ = q1 ε. ˙ µ

Z počáteční podmínky získáme p1 σ(0+) = q1 ε(0+). Nyní počítejme, jak se chová Maxwellův prvek v creep testu, nechť tedy σ(t) = σ0 , pak ε(t) =

σ0 (p1 + p0 t) q1

a creep funkce je rovna J(t) =

ε(t) 1 = (p1 + p0 t). σ0 q1

Na obrázku ?? lze vidět, že odezva je nespojitá, ovšem lineární, to není uspokojivý výsledek pro Maxwellův prvek. 28

Dále počítejme, jak se chová Maxwellův prvek v testu napěťové relaxace, nechť tedy ε(t) = ε0 H(t), kde H(t) je Heavysidova funkce, pak dostáváme p0 q1 σ(t) = ε0 e− p1 t p1 a funkce napěťové relaxace je rovna G(t) =

σ(t) q1 p0 = e− p1 t . ε0 p1

Na obrázku ?? lze vidět, že odezva vypadá pěkně, a tedy Maxwellův prvek dává uspokojivý výsledek. Nyní vypočítáme, jak závisí napětí σ na deformaci ε. Upravujme 

˙ p0  q 1 p0 t σ(t)e p1 t = εe ˙ p1 p1

p

− p0 t

q1 + p1

Z

t

p0

− p (t−τ ) 1 ε(t)e ˙ dτ 0 Z p p0 q1 t poč. podmínka q1 − p0 (t−τ ) 1 ε(t)e ˙ ε(0+)e− p1 t + dτ = p1 p1 0 Z t q1 = ε(0+)G(t) + ε(t)G(t ˙ − τ ) dτ. p1 0

σ(t) = σ(0+)e

1

Ve vyšší dimenzi tento výsledek odpovídá tomu, že jsme dostali integrální vazbu mezi T a D. Počítejme nyní duální zápis a vypočítejme, jak ε závisí na σ Z t Z t q1 ε(t) = q1 ε(0+) + p1 σ(τ ˙ ) dτ + p0 σ(τ ) dτ 0 0 Z t = q1 ε(0+) + σ(τ ˙ )(p1 + p0 (t − τ )) dτ + p0 σ(0+)t. 0

Použitím počátečních podmínek pak získáváme     Z t p1 p0 p0 p1 ε(t) = + t σ(0+) + + (t − τ ) dτ σ(τ ˙ ) q1 q1 q1 q1 0 Z t = J(t)σ(0+) + σ(τ ˙ )J(t − τ ) dτ. 0

Dostali jsme model rychlostního typu (rate-type model) a dva modely integrálního typu, všechny ekvivalentní díky tomu, že jsme studovali lineární modely. Všechny modely platí v jedné dimenzi, není ale zřejmé, jak zobecnit tyto modely do tří dimenzí. 1. Lineární tlumič není aproximativní teorie, zatímco lineární pružnost je aproximativní teorie! 2. Záleží, odkud vyjdeme, zobecňujeme model rychlostního typu, nebo integrální model? Jde o zobecnění lineární, nebo nelineární? Různými přístupy dostáváme různé výsledky, které již nejsou ekvivalentní. 3. Co dělat, pokud tlumič či pružina budou záviset na deformaci nelineárním způsobem? 4. Parciální časová derivace, ani materiálová časová derivace nejsou objektivní, zobecnění už nejsou jednoznačné. 5. Mnoho nelineárních 3D modelů se mohou redukovat v 1D na stejnou rovnici. K zachycení jednorozměrného experimentu jich lze tedy použít mnoho. Kelvinův-Voigtův prvek Kelvinův-Voigtův prvek vznikne paralelním zapojením lineární pružiny (linear spring) a lineárního tlumiče (linear dashpot) Deformace pružiny i tlumiče je stejná, celková síla je rovna součtu sil v pružině i tlumiči F = FS + FD , ∆ = ∆S = ∆D . Použitím konstitutivních vztahů pro lineární pružinu a lineární tlumič dostáváme p0 σ = q0 ε + q1 ε˙ ε(0+) = 0.

29

Spring

Dashpot l Obrázek 32: Kelvinův-Voigtův prvek

Funkce napěťové relaxace je rovna G(t) =

q0 q1 + δ(t), p0 p0

nevýhoda pro Kelvin-Voigt. Funkce nelineárního creepu q0  p0  1 − e− q 1 t , J(t) = q0 výhoda pro Kelvin-Voigt. Domácí úkol č. 16 Pro Kelvinův-Voigtův prvek vypočítejte funkci napěťové relaxace a funkci creepu. Dále odvoďte, jak napětí σ závisí na deformaci ε, a naopak, jak deformace ε závisí na napětí σ. Domácí úkol č. 17 Uvažujte prvek se sériovým zapojením Kelvinova-Voigtova prvku a lineárního tlumiče jako na obrázku 33. Odvoďte konstitutivní vztah pro tento prvek včetně počátečních podmínek.

E1

µ1 µ2 Obrázek 33: Trojprvek 11. přednáška, 14. prosince 2011 Až doteď jsme studovali případ, kde pružiny a tlumiče splňují lineární závislost deformace na napětí. Předpokládejme, že závislost je nelineární, např. µ = µ(ε), ˙ pak dostáváme µ(ε) ˙

σ˙ + σ = µ(ε) ˙ ε, ˙ E

není snadné vypočítat, jak vypadá relaxační funkce. Všechny ukázáné modely, byly jednodimenzionální, ukážeme přechod do vyšší dimenze. Předpokládejme nestlačitelnou tekutinu tr D = 0. Uvažujme tenzor napětí ve tvaru T = −pI + S a řekněme, že S splňuje stejnou rovnici, co jsme odvodili dříve, zkusme tedy S˙ + λS = 2µD. Vezměme vztažnou soustavu x∗ = Q(x − x0 ) + c, 30

kde Q je ortogonální, chceme, aby se tenzor S transformoval S∗ = QSQT . Ve hvězdičkované soustavě materiál splňuje rovnici ∗ S∗ + λS˙ = 2µD∗ ,

což je rovno QSQT + λ


d QSQT = 2µQDQT , dt

kdyby platilo QSQT + λ

 
d QSQT = Q S + λS˙ QT , dt

je to v pořádku, ale to neplatí, protože
d T ˙ ˙ T + QSQ ˙T QSQT = QSQ + QSQ dt ale součet prvního a posledního členu nemůže být nikdy roven nule. Tedy návrh, který jsme zkusili, je špatný – ◦

není objektivní. Označme objektivní derivaci S znakem S splňující ◦

◦

S∗ = Q S QT . Definujme nyní dS − LS − SLT + (tr L)S. dt Ověřme, že tato derivace je objektivní. Chceme   ◦ ◦ dS T T ∗ S =QSQ =Q − LS − SL + (tr L)S QT . dt ◦

S=

Počítejme ◦ dS∗ − L∗ S∗ − S∗ (LT )∗ + (tr L∗ )S∗ =QSQT = dt

=


d QSQT − L∗ QSQT − QSQT (LT )∗ + (tr L∗)QSQT dt

Potřebujeme vědět, co je L∗ , už víme, že je rovno ˙ T + QLQT . L∗ = QQ Pokračujme ve výpočtu
d QSQT − L∗ QSQT − QSQT (LT )∗ + (tr L∗)QSQT = dt   T ˙ ˙ T + QSQ ˙ T − QQ ˙ T + QLQT QSQT − = QSQ + QSQ     ˙ T + QLQT + tr QQ ˙ T + QLQT QSQT = QSQT QQ   ˙ T − QLSQT − QSLT QT + tr QQ ˙ T + QLQT QSQT = = QSQ   = Q S˙ − LS − SLT QT + Q(tr L)SQT , neboť platí cykličnost stopy a ˙ T + QQ ˙ T. QQ Derivace je objektivní, i když v ní nebude poslední část se stopou O

S=

dS − LS − SLT . dt

31

Domácí úkol č. 18

Ukažte, že pokud je bilance celkové energie      ∂ 1 1 2 2 ρe + ρ|v| + div v ρe + ρ|v | = div(Tv) + div q + f · v ∂t 2 2

invariantní vůči transformaci x∗ = x + wt,

t∗ = t

pak systém automaticky splňuje bilanci hmoty, hybnosti a bilanci pro vnitřní energii. Hint: Ohvězdičkujte systém, dosaďte, zderivujte a porovnejte členy stejné mocniny. Pokračujme dále s objektivní derivací O

S=

dS − LS − SLT dt

a Maxwellovým modelem O

T = −pI + S,

S + λ S= 2µD.

Mohlo by se zdát, že v případě ustáleného proudění vypadne derivace a dostane se Newtonský model. Ustálené proudění ve tvaru v = (u(y, t), 0, 0) a S = S(y, t), dosadáme-li do Maxwellova modelu, dostaneme   ∂u ∂u ∂Syx Syx + λ − Syy = ∂t ∂y ∂y ∂Syy Syy + λ =0 ∂t Dostaneme tedy pro smykové napětí vztah, který známe z jednodimenzionálního problému Syx + λ

∂u ∂Syx = ∂t ∂y

a pro normálové složky Sxx = 2λSxy

∂u ∂y

Syy = 0 Szz = 0 a model nám umožňuje zachytit rozdíl normálových napětí v jednoduchém smykovém poli. Domácí úkol č. 19 Proveďte podrobně výpočet pro ustálené proudění ve tvaru v = (u(y, t), 0, 0) a S = S(y, t) s použitím Maxwellova modelu ve třech dimenzích T = −pI + S, O

S + λ S = 2µD.

Meření viskozity Viskozity různých tekutin se pohybují od 10−2 do 1025 , je tedy zřejmé, že viskozitu musíme měřit různými metodami pro různé viskozity. Pokud uvažujeme model s viskozitou závislou na rychlosti smyku, zaleží na tom, v jakém rozmezí experimentátoři viskozitu naměřili. Coutteův viskozimetr (19. století) Předpokládá se, že máme model, kterým chceme materiál popsat a podle toho viskozimetrem měříme. Pro Newtonovské modely máme analytický vztah, pro Nenewtonské tekutiny nemusí existovat analytické řešení. Pressure hole error Cup viscometer Plate cylinder Falling cylinder viscometer Earth’s mantle

32

12. přednáška, 21. prosince 2011

Nestlačitelné tekutiny rychlostního typu – odvození termodynamicky konzistentních modelů Motivace a cíl Kombinací dvou lineárních tlumičů a jedné lineární pružiny, dostaneme rovnici p1 T˙xy + p0 Txy = q1 D˙ xy + q0 Dxy . Jak zobecnit tento vztah do 3D? Řešení není triviální z následujících důvodů: i) V odvození byly využity jen lineární modely. ii) Pro objektivní tenzor platí, že ani obyčejná (parciální) časová derivace, ani materiálová časová derivace nejsou objektivní. Bylo zavedeno mnoho objektivních derivací. Není jasné, kterou objektivní derivaci zvolit. iii) Odvozený model je vhodný k zachycení napěťové relaxace a nelineárního creepu. Jak zahrnout další nenewtonské jevy do 3D modelu? iv) Yeleswerapu model pro popis vlastnosti krve div v = 0 ρv˙ = div T + ρb S + λ1



T = −pI + S  ˙ − LD − DLT ) S˙ − LS − SLT = µ(|D|)D + λ2 (D   1 + ln(1 + Λ|D|) µ(|D|) = µ∞ + (µ0 − µ∞ ) 1 + Λ|D|

Je tento model dobrý? Je termodynamicky konzistentní? Nástroje K odvozování bude použito: • implicitní konstitutivní teorie • maximalizace produkce entropie • koncept přirozené konfigurace přiřazené k současné konfiguraci Připomenutí

Mámé bilanční rovnice ρ˙ = −ρ div v ρv˙ = div T, T = TT ρE˙ = div(Tv + q)

Kelvin-Voigtův model Předpokládejme entropii tvaru η = η˜(e, ρ, BkR ) ⇔ e = e˜(η, ρ, BkR ) = eˆ(η, ρ, tr BkR ). Definujme kinematické veličiny. Deformační gradient FkR =

∂χkR ∂X

a levý a pravý Cauchy-Greenův tenzor napětí CkR = FT kR FkR .

BkR = FkR FT kR , Počítejme

∂ˆ e ∂ˆ e ∂ˆ e ˙k . I·B ρE˙ − ρv˙ · v = ρe˙ = ρ η˙ + ρ ρ˙ + ρ R ∂η ∂ρ ∂ (tr BkR ) 33

˙ k = LFk , z toho vyplývá, že Víme, že F R R ˙ k = LBk + Bk LT ⇒ I · B ˙ k = 2Bk · D. B R R R R R Pokračujme ve výpočtu ∂ˆ e Bk · D = ∂ (tr BkR ) R  d   ∂ˆ e 2 ∂ˆ e · Dd + m − ρ = T − 2ρ BkR tr BkR + p div v. ∂ (tr BkR ) 3 ∂ (tr BkR )

ρθη˙ − div q = T · D + p div v − 2ρ

Podělíme teplotou a máme ρη˙ − div

q θ

=

  1 q · ∇θ . Tddis · Dd + tdis div v + θ θ {z } | ξ

Nestlačitelnost a konstantní teplota dává ξ = Tddis · Dd =

 T − 2ρ

∂ˆ e Bk ∂ (tr BkR ) R

d

· Dd .

˜ dis , D) volme konstitutivní vztahy buď Dále obecně pro ξ = ξ(T 1 |Tdis |2 ξ˜ = 2ν nebo ξ˜ = 2ν|D|2 . V prvním případě maximalizujeme ξ vzhledem k Tdis bez další vazby, v druhém případě maximalizujeme vzhledem k D s vazbou nestlačitelnosti tr D = 0. Maximalizací dostáváme T = −mI + 2νD + 2ρ

∂ˆ e Bd . ∂ (tr BkR ) kR

Pokud pro volnou energii platí (neo-Hookeův materiál) ρψ =

µ (tr BkR − 3) , 2

pak platí ∂ˆ e µ = I ∂ (tr BkR ) 2 a dostáváme Kelvin-Voigtův model T = −mI + 2νD + µBdkR . Pokud je nulová produkce entropie, pak Tddis = 0 1 p + m + µ tr BkR = 0 3 q=0 a dostáváme T = Td + mI = Tddis + µBdkR + mI = −pI + µBkR . Nevyžadujeme-li nestlačitelnost a použijeme vztah jako v NSF Tddis = 2νDd 1 2ν + 3λ p + m + µ tr BkR = div v 3 3 q = k∇θ, dostaneme T = −pI + 2νD + λ(div v)I + µBkR q = k∇θ. 34

Viskoelastický model (K.R. Rajagopal, A.R. Srinivasa, JNFM 2000) Těleso tvořené viskoelastickým materiálem se nachází v čase t v konfiguraci κ(t), tato konfigurace se nazývá aktuální (současná) konfigurace. Tuto konfiguraci budeme porovnávat s konfigurací systému κ(r), což je konfigurace systému v klidu na počátku, tzv. referenční konfigurace. Zavedeme přirozenou konfiguraci κ(p(t)). Je to taková konfigurace systému odpovídající aktuální konfiguraci κ(t), při které je těleso oproštěno od všech vnějších podnětů. Pak se aktuální konfigurace a jí takto odpovídající přirozená konfigurace od sebe liší čistě elastickou deformací. Zavedením pojmu přirozená konfigurace tedy rozdělíme proces na část čistě elastickou a čistě viskózní. a)

b)

c)

Obrázek 34: Je dnoduchý příklad přirozené konfigurace Ukažme si na jednoduchém konkrétním příkladě, jak lze motivovat pojem přirozené konfiguraci. Jednoduchý lineární viskoelastický model lze například znázornit sériově zapojenou pružinou a tlumičem jako na obr. 34. Tlumič tvoří dva souosé válce vyplněné tekutinou. Pružina představuje elastickou část a tlumič viskózní část. Na obr. 34a) je na počátku materiál relaxován a pružina i tlumič nejsou natažené. Systém se nachází v referenční konfiguraci κ(r). Poté zapojení natáhneme (obr. 34b)), dojde k deformaci a systém se ocitne v aktuální konfiguraci κ(t). Pak soustavu necháme zrelaxovat, pružina se stáhne do svého klidového stavu, tlumič zůstane natažený (obr. 34c)). Konfiguraci tohoto systému nazveme přirozenou konfigurací κ(p(t)). Protože při relaxaci se deformace tlumiče nemění, budeme uvažovat, že celková volná energie systému je závislá jen na deformaci z přirozené konfigurace.

F

κ(r)

κ(t)

Fp(t)

G κ(p(t))

Obrázek 35: Přirozená konfigurace Zavedeme deformační gradient (na obr. 35 označeno F) FkR = ∇χ zobrazující infinitezimální element z κ(r) do κ(t). Ze znalosti deformačního gradientu můžeme určit následující veličiny. Levý a pravý Cauchyho-Greenův tenzor napětí BkR = FkR FT kR ,

CkR = FT kR FkR , 35

Dále definujme Fkp(t) jako zobrazení infinitezimálního elementu z κ(p(t)) do κ(t) (na obr. 35 označeno Fp(t) ) a G zobrazuje z κ(r) do κ(p(t)). Tyto dvě nové veličiny mají tedy vůči sobě vztah G = F−1 kp(t) FkR . Zaveďme dále levý a pravý Cauchyho-Greenův tenzor napětí Bkp(t) = Fkp(t) FT kp(t) ,

Ckp(t) = FT kp(t) Fkp(t)

a rychlost deformace (v analogii s obvyklým vztahem mezi gradientem rychlosti a deformačním gradientem), resp. její symetrickou část Lkp(t) + LT kp(t) ˙ −1 , Dk . Lkp(t) = GG = p(t) 2 Počítejme, čemu je rovna materiálová derivace ˙T = ˙k ˙ k FT + Fk F B =F kp(t) kp(t) p(t) p(t) p(t) ˙ −T ˙ T ˙ k G−1 FT + Fk G˙−1 FT + Fk G−T F FT FkR = kp(t) kp(t) kR + Fkp(t) G R R p(t) LBkp(t) + Bkp(t) LT − 2Fkp(t) Dkp(t) FT kp(t) , kde jsme použili vztahu ˙ ˙ −1 . A−1 = −A−1 AA Objektivní derivace se nám objeví přirozeným způsobem O

˙k − LBkp(t) − Bkp(t) LT = −2Fkp(t) Dkp(t) FT Bk(p(t)) = B kp(t) . p(t) O

Všimněme si, že I = −2D. 13. přednáška, 4. ledna 2012 Domácí úkol č. 20 Proveďte podrobně odvození pro nestlačitelný Kelvinův-Voigtův-Fourierův model pro tři varianty α) det BkR = 1 β) div v = 0 γ) det BkR = 1 a zároveň div v = 0. Platí ˙k tr B = 2Bkp(t) · D − 2Fkp(t) Dkp(t) · Fkp(t) p(t) Volme entropii η = η˜(e, ρ, Bkp(t) ), můžeme invertovat a používat vnitřní energii e = e˜(η, ρ, Bkp(t) ) a speciálně volme
µ tr Bkp(t) − 3 e = e0 (η, ρ) + 2ρ a jako u Kelvin-Voigtova modelu dostaneme q · ∇θ + µFkp(t) Dkp(t) · Fkp(t) θ Pro jednoduchost q = 0 a div v = 0 a v analogii s Rajagopal, Srinivasa (2000) předpokládáme, že det Bkp(t) = 1. Z toho a nestlačitelnosti FkR plyne, že det G = 1. Derivací det G získáváme tr Dkp(t) = 0. ξ = (T − µBkp(t) ) · D + (p + m) ˜ div v +

Domácí úkol č. 21 Ukažte, že pokud platí det G = 1, pak platí tr Dkp(t) = 0. Pro produkci entropie tedy platí (využijeme nulovosti stop D a Dkp(t) ) ξ = (T − µBkp(t) )d · D + µCdkp(t) · Dkp(t) = Tddis · D + µCdkp(t) · Dkp(t) . V analogii s předchozími zkušenostmi můžeme volit produkci entropie ˜ ξ = ξ(D, Dkp(t) ) = 2ν0 |D|2 + 2ν1 |Dkp(t) |2 nebo třeba taky ˆ d , S) = 1 |Td |2 + 1 |S|2 , ξ = ξ(T dis 2ν0 dis 2ν1 d kde S = µCkp(t) . Konstitutivní rovnici pro T dostaneme maximalizací produkce entropie. Získaný model není „pěknýÿ. 36

Poznámky: Pokud ν0 = 0, dostaneme zobecněný Maxwellův model, pokud ν0 > 0 dostaneme zobecněný Oldroydův model. Klasický Maxwellův model a Oldroyd-B model se dostanou linearizací elastické odezvy, předpokládáme, že Bkp(t) = I + A, kAk = O(ε). Je důležité, že zde elasticita je konečná pružnost. Dřívější přístupy dávaly jen lineární elasticitu. Rajagopal a Srinivasa použili jinou produkci entropie ˜ ξ = ξ(D, Dk(p(t) , Bkp(t) ) = 2ν0 |D|2 + 2ν1 Dkp(t) · Bkp(t) Dkp(t) . Pozorujeme, že Dkp(t) · Bkp(t) Dkp(t) = |Dkp(t) Fkp(t) |2 ≥ 0, a tedy druhý zákon termodynamiky je splněn. Definujme Lagrangeovu funkci ˜ . . ) + λ1 (ξ(. ˜ . . ) − Td · D − µCd L(D, Dkp(t) ) = ξ(. dis kp(t) · Dkp(t) ) + λ2 tr D + λ3 tr Dkp(t) . Nutné podmínky jsou ∂L = 0, ∂D

∂L = 0. ∂Dkp(t)

Dosadíme 1 + λ1 ∂ ξ˜ λ2 = Tddis + I λ1 ∂D λ1 ˜ λ3 1 + λ1 ∂ ξ = µCdkp(t) + I. λ1 ∂Dkp(t) λ1

(9) (10)

Eliminujme nyní Lagrangeovy multiplikátory, nejprve derivujme ∂ ξ˜ = 4ν0 D, ∂D

∂ ξ˜ = 4ν1 Bkp(t) Dkp(t) . ∂Dkp(t)

Nyní proveďme (9) · D + (10) · Dkp(t) , použijeme vazeb nestlačitelnosti a dostaneme Tddis · D + µCdkp(t) · Dkp(t) 1 ξ˜ 1 + λ1 = . = ˜ = ∂ ξ˜ ∂ξ ˜ λ1 2 2ξ ·D+ · Dkp(t) ∂D

∂Dkp(t)

Nyní provedeme tr (9) a získáme 0 = tr Tddis − 3

λ2 λ2 ⇒ = 0. λ1 λ1

Provedeme-li tr (10), dostaneme λ3 2 = ν1 Dkp(t) · Bkp(t) . λ1 3 Získali jsme Tddis = 2ν0 D,   1 2ν1 Bkp(t) Dkp(t) − (Dkp(t) · Bkp(t) )I = µCdkp(t) . 3

(11) (12)

Z (11) plyne, že T = −pI + 2ν0 D + µBdkp(t) . Před dalším výpočtem počítejme, připomeňme, že ˙k B = LBkp(t) + Bkp(t) LT − 2Fkp(t) Dkp(t) FT kp(t) p(t) a vypočítejme, čemu je rovna Oldroydova derivace deviatorické části Bkp(t)  O d Bkp(t) = Bkp(t)

O  O O 1 1 ˙ k )I − 1 (tr Bk ) I = − (tr Bkp(t) )I =Bkp(t) − (tr B p(t) p(t) 3 3 3 1 = −2Fkp(t) Dkp(t) FT tr(LBkp(t) + LT Bkp(t) )I+ kp(t) − 3 2 2 + tr(FT kp(t) Fkp(t) Dkp(t) )I + (tr Bkp(t) )D = 3 3 2 2 2 = −2Fkp(t) Dkp(t) FT kp(t) − (D · Bkp(t) )I + (Ckp(t) · Dkp(t) )I + (tr Bkp(t) )D. 3 3 3

37

Nyní provedeme polární rozklad deformačního gradientu Fkp(t) = Rkp(t) Ukp(t) = Vkp(t) Rkp(t) . Vzhledem k tomu, že uvažujeme izotropní materiál, můžeme předpokládat, že materiál je natočený správným směrem, a tedy Rkp(t) = I. Pak je Fkp(t) = Ukp(t) = Vkp(t) a Fkp(t) je tedy symetrický tenzor. Vynásobíme-li rovnici (12) zleva F−1 kp(t) a zprava Fkp(t) , dostaneme (víme, tedy že v tom případě Ckp(t) = Bkp(t) )   1 2ν1 Fkp(t) Dkp(t) Fkp(t) − (Bkp(t) · Dkp(t) )I = µBdkp(t) . 3 Ze vztahu pro Oldroydovu derivaci deviatorické části Bkp(t) pak dostáváme ! 1 1 1 Od d 2ν1 − Bkp(t) − (D · Bkp(t) )I + (tr Bkp(t) )D = µBdkp(t) . 2 3 3 Všimneme si nyní, že stopa této rovnice je nulová, ve třech rozměrech se tedy jedná o soustavu pěti skalárních rovnic pro šest neznámých. Soustavu uzavřeme rovnicí pro nestlačitelnost   1 det Bkp(t) + (tr Bkp(t) )I = 1. 3 Dostáváme soustavu rovnic pro nelineární viskoelastický model rychlostního typu div v = 0, ρv˙ = div T, T = −pI + 2ν0 D + µBdkp(t) , !

2ν1

−

1 1 1 Od Bkp(t) − (D · Bdkp(t) )I + (tr Bkp(t) )D = µBdkp(t) , 2 3 3   1 det Bkp(t) + (tr Bkp(t) )I = 1. 3

(13) (14)

Linearizace viskoelastického modelu Provedeme linearizaci elastické odezvy tohoto modelu a ukážeme, že pak tento model přechází na Oldroyd-B model. Uvažujme tedy levý Cauchyho-Greenův tenzor ve tvaru Bkp(t) = I + A, kde kAk = ε, 0 < ε  1. Linearizací myslíme to, že po dosazení do nelineárního viskoelastického modelu budeme zanedbávat členy ε2 a vyšší. Použijme poslední rovnici pro nestlačitelnost (14) a proveďme Taylorův rozvoj determinantu


1 (tr A)2 − tr A2 + O(ε3 ) ⇒ 1 = det Bkp(t) = det(I + A) = 1 + tr(A) + 2

1 tr(A) = tr A2 − (tr A)2 + O(ε3 ) . 2 | {z } O(ε2 )

Z toho vyplývá, že tr Bkp(t) = 3 +



1 tr A2 − (tr A)2 + O(ε3 ), 2

dosadíme a získáme vztah pro Bdkp(t) ,

1 1 Bdkp(t) = Bkp(t) − (tr Bkp(t) )I = A − tr A2 − (tr A)2 I + O(ε3 )I. 3 6 Než dosadíme do rovnice pro nelineární viskoelastický model (13) vypočítejme, čemu je rovno O

O

d d T ˙d Bdkp(t) = B kp(t) − LBkp(t) − Bkp(t) L =A −


1 ˙ ˙ I + O(ε2 ). A · A − (tr A)(tr A) | {z } 3 O(ε2 )

Dosadíme nyní do (13) a získáme  i O 1 h ˙ · A I + O(ε2 ). µA + ν1 A= 2ν1 D − ν1 2D − A 3 38

(15)

Tento model je velmi podobný Oldroyd-B modelu, liší se s ním ale v posledním členu. Ukážeme, že tento člen je velikostí O(ε2 ) a že ho tedy můžeme zanedbat. Násobme (15) skalárně s tenzorem A, dostáváme h  i ˙ · A tr A . ˙ · A − ν1 (LA · A + ALT · A) = 2ν1 (D · A) − 1 ν1 2D − A µ |A|2 +ν1 A |{z} |{z} {z } | 3 O(ε2 )

O(ε2 )

O(ε2 )

Z toho vyplývá, že10 

 ˙ · A = O(ε2 ). 2D − A

Tím jsme ukázali, že nelineární viskoelastický model přechází linearizací na Oldroyd-B model div v = 0, ρv˙ = div T, T = −pI + 2ν0 D + µA, O

µA + ν1 A = 2ν1 D a odvozený termodynamicky kompatibilní nelineární viskoelastický model je v jistém smyslu konzistentní s používanými viskoelastickými modely. Následující část textu nebyla odpřednesena a nebude zkoušena.

Nelineární dvousložkový viskoelastický model Uvažujme nyní materiál skládající se ze dvou složek, každé složce odpovídá jedna přirozená konfigurace kp1 (t) a kp2 (t) (viz obr. (36)). Volme entropii η = η˜(e, ρ, Bkp1 (t) , Bkp2 (t) ), můžeme invertovat a používat vnitřní energii F κ(t)

κ(r)

Fkp1 (t)

G1

κ(p1 (t))

Fkp2 (t)

G2 κ(p2 (t))

Obrázek 36: Dvě přirozené konfigurace e = e˜(η, ρ, Bkp1 (t) , Bkp2 (t) ) a speciálně volme e = e0 (η, ρ) +

 µ   µ1  2 tr Bkp1 (t) − 3 + tr Bkp2 (t) − 3 , 2ρ 2ρ

uvažujeme izotermální nestlačitelný model, pak platí, že tr D = tr Dkp1 (t) = tr Dkp2 (t) = 0. Obdobně jako pro případ s jednou přirozenou konfigurací, dostaneme vztah pro rychlost produkce entropie ξ = (T − µ1 Bkp1 (t) − µ2 Bkp2 (t) ) · D + µ1 Cdkp = Tddis · D + µ1 Cdkp 10 Všimněte

1 (t)

· Dkp1 (t) + µ2 Cdkp

2 (t)

1 (t)

· Dkp2 (t) .

si, že z daného lze vyčíst, že ˙ = 2D + O(ε), A

což u malých deformací očekáváme.

39

· Dkp1 (t) + µ2 Cdkp

2 (t)

· Dkp2 (t)

Volme následující vztah pro rychlost produkce entropie ˜ ξ = ξ(D, Dkp1 (t) , Dkp2 (t) , Bkp1 (t) , Bkp2 (t) ) = 2ν0 |D|2 + 2ν1 Dkp1 (t) · Bkp1 (t) Dkp1 (t) + 2ν2 Dkp2 (t) · Bkp2 (t) Dkp2 (t) . Opět pozorujeme, že rychlost produkce entropie je nezáporná a splňuje tedy druhý zákon termodynamiky. Definujme Lagrangeovu funkci ˜ . . ) + λ1 (ξ(. ˜ . . ) − Td · D − µ1 Cd L(D, Dkp1 (t) , Dkp2 (t) ) = ξ(. dis kp

1 (t)

· Dkp1 (t) − µ2 Cdkp

2 (t)

· Dkp2 (t) )+

λ2 tr D + λ3 tr Dkp1 (t) + λ4 tr Dkp2 (t) . Nutné podmínky pro maximum jsou ∂L = 0, ∂D

∂L = 0. ∂Dkp2 (t)

∂L = 0, ∂Dkp1 (t)

Dosadíme 1 + λ1 ∂ ξ˜ λ2 = Tddis + I, λ1 ∂D λ1 λ3 1 + λ1 ∂ ξ˜ = µ1 Cdkp (t) + I, 1 λ1 ∂Dkp1 (t) λ1

(16) (17)

λ4 1 + λ1 ∂ ξ˜ = µ1 Cdkp (t) + I. 2 λ1 ∂Dkp2 (t) λ1

(18)

Eliminujme nyní Lagrangeovy multiplikátory, nejprve derivujme ∂ ξ˜ = 4ν1 Bkp1 (t) Dkp1 (t) , ∂Dkp1 (t)

∂ ξ˜ = 4ν0 D, ∂D

∂ ξ˜ = 4ν2 Bkp2 (t) Dkp2 (t) . ∂Dkp2 (t)

Nyní proveďme (16) · D + (17) · Dkp1 (t) + (18) · Dkp2 (t) , použijeme vazeb nestlačitelnosti a dostaneme Tddis · D + µ1 Cdkp (t) · Dkp1 (t) + µ2 Cdkp (t) · Dkp2 (t) 1 ξ˜ 1 + λ1 1 2 = . = = ˜ ˜ ˜ ∂ξ ∂ξ ∂ξ ˜ λ1 2 2 ξ ·D+ · Dkp (t) + · Dkp (t) ∂D

∂Dkp

∂Dkp

1

1 (t)

2

2 (t)

Nyní provedeme tr (16) a získáme 0 = tr Tddis − 3

λ2 λ2 ⇒ = 0. λ1 λ1

Provedeme-li tr (17) a tr (18), dostaneme 2 λ3 = ν1 Dkp1 (t) · Bkp1 (t) , λ1 3 2 λ4 = ν2 Dkp2 (t) · Bkp2 (t) . λ1 3 Získali jsme

 2ν1 Bkp1 (t) Dkp1 (t)  2ν2 Bkp2 (t) Dkp2 (t)

1 − (Dkp1 (t) 3 1 − (Dkp2 (t) 3

Tddis = 2ν0 D,  · Bkp1 (t) )I = µ1 Cdkp (t) , 1  · Bkp2 (t) )I = µ2 Cdkp (t) .

(19) (20) (21)

2

Z (19) plyne, že T = −pI + 2ν0 D + µ1 Bdkp

1 (t)

+ µ2 Bdkp

2 (t)

.

Jako v případě jedné přirozené konfigurace platí O

Bdkp

1 (t)

= −2Fkp1 (t) Dkp1 (t) FT kp

1 (t)

2 (t)

= −2Fkp2 (t) Dkp2 (t) FT kp

2 (t)

O

Bdkp

2 − (D · Bkp1 (t) )I + 3 2 − (D · Bkp2 (t) )I + 3 40

2 (Ckp1 (t) · Dkp1 (t) )I + 3 2 (Ckp2 (t) · Dkp2 (t) )I + 3

2 (tr Bkp1 (t) )D, 3 2 (tr Bkp2 (t) )D 3

a díky izotropii materiálu víme, že Fkp1 (t) = Vkp1 (t) a Fkp2 (t) = Vkp2 (t) jsou symetrické tenzory. Vynásobíme-li −1 rovnici (20) zleva F−1 kp (t) a zprava Fkp1 (t) , a rovnici (21) zleva Fkp (t) a zprava Fkp2 (t) , dostaneme 1

2



 1 2ν1 Fkp1 (t) Dkp1 (t) Fkp1 (t) − (Bkp1 (t) · Dkp1 (t) )I = µ1 Bdkp (t) , 1 3   1 2ν2 Fkp2 (t) Dkp2 (t) Fkp2 (t) − (Bkp2 (t) · Dkp2 (t) )I = µ2 Bdkp (t) . 2 3 Ze vztahů pro Oldroydovu derivaci Bdkp

1 (t)

a Bdkp

2 (t)

pak dostáváme !

2ν1

1 1 O 1 − Bdkp (t) − (D · Bdkp (t) )I + (tr Bkp1 (t) )D 1 1 2 3 3

!

2ν2

1 1 1 O − Bdkp (t) − (D · Bdkp (t) )I + (tr Bkp2 (t) )D 2 2 2 3 3

= µ1 Bdkp

,

= µ2 Bdkp

.

1 (t)

2 (t)

Všimneme si nyní, že stopa obou rovnic je nulová, ve třech rozměrech se tedy jedná o soustavu dvakrát pět skalárních rovnic pro dvakrát šest neznámých. Soustavu uzavřeme rovnicemi pro nestlačitelnost   1 det Bkp1 (t) + (tr Bkp1 (t) )I = 1, 3   1 det Bkp2 (t) + (tr Bkp2 (t) )I = 1. 3 Dostáváme soustavu rovnic pro dvousložkový nelineární viskoelastický model rychlostního typu div v = 0, ρv˙ = div T, T = −pI + 2ν0 D + µ1 Bdkp (t) + µ2 Bdkp (t) , 1 2 !

2ν1 2ν2

1 1 1 Od Bkp (t) − (D · Bdkp (t) )I + (tr Bkp1 (t) )D = µ1 Bdkp (t) , 1 1 1 2 3 3 ! 1 O 1 1 − Bdkp (t) − (D · Bdkp (t) )I + (tr Bkp2 (t) )D = µ2 Bdkp (t) , 2 2 2 2 3 3   1 det Bkp1 (t) + (tr Bkp1 (t) )I = 1, 3   1 det Bkp2 (t) + (tr Bkp2 (t) )I = 1. 3 −

41