Równania Naviera-Stokesa

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacja, szukaj

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 pochodnej ciśnienia.

Rozwiązania równań dla danego problemu fizycznego muszą być znalezione na drodze rachunku różniczkowego i całkowego. W praktyce, jedynie najprostsze przypadki mogą być rozwiązane dokładnie na tej drodze. To znaczy przypadki nieturbulentnego, spokojnego przepływu (nie zmieniającego się w czasie), w których liczba Reynoldsa ma małą wartość.

W bardziej złożonych przypadkach, takich jak systemy badania pogody na Ziemi, takie jak El Niño lub przy obliczeniach 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]

Ogólna forma równań[edytuj | edytuj kod]

Ogólna forma równań Naviera-Stokes’a[2] (dwa ostatnie człony po prawej stronie równania wynikają z uwzględnienia niestałości lepkości w obrębie przepływu - bardziej uproszczone formy równań Naviera-Stoeksa znajdują sie w sekcji #Dodatkowe założenia):

\rho\left(\frac{\partial\vec v}{\partial t}+\vec v\cdot\nabla\vec v\right) =\rho\vec f -\nabla p + \mu\triangle\vec v + (\lambda+\mu)\nabla(\nabla\cdot v) + (\nabla\cdot v)(\nabla\lambda) + (\nabla\mu)\cdot\big(\nabla\vec v + (\nabla\vec v)^T\big)

gdzie:

  • \rho\ [kg/m^3]gęstość płynu,
  • \vec v\ [m/s] – wektor prędkości,
  • \nabla\ [1/m]Operator nabla, poniżej przykłady użycia w kontekście gradientu i dywergencji:
\nabla\vec v=\begin{bmatrix}
\partial_xu & \partial_xv & \partial_xw \\ 
\partial_yu & \partial_yv & \partial_yw \\ 
\partial_zu & \partial_zv & \partial_zw
\end{bmatrix}\ [1/s]gradient prędkości (rozwinięty w kartezjańskim układzie współrzednych gdzie np.\partial_xu=\frac{\partial u}{\partial x}),
\begin{align}
\nabla(\nabla\cdot\vec v) &= \begin{bmatrix}\partial_x(\nabla\cdot\vec v) \\ \partial_y(\nabla\cdot\vec v) \\ \partial_z(\nabla\cdot\vec v)\end{bmatrix} = \begin{bmatrix}\partial_{xx}u+\partial_{xy}v+\partial_{xz}w \\ \partial_{yx}u+\partial_{yy}v+\partial_{yz}w \\ \partial_{zx}u+\partial_{zy}v+\partial_{zz}w)\end{bmatrix}\ [1/ms^2]
\end{align} wyrażenie takiej postaci to gradient z dywergencji predkoś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. \partial_{xy}u=\frac{\partial^2u}{\partial x\partial y}
\nabla\cdot\vec v=\partial_xu+\partial_yv+\partial_zw\ [1/s]dywergencja prędkości,
\nabla p = [\partial_xp,\partial_yp,\partial_zp]^T\ [kg/m^2s^2]gradient cisnienia,
  • \vec f\ [m/s^2] – wektor przyspieszenia płynu (sił masowych np. przyspieszenie ziemskie)
  • p\ [Pa=kg/m\cdot s^2] – ciśnienie,
  • \mu\ [kg/m\cdot s]lepkość kinematyczna (powyższe równania uwzględniają że nie jest ona stałą tylko jest rózna dla różnych punktów przepływu)
  • \lambda\ [kg/m\cdot s]lepkość objętościowa (powyższe równania uwzględniają że nie jest ona stałą tylko jest rózna dla różnych punktów przepływu)
  • \triangle\vec v = \nabla(\nabla \cdot \vec v) - \nabla \times (\nabla \times \vec v)\ [1/m\cdot s] - wektorowy operator Laplace'a (nie mylić z skalarnym operatorem Laplace'a, który akurat w układzie kartezjańskim, wykonany na wektorze, jest mu równoważny (ale już nie w np. cylindrycznym)). W kartezjańskim układzie dostajemy taki wektor:
