Seltsame Schwingung beim Lösen der Advektionsgleichung durch Finite-Differenzen mit vollständig geschlossenen Neumann-Randbedingungen (Reflexion an Grenzen)

32

Ich versuche, die Advektionsgleichung zu lösen, aber es erscheint eine seltsame Schwingung in der Lösung, wenn die Welle von den Grenzen reflektiert wird. Wenn jemand dieses Artefakt schon einmal gesehen hat, wäre ich daran interessiert, die Ursache zu kennen und wie man sie vermeidet!

Dies ist ein animiertes GIF, das in einem separaten Fenster geöffnet wird, um die Animation anzuzeigen (es wird nur einmal abgespielt oder nicht, sobald es zwischengespeichert wurde!) Ausbreitung eines Gaußschen Pulses

Beachten Sie, dass die Ausbreitung sehr stabil zu sein scheint, bis die Welle an der ersten Grenze zu reflektieren beginnt. Was denkst du könnte hier passieren? Ich habe einige Tage damit verbracht, meinen Code zu überprüfen, und kann keine Fehler finden. Es ist seltsam, weil es zwei sich ausbreitende Lösungen zu geben scheint: eine positive und eine negative; nach der Überlegung von der ersten Grenze. Die Lösungen scheinen sich entlang benachbarter Maschenpunkte zu bewegen.

Die Implementierungsdetails folgen.

Die Advektionsgleichung,

ut=vux

Dabei ist die Ausbreitungsgeschwindigkeit.v

Der Crank-Nicolson ist eine bedingungslose (pdf-Link) stabile Diskretisierung für die Advektionsgleichung, vorausgesetzt, dass sich im Raum langsam ändert (enthält bei Fouriertransformation nur niederfrequente Komponenten).u(x)

Die Diskretisierung, die ich angewendet habe, ist,

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

Wenn Sie die Unbekannten auf die rechte Seite setzen, können Sie dies in linearer Form schreiben.

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

Dabei gilt: (um den Zeitmittelwert zu berechnen, der zwischen dem gegenwärtigen und dem zukünftigen Punkt gleichmäßig gewichtet ist) und .r = v Δ tβ=0.5r=vΔt2Δx

Diese Mengen von Gleichungen haben die Matrixform , wobeiAun+1=Mun

A=(1βr0βr1βrβr1βr0βr1)

M=(1(1β)r0(1β)r1(1β)r(1β)r1(1β)r0(1β)r1)

Die Vektoren und sind die bekannten und unbekannten Größen, nach denen wir auflösen wollen.u n + 1unun+1

Ich wende dann geschlossene Neumann-Randbedingungen an der linken und rechten Grenze an. Mit geschlossenen Grenzen meine ich an beiden Schnittstellen. Für geschlossene Grenzen stellt sich heraus, dass (ich werde meine Arbeit hier nicht zeigen) wir nur die obige Matrixgleichung lösen müssen. Wie von @DavidKetcheson ausgeführt, beschreiben die obigen Matrixgleichungen tatsächlich Dirichlet-Randbedingungen . Für Neumann-Randbedingungenux=0

A=(100βr1βrβr1βr001)

Aktualisieren

Das Verhalten scheint ziemlich unabhängig von der Auswahl der von mir verwendeten Konstanten zu sein, aber dies sind die Werte für das Diagramm, das Sie oben sehen:

  • v = 2
  • dx = 0,2
  • dt = 0,005
  • σ = 2 (Gaußsches hwhm)
  • β = 0,5

Update II

Bei einer Simulation mit einem Diffusionskoeffizienten ungleich Null, (siehe Kommentare unten), geht die Schwingung weg, aber die Welle reflektiert nicht mehr !? Ich verstehe nicht warum?D=1

Verbreitung und Beratung

