Übungsblatt 2
Grundlagen der Methode der finiten Elemente

2.1 Fourier-Basisfunktionen

2.1.1 Orthogonalität der 1D-Fourier-Basisfunktionen

3 Punkte
Auf Übungsblatt 1 haben Sie eine analytische Lösung der Diffusionsgleichung mithilfe der eindimensionalen Fourier-Reihe bestimmt. Weisen Sie nach, dass die Basisfunktionen dieser Fourier-Reihe eine orthogonale Basis bilden. Das Skalarprodukt zweier \(L\)-periodischen Funktion sei dabei definiert als: \begin {equation} \left ( f, g \right ) = \frac {1}{L} \int _0^L \dif x \, \left ( f^*(x) g(x) \right ) \end {equation} wobei \(f^*(x)\) für das komplex-konjugierte von \(f(x)\) steht.

2.1.2 Orthogonalität der 2D-Fourier-Basisfunktionen

3 Punkte
Zeigen Sie, dass auch die Basisfunktionen der zweidimensionalen Fourier-Reihe eine orthogonale Basis bilden. Das Skalarprodukt zweier Funktionen, die in x-Richtung \(L_x\)-periodisch und in y-Richtung \(L_y\)-periodisch sind, sei dabei gegeben als \begin {equation} \left (f, g \right ) = \frac {1}{L_x L_y} \int _0^{L_y} \int _0^{L_x} \dif x \dif y \, \left ( f^*(x,y) g(x,y) \right ) \end {equation}

2.1.3 Fourier-Koeffizienten

2 Punkte
Leiten Sie die Formel für die Koeffizienten der zweidimensionalen Fourierreihe her, indem Sie die Funktion \(f(x,y)\) auf die Fourierreihen-Basis projizieren.

2.2 Finite-Elemente

In dieser Aufgabe betrachten wir die Grundlagen der Finite Elemente Methode. Als Beispielproblem nehmen wir die stationäre Diffusionsgleichung in 1D (Vgl. Übungsblatt 1, Aufgabe 4): \begin {equation} D \frac {\partial ^2 c}{\partial x^2} = f(x) \label {eq:Diffusion} \end {equation} Wie in Übungsblatt 1 ist dabei \(D\) die Diffusionskonstante, \(c(x)\) eine Stoffmengenkonzentration und \(f(x)\) ein Quellterm. Die Gleichung gilt auf einem endlichen Gebiet \([0; L]\).

2.2.1 Diskretisierung für allgemeine Quellterme

5 Punkte
Diskretisieren Sie Gleichung \eqref{eq:Diffusion} mithilfe von Finiten Elementen, indem Sie zuerst die schwache Form herleiten und die Differenzierbarkeitsanforderung verringern, dann den Galerkin-Ansatz wählen und zuletzt die gewählten Basisfunktionen einsetzen. Die Gitterpunkte am Rand des Gebietes können Sie dabei zunächst vernachlässigen.

Als Diskretisierung wollen wir \(N\) gleichmäßig verteilte Gitterpunkte verwenden. Als Basisfunktionen wählen wir die stückweise linearen Zeltfunktionen, wie sie in Abb. 2.1 dargestellt sind.

2.2.2 Diskretisierung der Quellterme

2 Punkte
Berechnen Sie nun die diskretisierten Terme der Quellterme aus der vorherigen Aufgabe für die Quellterme:

1. \(f(x) = -\delta \left (x - \frac {L}{3} \right ) / \mathrm {m}^3 \mathrm {s} + \delta \left (x - \frac {2L}{3} \right ) \mathrm {m}^{-3} \mathrm {s}^{-1}\)
2. \(f(x) = f_0\)

2.3 Finite-Element im periodischen Raum

In dieser Aufgabe wollen wir die Diffusionsgleichung mit periodischen Randbedingungen auf dem Gebiet \([0;L]\), die in Übungsblatt 1, Aufgabe 4 analytisch gelöst wurde mit Finiten Elementen lösen.