\triangle\vec v = \begin{bmatrix}\partial_{xx}u+\partial_{yy}u+\partial_{zz}u \\ \partial_{xx}v+\partial_{yy}v+\partial_{zz}v\\ \partial_{xx}w+\partial_{yy}w+\partial_{zz}w\end{bmatrix}

Postać w kartezjańskim układzie współrzędnych[edytuj | edytuj kod]

\begin{align}
 \rho (\overbrace{\partial_tu + u\partial_xu + v\partial_yu + w\partial_zu}^{\frac{d\vec v}{dt}})
    &= \overbrace{\rho f_x}^{\rho\vec f}-\overbrace{\partial_xp}^{\nabla p} + 
\overbrace{\mu(\partial_{xx}u+\partial_{yy}u+\partial_{zz}u)}^{\mu\triangle\vec v} + 
\overbrace{(\lambda+\mu)(\partial_{xx}u+\partial_{xy}v+\partial_{xz}w)}^{(\lambda+\mu)(\nabla(\nabla\cdot\vec v))} + 
\overbrace{(\partial_xu+\partial_yv+\partial_zw)(\partial_x\lambda)}^{(\nabla\cdot\vec v)(\nabla\lambda)} + 
\overbrace{\big((\partial_x\mu)(\partial_xu) + (\partial_y\mu)(\partial_yu) + (\partial_z\mu)(\partial_zu)\big)}^{(\nabla\mu)\cdot(\nabla\vec v)} + 
\overbrace{\big((\partial_x\mu)(\partial_xu) +(\partial_y\mu)(\partial_xv) + (\partial_z\mu)(\partial_xw)\big)}^{(\nabla\mu)\cdot(\nabla\vec v)^T}
\\
  \rho \left(\partial_tv + u\partial_xv + v\partial_yv + w\partial_zv\right)
    &=\ \rho f_y-\partial_yp\ +\mu(\partial_{xx}v+\partial_{yy}v+\partial_{zz}v) + 
(\lambda+\mu)(\partial_{yx}u+\partial_{yy}v+\partial_{yz}w) +
(\partial_xu+\partial_yv+\partial_zw)(\partial_y\lambda) + 
\big((\partial_x\mu)(\partial_xv)+(\partial_y\mu)(\partial_yv)+(\partial_z\mu)(\partial_zv)\big) + 
\big((\partial_x\mu)(\partial_yu)+(\partial_y\mu)(\partial_yv)+(\partial_z\mu)(\partial_yw)\big) \\

  \rho \left(\partial_tw + u\partial_xw + v\partial_yw + w\partial_zw\right)
    &=\rho f_z-\partial_zp +\mu(\partial_{xx}w+\partial_{yy}w+\partial_{zz}w) + 
(\lambda+\mu)(\partial_{zx}u+\partial_{zy}v+\partial_{zz}w) +
(\partial_xu+\partial_yv+\partial_zw)(\partial_z\lambda) + 
\big((\partial_x\mu)(\partial_xw)+(\partial_y\mu)(\partial_yw)+(\partial_z\mu)(\partial_zw)\big) + 
\big((\partial_x\mu)(\partial_zu)+(\partial_y\mu)(\partial_zv)+(\partial_z\mu)(\partial_zw)\big)
\end{align}

gdzie uzyto oznaczeń pochodnych cząstkowych \partial_xf=\frac{\partial f}{\partial x} oraz \partial_{xy}f=\frac{\partial^2 f}{\partial x\partial y}

Wyprowadzenie[3][edytuj | edytuj kod]

Wychodzimy od równania pędu Cauchy'ego:

\rho \frac{d\vec v}{dt}=\nabla\cdot\boldsymbol\sigma + \rho\vec f

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

