Dobra skończona różnica dla równania ciągłości

Jaka byłaby dobra dyskretyzacja różnic skończonych dla następującego równania:

$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho u\right)=0$ ?

Możemy wziąć przypadek 1D:

$\frac{\partial \rho}{\partial t} + \frac{d}{dx}\left(\rho u\right)=0$

Z jakiegoś powodu wszystkie schematy, które mogę znaleźć, dotyczą sformułowania we współrzędnych Lagrangiana. Na razie wymyśliłem ten schemat (zignoruj indeks j ):

$\frac{\rho^{n+1}_{i,j}-\rho^{n}_{i,j}}{\tau} + \frac{1}{h_x}\left(\frac{\rho^{n+1}_{i+1,j}+\rho^{n+1}_{i,j}}{2}u^{n}_{_xi+1/2,j}- \frac{\rho^{n+1}_{i,j}+\rho^{n+1}_{i-1,j}}{2}u^{n}_{_xi-1/2}\right)=0$

Ale wydaje się być naprawdę niestabilny lub mieć okropny warunek stabilności. Czy tak jest

Prędkość jest faktycznie obliczana na podstawie prawa darcy . Dodatkowo mamy równanie stanu. Pełny układ składa się również z równania energii i równania stanu gazu doskonałego. Prędkości mogą stać się ujemne .$u=-\frac{k}{\mu}\nabla p$

Patrzysz na równanie zachowania masy:

$\dfrac{dm}{dt}=0$

Rozważając ewolucję masy na jednostkę objętości, sprowadza się to do równania doradztwa gęstości w postaci strumienia:

$\dfrac{\partial \rho}{\partial t} = -\nabla \cdot (\rho u)$

Dobrą rzeczą jest to, że jest to tylko równanie doradcze dowolnego pola skalarnego (w naszym przypadku jest to gęstość ) i jest (względnie) łatwe do rozwiązania, pod warunkiem odpowiednich schematów różnicowania czasu i przestrzeni oraz początkowego i warunki brzegowe.$\rho$

Projektując skończony schemat różnicowania, martwimy się o zbieżność, stabilność i dokładność. Schemat jest zbieżny, jeśli gdyΔt→0. Stabilność schematów zapewnia, że ​​ilośćApozostaje skończona, gdy. Formalna dokładność schematu wskazuje, gdzie leży błąd skracania w szeregu rozszerzalności Taylora pochodnej cząstkowej. Zajrzyj do podręcznika CFD, aby uzyskać więcej informacji na temat tych podstawowych właściwości schematu różnicowania.$\dfrac{\Delta A}{\Delta t} \rightarrow \dfrac{\partial A}{\partial t}$$\Delta t \rightarrow 0$$A$$t \rightarrow \infty$

Teraz najprostszym podejściem jest przejście bezpośrednio do pierwszego rzędu różnicowania. Ten schemat jest pozytywnie określony, konserwatywny i wydajny obliczeniowo. Dwie pierwsze właściwości są szczególnie ważne, gdy modelujemy ewolucję wielkości, która jest zawsze dodatnia (tj. Masa lub gęstość).

Dla uproszczenia spójrzmy na przypadek 1-D:

$\dfrac{\partial \rho}{\partial t} = -\dfrac{\partial(\rho u)}{\partial x}$

Teraz wygodnie jest zdefiniować strumień , aby:$\Phi = \rho u$

$\dfrac{\partial(\rho u)}{\partial x} = \dfrac{\partial \Phi}{\partial x} \approx \dfrac{\Delta \Phi}{\Delta x} \approx \dfrac{\Phi_{i+1/2}-\Phi_{i-1/2}}{\Delta x}$

Oto schemat tego, co symulujemy:

            u           u
|          -->         -->          |
|    rho    |    rho    |    rho    |
x-----o-----x-----o-----x-----o-----x
     i-1  i-1/2   i   i+1/2  i+1

