Zachowanie wielkości fizycznej przy zastosowaniu warunków brzegowych Neumanna zastosowanych do równania dyfuzyjno-doradczego

Nie rozumiem różnych zachowań równania rada-dyfuzja, kiedy stosuję różne warunki brzegowe. Moją motywacją jest symulacja rzeczywistej wielkości fizycznej (gęstości cząstek) w warunkach dyfuzji i doradztwa. Gęstość cząstek należy zachować we wnętrzu, chyba że wypłynie ona z krawędzi. Zgodnie z tą logiką, jeśli wymuszę warunki brzegowe Neumanna, końce systemu, takie jak (po lewej i prawej stronie), to system powinien być „zamknięty”, tj. jeśli strumień na granicy wynosi zero, wówczas żadne cząstki nie mogą uciec.$\frac{\partial \phi}{\partial x}=0$

Dla wszystkich poniższych symulacji zastosowałem dyskretyzację Cranka-Nicolsona do równania rada-dyfuzja i wszystkie symulacje mają warunków brzegowych. Jednak dla pierwszego i ostatniego wiersza macierzy (wiersze warunków brzegowych) zezwalam na zmianę niezależnie od wartości wewnętrznej. Pozwala to na pełne domniemanie punktów końcowych.$\frac{\partial \phi}{\partial x}=0$$\beta$

Poniżej omawiam 4 różne konfiguracje, tylko jedna z nich jest tym, czego się spodziewałem. Na koniec omawiam moje wdrożenie.

Limit tylko dyfuzji

Tutaj warunki porady są wyłączane przez ustawienie prędkości na zero.

Tylko dyfuzja, z $\boldsymbol{\beta}$ = 0,5 (Crank-Niscolson) we wszystkich punktach

Ilość nie jest zachowana, co widać po zmniejszeniu obszaru impulsu.

Tylko dyfuzja, przy czym = 0,5 (Crank-Niscolson) w punktach wewnętrznych, a = 1 (pełne niejawne) na granicach$\boldsymbol{\beta}$$\boldsymbol{\beta}$

Stosując w pełni niejawne równanie na granicach, osiągam to, czego oczekuję: żadne cząstki nie uciekają . Widać to po zachowaniu obszaru jako rozproszenia cząstek. Dlaczego wybór w punktach granicznych powinien wpływać na fizykę sytuacji? Czy to błąd, czy oczekiwany?$\beta$

Dyfuzja i porady

Po uwzględnieniu terminu porady wydaje się, że wartość na granicach nie wpływa na rozwiązanie. Jednak we wszystkich przypadkach, gdy granice wydają się być „otwarte”, tzn. Cząstki mogą uciec z granic. Dlaczego tak jest?$\beta$

Porady i dyfuzja z = 0,5 (Crank-Niscolson) we wszystkich punktach$\boldsymbol{\beta}$

Porady i dyfuzja z = 0,5 (Crank-Niscolson) w punktach wewnętrznych, a = 1 (pełny dorozumiany) na granicach$\boldsymbol{\beta}$$\boldsymbol{\beta}$

Implementacja równania rada-dyfuzja

Zaczynając od równania rada-dyfuzja,

$\frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial x^2} + \boldsymbol{v}\frac{\partial \phi}{\partial x}$

Pisanie przy użyciu Crank-Nicolson daje:

$\frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = D \left[ \frac{1 - \beta}{(\Delta x)^2} \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + \frac{\beta}{(\Delta x)^2} \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) \right] + \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right]$

Zauważ, że = 0,5 dla Crank-Nicolson, = 1 dla pełnego niejawnego i = 0 dla pełnego jawnego.β $\beta$$\beta$$\beta$

Aby uprościć notację, dokonajmy podstawienia,

$s = D\frac{\Delta t}{(\Delta x)^2} \\ r = \boldsymbol{v}\frac{\Delta t}{2 \Delta x}$

i przenieś znaną wartość pochodnej czasu na prawą stronę,$\phi_{j}^{n}$

$\phi_{j}^{n+1} = \phi_{j}^{n} + s \left( 1-\beta \right) \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + s \beta \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) + r \left( 1 - \beta \right) \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + r \beta \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right)$

Faktoryzacja warunków daje,$\phi$

$\underbrace{\beta(r - s)\phi_{j-1}^{n+1} + (1 + 2s\beta)\phi_{j}^{n+1} -\beta(s + r)\phi_{j+1}^{n+1}}_{\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}}} = \underbrace{ (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n}}_{\boldsymbol{M\cdot}\boldsymbol{\phi^n}}$

