9 Fallstudie: 2D-Kontinuitätsgleichung

Betrachten verschiedene (auch ungeeignete) Verfahren zum Lösen der Kontinuitätsgleichung in zwei Raumdimensionen.

9.1Problemstellung¶

9.1.1Gleichung¶

Die Kontinuitätsgleichung

beschreibt z.B. die Masseerhaltung beim Materialtransport in einer Strömung. Dabei ist $u(x,y,t)$ die Dichte des transportierten Materials (z.B. eine Ölschicht auf einer Wasseroberfläche) im Punkt $(x,y)$ zur Zeit $t$. Die Werte $v(x,y)$ und $w(x,y)$ beschreiben gemeinsam als Vektor Strömungsrichtung und Strömungsgeschwindigkeit (z.B. einer Wasseroberfläche). Zur durch die Funktionen $v$ und $w$ gegebenen Strömung ist die Funktion $u$ gesucht.

Zu beachten ist, dass die Kontinuitätsgleichung keine Diffussionsprozesse abbildet, sondern lediglich den Transport durch Strömung. Ist z.B. die Wasseroberfläche strömungsfrei, so wird die anfängliche Ölverteilung in der Lösung der Gleichung über die Zeit erhalten bleiben. Praktisch ist allerdings davon auszugehen, dass der Ölteppich “breitfließt”, also von sich aus die Dichte verringert und die beanspruchte Fläche vergrößert. Für die kombinierte Betrachtung von Transport- und Diffusionsprozessen kann man auch eine PDE herleiten und lösen. Der Einfachheit halber verzichten wir hier jedoch auf den Diffusionsanteil. Praktisch bedeutet dies, dass unser Ölteppich bereits eine Dichte erreicht hat, bei der die Ölschicht ohne externe Kräfte (Strömung) sich nicht mehr verändert.

Ebenfalls zur Vereinfachung beschränken wir uns auf eine im Zeitverlauf unveränderliche Strömung. Die Funktionen $v$ und $w$ hängen also nur vom Ort, aber nicht von der Zeit ab. Ein Teil der unten betrachteten Lösungsverfahren ist auch für zeitlich veränderliche Strömungen einsetzbar.

9.1.2Anfangsbedingung¶

Die Anfangskonzentration des durch die Strömung transportierten Stoffs kann durch eine vom Ort abhängige Funktion $f$ ausgedrückt werden:

9.1.3Randbedingungen¶

Das Formulieren von sinnvollen Randbedingungen ist schwierig. Folgende Bedinungungen sollen am Rand erfüllt sein:

• Kein Material strömt vom Rand des betrachteten Gebietes $B\subseteq\bbR^2$ ein.

• Material kann frei über den Rand mit der Strömung abtransportiert werden.

Je nach Strömungsrichtung in Bezug auf den Rand, ist die Randbedingung also “kein Fluss über den Rand” oder “Fluss über den Rand in unbekannter Höhe”. Diese Randbedingungen lassen sich nur mühsam in Formeln ausdrücken. Wir werden Sie deshalb je nach Lösungsverfahren individuell in die Algorithmen integrieren.

Wir interpretieren hier $B$ als Ausschnitt aus einem räumlich größeren Strömungsgeschehen. Durch einen Unfall auf einer Ölbohrplattform könnte Öl freigesetzt worden sein und wir wollen die Ausbreitungs des Ölteppichs in einem gewissen Bereich unter gegebenen Strömungsverhältnissen vorhersagen. Der Bereich soll dabei so groß sein, dass der anfänglich der komplette Ölteppich darin enthalten ist. Im Laufe der Zeit wird das Öl den Rand erreichen und unseren Vorhersagebereich verlassen.

Wählen

9.1.4Inkompressible Strömung¶

Zur Vereinfachung beschränken wir uns auf inkompressible Strömungen, also solche, bei denen Zufluss und Abfluss in jedem Punkt identisch sind. In der Sprache der Differentialoperatoren soll die Divergenz des Strömungsfeldes überall Null sein:

Auf eine Herleitung dieses Zusammenhangs verzichten wir hier (siehe z.B. Wikipedia zu “Incompressible flow” für eine Variante).

Bezogen auf das Beispiel eines Ölteppichs auf einer Wasserobfläche bedeutet Inkompressibilität insbesondere, dass keine Strömung in vertikale Richtung stattfindet. Es kann also kein Wasser (und damit Öl) nach unten “verschwinden” oder von unten “auftauchen”.

Durch die Annahme einer inkompressiblen Strömung kann die Kontinuitätsgleichung zu

vereinfacht werden (IDVID 910).

9.2Diskretisierung¶

Setzen

• $\Delta x:=\frac{a}{n_x}$,

• $\Delta y:=\frac{b}{n_y}$,

• $\Delta t:=\frac{T}{n_t}$

und wählen in den Variablen $x$, $y$, $z$ gleichmäßig verteilte Stützstellen

• $x_i:=i\,\Delta x$, $\quad i=0,1,\ldots,n_x$,

• $y_j:=i\,\Delta y$, $\quad j=0,1,\ldots,n_y$,

• $t_k:=i\,\Delta t$, $\quad k=0,1,\ldots,n_t$.

Bezeichnen die Funktionswerte an den Stützstellen entsprechend mit