Oceniamy ewolucję w komórce i . Zysk lub strata netto wynika z różnicy co przychodzi, cp I - 1 / 2 i co wychodzi, Φ i + 1 /$\rho$$i$$\Phi_{i-1/2}$$\Phi_{i+1/2}$. Tutaj zaczynamy odbiegać od odpowiedzi Pawła. W prawdziwym zachowawczym różnicowaniu w górę, ilość w centrum komórki jest przenoszona przez prędkość na krawędzi komórki, w kierunku jej ruchu. Innymi słowy, jeśli wyobrażasz sobie, że jesteś zalecaną ilością i siedzisz w centrum komórki, jesteś przenoszony do komórki przed sobą z prędkością na krawędzi komórki. Ocena strumienia na krawędzi komórki jako iloczynu gęstości i prędkości, zarówno na krawędzi komórki, jest nieprawidłowa i nie chroni zalecanej ilości.

Strumienie przychodzące i wychodzące są oceniane jako:

$\Phi_{i+1/2} = \dfrac{u_{i+1/2}+|u_{i+1/2}|}{2}\rho_{i} + \dfrac{u_{i+1/2}-|u_{i+1/2}|}{2}\rho_{i+1}$

$\Phi_{i-1/2} = \dfrac{u_{i-1/2}+|u_{i-1/2}|}{2}\rho_{i-1} + \dfrac{u_{i-1/2}-|u_{i-1/2}|}{2}\rho_{i}$

Powyższe podejście do różnicowania strumienia zapewnia ostateczność. Innymi słowy, dostosowuje kierunek różnicowania zgodnie ze znakiem prędkości.

Kryterium stabilności Couranta-Friedrichsa-Lewy'ego (CFL), gdy różnicowanie czasu odbywa się za pomocą prostego pierwszego rzędu, różnicowanie w przód Eulera jest podawane jako:

$\mu = \dfrac{u\Delta t}{\Delta x} \leq 1$

Zauważ, że w 2 wymiarach kryterium stabilności CFL jest bardziej rygorystyczne:

$\mu = \dfrac{c\Delta t}{\Delta x} \leq \dfrac{1}{\sqrt{2}}$

gdzie jest wielkością prędkości, $c$ .$\sqrt{u^{2}+v^{2}}$

Kilka rzeczy do rozważenia. Ten schemat może, ale nie musi być odpowiedni dla twojej aplikacji, w zależności od rodzaju symulowanego procesu. Ten schemat jest wysoce dyfuzyjny i jest odpowiedni dla bardzo płynnych przepływów bez ostrych gradientów. Jest również bardziej dyfuzyjny w przypadku krótszych przedziałów czasowych. W przypadku 1-D uzyskasz prawie dokładne rozwiązanie, jeśli gradienty są bardzo małe, a jeśli . W przypadku 2D nie jest to możliwe, a dyfuzja jest anizotropowa.$\mu = 1$

Jeśli twój system fizyczny rozważa fale uderzeniowe lub innego rodzaju wysokie gradienty, powinieneś rozważyć różnicowanie wyższego rzędu wyższego rzędu (np. 3 lub 5 rzędu). Warto również spojrzeć na rodzinę systemów z poprawkami Fluxed Transport (Zalesak, 1979, JCP); korekta antydyfuzyjna dla powyższego schematu autorstwa Smolarkiewicza (1984, JCP); Rodzina schematów MPDATA autorstwa Smolarkiewicza (1998, JCP).

W przypadku różnicowania czasu, różnicowanie Eulera pierwszego rzędu może być satysfakcjonujące dla twoich potrzeb. W przeciwnym razie spójrz na metody wyższego rzędu, takie jak Runge-Kutta (iteracyjny) lub Adams-Bashforth i Adams-Moulton (wielopoziomowy).

Warto byłoby zajrzeć do podręcznika dla absolwentów CFD, który zawiera podsumowanie wyżej wymienionych programów i wiele innych.

$(\frac{d}{dx})$

$\frac{\rho_i^{k+1}-\rho_i^k}{\Delta t} + \frac{\rho_i^k U_i^k-\rho_{i-1}^k U_{i-1}^k}{\Delta x} = 0$

$\Delta x$$\Delta t$