Równania Naviera-Stokesa

Z testwiki
Przejdź do nawigacji Przejdź do wyszukiwania

Równania Naviera-Stokesa (nazwane na cześć Claude’a-Louis Naviera i George’a Gabriela Stokesa) – zestaw równań opisujących zasadę zachowania pędu dla poruszającego się płynu. Według nich zmiany pędu elementu płynu zależą jedynie od sił masowych, zewnętrznego ciśnienia i wewnętrznych sił lepkości w płynie.

Dla płynu idealnego o zerowej lepkości równania mówią, że przyspieszenie jest proporcjonalne do gradientu ciśnienia.

Równania są wyrażone w postaci różniczkowej, dla danego problemu fizycznego muszą być znalezione na drodze rachunku różniczkowego i całkowego. W praktyce jedynie najprostsze przypadki mogą być rozwiązane analitycznie, na przykład nieturbulentnego (laminarnego), stacjonarnego przepływu (niezmieniającego się w czasie), w których liczba Reynoldsa ma małą wartośćSzablon:R.

W bardziej złożonych przypadkach, jak prognozowanie pogody na Ziemi, analiza El Niño lub obliczenia siły nośnej skrzydeł samolotów, rozwiązania równań Naviera-Stokesa mogą być znalezione jedynie metodami numerycznymi przy pomocy komputerów. Jest to oddzielna dziedzina nauki zwana obliczeniową mechaniką płynów.

W 2000 roku Instytut Matematyczny Claya ogłosił równania Naviera-Stokesa jednym z siedmiu problemów milenijnych matematyki i zaoferował 1 000 000 dolarów nagrody za podanie rozwiązania lub kontrprzykładu[1]. Oficjalny opis problemu przedstawił Charles Fefferman.[2]

Ogólna forma równań

Ogólna forma równań Naviera-Stokesa[3]. Dwa ostatnie człony po prawej stronie równania wynikają z uwzględnienia niestałości lepkości w obrębie przepływu – uproszczone formy równań Naviera-Stokesa znajdują się w sekcji #Dodatkowe założenia:

ρ(vt+(v)v)=ρfp+μv+(λ+μ)(v)+(v)(λ)+(v+(v)T)(μ),

gdzie:

v=[xxu+yyu+zzuxxv+yyv+zzvxxw+yyw+zzw].
v=[xuxvxwyuyvywzuzvzw] [1/s]gradient prędkości (rozwinięty w kartezjańskim układzie współrzędnych, gdzie np. xu=ux)
(v)=[x(v)y(v)z(v)]=[xxu+xyv+xzwyxu+yyv+yzwzxu+zyv+zzw)] [1/ms] wyrażenie takiej postaci to gradient z dywergencji prędkości w kartezjańskim układzie współrzędnych (jest to wektor, gdyż dywergencja prędkości to skalar, a gradient ze skalaru to wektor), gdzie oznaczenie np. xyu=2uxy
v=xu+yv+zw [1/s]dywergencja prędkości
p=[xp,yp,zp]T [kg/m2s2]gradient ciśnienia

Powyższe równania Naviera-Stokesa zostały podane w formie operatorowej, która dla wszystkich układów współrzędnych przestrzennych ma taką samą postać (ale już nie w sytuacji gdy zmieniamy współrzędne w czasie np. wprowadzenie wirującego układu współrzędnych spowoduje pojawienie się dodatkowych elementów w równaniach związanych z efektem Coriolisa). Ponadto, przy wyznaczaniu jawnej postaci różniczkowej równań dla konkretnego układu współrzędnych przestrzennych musimy znaleźć postać wszystkich elementów równania (w tym operatorów) dla docelowego układu współrzędnych – przykładowo: element μv zawierający wektorowy operator Laplace’a dla cylindrycznego układu współrzędnych będzie miał postać:

μv=μ[uu/r22φv/r2vv/r2+2φu/r2w]=μ[ru/r+rru+φφu/r2+zzuu/r22φv/r2rv/r+rrv+φφv/r2+zzvv/r2+2φu/r2rw/r+rrw+φφw/r2+zzw]

gdzie symbol pochodnej cząstkowej z dolnym indeksem/ami obejmuje tylko bezpośrednio następujący po nim symbol np. φφv/r2=(φφv)/r2=1r22vφφ

Postać w kartezjańskim układzie współrzędnych

Równanie wektorowe NS dla poszczególnych współrzędnych przyjmuje postać

x:ρ(tu+uxu+vyu+wzudvdt)=ρfxρfxpp+μ(xxu+yyu+zzu)μv+(λ+μ)(xxu+xyv+xzw)(λ+μ)((v))+(xu+yv+zw)(xλ)(v)(λ)+((xμ)(xu)+(yμ)(yu)+(zμ)(zu))(v)T(μ)+((xμ)(xu)+(yμ)(xv)+(zμ)(xw))(v)(μ)y:ρ(tv+uxv+vyv+wzv)=ρfyyp +μ(xxv+yyv+zzv)+(λ+μ)(yxu+yyv+yzw)+(xu+yv+zw)(yλ)+((xμ)(xv)+(yμ)(yv)+(zμ)(zv))+((xμ)(yu)+(yμ)(yv)+(zμ)(yw))z:ρ(tw+uxw+vyw+wzw)=ρfzzp+μ(xxw+yyw+zzw)+(λ+μ)(zxu+zyv+zzw)+(xu+yv+zw)(zλ)+((xμ)(xw)+(yμ)(yw)+(zμ)(zw))+((xμ)(zu)+(yμ)(zv)+(zμ)(zw))

gdzie użyto oznaczeń pochodnych cząstkowych xf=fx oraz xyf=2fxy.

Wyprowadzenie

Wyprowadzenie rozpoczyna się od równania pędu Cauchy’ego zapisanego w kartezjańskim układzie współrzędnych (upraszczający obliczenia i stanowiącym punkt wyjścia dla innych układów wsp.), dla którego σdiv(σT)=(Tσ)T gdzie w ostatniej równości po prawej działania możemy traktować jak mnożenie/transpozycje wektorów i macierzy (co również będziemy wykorzystywali w dalszej części wyprowadzenia)[4]:

ρdvdt=(Tσ)T+ρf.

Aby otrzymać równania Naviera-Stokesa, należy poczynić założenia w celu zamodelowania tensora naprężeń σ. Zwyczajowo tensor naprężenia dzieli się na dwie części:

σ=p𝑰+τ,

gdzie 𝑰 to macierz jednostkowa. Gdy płyn jest nieruchomy w polu sił, to tensor τ zeruje się natomiast ciśnienie p niekoniecznie.

W oznaczeniach poniżej przyjmuje się, że v=[u,v,w]=[v1,v2,v3] oraz x=[x,y,z]=[x1,x2,x3]. Ponadto będziemy stosować konwencję sumacyjną Einsteina.

Założenie 1

Równania opisują płyn newtonowski, tzn. tensor τ jest liniową funkcją gradientu prędkości, tzn. ma postać:

τij=Lijklvkxl,

gdzie: 𝑳=[Lijkl] jest tensorem czwartego rzędu (ma 81 składowych). Ponieważ stosujemy konwencję sumacyjną Einsteina, więc przed prawą stroną powyższego równania stoją dwa znaki sumy: po indeksie k oraz l.

Założenie 2

Płyn jest izotropowy i jednorodny. Implikuje to że tensor 𝑳 powinien być niezależny od kierunku. Innym znanym tensorem niezależnym od kierunku jest δij (delta Kroneckera). Delta Kroneckera jest rzędu 2 więc my musimy znaleźć jej odpowiednik, ale rzędu 4. W celu zbadania niezależności tensora od kierunku wprowadźmy 4-liniową formę S (funkcję dającą skalar w wyniku) będącą iloczynem wewnętrznym tensora 𝑳 i czterech wektorów a=[ai],b=[bi],c=[ci],d=[di] (ponieważ stosujemy konwencję sumacyjną Einsteina, więc przed prawą stroną równania poniżej stoją cztery znaki sumy: po indeksie i,j,k oraz l.):

