Równanie pędu Cauchy’ego

Wprowadzenie

Równanie pędu Cauchy’ego – wektorowe równanie różniczkowe cząstkowe zaproponowane przez Cauchy’ego, które opisuje nierelatywistyczny transport pędu w każdym ośrodku ciągłym. Dane jest ono następująco:

$\rho \frac{d\overrightarrow{v}}{dt}=\mathrm{\nabla }\cdot \mathit{\sigma }+\rho \overrightarrow{f}$

gdzie:

$\rho \mathrm{[}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}\mathrm{]}$ – gęstość,

$d\overrightarrow{v}/dt={\partial }_t\overrightarrow{v}+\overrightarrow{v}\cdot \mathrm{\nabla }\overrightarrow{v}\mathrm{[}\mathrm{m}/{\mathrm{s}}^{\mathrm{2}}\mathrm{]}$ – pochodna substancjalna prędkości,

$\mathrm{\nabla }\mathrm{[}\mathrm{1}/\mathrm{m}\mathrm{]}$ – operator nabla,

$\mathit{\sigma }\mathrm{[}\mathrm{P}\mathrm{a}\mathrm{=}\mathrm{N}/{\mathrm{m}}^{\mathrm{2}}\mathrm{]}$ – tensor naprężenia,

$\overrightarrow{f}\mathrm{[}\mathrm{m}/{\mathrm{s}}^{\mathrm{2}}\mathrm{]}$ – przyspieszenie związane z siłami masowymi.

Jest to równanie wektorowe (tzn. jego rozwiązaniem jest pole wektorowe) które po rozwinięciu w układzie kartezjańskim ma postać trzech równań – po jednym dla każdej składowej wynikowego pola wektorowego:

$\begin{array}{cccc}x& :& \rho (\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z})& =\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{yx}}{\partial y}+\frac{\partial {\sigma }_{zx}}{\partial z}+\rho f_x\\ y& :& \rho (\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z})& =\frac{\partial {\sigma }_{xy}}{\partial x}+\frac{\partial {\sigma }_{yy}}{\partial y}+\frac{\partial {\sigma }_{zy}}{\partial z}+\rho f_y\\ z& :& \rho (\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z})& =\frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{yz}}{\partial y}+\frac{\partial {\sigma }_{zz}}{\partial z}+\rho f_z\end{array}$

Jak widać układ nie jest zamknięty, gdyż mamy tylko 3 równania a 13 niewiadomych tj. $\rho$ (skalar – 1 niewiadoma), $\overrightarrow{v}$ (vektor – 3 niewiadome), $\mathit{\sigma }$ (macierz – 9 niewiadomych). Poza tym $t,x,y,z,\overrightarrow{f}$ są wiadome w ramach warunków początkowych/brzegowych.

Detale dotyczące $\mathrm{\nabla }\cdot \mathit{\sigma }$