\boldsymbol\sigma=-p\boldsymbol I + \boldsymbol\tau

gdzie \boldsymbol I to macierz jednostkowa. Gdy płyn stoi nieruchomo w polu sił to tensor \boldsymbol\tau zeruje się natomiast ciśnienie p niekoniecznie.

W oznaczeniach poniżej przyjmiemy że \vec v=[u,v,w]=[v_1,v_2,v_3] oraz \vec x=[x,y,z]=[x_1,x_2,x_3]. Ponadto będziemy stosować konwencję sumacyjną Einsteina.

  • Założenie 1: równania opisują płyn newtonowski tzn. tensor \boldsymbol\tau jest liniową funkcją gradientu prędkości tzn ma postać:
\tau_{ij}=L_{ijkl}\frac{\partial v_k}{\partial x_l}

gdzie: \boldsymbol L=[L_{ijkl}] jest tensorem czwartego rzędu (ma 81 składowych). Ponieważ sotsujemy 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 \boldsymbol L powinien być niezależny od kierunku. Innym znanym tenosrem niezależnym od kierunku jest \delta_{ij} (delta Kroneckera). Delta Kroneckera jest rzędu 2 więc my musimy znaleźć jej odpowiednik ale rzędu 4. Wprowadźmy więc skalar S będący iloczynem wewnętrznym tensora \boldsymbol L i czterech wektorów \vec a=[a_i], \vec b=[b_i], \vec c=[c_i], \vec d=[d_i] (ponieważ sotsujemy konwencję sumacyjną Einsteina więc przed prawą stroną równania ponizej stoją cztery znaki sumy: po indeksie i,j,k oraz l.):
S=L_{ijkl}a_ib_jc_kd_l

Aby skalar S był niezależny od kierunku, nie powinien być zależny od bezwzględnego położenia wektorów \vec a,\vec b,\vec c,\vec d, ale powinien zmieniać wartość gdy wektory zmieniają swoje długości lub położenia wzgledem siebie. Zatem skalar S powinien zmieniać wartość gdy zmieniają sie cosinusy kątów między poszczególnymi wektorami (i gdy zmienia sie długość wektorów). Funkcją która to umożliwia jest iloczyn skalarny, wówczas:

S=\alpha(\vec a\cdot\vec b)(\vec c\cdot\vec d) + \beta(\vec a\cdot\vec c)(\vec b\cdot\vec d) + \gamma(\vec a\cdot\vec d)(\vec b\cdot\vec c)

Inne dodatkowe funkcje wektorów ponad powyższą, niczego juz nie wniosą. Zatem rozpisując iloczyny skalarne mamy:

S=\alpha a_ib_ic_jd_j + \beta a_ic_ib_jd_j + \gamma a_id_ib_jc_j

przekształcając dalej:

S=a_ib_jc_kd_l(\alpha\delta_{ij}\delta_{kl} + \beta\delta_{ik}\delta_{jl} + \gamma\delta_{il}\delta_{jk})=L_{ijkl}a_ib_jc_kd_l

zatem tensor \boldsymbol L, aby był izotropowy, musi mieć postać:

L_{ijkl}=\alpha\delta_{ij}\delta_{kl} + \beta\delta_{ik}\delta_{jl} + \gamma\delta_{il}\delta_{jk}


  • Założenie 3. Symetria tensora naprężeń tj. \sigma_{ij}=\sigma_{ji}. Badając równowagę elemenatrnego sześcianu, zakładając że nie wsytępują naprężenia momentowe (dla których uogulnioną teorię sformułowali bracia Cosserat, 1909)[4], 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. 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:
\sigma_{ij}=\sigma_{ji} \qquad\Rightarrow\qquad L_{ijkl}=L_{jikl}

Podstawiając wyprowadzoną postać tensora \boldsymbol L pod lewą i prawą stronę równania otrzymamy:

\alpha\delta_{ij}\delta_{kl} + \beta\delta_{ik}\delta_{jl} + \gamma\delta_{il}\delta_{jk}= 
        \alpha\delta_{ji}\delta_{kl} + \beta\delta_{jk}\delta_{il} + \gamma\delta_{jl}\delta_{ik}

i dalej:

\beta\delta_{ik}\delta_{jl} + \gamma\delta_{il}\delta_{jk}= 
        \beta\delta_{jk}\delta_{il} + \gamma\delta_{jl}\delta_{ik}

co prowadzi do:

(\beta-\gamma)\delta_{ik}\delta_{jl} = 
        (\beta-\gamma)\delta_{jk}\delta_{il}

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

(\beta-\gamma)\delta_{11}\delta_{22} = 
        (\beta-\gamma)\delta_{21}\delta_{12}\qquad \Rightarrow \qquad
(\beta-\gamma)\cdot 1 \cdot 1 = (\beta-\gamma)\cdot 0\cdot 0

Zatem

\beta=\gamma

Podstawiając to do tenosra \boldsymbol L otrzymujemy:

L_{ijkl}=\alpha\delta_{ij}\delta_{kl} + \beta(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})

Zamieniajac nazwy parametrów na powszechnie używane tj.: \mu=\beta oraz \lambda=\alpha do tensora \boldsymbol L w tensorze \boldsymbol\taumamy:

\tau_{ij}=\left(\lambda\delta_{ij}\delta_{kl} + \mu(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})\right)\frac{\partial v_k}{\partial x_l}

Po uwzględnieniu działania delty Kroneckera otrzymamy:

\tau_{ij}=\lambda\delta_{ij}\frac{\partial v_k}{\partial x_k} + \mu\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)

Ostatecznie tensor naprężenia \boldsymbol\sigma można zapisać używając operatorów następuąco (literka T w górnym indeksie oznacza transpozycję macierzy):

\boldsymbol\sigma= -p\boldsymbol I + \lambda(\nabla\cdot v)\boldsymbol I +\mu(\nabla\vec v + (\nabla\vec v)^T)

W ten sposób doszlismy do końca definicji modelu i jest to właściwie koniec wyprowadzenia. Podstawiając taką formę tensora naprężeń \boldsymbol\sigma do równania pedu Cauchy'ego otrzymamy równania Naviera-Stokesa.

Jednak aby wszystko było w pełni jasne dokonamy ten krok podstawienia. Wyliczmy dywergencję z \boldsymbol\sigma:

\nabla\cdot\boldsymbol\sigma =  -\nabla\cdot p\boldsymbol I + \nabla\cdot\lambda(\nabla\cdot v)\boldsymbol I +\nabla\cdot\mu(\nabla\vec v + (\nabla\vec v)^T) = -(\nabla\cdot \boldsymbol I) p +(\nabla\cdot \boldsymbol I)(\lambda\nabla\cdot v) + \nabla\cdot(\mu\nabla\vec v) + \nabla\cdot(\mu\nabla\vec 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 zaleznymi od połozenia a nie stałymi):

\begin{align}
\nabla\cdot\boldsymbol\sigma &= -\nabla p +\nabla(\lambda\nabla\cdot v)+ \nabla\cdot(\mu\nabla\vec v) + \nabla\cdot(\mu\nabla\vec v)^T \\
&= -\nabla p +\overbrace{\lambda\nabla(\nabla\cdot\vec v) + (\nabla\lambda)(\nabla\cdot v)}^{\nabla(\lambda\nabla\cdot v)} + \overbrace{\underbrace{\mu\nabla\cdot(\nabla\vec v)}_{\mu\triangle\vec v} + (\nabla\mu)\cdot(\nabla\vec v)}^{\nabla\cdot(\mu\nabla\vec v)} + \overbrace{\underbrace{\mu\nabla\cdot(\nabla\vec v)^T}_{\mu\nabla(\nabla\cdot v)} + (\nabla\mu)\cdot(\nabla\vec v)^T}^{\nabla\cdot(\mu\nabla\vec v)^T} \\
&= -\nabla p + \mu\triangle\vec v + (\lambda+\mu)\nabla(\nabla\cdot v) + (\nabla\cdot v)(\nabla\lambda) + (\nabla\mu)\cdot\big(\nabla\vec v + (\nabla\vec v)^T\big)
\end{align}

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