który możemy zapisać w postaci macierzy jako gdzie,$\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}} = \boldsymbol{M}\cdot\boldsymbol{\phi^{n}}$

$\boldsymbol{A} = \left( \begin{matrix} 1+2s\beta & -\beta(s + r) & & 0 \\ \beta(r-s) & 1+2s\beta & -\beta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \beta(r-s) & 1+2s\beta & -\beta (s + r) \\ 0 & & \beta(r-s) & 1+2s\beta \\ \end{matrix} \right)$

$\boldsymbol{M} = \left( \begin{matrix} 1-2s(1-\beta) & (1 - \beta)(s + r) & & 0 \\ (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\ 0 & & (1 - \beta)(s - r) & 1-2s(1-\beta) \\ \end{matrix} \right)$

Stosowanie warunków brzegowych Neumanna

NB znowu pracuje nad wyprowadzeniem. Myślę, że zauważyłem błąd. Przy zapisywaniu skończonej różnicy warunku brzegowego założyłem w pełni niejawny schemat ( = 1). Jeśli przyjmiesz tutaj schemat Crank-Niscolson, złożoność staje się zbyt duża i nie mogłem rozwiązać wynikowych równań, aby wyeliminować węzły znajdujące się poza domeną. Wydaje się jednak możliwe, istnieją dwa równania z dwiema niewiadomymi, ale nie mogłem sobie z tym poradzić. To prawdopodobnie tłumaczy różnicę między pierwszą a drugą fabułą powyżej. Myślę, że możemy stwierdzić, że prawidłowe są tylko wykresy z = 0,5 w punktach granicznych.$\beta$$\beta$

Zakładając, że strumień po lewej stronie jest znany (przy założeniu całkowicie niejawnej postaci),

$\frac{\partial\phi_1^{n+1}}{\partial x} = \sigma_L$

Pisanie tego jako centralnej różnicy daje,

$\frac{\partial\phi_1^{n+1}}{\partial x} \approx \frac{\phi_2^{n+1} - \phi_0^{n+1}}{2\Delta x} = \sigma_L$

dlatego $\phi_0^{n+1} = \phi_{2}^{n+1} - 2 \Delta x\sigma_L$

Zauważ, że wprowadza to węzeł który jest poza domeną problemu. Ten węzeł można wyeliminować za pomocą drugiego równania. Możemy napisać węzeł jako,$\phi_0^{n+1}$$j=1$

$\beta(r - s)\phi_0^{n+1} + (1+2s\beta)\phi_1^{n+1} - \beta(s+r)\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n}$

Podstawienie wartości znalezionej w warunku brzegowym daje następujący wynik dla wiersza = 1,$\phi_0^{n+1}$$j$

$(1+2s\beta)\phi_1^{n+1} - 2s\beta\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n} + 2\beta(r-s)\Delta x\sigma_L$

Wykonanie tej samej procedury dla ostatniego wiersza (przy = ) daje,$j$$J$

$-2s\beta\phi_{J-1}^{n+1} + (1+2s\beta)\phi_J^{n+1} = (1-\beta)(s - r)\phi_{J-1}^{n} + (1 - 2s(1-\beta))\phi_{J}^{n} + 2\beta(s+r)\Delta x\sigma_R$

Ostatecznie ukrywanie wierszy granic (ustawienie = 1) daje,$\beta$

$(1+2s)\phi_1^{n+1} - 2s\phi_2^{n+1} = \phi_{j-1}^{n} + 1\phi_{j}^{n} + 2(r-s)\Delta x\sigma_L$

$-2s\phi_{J-1}^{n+1} + (1+2s)\phi_J^{n+1} = \phi_{J}^{n} + 2(s+r)\Delta x\sigma_R$

Dlatego w warunkach brzegowych Neumanna możemy zapisać równanie macierzowe, ,$\boldsymbol{A}\cdot\phi^{n+1} = \boldsymbol{M}\cdot\phi^{n} + \boldsymbol{b_N}$

gdzie,

$\boldsymbol{A} = \left( \begin{matrix} 1+2s & -2s & & 0 \\ \beta(r-s) & 1+2s\beta & -\beta (s + r) & \\ & \ddots & \ddots & \ddots \\ & \beta(r-s) & 1+2s\beta & -\beta (s + r) \\ 0 & & -2s & 1+2s \\ \end{matrix} \right)$