Dabei wählen wir die gleiche Diskretisierung und die gleichen Basisfunktionen wie in Aufgabe 2.3. Genau wie bei der analytischen Lösung ist der Quellterm \(f(x) = -\delta \left (x - \frac {L}{3} \right ) \mathrm {m}^{-3} \mathrm {s}^{-1} + \delta \left (x - \frac {2L}{3} \right ) \mathrm {m}^{-3} \mathrm {s}^{-1}\).

2.3.1 Mittelwertbedingung

1 Punkte
Wie wir auf dem Übungsblatt 1 gesehen haben, ist der Mittelwert der Lösungsfunktion ohne weitere Angaben nicht bestimmbar. Anders gesagt, wenn \(c_1(x)\) eine periodische Lösung der Diffusionsgleichung ist, dann ist \(c_2(x) = c_1(x) + \text {konst.}\) ebenfalls eine periodische Lösung der Diffusionsgleichung. Diese Unendlichkeit an Lösungen bedeutet, dass das Gleichungssystem, das wir mit den Finiten Elementen aufstellen, nicht eindeutig lösbar ist.

Damit wir ein lösbares Gleichungssystem bekommen, brauchen wir ein Problem mit einer eindeutigen Lösung. Wir müssen deshalb die Forderung ‘periodische Randbedingungen’ mit einer weiteren Bedingung ergänzen. Wir wählen dafür einen vorgegebenen Mittelwert von \(c_0\), d.h. \begin {equation} \frac {1}{L}\int _0^L \dif x \, c(x) = c_0 \end {equation}

Dadurch ergibt sich eine weitere Gleichung für die Koeffizienten \(a_i\). Stellen Sie diese Gleichung auf.

2.3.2 Systemmatrix

1 Punkte
Das Gleichungssystem aus diskretisierten Gleichungen für die Koeffizienten \(a_i\) wird üblicherweise in Matrix-Form geschrieben: \begin {equation} \underline {K} \cdot \overrightarrow {a} = \overrightarrow {f} \end {equation} \(\underline {K}\) nennt man dann die Systemmatrix.

Stellen Sie die Systemmatrix für die Lösung der 1D-Diffusionsgleichung mit periodischen Randbedingungen und einem vorgegebenen Mittelwert mithilfe von lineare Finiten Elementen für 4 Gitterpunkte auf.

2.3.3 Numerische Lösung

8 Punkte
Schreiben Sie eine python-Funktion, um die Diffusionsgleichung im periodischen Raum mit linearen Finiten Elementen zu lösen. Diese Funktion sollte als Argumente die Anzahl der Gitterpunkte \(N\), den Abstand zwischen zwei Gitterpunkten \(dx\) und einen Vektor mit der bereits diskretisierten rechten Seite des Problems erwarten. Verwenden Sie also bitte die Schnittstelle:

Nutzen Sie Ihre Funktion, um einen Plot, zu erstellen, auf dem die FEM-Lösung mit der analytischen Lösung verglichen wird. Auf dem Plot muss zu sehen sein:

• Analytische Lösung aus Übungsblatt 1, Aufgabe 4: \(c(x) = -\frac {1}{D}\max (0\mathrm {m}, x-L/3) \mathrm {m}^{-2} \mathrm {s}^{-1} + \frac {1}{D}\max (0\mathrm {m}, x-2L/3) \mathrm {m}^{-2} \mathrm {s}^{-1} + \frac {1}{3D} \mathrm {m}^{-2} \mathrm {s}^{-1} x + c_0\)
• Lösung mit FEM für N=4
• Lösung mit FEM für N=7

Kommentieren Sie kurz, was man auf dem Plot beobachten kann. Was würden Sie erwarten, wenn der Dirac-Impuls nicht direkt auf einem Gitterpunkt liegt? (Eine Berechnung ist nicht nötig, ein kurzer Kommentar genügt.)

Verwenden Sie für den Plot folgende Werte für die Parameter:

• \(D = 0.8 \mathrm {m}^{2} \mathrm {s}^{-1}\)
• \(L = 1.5\mathrm {m}\)
• \(1 / L \int _{0}^{L} \text {d}x\, c(x) = 0.3 \mathrm {m}^{-3}\)

2.4 Finite-Elemente mit Dirichlet- und Neumann-Randbedingungen

In dieser Aufgabe betrachten wir die Diffusionsgleichung mit einer Dirichlet-Randbedingung bei \(x=0\mathrm {m}\) und einer Neumann-Randbedingung bei \(x=L\): \begin {align} c(x=0) &= c_0 \quad \text {(Dirichlet)}\\ \left .\frac {\partial c}{\partial x}\right \vert _{x=L} &= c'_L \quad \text {(Neumann)} \end {align}

Physikalisch entspricht dies einem System, dass auf der einen Seite an ein unendliches Reservoir, in dem die Konzentration immer gleich bleibt, angeschlossen ist, während auf der anderen Seite ein konstanter Teilchenstrom in bzw. aus dem System fließt (vgl. Abb. 2.2).

Zur Lösung mit Finiten Elementen wählen wir die gleiche Diskretisierung und die gleichen Basisfunktionen wie in Aufgabe 2.3.

2.4.1 Randbedingungen

2 Punkte
Stellen Sie die Gleichungen für die Koeffizienten \(a_i\) auf, die die Dirichlet und die Neumann Randbedingungen beschreiben. Setzen Sie dabei für den Quellterm:
1. \(f(x) = -\delta \left (x - \frac {L}{3} \right ) \mathrm {m}^{-2} \mathrm {s}^{-1} + \delta \left (x - \frac {2L}{3} \right ) \mathrm {m}^{-2} \mathrm {s}^{-1}\)
2. \(f(x) = f_0\)

2.4.2 Numerische Lösung

10 Punkte
Schreiben Sie eine python-Funktion, um die Diffusionsgleichung mit einer Dirichlet- und einer Neumann-Randbedingung mit linearen Finiten Elementen zu lösen. Diese Funktion sollte als Argumente die Anzahl der Gitterpunkte \(N\), den Abstand zwischen zwei Gitterpunkten \(dx\) und einen Vektor mit der bereits diskretisierten rechten Seite des Problems erwarten. Ein zusätzliches Argument sollte angeben, ob die Systemmatrix zurückgegeben wird oder nicht. Verwenden Sie also bitte die Schittstelle:

Nutzen Sie Ihre Funktion, um 2 Plots zu erstellen, auf denen Folgendes dargestellt wird:

• Lösung der Diffusionsgleichung mit zwei Delta-Distributionen als Quellterm: FEM-Lösung für \(N=4\), \(N=7\) und \(N=10\) darstellen.
• Lösung der Diffusionsgleichung mit konstantem Quellterm: FEM-Lösung für \(N=4\), \(N=7\) und \(N=10\) sowie die analytische Lösung

Kommentieren Sie jeweils kurz, was auf den Plots zu beobachten ist.

Verwenden Sie folgende Werte für die Parameter:

• \(D = 0.8\frac {\mathrm {m}^2}{\mathrm {s}}\)
• \(L = 3\mathrm {m}\)
• \(c_0 = 2\mathrm {m}^{-3}\)
• \(c'_L = 0.7\mathrm {m}^{-4}\)
• \(f_0 = 0.5\mathrm {m}^{-3}\mathrm {s}^{-1}\)

2.5 Finite-Elemente Lösung der linearisierten Poisson-Boltzmann Gleichung

Nachdem wir die Grundlagen der FEM betrachtet haben, wollen wir jetzt zu unserem Modell-Problem zurückkehren und die linearisierte Poisson-Boltzmann-Gleichung für den symmetrischen Elektrolyten in 1D lösen, die in Übungsblatt 1, Aufgabe 7 hergeleitet wurde. In 1D lautet sie \begin {equation} \frac {\partial ^2 \Phi (x)}{\partial x^2} = \frac {2 c^{\infty } q^2}{\varepsilon k_B T} \Phi (x) = \lambda ^{-2} \Phi (x) \, \text {,} \label {eq:A5_Poisson_Boltzmann} \end {equation} wobei \(\Phi (x)\) das elektrostatische Potential, \(\varepsilon \) die Permittivität, \(k_B\) die Boltzmann-Konstante, \(q\) der Ladungsbetrag der betrachteten Ionen, \(T\) die Temperatur, \(c^{\infty }\) die Referenzkonzentration der Ionen und \(\lambda \) die Debye-Länge ist.

Wir betrachten die Gleichung auf einem endlichen Gebiet \([0; L]\). An den Rändern von diesem Gebiet befinden sich zwei inerte Elektroden mit elektrostatischem Potential \(\Phi _0\) und \(\Phi _1\), bei dem System könnte es sich also um einen Plattenkondensator handeln (vgl. Abbildung)

2.5.1 Diskretisierung

6 Punkte
Diskretisieren Sie Gleichung \eqref{eq:A5˙Poisson˙Boltzmann} mithilfe von Finiten Elementen, indem Sie zuerst die schwache Form herleiten und die Differenzierbarkeitsanforderung verringern, dann den Galerkin-Ansatz wählen und zuletzt die gewählten Basisfunktionen einsetzen.

Verwenden Sie die gleichen Basisfunktionen wie in Aufgabe 2 - 4.

2.5.2 Numerische Lösung

7 Punkte
Schreiben Sie eine python-Funktion, um die 1D-Poisson-Boltzmann Gleichung mit zwei Dirichlet-Randbedingungen mit linearen Finiten Elementen zu lösen. Diese Funktion sollte als Argumente die Anzahl der Gitterpunkte \(N\), den Abstand zwischen zwei Gitterpunkten \(l\), die Debye-Länge \(\lambda = \sqrt \frac {\varepsilon k_B T}{2 c^\infty q^2}\) und einen Vektor mit der bereits diskretisierten rechten Seite des Problems nehmen. Verwenden Sie also die Schnittstelle:

Nutzen Sie Ihre Funktion, um einen Plot mit der Lösung der 1D-Poisson-Boltzmann Gleichung zu erstellen. Der Plot sollte enthalten:

• Analytische Lösung aus Übungsblatt 1, Aufgabe 7: \begin {equation} \Phi (x) = K_1 \exp {\left ( \frac {1}{\lambda }x\right )} + K_2\exp {\left ( -\frac {1}{\lambda }x\right )} \end {equation} mit \(K_1 = \left ( \Phi _1 - \Phi _0 \exp {\left ( -\frac {L}{\lambda }\right )}\right )\left ( \exp {\left ( \frac {L}{\lambda }\right )} - \exp {\left ( -\frac {L}{\lambda }\right )} \right )^{-1}\) und \(K_2 = \Phi _0 - K_1\)
• FEM Lösung für N=4
• FEM Lösung für N=8
• FEM Lösung für N=16

Verwenden Sie folgende Werte für die Parameter:

• \(L = 3 \mathrm {nm}\)
• \(\Phi _0 = -0.01 \mathrm {V}\)
• \(\Phi _1 = 0.04 \mathrm {V}\)
• \(c^{\infty } = 10^{3}\mathrm {mol}\, \mathrm {m}^{-3}\)
• \(q = e = 1.602\cdot 10^{-19}\mathrm {A} \mathrm {s}\)
• \(\varepsilon = 80 \cdot 8.85\cdot 10^{-12}\mathrm {A} \, \mathrm {s} \, \mathrm {V}^{-1} \, \mathrm {m}^{-1}\)
• \(k_B = 1.380649\cdot 10^{-23} \mathrm {J} \, \mathrm {K}^{-1}\)
• \(T = 293.15 \mathrm {K}\)
• Avogrado-Zahl \(N_A = 6.022\cdot 10^{23}\)

---