• $u_{i,j,k}:=u(x_i,y_j,t_k)$,

• $v_{i,j}:=v(x_i,y_j)$,

• $w_{i,j}:=w(x_i,y_j)$,

• $f_{i,j}:=f(x_i,y_j)$.

9.3Einfaches Upwind-Verfahren¶

9.3.1Idee¶

Ersetzt man $u_t$ durch eine Vorwärtsdifferenz und $u_x$ und $u_y$ durch Rückwärtsdifferenzen, so erhält man

(IDVID 920). Die Zahlen

heißen Courant-Zahlen.

Zur Interpretation betrachten wir eine Strömung parallel zur $x$-Achse, also $w\equiv 0$:

Wie sehen hier:

• Sind alle $v_{i,j}$ positiv, so ändert sich die Materialmenge an der Stelle $(x_i,y_j)$ pro Zeitschritt indem der Anteil $v_{i,j}\,\frac{\Delta t}{\Delta x}\,u_{i,j,k}$ der vorhandenen Menge $u_{i,j,k}$ abfließt und der Anteil $v_{i,j}\,\frac{\Delta t}{\Delta x}\,u_{i-1,j,k}$ aus Strömungsrichtung hinzukommt.

• Es sollte $C_x\leq 1$ gelten, da sonst mehr Material pro Zeitschritt abfließen würde als vorhanden ist. Der Abstand $\Delta t$ der Zeitstützstellen muss also klein genug sein. Wir können $\frac{\Delta x}{\Delta t}$ als höchste mit unserer Diskretisierung simulierbare Strömungsgeschwindigkeit interpretieren, denn die Bedingung $C_x\leq 1$ besagt

  $$|v_{i,j}|\leq\frac{\Delta x}{\Delta t}\quad\text{für alle }i,j.$$

  (9)

• Sind die $v_{i,j}$ negativ, so wird der Materialzuwachs pro Zeitschritt aus der stromabwärts (!) gelegenen Materialmenge $u_{i-1,j,k}$ berechnet. Dies ist zwar theoretisch eine valide Approximation der Ableitung in $x$-Richtung, praktisch aber wenig sinnvoll. Deshalb sollte bei negativem $v_{i,j}$ eine Vorwärtsdifferenz verwendet werden.

9.3.2Umsetzung¶

Setzen wir

und

so erhalten wir die diskretisierte Gleichung

(IDVID 625). Dies ist das sogenannte Upwind-Verfahren (Upwind = gegen den Wind bzw. die Strömung)

Die Randbedingungen können hier leicht integriert werden. Liegt $(x_i,y_j)$ an einem Rand mit Zustrom in das Gebiet, so setzt man $u_{i,j,k}^{\mathrm{in},x}:=0$ und $u_{i,j,k}^{\mathrm{in},y}:=0$. An Rändern mit Abfluss aus dem Gebiet ist nichts zu tun.

9.3.3Kondition und Konvergenz¶

Man kann zeigen, dass die Berechnungen im Upwind-Verfahren gut konditioniert sind, wenn

gilt.

Allerdings liegt keine Konvergenz vor, da die Materialverteilung im Laufe der Zeit “breitläuft”. Hintergrund ist, dass man die gewählte Diskretisierung der Kontinuitätsgleichung auch als Diskretisierung (mittels zentraler Differenzen) einer kombinierten Kontinuitäts-Diffusions-Gleichung interpretieren kann. Zwei verschiedene PDE führen also zum selben diskretisierten Problem. Konvergenz kann aber höchstens bezüglich einer dieser beiden PDE vorliegen.

9.3.4Implementierung¶

Schwierigster Punkt bei der Implementierung ist die Fallunterscheidung je nach Strömungsrichtung. Eine effiziente Lösung in Python (also ohne explizite Schleifen) ist mittels Bool’scher Indizierung von NumPy-Arrays möglich.

import numpy as np
import matplotlib.pyplot as plt

def cont2d_upwind(a, b, T, v, w, f, nx, ny, nt):

    delta_x = a / nx
    delta_y = b / ny
    delta_t = T / nt
    
    x = np.linspace(0, a, nx + 1)
    y = np.linspace(0, b, ny + 1)
    grid_x, grid_y = np.meshgrid(x, y, indexing='ij')  # x ist erste Dimension

    v = v(grid_x, grid_y)
    w = w(grid_x, grid_y)
    Cx = np.abs(v).max() * delta_t / delta_x
    Cy = np.abs(w).max() * delta_t / delta_y
    print(f'Summe der Courant-Zahlen: {Cx + Cy}')

    u = np.empty((nt + 1, nx + 1, ny + 1), dtype=float)

    # Anfangsbedingungen
    u[0, :, :] = f(grid_x, grid_y)

    # Zeitschritte
    v_pos_mask = v >= 0
    v_pos_mask[0, :] = False  # kein Zufluss von linkem Rand
    v_neg_mask = np.logical_not(v_pos_mask)
    v_neg_mask[-1, :] = False  # kein Zufluss von rechtem Rand
    w_pos_mask = w >= 0
    w_pos_mask[:, 0] = False  # kein Zufluss vom oberen Rand
    w_neg_mask = np.logical_not(w_pos_mask)
    w_neg_mask[:, -1] = False  # kein Zufluss vom unteren Rand
    
    u_inx = np.zeros((nx + 1, ny + 1), dtype=float)
    u_iny = np.zeros((nx + 1, ny + 1), dtype=float)
    for k in range(0, nt):

        u_inx[0, :] = 0
        u_inx[-1, :] = 0
        u_inx[v_pos_mask] = u[k, :-1, :][v_pos_mask[1:, :]]
        u_inx[v_neg_mask] = u[k, 1:, :][v_neg_mask[:-1, :]]
        u_iny[:, 0] = 0
        u_iny[:, -1] = 0
        u_iny[w_pos_mask] = u[k, :, :-1][w_pos_mask[:, 1:]]
        u_iny[w_neg_mask] = u[k, :, 1:][w_neg_mask[:, :-1]]

        u[k + 1, :, :] = u[k, :, :] \
                         - np.abs(v) * delta_t / delta_x * (u[k, :, :] - u_inx) \
                         - np.abs(w) * delta_t / delta_y * (u[k, :, :] - u_iny)
        
    return u

