Erhaltung einer physikalischen Größe unter Verwendung von Neumann-Randbedingungen, die auf die Advektions-Diffusions-Gleichung angewendet werden

25

Ich verstehe das unterschiedliche Verhalten der Advektions-Diffusions-Gleichung nicht, wenn ich unterschiedliche Randbedingungen anwende. Meine Motivation ist die Simulation einer realen physikalischen Größe (Teilchendichte) unter Diffusion und Advektion. Die Teilchendichte sollte im Inneren erhalten bleiben, es sei denn, sie fließt aus den Rändern heraus. Wenn ich nach dieser Logik Neumann-Randbedingungen wie (auf der linken und rechten Seite) erzwinge , sollte das System "geschlossen" sein, d. Wenn der Fluss an der Grenze Null ist, können keine Partikel entweichen.ϕx=0

Für alle folgenden Simulationen habe ich die Crank-Nicolson-Diskretisierung auf die Advektions-Diffusions-Gleichung angewendet und alle Simulationen haben Randbedingungen. Für die erste und letzte Zeile der Matrix (die Randbedingungszeilen) erlaube ich jedoch, dass unabhängig vom inneren Wert geändert wird. Dadurch können die Endpunkte vollständig implizit sein.βϕx=0β

Im Folgenden werden 4 verschiedene Konfigurationen besprochen, von denen nur eine den Erwartungen entspricht. Am Ende diskutiere ich meine Implementierung.

Diffusionslimit

Hier werden die Advektionsterme ausgeschaltet, indem die Geschwindigkeit auf Null gesetzt wird.

Nur Diffusion mit = 0,5 (Crank-Niscolson) an allen Punktenβ

Nur Diffusion (Neumann-Grenzen mit Beta = 0,5)

Die Menge bleibt nicht erhalten, wie die Verringerung der Impulsfläche zeigt.

Nur Diffusion mit = 0,5 (Crank-Niscolson) an den inneren Punkten und = 1 (vollständig implizit) an den Grenzenββ

Nur die Diffusion (Neumann-Grenzen mit Beta = 0,5 für das Innere, Beta = 1 vollständig implizit) die Grenzen

Durch die Verwendung einer vollständig impliziten Gleichung an den Grenzen erreiche ich das, was ich erwarte: Keine Partikel entweichen . Sie können dies daran erkennen, dass die Fläche als Partikel diffus konserviert wird. Warum sollte die Wahl von an den Grenzpunkten die Physik der Situation beeinflussen? Ist das ein Bug oder erwartet?β

Verbreitung und Beratung

Wenn der Advektionsbegriff eingeschlossen ist, scheint der Wert von an den Grenzen die Lösung nicht zu beeinflussen. In allen Fällen, in denen die Grenzen offen zu sein scheinen, können Partikel den Grenzen entkommen. Warum ist das so?β

Advektion und Diffusion mit = 0,5 (Crank-Niscolson) an allen Punktenβ

Advektionsdiffusion (Neumann-Grenzen mit Beta = 0,5)

Advektion und Diffusion mit = 0,5 (Crank-Niscolson) an inneren Punkten und = 1 (vollständig implizit) an den Grenzenββ

Advektion und Diffusion (Neumann-Grenzen mit Beta = 0,5 für das Innere, Beta = 1 vollständig implizit) die Grenzen

Implementierung der Advektions-Diffusions-Gleichung

Ausgehend von der Advektions-Diffusions-Gleichung

ϕt=D2ϕx2+vϕx

Schreiben mit Crank-Nicolson gibt,

ϕ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)]

Beachten Sie, dass = 0,5 für Crank-Nicolson, = 1 für vollständig implizite und = 0 für vollständig explizite.βββ

Um die Notation zu vereinfachen, nehmen wir die Ersetzung vor,

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

und verschiebe den bekannten Wert der Zeitableitung nach rechts,ϕ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)

Die Berücksichtigung der Terme ergibt,ϕ

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

die wir in Matrixform wie schreiben kann , wo,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β))

