Ein guter endlicher Unterschied für die Kontinuitätsgleichung

22

Was wäre eine gute Finite-Differenzen-Diskretisierung für die folgende Gleichung:

ρt+(ρu)=0?

Wir können den 1D Fall nehmen:

ρt+ddx(ρu)=0

Aus irgendeinem Grund sind alle Schemata, die ich finden kann, für die Formulierung in Lagrange-Koordinaten. Ich habe mir vorerst dieses Schema ausgedacht (ohne Berücksichtigung des j- Index):

ρi,jn+1ρi,jnτ+1hx(ρi+1,jn+1+ρi,jn+12uxi+1/2,jnρi,jn+1+ρi1,jn+12uxi1/2n)=0

Scheint aber wirklich instabil zu sein oder einen schrecklichen Stabilitätszustand zu haben. Ist das so?

Die Geschwindigkeit wird tatsächlich durch das Darcy-Gesetz u = - k berechnet. Außerdem haben wir die Zustandsgleichung. Das Gesamtsystem besteht auch aus einer Energiegleichung und der Zustandsgleichung für das ideale Gas. Die Geschwindigkeiten können negativ werden.u=kμp

tiam
quelle
Im 1D-Fall ist das Problem im Wesentlichen ein hyperbolisches pde 1. Ordnung. Haben Sie versucht, ein Finite-Differenzen-Schema erster Ordnung gegen den Wind zu verwenden?
Paul
Bisher laufe ich mit dem, was ich in der Frage geschrieben habe. Mein Fall ist eigentlich 2d. Aber da dies eine klassische Gleichung ist, dachte ich, dass es auch eine klassische Diskretisierung geben würde.
Tiam
Könnten Sie zeigen, wie ein Gegenwind-Schema dafür aussehen würde? Ich kenne das Konzept aus der Methode des endlichen Volumens, wenn Sie es konvektiv verwenden, aber dort haben Sie kein räumliches Derivat eines Produkts mehr.
Tiam
Ist das Geschwindigkeitsfeld gegeben oder erfüllt es auch eine Evolutionsgleichung?
David Ketcheson
Die Geschwindigkeit wird tatsächlich durch das Darcy-Gesetz u = - k berechnet. Das Gesamtsystem besteht auch aus einer Energiegleichung und der Zustandsgleichung für das ideale Gas. Die Geschwindigkeiten können negativ werden. u=kμp
Tiam

Antworten:

21

Sie betrachten die Massenerhaltungsgleichung:

dmdt=0

Wenn man die Massenentwicklung pro Volumeneinheit betrachtet, läuft dies auf die Dichte-Advektionsgleichung in Flussform hinaus:

ρt=(ρu)

Das Gute daran ist, dass es sich nur um die Advektionsgleichung eines beliebigen Skalarfelds handelt (in unserem Fall ist dies die Dichte ) und dass es (relativ) leicht zu lösen ist, vorausgesetzt, dass geeignete zeitliche und räumliche Differenzierungsschemata sowie initiale und Randbedingungen.ρ

Bei der Entwicklung eines Schemas mit endlichen Differenzen achten wir auf Konvergenz, Stabilität und Genauigkeit. Ein Schema konvergiert, wenn wennΔt0. Die Stabilität der Schemata stellt sicher, dass die MengeAendlich bleibt, wennt. Die formale Genauigkeit des Schemas gibt an, wo der Kürzungsfehler in der Taylor-Expansionsreihe der partiellen Ableitung liegt. In einem CFD-Lehrbuch finden Sie weitere Einzelheiten zu diesen grundlegenden Eigenschaften eines Differenzierungsschemas.ΔAΔtAtΔt0At

Der einfachste Ansatz besteht nun darin, direkt zum Upstream-Differenzieren erster Ordnung überzugehen. Dieses Schema ist eindeutig positiv, konservativ und rechnerisch effizient. Die ersten beiden Eigenschaften sind besonders wichtig, wenn wir die Entwicklung einer immer positiven Größe (dh Masse oder Dichte) modellieren.

Schauen wir uns der Einfachheit halber den 1-D-Fall an:

ρt=(ρu)x

Es ist jetzt zweckmäßig, den Fluss zu definieren , so dass:Φ=ρu

(ρu)x=ΦxΔΦΔxΦi+1/2Φi1/2Δx

Hier ist eine schematische Darstellung dessen, was wir simulieren:

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

Wir bewerten die Entwicklung von in Zelle i . Der Nettogewinn oder -verlust kommt aus der Differenz von dem, was kommt, Φ i - 1 / 2 und was geht, Φ i + 1 / 2ρiΦi1/2Φi+1/2. Hier beginnen wir, von der Antwort des Paulus abzuweichen. Bei der echten konservativen Differenzierung in Aufwärtsrichtung wird die Größe im Zellzentrum durch die Geschwindigkeit an der Zellkante in Bewegungsrichtung übertragen. Mit anderen Worten, wenn Sie sich vorstellen, dass Sie die empfohlene Menge sind und in der Mitte der Zelle sitzen, werden Sie durch die Geschwindigkeit am Zellenrand in die vor Ihnen liegende Zelle befördert. Die Auswertung des Flusses am Zellrand als Produkt aus Dichte und Geschwindigkeit am Zellrand ist nicht korrekt und konserviert die empfohlene Menge nicht.