Neben der Lösung $u$ visualisieren wir auch das Strömungsfeld.

import matplotlib.animation as anim
from IPython.display import HTML

def show_anim(u, a, b, T, v, w, quiver_step=1):

    nt = u.shape[0] - 1
    nx = u.shape[1] - 1
    ny = u.shape[2] - 1
    x = np.linspace(0, a, nx + 1)
    y = np.linspace(0, b, ny + 1)
    grid_x, grid_y = np.meshgrid(x, y, indexing='ij')  # x ist erste Dimension

    v = v(grid_x, grid_y)
    w = w(grid_x, grid_y)

    # Bildrate auf 25 fps begrenzen
    fps_real = nt / T
    fps_show = 25
    skip_frames = int(np.floor(fps_real / fps_show))
    n_frames = (nt + 1) // (skip_frames + 1)
    
    # Figure/Axes erstellen
    fig, ax = plt.subplots()
    ax.set_xlim(0, a)
    ax.set_ylim(0, b)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    
    # Update-Funktion
    def update(frame):
        ax.clear()
        ax.contourf(
            grid_x, grid_y, u[(skip_frames + 1) * frame, :, :],
            levels=100, cmap='jet', vmin=0, vmax=1
        )
        ax.quiver(
            grid_x[::quiver_step, ::quiver_step],
            grid_y[::quiver_step, ::quiver_step],
            v[::quiver_step, ::quiver_step],
            w[::quiver_step, ::quiver_step],
            angles='xy', pivot='mid', color='#ffffff'
        )
        ax.axis('equal')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
    
    # Animation erstellen und anzeigen
    fa = anim.FuncAnimation(fig, update, frames=n_frames, interval=1000/fps_show)
    display(HTML(fa.to_jshtml()))
    plt.close()  # verhindert automatischen Aufruf von plt.show() durch Jupyter

Als Anfangsverteilung $f$ verwenden wir einen kreisförmigen Bereich mit zum Rand linear abfallender Konzentration.

Eine einfache Strömung mit konstanter Richtung und konstanter Geschwindigkeit:

a = 3
b = 2
T = 2

nx = 40
ny = int(np.ceil(nx / a * b))  # \delta y = \delta x
grid_x, grid_y = np.meshgrid(
    np.linspace(0, a, nx + 1),
    np.linspace(0, b, ny + 1),
    indexing='ij'
)

v = lambda x, y: np.ones_like(x)
w = lambda x, y: np.ones_like(x)

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.3) ** 2 + (y/b - 0.3) ** 2))

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(T * (v_max * nx / a + w_max * ny / b)))  # Courant <= 1

u = cont2d_upwind(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Summe der Courant-Zahlen: 0.9938271604938271

Wir sehen, dass der Ölteppich breitläuft wie wir oben schon festgestellt hatten, obwohl das kontinuierliche Problem gar keine Diffusion modelliert. Zusätzlich verändert sich auch die Form: quer zur Strömung läuft der Ölteppich stärker breit als in Strömungsrichtung. Mehr dazu beim nächsten Verfahren.

Betrachten noch eine kreisförmige Strömung. Nach Ablauf des Zeitintervalls sollte bei exakter Rechnung wieder die Anfangssituation vorliegen, was offensichtlich nicht der Falls ist.

T = 2 * np.pi

def r_phi(x, y):
    x = x - 0.5 * a
    y = y - 0.5 * b
    r = np.sqrt(x ** 2 + y ** 2)
    phi = np.arctan2(y, x)
    return r, phi
def v(x, y):
    r, phi = r_phi(x, y)
    return -r * np.sin(phi)
def w(x, y):
    r, phi = r_phi(x, y)
    return r * np.cos(phi)

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(T * (v_max * nx / a + w_max * ny / b)))  # Courant <= 1

