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

25

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.ϕ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.ϕx=0β

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 β = 0,5 (Crank-Niscolson) we wszystkich punktach

Tylko dyfuzja (granice Neumanna z beta = 0,5)

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βββ

Tylko dyfuzja (granice Neumanna z beta = 0,5 dla wnętrza, beta = 1 całkowicie niejawnie) granice

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?β

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?β

Porady i dyfuzja z = 0,5 (Crank-Niscolson) we wszystkich punktachβ

Porada-dyfuzja (granice Neumanna z beta = 0,5)

Porady i dyfuzja z = 0,5 (Crank-Niscolson) w punktach wewnętrznych, a = 1 (pełny dorozumiany) na granicachβββ

Porady i dyfuzja (granice Neumanna z beta = 0,5 dla wnętrza, beta = 1 całkowicie niejawnie) granice

Implementacja równania rada-dyfuzja

Zaczynając od równania rada-dyfuzja,

ϕt=D2ϕx2+vϕx

Pisanie przy użyciu Crank-Nicolson daje:

ϕjn+1ϕjnΔt=D[1β(Δx)2(ϕj1n2ϕjn+ϕj+1n)+β(Δx)2(ϕj1n+12ϕjn+1+ϕj+1n+1)]+v[1β2Δx(ϕj+1nϕj1n)+β2Δx(ϕj+1n+1ϕj1n+1)]

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

Aby uprościć notację, dokonajmy podstawienia,

s=DΔt(Δx)2r=vΔt2Δx

i przenieś znaną wartość pochodnej czasu na prawą stronę,ϕjn

ϕjn+1=ϕjn+s(1β)(ϕj1n2ϕjn+ϕj+1n)+sβ(ϕj1n+12ϕjn+1+ϕj+1n+1)+r(1β)(ϕj+1nϕj1n)+rβ(ϕj+1n+1ϕj1n+1)

Faktoryzacja warunków daje,ϕ

β(rs)ϕj1n+1+(1+2sβ)ϕjn+1β(s+r)ϕj+1n+1Aϕn+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1nMϕn

który możemy zapisać w postaci macierzy jako gdzie,Aϕn+1=Mϕn

A=(1+2sββ(s+r)0β(rs)1+2sββ(s+r)β(rs)1+2sββ(s+r)0β(rs)1+2sβ)

M=(12s(1β)(1β)(s+r)0(1β)(sr)12s(1β)(1β)(s+r)(1β)(sr)12s(1β)(1β)(s+r)0(1β)(sr)12s(1β))

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.ββ

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

ϕ1n+1x=σL

Pisanie tego jako centralnej różnicy daje,

ϕ1n+1xϕ2n+1ϕ0n+12Δx=σL

dlatego ϕ0n+1=ϕ2n+12Δxσ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,ϕ0n+1j=1

β(rs)ϕ0n+1+(1+2sβ)ϕ1n+1β(s+r)ϕ2n+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1n

Podstawienie wartości znalezionej w warunku brzegowym daje następujący wynik dla wiersza = 1,ϕ0n+1j

(1+2sβ)ϕ1n+12sβϕ2n+1=(1β)(sr)ϕj1n+(12s[1β])ϕjn+(1β)(s+r)ϕj+1n+2β(rs)ΔxσL

Wykonanie tej samej procedury dla ostatniego wiersza (przy = ) daje,jJ

2sβϕJ1n+1+(1+2sβ)ϕJn+1=(1β)(sr)ϕJ1n+(12s(1β))ϕJn+2β(s+r)ΔxσR

Ostatecznie ukrywanie wierszy granic (ustawienie = 1) daje,β

(1+2s)ϕ1n+12sϕ2n+1=ϕj1n+1ϕjn+2(rs)ΔxσL

2sϕJ1n+1+(1+2s)ϕJn+1=ϕJn+2(s+r)ΔxσR

Dlatego w warunkach brzegowych Neumanna możemy zapisać równanie macierzowe, ,Aϕn+1=Mϕn+bN

gdzie,

A=(1+2s2s0β(rs)1+2sββ(s+r)β(rs)1+2sββ(s+r)02s1+2s)

M=(100(1β)(sr)12s(1β)(1β)(s+r)(1β)(sr)12s(1β)(1β)(s+r)001)

bN=(2(rs)ΔxσL002(s+r)ΔxσR)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ϕx=0j=Dϕx+vϕ=0.

Zgodziłbyś się?

boyfarrell
źródło
Wydaje się, że warunki brzegowe nie są poprawnie realizowane. Czy możesz nam pokazać, w jaki sposób nałożyłeś warunki brzegowe?
David Ketcheson
OK Zaktualizowałem implementację i myślę, że zauważyłem błąd dotyczący zastosowania = 0,5 tylko w wierszach granicznych. Zaktualizowałem moje „bieżące rozumienie” u dołu pytania. Czy masz jakiś komentarz? β
boyfarrell
Więc ... jak wygląda dyskretyzacja na granicach w przypadku granic Robin? Pokazałeś to dla granic Neumanna, ale nie dla granic Robin.

Odpowiedzi:

15

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

ϕt=x(Dϕx+vϕ)+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ϕx+vϕ=0ϕS(x,t)(10,10)

Dϕx(10)+vϕ(10)=0

dla lewej strony i

Dϕx(10)+vϕ(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

ϕtdx=x(Dϕx+vϕ)dx+S(x,t)dx

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

ϕtdx=(Dϕx+vϕ)(10)(Dϕx+vϕ)(10)+S(x,t)dx

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

0Tϕtdxdt=0TS(x,t)dxdt

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

ϕ(x,T)dxϕ(x,0)dx=0TS(x,t)dx

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

Dr_Sam
źródło
Rozumiem teraz, dlaczego działało to tylko wtedy, gdy = 0; ponieważ oznaczałoby to ochronę zgodnie z powyższym podejściem. Jaka byłaby konsekwencja zastosowania powyższego warunku brzegowego, gdyby fala odbiła? Myślałem, że to nie będzie możliwe, ponieważ w równaniu nie ma nic, co dałoby mi ujemną prędkość? v
boyfarrell
Najlepszym sposobem, aby się dowiedzieć, jest prawdopodobnie spróbować! Ale jeśli to działa poprawnie (i robi to IMO), powinieneś zobaczyć pewną ilość która zaczyna się gromadzić po lewej stronie domeny: wskazówka popycha w tym kierunku, ale granica jest zamknięta. Akumulacja zatrzymuje się, gdy dyfuzja jest wystarczająco duża, aby ją zrównoważyć. Więc nie, nie powinno być fali odbitej. ϕϕ
Dr_Sam
@DrSam Tylko pytanie dotyczące wdrożenia. Rozumiem twój punkt widzenia, jak zrobić zero po lewej stronie. Ale kiedy mówisz „po prawej stronie tylko termin”, co to oznacza? Myślałem, że warunki brzegowe powinny być albo Neumann, albo Dirichlet (lub połączenie obu)?
boyfarrell
@boyfarrell Lewo / prawo w odpowiedzi odnosiło się do wyprowadzenia prawidłowych warunków brzegowych, a nie sposobu, w jaki jest ono implementowane (edytowane dla jasności). Warunki Robin są warunkami klasycznymi, choć są mniej znane niż Dirichlet i Neumann.
Dr_Sam
Więc jeśli chodzi o wdrożenie, czy uważasz, że powinienem wyprowadzić warunki brzegowe Robina dla obu granic? Ponadto, jeśli równanie zawiera warunek reakcji (np. ma warunek brzegowy potrzeby także ten termin?
ϕt=x(Dϕx+vϕ)+S(x,t)
boyfarrell