Eingehende und ausgehende Flüsse werden wie folgt bewertet:

Φi+1/2=ui+1/2+|ui+1/2|2ρi+ui+1/2|ui+1/2|2ρi+1

Φi1/2=ui1/2+|ui1/2|2ρi1+ui1/2|ui1/2|2ρi

Die obige Behandlung der Flussdifferenzierung stellt eine Upstream-Definitivität sicher. Mit anderen Worten, es stellt die Differenzierungsrichtung gemäß dem Vorzeichen der Geschwindigkeit ein.

Das Stabilitätskriterium Courant-Friedrichs-Lewy (CFL) für die Zeitdifferenzierung mit einfacher erster Ordnung lautet wie folgt:

μ=uΔtΔx1

Beachten Sie, dass das CFL-Stabilitätskriterium in zwei Dimensionen strenger ist:

μ=cΔtΔx12

wo cu2+v2.

Einige Dinge zu beachten. Dieses Schema ist möglicherweise für Ihre Anwendung geeignet oder nicht, je nachdem, welche Art von Prozess Sie simulieren. Dieses Schema ist sehr diffusiv und eignet sich für sehr gleichmäßige Strömungen ohne scharfe Gefälle. Es ist auch diffusiver für kürzere Zeitschritte. Im 1-D-Fall erhalten Sie eine fast exakte Lösung, wenn die Gradienten sehr klein sind und wennμ=1. Im 2-D-Fall ist dies nicht möglich und die Diffusion ist anisotrop.

Wenn Ihr physisches System Stoßwellen oder hohe Gradienten anderer Art berücksichtigt, sollten Sie die Differenzierung höherer Ordnung (z. B. 3. oder 5. Ordnung) vorab untersuchen. Es kann sich auch lohnen, einen Blick auf die Familie der flusskorrigierten Transportsysteme zu werfen (Zalesak, 1979, JCP). Antidiffusionskorrektur für das obige Schema von Smolarkiewicz (1984, JCP); MPDATA-Schemafamilie von Smolarkiewicz (1998, JCP).

Für die Zeitdifferenzierung kann die Euler-Differenzierung erster Ordnung für Ihre Anforderungen zufriedenstellend sein. Andernfalls sollten Sie Methoden höherer Ordnung wie Runge-Kutta (iterativ) oder Adams-Bashforth und Adams-Moulton (mehrstufig) untersuchen.

Es lohnt sich, in einem Lehrbuch für CFD-Absolventen nach einer Zusammenfassung der oben genannten und vieler weiterer Programme zu suchen.

milancurcic
quelle
Danke für die Antwort. Jetzt sehe ich deutlich den Aufwind :). Ich werde versuchen, das jetzt umzusetzen! Ich frage mich, kann die Tatsache, dassuÄnderungen bei jedem Zeitschritt wirken sich auf die Stabilität aus?
12.
1
Nein, solange Sie die CFL-Bedingung erfüllen. Sie können entweder adaptive Zeitschritte durchführen, dΔt=Δxmeinx(u)oder setzen Δtkonstant entsprechend der maximal zu erwartenden Geschwindigkeit in Ihrem Problem. Denken Sie daran, dass verschiedene Kombinationen von Zeit- und Raumdifferenzierungsmethoden zu unterschiedlichen CFL-Einschränkungen führen.
milancurcic
Es ist ein bisschen seltsam, ich habe das Schema implementiert und es geschafft, einen Impuls von einer Grenze zur anderen und wieder zurück zu senden (durch Umkehren der Geschwindigkeit). Aber sobald ich das sageu=-Cρ Es beginnt, einen extrem kleinen Zeitschritt zu erfordern, obwohl die Geschwindigkeit unter 1 liegt. Das Festlegen des dynamischen Zeitschritts, wie oben definiert, hat auch nicht geholfen.
14.
Oder vielleicht ist es gar nicht so seltsam, vielleicht hat Ihr Kommentar oben den Fall nicht geregelt, in dem u und ρ are coupled.
tiam
The stability constraints and order of accuracy are formally derived and valid for linearized advection equation - where u does not depend on ρ. In the past, I have successfully coupled this equation with non-linear Navier-Stokes equations for u,v. Formal stability constraints are not satisfied in that case, but keeping your increments reasonably low works. When setting u=Cρ, your equation becomes ρt=C[(ρ)2+ρ2ρ]. You should investigate (if possible) what is the stability criterion for your equation.
milancurcic
13

In the 1D case, you don't want to use a forward or central difference scheme for the spatial derivative term (ddx) because they are numerically unstable. Instead, it is better to discretize the equation with an explicit backwards (upwind) finite difference for the spatial derivative:

ρik+1ρikΔt+ρikUikρi1kUi1kΔx=0.

If the velocities are positive, then this backward scheme is stable. If they are negative, then a forward difference will work. Regardless, there is always a constraint on your choice of Δx and Δt (courant number) to make the scheme stable.

Paul
quelle
Would evaluating ρ at k+1 instead remove the Δt constraint?
tiam
I'm not entirely sure... I think you would have to check the truncation error to make sure that it approximates the PDE correctly. You may want to consider the other implict schemes on this website: web.mit.edu/dongs/www/publications/projects/…
Paul