u = cont2d_upwind(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Summe der Courant-Zahlen: 0.9953316347458231

9.4Verbessertes Upwind-Verfahren¶

9.4.1Idee¶

Hatten beim Upwind-Verfahren beobachtet, dass der Ölteppich quer zur Strömung stärker breitläuft als in Strömungsrichtung. Dies liegt na der gewählten Diskretisierung: Die Strömung wird separat in $x$- und in $y$-Richtung simuliert. Ist in einem Punkt $(x_i,y_j)$ die Strömungsrichtung $(1, 1)$ (also 45°), so müsste eigentlich Material vom Punkt $(x_{i-1},y_{j-1})$ (links unten) einströmen. In der Simulation kommt aber sämtlicher Zufluss nur von $(x_{i-1},y_j)$ (links) und $(x_i,y_{j-1})$ (unten). Entsprechend wirkt der Breitlaufen stärker nach links und unten (also nicht in Strömungsrichtung) als direkt nach links unten (also in Strömungsrichtung).

Um diesen Effekt mit einfachen Mitteln zu veringern, kann man die Berechnungen in $x$- und $y$-Richtung zeitlich nacheinander durchführen. Es wird also erst der Zufluss in $x$-Richtung in jedem Punkt berechnet und anschließend im nächsten Zeitschritt der Zufluss in $y$-Richtung. Dadurch fließt bei Betrachtung eines Punktes $(x_i,y_j)$ zunächst Material von $(x_{i-1},y_{j-1})$ (links unten) nach $(x_i,y_{j-1})$ (unten) und dann nach oben in den betrachteten Punkt. Der Punkt $(x_i,y_j)$ bekommt nun also Zufluss von links, links unten und unten.

Dieses Vorgehen wird auch als Corner-Transport-Upwind-Verfahren (CTU) bezeichnet.

9.4.2Umsetzung¶

Um die Anzahl der Zeitschritt nicht zu verdoppeln verwendet man bei der Darstellung dieses Vorgehens in Formeln gern halbe Zeitschritte:

9.4.3Implementierung¶

def cont2d_corner(a, b, T, v, w, f, nx, ny, nt):

    delta_x = a / nx
    delta_y = b / ny
    delta_t = T / nt
    
    x = np.linspace(0, a, nx + 1)
    y = np.linspace(0, b, ny + 1)
    grid_x, grid_y = np.meshgrid(x, y, indexing='ij')  # x ist erste Dimension

    v = v(grid_x, grid_y)
    w = w(grid_x, grid_y)
    Cx = np.abs(v).max() * delta_t / delta_x
    Cy = np.abs(w).max() * delta_t / delta_y
    print(f'Summe der Courant-Zahlen: {Cx + Cy}')

    u = np.empty((nt + 1, nx + 1, ny + 1), dtype=float)

    # Anfangsbedingungen
    u[0, :, :] = f(grid_x, grid_y)

    # Zeitschritte
    v_pos_mask = v >= 0
    v_pos_mask[0, :] = False  # kein Zufluss von linkem Rand
    v_neg_mask = np.logical_not(v_pos_mask)
    v_neg_mask[-1, :] = False  # kein Zufluss von rechtem Rand
    w_pos_mask = w >= 0
    w_pos_mask[:, 0] = False  # kein Zufluss vom oberen Rand
    w_neg_mask = np.logical_not(w_pos_mask)
    w_neg_mask[:, -1] = False  # kein Zufluss vom unteren Rand
    
    u_inx = np.zeros((nx + 1, ny + 1), dtype=float)
    u_iny = np.zeros((nx + 1, ny + 1), dtype=float)
    for k in range(0, nt):

        u_inx[0, :] = 0
        u_inx[-1, :] = 0
        u_inx[v_pos_mask] = u[k, :-1, :][v_pos_mask[1:, :]]
        u_inx[v_neg_mask] = u[k, 1:, :][v_neg_mask[:-1, :]]

        u_tmp = u[k, :, :] - np.abs(v) * delta_t / delta_x * (u[k, :, :] - u_inx)
        
        u_iny[:, 0] = 0
        u_iny[:, -1] = 0
        u_iny[w_pos_mask] = u_tmp[:, :-1][w_pos_mask[:, 1:]]
        u_iny[w_neg_mask] = u_tmp[:, 1:][w_neg_mask[:, :-1]]

        u[k + 1, :, :] = u_tmp - np.abs(w) * delta_t / delta_y * (u_tmp - u_iny)
        
    return u

Nun bleibt der Ölteppich kreisförmig:

a = 3
b = 2
T = 2

nx = 40
ny = int(np.ceil(nx / a * b))  # \delta y = \delta x
grid_x, grid_y = np.meshgrid(
    np.linspace(0, a, nx + 1),
    np.linspace(0, b, ny + 1),
    indexing='ij'
)

v = lambda x, y: np.ones_like(x)
w = lambda x, y: np.ones_like(x)

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.3) ** 2 + (y/b - 0.3) ** 2))

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(T * (v_max * nx / a + w_max * ny / b)))  # Courant <= 1