S(a,b,c,d)=Lijklaibjckdl

Tensor 𝑳 jest niezależny od kierunku wektorów na których działa tylko wtedy gdy również funkcja S jest niezależna. Aby funkcja S była niezależna od kierunku, nie powinna być zależna od bezwzględnego położenia wektorów a,b,c,d, ale powinna zmieniać wartość gdy wektory zmieniają swoje długości lub położenia względem siebie. Zatem S powinna zmieniać wartość, gdy zmieniają się cosinusy kątów między poszczególnymi wektorami (i gdy zmienia się długość wektorów). Funkcją, która to umożliwia jest iloczyn skalarny, wówczas S możemy skonstruować następująco (dla zwięzłości pomijamy listę argumentów S):

S=α(ab)(cd)+β(ac)(bd)+γ(ad)(bc).

Inne liniowe funkcje ponad powyższą (np. z wyrazami typu (a×b)(c×d)), niczego więcej nie wniosą. Zatem rozpisując iloczyny skalarne, mamy:

S=αaibicjdj+βaicibjdj+γaidibjcj,

przekształcając dalej:

S=aibjckdl(αδijδkl+βδikδjl+γδilδjk)=Lijklaibjckdl,

zatem tensor 𝑳, aby był izotropowy, musi mieć postać:

Lijkl=αδijδkl+βδikδjl+γδilδjk.
Założenie 3

Symetria tensora naprężeń, tj. σij=σji. Badając równowagę elementarnego sześcianu, zakładając, że nie występują naprężenia momentowe (dla których uogólnioną teorię sformułowali bracia Cosserat, 1909)[5], można dowieść, że tensor naprężenia jest symetryczny. Jednak w ogólności tak nie jest i np. dla płynów polarnych w polu elektromagnetycznym takie założenie jest błędne. Klasyczne równania Naviera-Stokesa nie uwzględniają tego typu płynów, ale za to uwzględniają szeroką klasę płynów powszechnie używanych. Zatem:

σij=σjiLijkl=Ljikl.

Podstawiając wyprowadzoną postać tensora 𝑳 pod lewą i prawą stronę równania, otrzymamy:

αδijδkl+βδikδjl+γδilδjk=αδjiδkl+βδjkδil+γδjlδik

i dalej:

βδikδjl+γδilδjk=βδjkδil+γδjlδik,

co prowadzi do:

(βγ)δikδjl=(βγ)δjkδil.

Ponieważ w szczególności ta równość musi zachodzić dla i=k,j=l,ij, np. i=1,j=2, to mamy:

(βγ)δ11δ22=(βγ)δ21δ12(βγ)11=(βγ)00.

Zatem

β=γ.

Podstawiając to do tensora 𝑳, otrzymujemy:

Lijkl=αδijδkl+β(δikδjl+δilδjk).
Ostatnie kroki wyprowadzenia

Zamieniając nazwy parametrów na powszechnie używane (a pierwotnie wywodzące się z rozważań na temat tensora sztywności dla ciał stałych, gdzie nazywane są stałymi Lamégo), tj.: μ=β oraz λ=α w tensorze 𝑳 i podstawiając go do tensora τ, mamy:

τij=(λδijδkl+μ(δikδjl+δilδjk))vkxl.

Po uwzględnieniu działania delty Kroneckera otrzymamy:

τij=λδijvkxk+μ(vixj+vjxi).

Ostatecznie tensor naprężenia σ=p𝑰+τ można zapisać, używając operatorów następująco (literka T w górnym indeksie oznacza transpozycję macierzy):

σ=p𝑰+λ(v)𝑰+μ(v+(v)T).