Dla $\mathrm{\nabla }\cdot \mathit{\sigma }$ operacja „$\cdot$” pomiędzy operatorem nabla $\mathrm{\nabla }$ a tensorem naprężenia $\mathit{\sigma }$ NIE JEST mnożeniem macierzy – nie można użyć lewostronnego mnożenia pomiędzy wektorem kontrawariantnym (kolumnowym; „pionowym”; co wynika z definicji operatora) a macierzą. Wyrażenie $\mathrm{\nabla }\cdot \mathit{\sigma }$ ma tylko symboliczne znaczenie. Standardowo dywergencja jest zdefiniowana jako operator na polu wektorowym (jako Iloczyn skalarny pomiędzy operatorem nabla i wektorem (więcej: Pochodna kowariantna)), ale może zostać rozszerzona na pola tensorowe 2 rzędu (macierze), jednak w tym przypadku mamy do czynienia z inną operacją: $\mathrm{\nabla }\cdot \mathit{\sigma }=\frac{\partial {\sigma }_{ki}}{\partial x_k}{𝐞}_i={\sigma }_{ki,k}{𝐞}_i$. W ogólności mamy $\mathrm{\nabla }\cdot \mathit{\sigma }\ne \mathrm{d}\mathrm{i}\mathrm{v}(\mathit{\sigma })=\mathrm{\nabla }\cdot {\mathit{\sigma }}^T.$ W kartezjańskim układzie współrzędnych możemy tę operację zapisać jako mnożenie macierzy (poniżej używamy „pusty” symbol mnożenia) w następujący sposób: $\mathrm{\nabla }\cdot \mathit{\sigma }=({\mathrm{\nabla }}^T\mathit{\sigma }{)}^T={\left({\left[\begin{array}{c}\frac{\partial }{\partial x}\\ \frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}\end{array}\right]}^T\left[\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{yx}& {\sigma }_{yy}& {\sigma }_{yz}\\ {\sigma }_{zx}& {\sigma }_{zy}& {\sigma }_{zz}\end{array}\right]\right)}^T={\left(\left[\frac{\partial }{\partial x}\frac{\partial }{\partial y}\frac{\partial }{\partial z}\right]\left[\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{yx}& {\sigma }_{yy}& {\sigma }_{yz}\\ {\sigma }_{zx}& {\sigma }_{zy}& {\sigma }_{zz}\end{array}\right]\right)}^T=\left[\begin{array}{c}\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{yx}}{\partial y}+\frac{\partial {\sigma }_{zx}}{\partial z}\\ \frac{\partial {\sigma }_{xy}}{\partial x}+\frac{\partial {\sigma }_{yy}}{\partial y}+\frac{\partial {\sigma }_{zy}}{\partial z}\\ \frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{yz}}{\partial y}+\frac{\partial {\sigma }_{zz}}{\partial z}\end{array}\right]$ gdzie rezultat (po prawej) jest kolumnowym wektorem, ${}^T$ oznacza transpozycję. Transpozycja operatora nabla w powyższych równaniach pozwala na użycie lewostronnego mnożenia macierzy (od drugiej równości), kolejna transpozycja (wykonana po mnożeniu macierzowym) zmienia rezultat z wektora wierszowego na kolumnowy (kontrawariantnym) co powzowli na dodawanie pozostałych składowych równania pędu (które są wektorami kolumnowymi) jak np. $\rho \overrightarrow{f}.$ Dla przypadku gdy sigma jest tensorem symetrycznym ${\sigma }_{ij}={\sigma }_{ji}$ (np. dla płynu Newtonowskiego w równaniach Naviera-Stokesa bazujących na równaniu pędu Cauchego) zachodzi $\mathrm{\nabla }\cdot \mathit{\sigma }=\mathrm{d}\mathrm{i}\mathrm{v}(\mathit{\sigma }).$ Jednak w ogólności nie powinniśmy mylić operacji $\mathrm{\nabla }\cdot \mathit{\sigma }$ z operacją $\mathrm{d}\mathrm{i}\mathrm{v}(\mathit{\sigma })$ która dla niesymetrycznego tensora sigma (używanym np. dla modelowania płynów nienewtonowskich jak polimery) daje inne rezultaty

Wyprowadzenie różniczkowe

Wychodzimy od uogólnionej zasady zachowania pędu, którą można zapisać następująco: „zmiana pędu układu jest proporcjonalna do siły wypadkowej działającej na ten układ” wyraża się ona wzorem:

$\overrightarrow{p}(t+\mathrm{\Delta }t)-\overrightarrow{p}(t)=\mathrm{\Delta }t\overrightarrow{\overline{F}}$

gdzie:

$\overrightarrow{p}(t)$ – pęd w chwili $t,$

$\overrightarrow{\overline{F}}$ – uśredniona w czasie $\mathrm{\Delta }t$ siła.

Po podzieleniu przez $\mathrm{\Delta }t$ i przejściu do granicy $\mathrm{\Delta }t\to 0$ otrzymujemy:

$\frac{d\overrightarrow{p}}{dt}=\overrightarrow{F}$

Skupimy się teraz kolejno na wyznaczeniu lewej, a następnie prawej strony powyższego równania dla stałej masy w różniczkowej sześciennej objętości kontrolnej (czyli elemencie ciała) której pęd i działające nań siły chcemy zbadać. Na koniec zestawimy lewą i prawą stronę powyższego równania, otrzymując równanie Pędu Cauchy’ego.

Zacznijmy od prawej strony równania

Składowa X sił działających na ścianki sześciennego elementu płynu (zielone dla górnej-dolnej ścianki; czerwone dla lewej-prawej; czarne dla przedniej-tylnej).

Na górnym wykresie widzimy przybliżenie funkcji $f(x)$ (niebieska linia) za pomocą różnicy skończonej (żółta linia). Na dolnym wykresie jest „nieskończenie wiele razy powiększone otoczenie punktu $x_1$” (fioletowy kwadracik z górnego wykresu). Na dolnym wykresie żółta linia nie została naniesiona. Na rysunku użyto dwóch równoważnych oznaczeń pochodnej $f^{\prime }(x_1)=\frac{df(x_1)}{d{x}_1}$ ponadto oznaczono $\mathrm{\Delta }f=f(x_1+\mathrm{\Delta }x)-f(x_1).$

Siły dzielimy na masowe i powierzchniowe

$\overrightarrow{F}={\overrightarrow{F}}_p+{\overrightarrow{F}}_m$

Na ścianki objętości kontrolnej działają siły powierzchniowe. Składowa X tych sił (w formie iloczynu naprężenia i pola powierzchni np. $-{\sigma }_{xx}dydz\mathrm{[}\mathrm{P}\mathrm{a}\mathrm{\cdot }\mathrm{m}\mathrm{\cdot }\mathrm{m}\mathrm{=}\frac{N}{m^2}\mathrm{\cdot }\mathrm{m}\mathrm{\cdot }\mathrm{m}\mathrm{=}\mathrm{N}\mathrm{]}$), dla każdej ścianki, została umieszczona na rysunku z elementem sześciennym.

Wyjaśnienie wartości sił (przybliżeń i znaków minus) na ściankach sześcianu

Wyjaśnienia wymaga dlaczego przy naprężeniach na ściankach leżących na osiach współrzędnych mamy znak minus (np. na lewej ściance mamy wartość $-{\sigma }_{xx}$). Dla uproszczenia skupmy się na lewej ściance z naprężeniem $-{\sigma }_{xx}.$ Znak minus spowodowany jest tym, że wektor normalny do tej ścianki $\overrightarrow{n}=[-1,0,0]=-{\overrightarrow{e}}_x$ jest ujemnym wektorem jednostkowym. Wówczas przy obliczaniu wektora naprężenia z definicji $\overrightarrow{s}=\overrightarrow{n}\cdot \mathit{\sigma }=[-{\sigma }_{xx},-{\sigma }_{xy},-{\sigma }_{xy}]$ zatem składowa x tego wektora to $s_x=-{\sigma }_{xx}$ (podobne rozumowanie używamy dla naprężeń na ściance dolnej i tylnej tj.: $-{\sigma }_{yx},-{\sigma }_{zx}$). Drugim elementem wymagającym wyjaśnienia jest przybliżenie wartości naprężeń na ściankach przeciwległych do ścianek leżących na osiach. Skupmy się na ściance prawej na której naprężenie jest przybliżeniem naprężenia ${\sigma }_{xx}$ ze ścianki lewej w punktach o wsp. $x+dx$ i wynosi ono ${\sigma }_{xx}+\frac{\partial {\sigma }_{xx}}{\partial x}dx.$ To przybliżenie wzięło się z zastosowania wzoru Taylora na przybliżenie funkcji tj.: ${\sigma }_{xx}(x+dx)={\sigma }_{xx}(x)+dx\frac{\partial {\sigma }_{xx}(x)}{\partial x}+\frac{d{x}^2}{2}\frac{{\partial }^2{\sigma }_{xx}(x)}{\partial x^2}+...$ Ponieważ wartość $d{x}^2$ jest nieskończenie mniejsza od wartości $dx$ więc wszystkie wyrazy z $dx$ w potęgach wyższych niż jeden możemy porzucić, uznając za nieistotne (argument Leibnitza - związany z iloczynem pochodnych). W ten sposób otrzymaliśmy szukane przybliżenie naprężenia na przeciwległej ściance. Bardziej intuicyjne przedstawienie przybliżenia wartości ${\sigma }_{xx}$ w punkcie $x+dx$ obrazuje rysunek pod sześcianem. Podobne rozumowanie przeprowadzamy dla przybliżeń naprężeń ${\sigma }_{yx},{\sigma }_{zx}$