u = cont2d_corner(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Summe der Courant-Zahlen: 0.9938271604938271

9.5Zentrale Differenzen¶

9.5.1Idee¶

Um eventuell besser Ergebnisse zu erzielen könnte man zentrale Differenzenquotienten verwenden. In der Zeit ist dies ungünstig, da wir dann das entstehende Gleichungssystem nicht mehr einfach durch Umstellen lösen können. In den Raumkoordinaten spricht aber erstmal nichts dagegen. Insbesondere ist damit die Fallunterscheidung je nach Strömungsrichtung nicht mehr nötig. Diesen Ansatz nennt man auch FTCS-Verfahren (für “forward time centered space”).

Es wird sicher allerdings zeigen, dass die Lösung des diskretisierten Problems schlecht konditioniert ist. Das Verfahren ist also praktisch nicht einsetzbar.

9.5.2Umsetzung¶

Einsetzen der Differenzenquotienten in die PDE und umstellen nach $u_{i,j,k+1}$ liefert

für $i=1,\ldots,n_x-1$ und $j=1,\ldots,n_y-1$ (IDVID 930).

An den Rändern des Gebiets stehen keine zentralen Differenzen zur Verfügung. Ist die Strömung in das Gebiet hinein gerichtet, so ist der Funktionswert außerhalb des Randes Null. Bei Abfluss aus dem Gebiet, ist der Funktionswert allerdings unklar. Hier können wir nur auf einseitige Differenzen ausweichen.

9.5.3Implementierung¶

def cont2d_ftcs(a, b, T, v, w, f, nx, ny, nt):

    delta_x = a / nx
    delta_y = b / ny
    delta_t = T / nt
    
    x = np.linspace(0, a, nx + 1)
    y = np.linspace(0, b, ny + 1)
    grid_x, grid_y = np.meshgrid(x, y, indexing='ij')  # x ist erste Dimension

    v = v(grid_x, grid_y)
    w = w(grid_x, grid_y)

    u = np.empty((nt + 1, nx + 1, ny + 1), dtype=float)

    # Anfangsbedingungen
    u[0, :, :] = f(grid_x, grid_y)

    # Zeitschritte
    v_pos_left_mask = v[0, :] >= 0
    v_neg_right_mask = v[-1, :] < 0
    w_pos_bottom_mask = w[:, 0] >= 0
    w_neg_top_mask = w[:, -1,] < 0
    
    for k in range(0, nt):

        # innere Punkte
        u[k + 1, :, :] = u[k, :, :]
        u[k + 1, 1:-1, :] -= 0.5 * v[1:-1, :] * delta_t / delta_x * (u[k, 2:, :] - u[k, :-2, :])
        u[k + 1, :, 1:-1] -= 0.5 * w[:, 1:-1] * delta_t / delta_y * (u[k, :, 2:] - u[k, :, :-2])

        # linker Rand
        u_ext = 2 * u[k, 0, :] - u[k, 1, :]
        u_ext[v_pos_left_mask] = 0
        u[k + 1, 0, :] -= 0.5 * v[0, :] * delta_t / delta_x * (u[k, 1, :] - u_ext)

        # rechter Rand
        u_ext = 2 * u[k, -2, :] - u[k, -1, :]
        u_ext[v_neg_right_mask] = 0
        u[k + 1, -1, :] -= 0.5 * v[-1, :] * delta_t / delta_x * (u_ext - u[k, -2, :])

        # unterer Rand
        u_ext = 2 * u[k, :, 0] - u[k, :, 1]
        u_ext[w_pos_bottom_mask] = 0
        u[k + 1, :, 0] -= 0.5 * w[:, 0] * delta_t / delta_y * (u[k, :, 1] - u_ext)

        # oberer Rand
        u_ext = 2 * u[k, :, -2] - u[k, :, -1]
        u_ext[w_neg_top_mask] = 0
        u[k + 1, :, -1] -= 0.5 * w[:, -1] * delta_t / delta_y * (u_ext - u[k, :, -2])
        
    return u

a = 3
b = 2
T = 2

nx = 40
ny = int(np.ceil(nx / a * b))  # \delta y = \delta x
grid_x, grid_y = np.meshgrid(
    np.linspace(0, a, nx + 1),
    np.linspace(0, b, ny + 1),
    indexing='ij'
)

v = lambda x, y: np.ones_like(x)
w = lambda x, y: np.ones_like(x)

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.3) ** 2 + (y/b - 0.3) ** 2))

nt = 200  # Wert viel höher als bei vorherigen Simulationen!