Boyfarrell
quelle
was hast du für ? v
Chris
v=2 in diesen Simulationen. Ich werde mit der Simulationseinstellung aktualisieren. Gute Idee.
Boyfarrell
Dann würde ich erwarten, dass der Anfangszustand nach rechts befördert wird und durch die rechte Grenze verschwindet. Alles, was mir einfällt, ist, dass das zentrale Schema störende Oszillationen erzeugen kann, sofern es nicht auf die Advektions-Diffusions-Gleichung mit der Zellen-Peclet-Nummer unter 2 angewendet wird.
Chris
Glauben Sie, es könnte ein Vorzeichenfehler bei der Gleichung vorliegen? Eigentlich ist es mein Endziel, dies mit der Advektions-Diffusions-Gleichung anzuwenden. Ich teste derzeit verschiedene Grenzfälle. Im obigen Beispiel wurde der Diffusionskoeffizient auf Null gesetzt. Ich habe oben eine neue Animation eingefügt. Ich verstehe nicht, warum der Peak nicht reflektiert wird, wenn der Diffusionskoeffizient ungleich Null ist. Es macht genau das, was Sie erwähnt haben (abgesehen von der Richtung).
Boyfarrell
Ich habe an gedacht , also ist das Vorzeichen in Ordnung. Die zweite Handlung sieht für mich in Ordnung aus. Warum sollten Sie erwarten, dass etwas reflektiert wird? Das könnte nur passieren, wenn sich irgendwie ändert. Versuchen Sie es mit dem Aufwindschema für Advektion anstelle des zentralen Schemas, dann sollten Sie für etwas Ähnliches sehen . v D = 0tu+vxu=0vD=0
Chris

Antworten:

28

Die Gleichung, die Sie lösen, lässt keine richtigen Lösungen zu, daher gibt es für diese Gleichung keine reflektierende Randbedingung . Wenn Sie die Merkmale berücksichtigen, werden Sie feststellen, dass Sie nur an der rechten Grenze eine Randbedingung auferlegen können. Sie versuchen, an der linken Grenze eine homogene Dirichlet-Randbedingung aufzuerlegen, die mathematisch ungültig ist.

Um es noch einmal zu wiederholen: Die Methode der Eigenschaften besagt, dass die Lösung entlang einer beliebigen Linie der Form für eine beliebige Konstante konstant sein muss . Somit wird die Lösung entlang der linken Grenze durch die Lösung zu früheren Zeitpunkten innerhalb Ihrer Problemdomäne bestimmt. Sie können dort keine Lösung aufzwingen.xνt=CC

Im Gegensatz zu der Gleichung, Ihr numerisches Schema nicht zugeben rechts laufende Lösungen. Die nach rechts gerichteten Moden werden als parasitäre Moden bezeichnet und umfassen sehr hohe Frequenzen. Beachten Sie, dass die rechtsgerichtete Welle ein Sägezahnwellenpaket ist, das den höchsten Frequenzen zugeordnet ist, die auf Ihrem Gitter dargestellt werden können. Diese Welle ist ein rein numerisches Artefakt, das durch Ihre Diskretisierung erzeugt wurde.

Zur Hervorhebung: Sie haben nicht das vollständige Problem mit den anfänglichen Grenzwerten aufgeschrieben, das Sie zu lösen versuchen. Wenn Sie dies tun, wird klar sein, dass es sich nicht um ein mathematisch gut gestelltes Problem handelt.

Ich bin froh, dass Sie dies hier gepostet haben, da es ein schönes Beispiel dafür ist, was passieren kann, wenn Sie ein Problem diskretisieren, das nicht gut gestellt ist, und das Phänomen der parasitären Modi. Ein großes +1 für deine Frage von mir.

David Ketcheson
quelle
danke für die diskussion und korrekturen. Ich hatte nicht gewürdigt, dass die Matrix nicht die gleichen Eigenschaften wie die Differentialgleichung haben wird. Weitere Kommentare folgen ...
Boyfarrell
Ja, ich sehe jetzt, wie das eigentlich Dirichlet-Randbedingungen sind, die ich oben zur Korrektur notiert habe. Dies ist das erste Mal, dass ich (im Ernst) versucht habe, den Prozess der Lösung dieser Gleichungen wirklich zu verstehen. Ich treffe immer wieder Unebenheiten auf der Straße. Ich freue mich, meine Fortschritte zu posten!
Boyfarrell
@ David Ketcheson: Ich stoße auf dasselbe Problem und habe mein Problem unter folgendem Link gepostet: scicomp.stackexchange.com/questions/30329/… Kannst du mir bitte erklären, warum du sagst, dass das Problem nicht "mathematisch gut gestellt" ist ? ? Vielen Dank.
Herman Jaramillo
@HermanJaramillo Sie versuchen, einen Lösungswert an der linken Grenze einzufügen, der bereits von der PDE bestimmt wird. In jedem Lehrbuch, in dem die Advektion besprochen wird, wird auch angegeben, welche Randbedingungen gültig sind und warum. Ich habe einen zweiten Absatz mit zusätzlicher Erklärung hinzugefügt. hoffentlich hilft das.
David Ketcheson
1
@HermanJaramillo: Nicht "mathematisch gut gestellt" bedeutet im Grunde, dass Sie zwei Gleichungen für einen Funktionswert an der Grenze, die Randbedingung sowie die PDE selbst haben. Im Allgemeinen widersprechen sich diese beiden Gleichungen. Allgemeiner kann man dies als Optimierungsproblem betrachten, bei dem beide Ziele so gut wie möglich erfüllt werden sollen.
Davidhigh
0