W ten sposób doszliśmy do końca definicji modelu i jest to właściwie koniec wyprowadzenia. Podstawiając taką formę tensora naprężeń σ do równania pędu Cauchy’ego, otrzymamy równania Naviera-Stokesa.

Końcowe podstawienie

Ponieważ jednak owo podstawienie nie jest trywialne (ze względów rachunkowych) dokonamy go poniżej by otrzymać jawną postać operatorową równań Naviera-Stokesa. Wyliczmy dywergencję transponowaną z σ, korzystając z faktu, że w układzie kartezjańskim możemy ją wyliczyć jako mnożenie wektora z macierzą:

Tσ=Tp𝑰+Tλ(v)𝑰+Tμ(v+(v)T)=(T𝑰)p+(T𝑰)(λv)+T(μv)+T(μv)T.

W przejściu do ostatniej równości skorzystaliśmy z przemienności mnożenia skalaru przez macierz. Rozpisując dokładnie wszystkie operatory (na postać pochodnych cząstkowych), można sprawdzić, że prawdziwy jest poniższy ciąg równości (uwaga! uwzględniono tu, że lepkości są funkcjami zależnymi od położenia, a nie stałymi):

(Tσ)T=p+T(λv)+(T(μv))T+(T(μv)T)T=p+λ(v)+(λ)(v)T(λv)+μ(Tv)Tμv+(v)T(μ)(T(μv))T+μ(T(v)T)Tμ(v)+(v)(μ)(T(μv)T)T=p+μv+(λ+μ)(v)+(v)(λ)+(v+(v)T)(μ).

Zatem podstawiając pod równanie Cauchy’ego, mamy:

ρdvdt=ρfp+μv+(λ+μ)(v)+(v)(λ)+(v+(v)T)(μ).

Rozpisując pochodną substancjalną z lewej strony, otrzymamy formę równania taką jak w sekcji #Ogólna forma równań.

Dodatkowe założenia

Możemy dodać kolejne założenie:

Założenie 4

Doświadczalnie wyznaczona zależność (tr to ślad macierzy)

tr(σ)=3p,

która jest dokładna dla gazów jednoatomowych, ale w praktyce znajduje zastosowanie do znacznie szerszej klasy płynów. Szczególnie gdy płyn jest prawie nieściśliwy, bo wówczas v0 i człon z λ w równaniach Naviera-Stokesa przestaje mieć znaczenie.

Jeżeli teraz wyliczymy ślad z tensora σ, to otrzymamy:

tr(σ)=tr(p𝑰+λ(v)𝑰+μ(v+(v)T)),

co po podstawieniu i rozwinięciu daje:

3p=3p+3λ(v)+2μ(v).

Zatem:

3λ(v)=2μ(v).

Więc:

λ=23μ.

Po podstawieniu tej zależności do równań Naviera-Stokesa przybiorą one postać:

ρ(vt+(v)v)=ρfp+μv+13μ(v)23(v)(μ)+(v+(v)T)(μ).
Założenie 5

Jeżeli przyjmiemy, że lepkość jest stała (lepkość silnie zależy od temperatury, więc gdy rozpatrujemy „płyn zimny” (bez uwzględniania równania energii), to wówczas taka sytuacja może mieć miejsce), to wówczas μ=0 i równania upraszczają się do postaci:

ρ(vt+(v)v)=ρf+(μ3vp)+μv.
Założenie 6

Przyjmując, że płyn jest nieściśliwy v=0, otrzymamy:

ρ(vt+(v)v)=ρfp+μv.

Właśnie ten przypadek ze stałą ρ=1 został przedstawiony jako problem milenijny[2].

Założenie 7

Jeżeli całkowicie pominiemy lepkość, tj. μ=0, to otrzymamy równania Eulera:

ρ(vt+(v)v)=ρfp.

Przypisy

Szablon:Przypisy

Linki zewnętrzne

Szablon:Równania różniczkowe Szablon:Problemy milenijne

Szablon:Kontrola autorytatywna