u = cont2d_ftcs(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Wir sehen, dass trotz recht feiner Zeitdiskretisierung schon nach wenigen Schritten Oszillationen entstehen, die sich immer weiter verstärken.

9.6Lax-Friedrichs-Verfahren¶

Mit etwas Aufwand kann man zeigen, dass das FTCS-Verfahren gut konditioniert wird, wenn man den Wert $u_{i,j,k}$ in der Formel zur Berechnung von $u_{i,j,k+1}$ durch den Mittelwert über die vier Nachbarpunkte ersetzt, also

mit

Gilt

mit $C_x$ und $C_y$ aus (7), so ist die Berechung gut konditioniert.

Das diskretisierte Problem kann gleichzeitig als Diskretisierung einer PDE mit Diffusionsterm interpretiert werden, sodass die anfängliche Materialdichte im Laufe der Zeit wieder breitlaufen wird.

Der Funktionswert der außerhalb des Gebietes liegenden Nachbarn von Randpunkten muss je nach Strömungsrichtung gewählt werden. Bei Strömung in das Gebiet hinein ist der Funktionswert Null. Bei Strömung aus dem Gebiet heraus ist er unbekannt. Eine Näherung können wir durch lineare Extrapolation ermitteln, vgl. Nebenbemerkung beim FTCS-Verfahren.

def cont2d_lax_friedrichs(a, b, T, v, w, f, nx, ny, nt):

    delta_x = a / nx
    delta_y = b / ny
    delta_t = T / nt
    
    x = np.linspace(0, a, nx + 1)
    y = np.linspace(0, b, ny + 1)
    grid_x, grid_y = np.meshgrid(x, y, indexing='ij')  # x ist erste Dimension

    v = v(grid_x, grid_y)
    w = w(grid_x, grid_y)

    Cx = np.abs(v).max() * delta_t / delta_x
    Cy = np.abs(w).max() * delta_t / delta_y
    print(f'Courant-Zahlen: {Cx}, {Cy}')
    
    u = np.empty((nt + 1, nx + 1, ny + 1), dtype=float)

    # Anfangsbedingungen
    u[0, :, :] = f(grid_x, grid_y)

    # Zeitschritte
    v_pos_left_mask = v[0, :] >= 0
    v_neg_right_mask = v[-1, :] < 0
    w_pos_bottom_mask = w[:, 0] >= 0
    w_neg_top_mask = w[:, -1,] < 0
    
    for k in range(0, nt):

        # innere Punkte
        u[k + 1, :, :] = 0
        u[k + 1, 1:, :] += 0.25 * u[k, :-1, :]  # linker Nachbar
        u[k + 1, :-1, :] += 0.25 * u[k, 1:, :]  # rechter Nachbar
        u[k + 1, :, 1:] += 0.25 * u[k, :, :-1]  # unterer Nachbar
        u[k + 1, :, :-1] += 0.25 * u[k, :, 1:]  # oberer Nachbar
        u[k + 1, 1:-1, :] -= 0.5 * v[1:-1, :] * delta_t / delta_x * (u[k, 2:, :] - u[k, :-2, :])
        u[k + 1, :, 1:-1] -= 0.5 * w[:, 1:-1] * delta_t / delta_y * (u[k, :, 2:] - u[k, :, :-2])

        # linker Rand
        u_ext = 2 * u[k, 0, :] - u[k, 1, :]
        u_ext[v_pos_left_mask] = 0
        u[k + 1, 0, :] += 0.25 * u_ext
        u[k + 1, 0, :] -= 0.5 * v[0, :] * delta_t / delta_x * (u[k, 1, :] - u_ext)

        # rechter Rand
        u_ext = 2 * u[k, -2, :] - u[k, -1, :]
        u_ext[v_neg_right_mask] = 0
        u[k + 1, -1, :] += 0.25 * u_ext
        u[k + 1, -1, :] -= 0.5 * v[-1, :] * delta_t / delta_x * (u_ext - u[k, -2, :])

        # unterer Rand
        u_ext = 2 * u[k, :, 0] - u[k, :, 1]
        u_ext[w_pos_bottom_mask] = 0
        u[k + 1, :, 0] += 0.25 * u_ext
        u[k + 1, :, 0] -= 0.5 * w[:, 0] * delta_t / delta_y * (u[k, :, 1] - u_ext)

        # oberer Rand
        u_ext = 2 * u[k, :, -2] - u[k, :, -1]
        u_ext[w_neg_top_mask] = 0
        u[k + 1, :, -1] += 0.25 * u_ext
        u[k + 1, :, -1] -= 0.5 * w[:, -1] * delta_t / delta_y * (u_ext - u[k, :, -2])
        
    return u

a = 3
b = 2
T = 4

nx = 40
ny = int(np.ceil(nx / a * b))  # \delta y = \delta x
grid_x, grid_y = np.meshgrid(
    np.linspace(0, a, nx + 1),
    np.linspace(0, b, ny + 1),
    indexing='ij'
)

v = lambda x, y: np.ones_like(x)
w = lambda x, y: np.ones_like(x)

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.3) ** 2 + (y/b - 0.3) ** 2))

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(2 * T * np.maximum(v_max * nx / a, w_max * ny / b)))  # Courant <= 0.5