Aus den obigen Antworten habe ich viel gelernt. Ich möchte diese Antwort einschließen, weil ich glaube, dass sie unterschiedliche Einblicke in das Problem bietet.

uxx+1cutt=0.
c

u(x,t)=f(xct)

u(x,t0)=p(x)x(,)p[xc(tt0)]

x[a,b]ab

t0

u(a,t0)=0,u(b,t0)=p[bc(tt0)].

Dies erzeugt einen Impuls, der nach rechts läuft, bis er am rechten Rand verschwindet.

Hier klicken für Animation auf Dirichlet am linken Rand

Ich bekomme immer noch ein Geräusch, das ich nicht verstehen kann (könnte hier bitte jemand helfen?)

Die andere Möglichkeit besteht darin, periodische Randbedingungen festzulegen. Das heißt, anstatt die Dirichlet-Randbedingung auf der linken Seite aufzuerlegen, können wir das Wellenpaket auferlegen, das mit der Grenze auf der linken Seite übereinstimmt. Das ist:

u(a,t0)=p[ac(tt0)],u(b,t0)=p[bc(tt0)].

Jedoch ac(tt0)<a für t>t0c>0[a,b]baac(tt0)+ba=a+b(tt0)u(a,t0)=p[bc(tt0]

Dieser Link zeigt, was ich als periodische Randbedingungen bezeichnen würde.

Ich habe die Animationen in Python erstellt und das Schema ist ein Crank-Nicholson-Schema, wie in der Frage hier angegeben.

Ich kämpfe immer noch mit dem Rauschmuster, nachdem die Welle die erste (rechte) Grenze erreicht hat.

Herman Jaramillo
quelle
1
Ich konnte die Animation auf meinem Handy nicht sehen, aber ich glaube, Ihr Rauschmuster ist auf die fehlende numerische Genauigkeit zurückzuführen. Die Absorption funktioniert nur, weil Sie der Grenze die exakte Lösung auferlegen. Bei dieser exakten Lösung unterscheidet sich die an der rechten Grenze ankommende numerische Lösung geringfügig in Frequenz und Phase. Dies führt wiederum zu Reflexionen und damit zu Störungen.
Davidhigh
@davidhigh: Danke für deine Informationen. Ich werde das überprüfen. Es tut mir leid, dass die Animation in Ihrem Telefon nicht funktioniert hat. Auch bei meinem Handy (Samsung) hat es nicht funktioniert. Dies könnte eine fehlende Software in den Telefonen sein. Es sollte in einem Computer funktionieren. Danke noch einmal.
Herman Jaramillo
Bitte. Die Animation selbst ist nicht so wichtig, ich wollte nur sehen, wie groß die Fehler sind. Übrigens: Indem Sie der Grenze die exakte Lösung auferlegen, vermeiden Sie die "schlechte Haltung" von vornherein. Das heißt, Sie haben immer noch zwei Gleichungen für einen Wert an der Grenze, aber wenn Sie das genaue Ergebnis verwenden, erzwingen Sie, dass sie konsistent sind. Dies funktioniert aber nur annähernd, da das numerische Ergebnis nicht ganz genau ist.
Davidhigh
Und noch ein Kommentar: Die allgemeine analytische Lösung der Wellengleichung ist eine Überlagerung eines nach links und eines nach rechts bewegten Impulses. In Ihrem Fall berücksichtigen Sie nur den rechtsdrehenden Impuls, sodass Sie bereits Anfangsbedingungen angewendet haben - im Gegensatz zu den Angaben, die Sie im Text machen.
Davidhigh
@davidhigh: Ich habe ein bisschen über deine Einsichten über das Rauschen nachgedacht, nachdem der Puls die Grenze erreicht hat. Ich glaube, Sie haben Recht und es gibt einen Unterschied zwischen der exakten analytischen Lösung und dem numerisch propagierten Puls. In der Grenze erzeugt dieser kleine Unterschied das kleine Rauschmuster, das dort gesehen wird. Das in dieser Diskussion gezeigte CN-Advektionssystem ist dispersiv und ich glaube, dass Dispersion zwar nicht bemerkt wird, bevor der Puls die Grenze erreicht, aber möglicherweise die kleine Störung (Unterschied zwischen analytischen und numerischen Lösungen) an der Grenze auslöst. Danke noch einmal.
Herman Jaramillo