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{\vec {v}}}{dt}}=\nabla \cdot {\boldsymbol {\sigma }}+\rho {\vec {f}}$

gdzie:

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

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

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

${\boldsymbol {\sigma }}\;\mathrm {[Pa=N/m^{2}]}$ – tensor naprężenia,

${\vec {f}}\;\mathrm {[m/s^{2}]}$ – 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{aligned}x&:&\rho \left({\frac {\partial u}{\partial t}}+u{\frac {\partial u}{\partial x}}+v{\frac {\partial u}{\partial y}}+w{\frac {\partial u}{\partial z}}\right)&={\frac {\partial \sigma _{xx}}{\partial x}}+{\frac {\partial \sigma _{yx}}{\partial y}}+{\frac {\partial \sigma _{zx}}{\partial z}}+\rho f_{x}\\[8pt]y&:&\rho \left({\frac {\partial v}{\partial t}}+u{\frac {\partial v}{\partial x}}+v{\frac {\partial v}{\partial y}}+w{\frac {\partial v}{\partial z}}\right)&={\frac {\partial \sigma _{xy}}{\partial x}}+{\frac {\partial \sigma _{yy}}{\partial y}}+{\frac {\partial \sigma _{zy}}{\partial z}}+\rho f_{y}\\[8pt]z&:&\rho \left({\frac {\partial w}{\partial t}}+u{\frac {\partial w}{\partial x}}+v{\frac {\partial w}{\partial y}}+w{\frac {\partial w}{\partial z}}\right)&={\frac {\partial \sigma _{xz}}{\partial x}}+{\frac {\partial \sigma _{yz}}{\partial y}}+{\frac {\partial \sigma _{zz}}{\partial z}}+\rho f_{z}\end{aligned}}$

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

Detale dotyczące $\nabla \cdot {\boldsymbol {\sigma }}$

Dla $\nabla \cdot {\boldsymbol {\sigma }}$ operacja „$\cdot$” pomiędzy operatorem nabla $\nabla$ a tensorem naprężenia ${\boldsymbol {\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 $\nabla \cdot {\boldsymbol {\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ą: $\nabla \cdot {\boldsymbol {\sigma }}={\cfrac {\partial \sigma _{ki}}{\partial x_{k}}}~\mathbf {e} _{i}=\sigma _{ki,k}~\mathbf {e} _{i}$. W ogólności mamy $\nabla \cdot {\boldsymbol {\sigma }}\neq \mathrm {div} ({\boldsymbol {\sigma }})=\nabla \cdot {\boldsymbol {\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: $\nabla \cdot {\boldsymbol {\sigma }}=(\nabla ^{T}{\boldsymbol {\sigma }})^{T}=\left({\begin{bmatrix}{\frac {\partial }{\partial x}}\\{\frac {\partial }{\partial y}}\\{\frac {\partial }{\partial z}}\end{bmatrix}}^{T}{\begin{bmatrix}\sigma _{xx}&\sigma _{xy}&\sigma _{xz}\\\sigma _{yx}&\sigma _{yy}&\sigma _{yz}\\\sigma _{zx}&\sigma _{zy}&\sigma _{zz}\end{bmatrix}}\right)^{T}=\left(\left[{\tfrac {\partial }{\partial x}}\ \ {\tfrac {\partial }{\partial y}}\ \ {\tfrac {\partial }{\partial z}}\right]{\begin{bmatrix}\sigma _{xx}&\sigma _{xy}&\sigma _{xz}\\\sigma _{yx}&\sigma _{yy}&\sigma _{yz}\\\sigma _{zx}&\sigma _{zy}&\sigma _{zz}\end{bmatrix}}\right)^{T}={\begin{bmatrix}{\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{bmatrix}}$ gdzie rezultat (po prawej) jest kolumnowym wektorem, $\mathbb {\ } ^{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 {\vec {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 $\nabla \cdot {\boldsymbol {\sigma }}=\mathrm {div} ({\boldsymbol {\sigma }}).$ Jednak w ogólności nie powinniśmy mylić operacji $\nabla \cdot {\boldsymbol {\sigma }}$ z operacją $\mathrm {div} ({\boldsymbol {\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:

${\vec {p}}(t+\Delta t)-{\vec {p}}(t)=\Delta t{\vec {\bar {F}}}$

gdzie:

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

${\vec {\bar {F}}}$ – uśredniona w czasie $\Delta t$ siła.

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

${\frac {d{\vec {p}}}{dt}}={\vec {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'(x_{1})={\frac {df(x_{1})}{dx_{1}}}$ ponadto oznaczono $\Delta f=f(x_{1}+\Delta x)-f(x_{1}).$

Siły dzielimy na masowe i powierzchniowe

${\vec {F}}={\vec {F}}_{p}+{\vec {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 {[Pa\cdot m\cdot m={\frac {N}{m^{2}}}\cdot m\cdot m=N]}$), 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 ${\vec {n}}=[-1,0,0]=-{\vec {e}}_{x}$ jest ujemnym wektorem jednostkowym. Wówczas przy obliczaniu wektora naprężenia z definicji ${\vec {s}}={\vec {n}}\cdot {\boldsymbol {\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 {dx^{2}}{2}}{\frac {\partial ^{2}\sigma _{xx}(x)}{\partial x^{2}}}+...$ Ponieważ wartość $dx^{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ć:

${\vec {F}}_{p}=(\nabla \cdot {\boldsymbol {\sigma }})dxdydz$

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

${\vec {F}}_{m}=\rho {\vec {f}}dxdydz$

Lewa strona równania

Wyznaczmy pęd:

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

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

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

Porównanie lewej i prawej strony równania

Otrzymamy:

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

Dzieląc przez $dxdydz,$ dostaniemy:

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

co kończy wywód.