u = cont2d_lax_friedrichs(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Courant-Zahlen: 0.49382716049382713, 0.5

Wie beim Upwind-Verfahren läuft die Anfangsverteilung auch beim Lax-Friedrichs-Verfahren quer zur Strömungsrichtung stärker breit als in Strömungsrichtung.s

9.7Semi-Lagrange-Verfahren¶

Gut funktionierende Finite-Differenzen-Ansätze für die Kontinuitätsgleichung sind schwierig zu finden. Das Crank-Nicolson-Verfahren liefert gute Ergebnisse, erfordert aber deutlich mehr Rechenaufwand, da die Lösung für den nächsten Zeitschritt nur implizit gegeben ist, also als Lösung eines linearen Gleichungssystems bestimmt werden muss.

Ein anderer, nicht auf finiten Differenzen beruhender Ansatz, ist das sogenannte Semi-Langrange-Verfahren. Dieses beruht wie das nachfolgend betrachtete Partikelverfahren auf der Idee, die Bewegung einzelner Partikel in der Strömung zu verfolgen.

9.7.1Idee¶

Um die Stoffkonzentration $u_{i,j,k+1}$ an einem Punkt $(x_i,y_j)$ zur Zeit $t_{k+1}$ zu ermitteln, berechnen wir, woher ein Partikel zur Zeit $t_k$ gewesen sein muss, damit es zur Zeit $t_{k+1}$ durch die Strömung zum betrachteten Punkt getragen wird. Der gesuchte Wert $u_{i,j,k+1}$ wird dann etwa so groß sein wie $u$ am Ausgangspunkt des Partikels zur Zeit $t_k$.

Führen wir diese “Rückwärtsrechnung” mit jedem Gitterpunkt durch, so erhalten wieder Zeitschritt für Zeitschritt die Stoffkonzentration im gesamten Gebiet.

9.7.2Umsetzung¶

Berechnen den Herkunftsort des in $(x_i,y_j)$ ankommenden Partikels näherungsweise als

und setzen

So einfach die Idee klingt, gibt es doch einige Probleme zu lösen:

• Die Bestimmung des Herkunftsortes eines Partikels ist nur einigermaßen genau, wenn das Strömungsfeld entlang der zurückgelegten Strecke nahezu konstant ist. Je nach Strömungsfeld, müssen die Zeitschritte also klein genug sein.

• Das Auswerten von $u(\cdot,\cdot,t_k)$ ist nur an den Gitterpunkten möglich. Die Punkte $(\tilde{x}_i,\tilde{y}_j)$ werden aber nicht mit den Gitterpunkten zusammenfallen. Wir müssen also interpolieren.

• Liegt ein Punkt $(\tilde{x}_i,\tilde{y}_j)$ außerhalb des betrachteten Gebiets, so ist klar, dass die Strömung bei $(x_i,y_j)$ in das Gebiet hinein gerichtet ist. Wir können in diesem Fall also $u_{i,j,k+1}:=0$ setzen entsprechend den Randbedingungen.

Wird stückweise lineare Interpolation genutz, so kann man zeigen, dass das Verfahren im Wesentlichen äquivalent zum Upwind-Verfahren ist. Um bessere Ergebnisse (weniger Diffusion) zu erzielen, kommen also nur glattere Interpolationen in Frage.

9.7.3Implementierung¶

Die Interpolation überlassen wir scipy.interpolate.RegularGridInterpolator. Das vereinfacht die Implementierung erheblich.

from scipy.interpolate import RegularGridInterpolator

def cont2d_semi_lagrange(a, b, T, v, w, f, nx, ny, nt):

    delta_x = a / nx
    delta_y = b / ny
    delta_t = T / nt
    
    x = np.linspace(0, a, nx + 1)
    y = np.linspace(0, b, ny + 1)
    grid_x, grid_y = np.meshgrid(x, y, indexing='ij')  # x ist erste Dimension

    v = v(grid_x, grid_y)
    w = w(grid_x, grid_y)

    u = np.empty((nt + 1, nx + 1, ny + 1), dtype=float)

    # Anfangsbedingungen
    u[0, :, :] = f(grid_x, grid_y)

    
    # Zeitschritte
    grid_x_tilde = grid_x - delta_t * v
    grid_y_tilde = grid_y - delta_t * w
    
    for k in range(0, nt):

        interpolator = RegularGridInterpolator(
            (x, y), u[k, :, :],
            method='cubic',
            bounds_error=False, fill_value=0  # nicht extrapolieren, sondern 0
        )
        u[k + 1, :, :] = interpolator((grid_x_tilde, grid_y_tilde))
        
    return u

a = 3
b = 2
T = 2

nx = 40
ny = int(np.ceil(nx / a * b))  # \delta y = \delta x
grid_x, grid_y = np.meshgrid(
    np.linspace(0, a, nx + 1),
    np.linspace(0, b, ny + 1),
    indexing='ij'
)

v = lambda x, y: np.ones_like(x)
w = lambda x, y: np.ones_like(x)

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.3) ** 2 + (y/b - 0.3) ** 2))

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(2 * T * np.maximum(v_max * nx / a, w_max * ny / b)))  # Courant <= 0.5

u = cont2d_semi_lagrange(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Die Diffusion ist nun deutlich geringer, aber nicht ganz verschwunden.

T = 2 * np.pi

def r_phi(x, y):
    x = x - 0.5 * a
    y = y - 0.5 * b
    r = np.sqrt(x ** 2 + y ** 2)
    phi = np.arctan2(y, x)
    return r, phi
def v(x, y):
    r, phi = r_phi(x, y)
    return -r * np.sin(phi)
def w(x, y):
    r, phi = r_phi(x, y)
    return r * np.cos(phi)

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(T * (v_max * nx / a + w_max * ny / b)))  # Courant <= 1