\rho \frac{d\vec v}{dt}=\rho\vec f -\nabla p + \mu\triangle\vec v + (\lambda+\mu)\nabla(\nabla\cdot v) + (\nabla\cdot v)(\nabla\lambda) + (\nabla\mu)\cdot\big(\nabla\vec v + (\nabla\vec v)^T\big)

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

Dodatkowe założenia[edytuj | edytuj kod]

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

  • Założenie 4: Doswiadczalnie wyznaczona zależność (tr to ślad macierzy)
tr(\boldsymbol\sigma)=-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 \nabla\cdot\vec v \rightarrow 0 i człon z \lambda w równaniach Naviera-Stokesa przestaje mieć znaczenie. Jeżeli teraz wyliczymy ślad z tensora \boldsymbol\sigma to otrzymay:

tr(\boldsymbol\sigma)=tr\Big(-p\boldsymbol I + \lambda(\nabla\cdot v)\boldsymbol I +\mu(\nabla\vec v + (\nabla\vec v)^T)\Big)

co po podstawieniu i rozwinięciu daje:

-3p=-3p + 3\lambda(\nabla\cdot v) +2\mu(\nabla\cdot v)

Zatem:

 3\lambda(\nabla\cdot v)=-2\mu(\nabla\cdot v)

Więc:

\lambda=-\frac{2}{3}\mu

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

\rho\left(\frac{\partial\vec v}{\partial t}+\vec v\cdot\nabla\vec v\right) =\rho\vec f -\nabla p + \mu\triangle\vec v + \frac{1}{3}\mu\nabla(\nabla\cdot v) - \frac{2}{3}(\nabla\cdot v)(\nabla\mu) + (\nabla\mu)\cdot\big(\nabla\vec v + (\nabla\vec v)^T\big)
  • Założenie 5: jeżeli przyjmiemy że lepkość jest stała (lepkość silnie zalezy 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 \nabla\mu=0 i równania upraszczają się do postaci:
\rho\left(\frac{\partial\vec v}{\partial t}+\vec v\cdot\nabla\vec v\right) =\rho\vec f + \nabla\left(\frac{\mu}{3}\nabla\cdot\vec v-p\right) + \mu\triangle\vec v
  • Założenie 6: przyjmując że płyn jest nieściśliwy \nabla\cdot v=0 otrzymamy:
\rho\left(\frac{\partial\vec v}{\partial t}+\vec v\cdot\nabla\vec v\right) =\rho\vec f - \nabla p + \mu\triangle\vec v
  • Założenie 7: jeżeli całkowicie pominiemy lepkość tj. \mu=0 to otrzymamy równania Eulera:
\rho\left(\frac{\partial\vec v}{\partial t}+\vec v\cdot\nabla\vec v\right) =\rho\vec f - \nabla p

Przypisy

  1. Millennium Prize Problems (ang.). [dostęp 2011-01-07].
  2. Joe Pedlosky: Fluid Dynamics of the Atmosphere and Ocean. Woods Hole Oceanographic Institution, 2014, s. rozdział 3 s.27.
  3. Joe Pedlosky: Fluid Dynamics of the Atmosphere and Ocean. Woods Hole Oceanographic Institution, 2014, s. rozdział 3 s.16-21.
  4. Andrzej Gawęcki: Mechanika materiałów i konstrukcji prętowych. Alma Mater, 2003, s. część1: s.3, 10.