$\boldsymbol{M} = \left( \begin{matrix} 1 & 0 & & 0 \\ (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\ & \ddots & \ddots & \ddots \\ & (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\ 0 & & 0 & 1 \\ \end{matrix} \right)$

$\boldsymbol{b_N} = \left( \begin{matrix} 2 (r - s) \Delta x \sigma_L & 0 & \ldots & 0 & 2 (s + r) \Delta x \sigma_R \end{matrix} \right)^{T}$

Moje obecne zrozumienie

• Myślę, że różnicę między pierwszym a drugim polem można wyjaśnić, zwracając uwagę na błąd opisany powyżej.

• W odniesieniu do zachowania ilości fizycznej. Uważam, że przyczyną jest to, że, jak wskazano tutaj , równanie doradcze w postaci, którą napisałem, nie pozwala na propagację w przeciwnym kierunku, więc fala po prostu przechodzi, nawet w warunkach brzegowych o zerowym strumieniu . Moja początkowa intuicja dotycząca ochrony była stosowana tylko wtedy, gdy termin porady wynosi zero (jest to rozwiązanie na działce nr 2, gdzie obszar jest zachowany).

• Nawet przy warunkach brzegowych zerowego strumienia Neumanna masa nadal może opuścić układ, ponieważ prawidłowe warunki brzegowe w tym przypadku to warunki brzegowe Robina , w których całkowity strumień jest określony . Ponadto warunek Neunmanna określa, że ​​masa nie może opuścić dziedziny przez dyfuzję , nie mówi nic o radach. Zasadniczo to, co słyszymy, to zamknięte warunki brzegowe dla dyfuzji i otwarte warunki brzegowe dla porady. Aby uzyskać więcej informacji, zobacz odpowiedź tutaj: Implementacja warunku granicznego zerowego gradientu w równaniu doradczym-dyfuzyjnym$\frac{\partial \phi}{\partial x} = 0$$j = D\frac{\partial \phi}{\partial x} + \boldsymbol{v}\phi = 0$.

Zgodziłbyś się?

• Samouczek dotyczący sposobu wdrażania warunków brzegowych Robina.

Myślę, że jednym z twoich problemów jest to, że (jak zauważyłeś w swoich komentarzach) warunki Neumanna nie są warunkami, których szukasz , w tym sensie, że nie oznaczają zachowania twojej ilości. Aby znaleźć odpowiedni warunek, przepisz swój PDE jako

$$\frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x}\left( D\frac{\partial \phi}{\partial x} + v \phi \right) + S(x,t) .$$

Teraz termin, który pojawia się w nawiasach, jest całkowitym strumieniem i jest to ilość, którą należy wyzerować na granicach, aby zachować . (Dodałem ze względu na ogólność i dla twoich komentarzy.) Warunki brzegowe, które musisz nałożyć, to (zakładając, że twoja domena kosmiczna to )$D\frac{\partial \phi}{\partial x} + v \phi = 0$$\phi$$S(x,t)$$(-10,10)$

$$D\frac{\partial \phi}{\partial x}(-10) + v \phi(-10) = 0$$

dla lewej strony i

$$D\frac{\partial \phi}{\partial x}(10) + v \phi(10) = 0$$

po prawej stronie. Są to tak zwane warunki brzegowe Robina (zauważ, że Wikipedia wyraźnie stwierdza, że ​​są to warunki izolacyjne dla równań doradztwa i dyfuzji).

Jeśli skonfigurujesz te warunki brzegowe, otrzymasz właściwości ochrony, których szukałeś. Rzeczywiście, mamy integrację w przestrzeni kosmicznej

$$\int \frac{\partial \phi}{\partial t} dx = \int \frac{\partial}{\partial x} \left( D \frac{\partial \phi}{\partial x} + v \phi \right) dx + \int S(x,t) dx$$

Używamy integracji z częściami po prawej stronie

$$\int \frac{\partial \phi}{\partial t} dx = \left( D \frac{\partial \phi}{\partial x} + v \phi \right)(10) - \left( D \frac{\partial \phi}{\partial x} + v \phi \right)(-10) + \int S(x,t) dx$$

Teraz dwa centralne warunki znikają dzięki warunkom brzegowym. Otrzymujemy integrując się w czasie

$$\int_0^T \int \frac{\partial \phi}{\partial t} dx dt = \int_0^T \int S(x,t) dx dt$$

a jeśli wolno nam zamienić dwie pierwsze całki ,

$$\int \phi(x,T) dx - \int \phi(x,0) dx = \int_0^T \int S(x,t) dx$$

To pokazuje, że domena jest izolowana od zewnątrz. W szczególności, jeśli , otrzymamy zachowanie .$S=0$$\phi$