u = cont2d_semi_lagrange(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

9.8Partikelverfahren¶

Beim Semi-Lagrange-Verfahren wird in jedem Zeitschritt ein neuer “Satz” Partikel betrachtet und über einen Zeitschritt rückwärts verfolgt. Man kann die Idee auch umgekehrt verwenden: Man startet zur Zeit $t_0$ mit einer zufälligen oder regelmäßigen Partikelverteilung und verfolgt die Bewegungen dieser festen Partikelmenge Zeitschritt für Zeitschritt. Die Werte $u_{i,j,k+1}$ erhält man dann aus Interpolation der Werte $u(\tilde{x}_{i,0},\tilde{y}_{j,0},t_0)$, wobei $(\tilde{x}_{i,0},\tilde{y}_{j,0})$ die anfänglichen Partikelpositionen sind.

Auch hier sind wieder einige Detailprobleme zu lösen:

• Die Interpolation ist aufwendiger als beim Semi-Lagrange-Verfahren, da die Stützstellen kein regelmäßiges Gitter bilden.

• Es können Teilgebiete entstehen, die keinerlei Partikel enthalten. Je nach Situation kann dies bedeuten, dass dort die Stoffkonzentration Null ist oder dass die Partikeldichte sehr gering ist. Ein klare Abgrenzung zwischen “keine Partikel” und “wenige Partikel” ist schwierig.

• In der nähe eines Gitterpunktes können sich viele Partikel aufhalten, was allerdings keinen nennenswerten Einfluss auf die Interpolation hat. Dadurch gehen Partikel “verloren” und hohe Stoffkonzentrationen werden nicht erkannt.

Die letzten beiden Problempunkte kann man lösen indem man das Gebiet in rechteckige Zellen teilt; eine Zelle pro Gitterpunkt mit dem Gitterpunkt als Mittelpunkt. Die Stoffdichte in jeder Zelle ergibt sich dann als Summe aller von den Partikeln in der Zelle getragenen Stoffkonzentrationen.

9.9Strömung um Hindernisse¶

Für komplexere Gebiete $B\subseteq\bbR^2$, die beispielsweise auch Hindernisse enthalten können, muss das Strömungsfeld separat simuliert werden. Für inkompressible Strömungen kann das Geschwindigkeitsfeld stets als Gradient eines sogenannten Potentials ausgedrückt werden:

mit einer Funktion $p:B\to\bbR$. Das Potential ist Lösung der Potentialgleichung

Mittels Neumann-Randbedingungen wird festgelegt über welche Teile des Randes Zu- bzw. Abfluss stattfinden (positive bzw. negative Ableitung in Normalenrichtung) und welche Teile undurchlässig sein (Ableitung in Normalenrichtung ist Null).

Beachtet die Anfangsverteilung die Form des Gebietes $B$ (Dicht ist Null innerhalb von Hindernissen), so kann das Gebiet $B$ für das Lösen der Kontinuitätsgleichung vereinfacht werden, indem Hindernisse nicht mit modelliert werden. Allerdings kann dieser Ansatz zu unerwartetem verhalten der Dichte am Rand von Hindernissen führen. Insbesondere beim Semi-Lagrange-Verfahren kann Masse verloren gehen. Partikelverfahren liefern hier bessere Ergebnisse.

Als Beispiel betrachten wir die Strömung um eine Kreisscheibe mit Radius $R$ und Mittelpunkt $(x_\mathrm{M},y_\mathrm{M})$. Für dieses Fall ist die Lösung der Potentialgleichung und damit das Strömungsfeld bekannt (siehe Potential flow around a circular cylinder für eine Herleitung):

für $\varphi\in[0,2\,\pi)$ und $r\geq R$, wobei

Innerhalb des Kreises ($r<R$) setzen wir das Strömungsfeld auf Null.

a = 3
b = 2
T = 3

nx = 40
ny = int(np.ceil(nx / a * b))  # \delta y = \delta x
grid_x, grid_y = np.meshgrid(
    np.linspace(0, a, nx + 1),
    np.linspace(0, b, ny + 1),
    indexing='ij'
)

def flow_field(x, y, x0, y0):
    R = 0.3
    x = x - x0
    y = y - y0
    r = np.sqrt(x ** 2 + y ** 2)
    phi = np.arctan2(y, x)
    r_mask = r > R
    s = np.zeros_like(x)
    psi = np.zeros_like(x)
    s[r_mask] = (1 - R ** 2 / r[r_mask] ** 2) * np.cos(phi[r_mask])
    psi[r_mask] = -(1 + R ** 2 / r[r_mask] ** 2) * np.sin(phi[r_mask])
    return s * np.cos(phi) - psi * np.sin(phi), s * np.sin(phi) + psi * np.cos(phi)

def v(x, y):
    return flow_field(x, y, 0.5 * a, 0.5 * b)[0]
def w(x, y):
    return flow_field(x, y, 0.5 * a, 0.5 * b)[1]

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.2) ** 2 + (y/b - 0.3) ** 2))

v_max = np.abs(v(grid_x, grid_y)).max()
w_max = np.abs(w(grid_x, grid_y)).max()
nt = int(np.ceil(2 * T * np.maximum(v_max * nx / a, w_max * ny / b)))  # Courant <= 0.5

u = cont2d_semi_lagrange(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)

Direkt vor und hinter der Kreisscheibe ändert sich die Richtung der Strömung stark, sodass die für das Semi-Langrange-Verfahren notwendige Bedingung einer lokal konstanten Strömung nicht mehr gegeben ist. Entsprechend groß sind die entstehenden Fehler in der Simulation. Dem könnte man mit einer Verkleinerung von $Delta t$ entgegenwirken, was aber wiederum die Diffusionseffekte verstärkt.

Am folgenden Beispiel sind die Fehler noch deutlicher. Die Artefakte in der linken Bildhälfte entstehen durch die Interpolation, welche von SciPy nur näherungsweise vorgenommen wird (obwohl alle umgebenden Punkte den Funktionswert 0 haben, sind die interpolierten Werte im Bereich 10^-7).

f = lambda x, y: np.maximum(0, 1 - 10 * np.sqrt((x/a - 0.2) ** 2 + (y/b - 0.5) ** 2))

u = cont2d_semi_lagrange(a, b, T, v, w, f, nx, ny, nt)
show_anim(u, a, b, T, v, w, quiver_step=2)