Sumując siły (ich składowe X) działające na każdą ze ścian, otrzymujemy:

$F_p^x=({\sigma }_{xx}+\frac{\partial {\sigma }_{xx}}{\partial x}dx)dydz-{\sigma }_{xx}dydz+({\sigma }_{yx}+\frac{\partial {\sigma }_{yx}}{\partial y}dy)dxdz-{\sigma }_{yx}dxdz+({\sigma }_{zx}+\frac{\partial {\sigma }_{zx}}{\partial z}dz)dxdy-{\sigma }_{zx}dxdy$

Po uporządkowaniu $F_p^x$ oraz po wykonaniu podobnego rozumowania dla składowych $F_p^y,F_p^z$ (nie ma ich na rysunku – będą to wektory równoległe odpowiednio do osi Y i Z) otrzymamy:

$F_p^x=\frac{\partial {\sigma }_{xx}}{\partial x}dxdydz+\frac{\partial {\sigma }_{yx}}{\partial y}dydxdz+\frac{\partial {\sigma }_{zx}}{\partial z}dzdxdy$

$F_p^y=\frac{\partial {\sigma }_{xy}}{\partial x}dxdydz+\frac{\partial {\sigma }_{yy}}{\partial y}dydxdz+\frac{\partial {\sigma }_{zy}}{\partial z}dzdxdy$

$F_p^z=\frac{\partial {\sigma }_{xz}}{\partial x}dxdydz+\frac{\partial {\sigma }_{yz}}{\partial y}dydxdz+\frac{\partial {\sigma }_{zz}}{\partial z}dzdxdy$

W zapisie operatorowym możemy to wówczas zapisać:

${\overrightarrow{F}}_p=(\mathrm{\nabla }\cdot \mathit{\sigma })dxdydz$

Na wnętrze objętości kontrolnej działają siły masowe, które zapiszemy z wykorzystaniem pola przyspieszenia $\overrightarrow{f}$ (którym może być np. przyspieszenie ziemskie $\overrightarrow{g}$):

${\overrightarrow{F}}_m=\rho \overrightarrow{f}dxdydz$

Lewa strona równania

Wyznaczmy pęd:

$\overrightarrow{p}=m\overrightarrow{v}=\rho \overrightarrow{v}dxdydz$

Ponieważ założyliśmy, że badana masa jest stała więc

$\frac{d\overrightarrow{p}}{dt}=\rho \frac{d\overrightarrow{v}}{dt}dxdydz$

Porównanie lewej i prawej strony równania

Otrzymamy:

$\rho \frac{d\overrightarrow{v}}{dt}dxdydz=(\mathrm{\nabla }\cdot \mathit{\sigma })dxdydz+\rho \overrightarrow{f}dxdydz$

Dzieląc przez $dxdydz,$ dostaniemy:

$\rho \frac{d\overrightarrow{v}}{dt}=\mathrm{\nabla }\cdot \mathit{\sigma }+\rho \overrightarrow{f}$

co kończy wywód.