Neumann-Randbedingungen anwenden

NB arbeitet die Ableitung erneut ab. Ich glaube, ich habe den Fehler entdeckt. Ich habe beim Schreiben der endlichen Differenz der Randbedingung ein vollständig implizites Schema ( = 1) angenommen. Wenn Sie hier ein Crank-Niscolson-Schema annehmen, wird die Komplexität zu groß und ich könnte die resultierenden Gleichungen nicht lösen, um die Knoten zu beseitigen, die außerhalb der Domäne liegen. Es scheint jedoch möglich, dass es zwei Gleichungen mit zwei Unbekannten gibt, aber ich konnte es nicht schaffen. Dies erklärt wahrscheinlich den Unterschied zwischen dem ersten und dem zweiten Diagramm oben. Ich denke, wir können daraus schließen, dass nur die Diagramme mit = 0.5 an den Grenzpunkten gültig sind.ββ

Unter der Annahme, dass der Fluss auf der linken Seite bekannt ist (unter Annahme einer vollständig impliziten Form),

ϕ1n+1x=σL

Schreiben Sie dies als zentrierte Differenz gibt,

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

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

Beachten Sie, dass dadurch ein Knoten der sich außerhalb der Domäne des Problems befindet. Dieser Knoten kann mit einer zweiten Gleichung eliminiert werden. Wir können den Knoten schreiben alsϕ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

Einsetzen des Wertes von aus der Randbedingung ergibt folgendes Ergebnis für die Zeile = 1,ϕ0n+1j

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

Durch Ausführen des gleichen Verfahrens für die letzte Reihe (bei = ) erhält manjJ

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

Wenn Sie die Randzeilen implizit machen (setze = 1), erhalten Sie:β

(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

Daher mit Neumann - Randbedingungen können wir die Matrixgleichung schreiben, ,Aϕn+1=Mϕn+bN

woher,

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

Mein aktuelles Verständnis

  • Ich denke, der Unterschied zwischen dem ersten und dem zweiten Plot erklärt sich aus dem oben beschriebenen Fehler.

  • In Bezug auf die Erhaltung der physikalischen Größe. Ich glaube, die Ursache liegt darin, dass, wie hier ausgeführt , die Advektionsgleichung in der von mir geschriebenen Form keine Ausbreitung in der umgekehrten Richtung zulässt, sodass die Welle auch bei Bedingungen mit null Fluss nur durchläuft . Meine anfängliche Vorstellung von der Erhaltung galt nur, wenn der Advektionszeitraum Null ist (dies ist die Lösung in Parzelle Nr. 2, in der die Fläche erhalten bleibt).

  • Auch bei Neumann-Nullfluss- Randbedingungen die Masse das System noch verlassen, weil die korrekten Randbedingungen in diesem Fall Robin- Randbedingungen sind, bei denen der Gesamtfluss ist angegeben ist . Darüber hinaus gibt die Neunmann-Bedingung an, dass die Masse die Domäne nicht durch Diffusion verlassen kann , und sagt nichts über die Advektion aus. Wir hören im Wesentlichen geschlossene Randbedingungen für die Diffusion und offene Randbedingungen für die Advektion. Weitere Informationen finden Sie in der Antwort hier, Implementierung der Randbedingung für den Gradienten Null in der Advektions-Diffusions-Gleichungϕx=0j=Dϕx+vϕ=0.

Würdest du zustimmen?

Boyfarrell
quelle
Es scheint, dass die Randbedingungen nicht korrekt implementiert sind. Können Sie uns zeigen, wie Sie die Randbedingungen auferlegt haben?
David Ketcheson
OK, ich habe die Implementierung aktualisiert und glaube, ich habe den Fehler beim Anwenden von = 0.5 nur an den Begrenzungszeilen festgestellt. Ich habe mein "aktuelles Verständnis" am Ende der Frage aktualisiert. Hast du einen Kommentar? β
Boyfarrell
Also ... wie sieht die Diskretisierung an den Grenzen im Falle der Robin-Grenzen aus? Sie haben es für die Neumann-Grenzen gezeigt, aber nicht für die Robin-Grenzen.

Antworten:

15

Ich denke, eines Ihrer Probleme ist, dass (wie Sie in Ihren Kommentaren bemerkt haben) Neumann-Bedingungen nicht die Bedingungen sind, nach denen Sie suchen , in dem Sinne, dass sie nicht die Erhaltung Ihrer Quantität bedeuten. Um den richtigen Zustand zu finden, schreiben Sie Ihre PDE wie folgt um

ϕt=x(Dϕx+vϕ)+S(x,t).

Der in Klammern gesetzte Ausdruck ist der Gesamtfluss, und dies ist die Menge, die Sie an den Grenzen auf Null setzen müssen, um . (Ich habe der Allgemeinheit halber und für Ihre Kommentare hinzugefügt .) Die Randbedingungen, die Sie auferlegen müssen, sind dann (angenommen, Ihre ist )Dϕx+vϕ=0ϕS(x,t)(10,10)

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

für die linke Seite und

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

für die rechte Seite. Dies sind die sogenannten Robin-Randbedingungen (Beachten Sie, dass Wikipedia ausdrücklich sagt, dass dies die isolierenden Bedingungen für Advektions-Diffusions-Gleichungen sind).

Wenn Sie diese Randbedingungen festlegen, erhalten Sie die gewünschten Konservierungseigenschaften. In der Tat haben wir die Integration über den Raumbereich

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

Mit der Integration von Teilen auf der rechten Seite haben wir

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

Nun verschwinden die beiden zentralen Begriffe dank der Randbedingungen. Durch die Integration in die Zeit erhalten wir

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

und wenn wir die ersten beiden Integrale wechseln dürfen ,

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

Dies zeigt, dass die Domäne von außen isoliert ist. Insbesondere wenn , erhalten wir die Erhaltung von .S=0ϕ

Dr_Sam
quelle
Mir ist jetzt klar, warum das nur funktioniert hat, wenn = 0; denn dies würde eine Erhaltung nach Ihrem obigen Ansatz bedeuten. Was wäre die Konsequenz dieser Randbedingung, würde die Welle reflektieren? Ich dachte, das wäre nicht möglich, weil die Gleichung nichts enthält, was mir eine negative Geschwindigkeit verleiht. v
Boyfarrell
Der beste Weg, es zu wissen, ist wahrscheinlich, es zu versuchen! Wenn sich dies jedoch korrekt verhält (und IMO tut es dies), sollten Sie eine bestimmte Menge von , die sich auf der linken Seite der Domäne ansammelt: Die Advektion drückt in diese Richtung, aber die Grenze ist geschlossen. Die Akkumulation stoppt, wenn die Diffusion groß genug ist, um sie auszugleichen. Also nein, es sollte keine reflektierte Welle geben. ϕϕ
Dr_Sam
@DrSam Nur eine Frage zur Implementierung. Ich verstehe Ihren Punkt darüber, wie man die Menge auf der linken Seite zu Null macht. Aber wenn Sie "rechts nur einen Grenzbegriff" sagen, was bedeutet das? Ich dachte, dass die Randbedingungen entweder Neumann oder Dirichlet sein sollten (oder eine Mischung aus beiden)?
Boyfarrell
@boyfarrell Das Links / Rechts in der Antwort bezog sich auf eine Ableitung der korrekten Randbedingungen, nicht auf die Art und Weise, wie sie implementiert sind (der Übersichtlichkeit halber bearbeitet). Robin Bedingungen sind klassische Bedingungen, obwohl es weniger bekannt als Dirichlet und Neumann sind.
Dr_Sam
Denken Sie, dass ich bezüglich der Implementierung Robin-Randbedingungen für beide Grenzen ableiten sollte? Auch wenn die Gleichung einen Reaktionsterm hat (zB
ϕt=x(Dϕx+vϕ)+S(x,t)
Muss