Übungsblatt 3
Finite-Elemente in 2D
Anmerkung: Die Abgabe von Arbeitsblatt 1 bis 4 ist verpflichtend und konstituiert die Studienleistung der Veranstaltung Simulationstechniken. Die Arbeitsblätter führen von der mathematischen Formulierung eines Modellproblems hin zur numerischen Lösung dieses Problems und bauen aufeinander auf. Zum Bestehen der Veranstaltung müssen auf jedem Blatt mindestens 50% der erzielbaren Punkte erreicht werden.
Geben Sie bei allen Aufgaben die Lösungswege und Zwischenergebnisse mit an. Das Endergebnis alleine ist nicht ausreichend! Wir empfehlen Ihnen die Nutzung von Python. Sollten Sie ein Jupyter-Notebook verwenden, dann können Sie dieses einfach direkt als Lösung bei uns einreichen. In allen anderen Fällen erzeugen Sie bitte ein PDF und legen die numerischen Codes als separate Datei dazu.
Sie werden durch die einzelnen Schritte der Modellimplementierung geleitet, und wir geben Hinweise zur Implementierung. Es ist nicht zwingend notwendig, diese 1-zu-1 zu verfolgen. Im Rahmen dieser Hinweise finden Sie Codeabschnitte, die Sie verwenden können. Sie dürfen natürlich auch die Codebeispiele aus dem Vorlesungsmaterial hier verwenden.
3.1 Entdimensionalisierung der linearisierten Poisson-Boltzmann-Gleichung und der Poisson-Nernst-Planck-Gleichungen
6 Punkte
Kontext:
- Die Entdimensionalisierung rückt alle Werte näher an die Größenordnung 1 heran. Dies kann bei der Konvergenz der Lösung helfen und vereinfacht die Wahl der Konvergenzkriterien.
- Durch Entdimensionalisierung identifizieren wir alle echten Parameter der Gleichung, also Parameter die nicht durch eine einfache Koordinatentransformation eingeführt werden können.
Winzige oder riesige Zahlen sind in der Regel für den Computer kein großes Problem. Über die interne Darstellung als Fließkommazahlen in der Form \(1.564912\,10^4\) mit Mantisse, Basis und Exponent kann der Computer mit Werte ganz verschiedener Größendordnungen umgehen. (Mehr Informationen dazu findet man z.B. auf Wikipedia.)
numpy arbeitet standardmäßig mit 64-bit Fließkommazahlen. Damit lässt sich etwa eine Präzision von \(10^{-15}\) erreichen. Vergleichen Sie dazu die Ausgabe von
2>>> np.finfo(float)
3finfo(resolution=1e-15, min=-1.7976931348623157e+308, max=1.7976931348623157e+308, dtype=float64)
und
2finfo(resolution=1e-15, min=-1.7976931348623157e+308, max=1.7976931348623157e+308, dtype=float64)
Problematisch wird es, wenn mit verschiedenen Variablen gearbeitet wird, die sich um mehr als 15 Größenordnungen unterscheiden. Wenn man solche Werte addiert, macht man also numerische Fehler. Vergleichen Sie dazu die Ausgabe von
mit
Bei unseren Differentialgleichungen, in denen physikalische Größen verschiedener Natur auftreten, können die Unterschiede in der Größenordnung gerade in der elektrochemischen Doppelschicht mit stark erhöhten oder verdünnten Konzentrationen schnell an diese Grenze stoßen, wenn in SI-Einheiten gerechnet wird. Zudem ist es nicht möglich, sinnvolle Konvergenzkriterien für iterative Löser zu definieren, weil unterschiedliche absolute Fehler für Spannungen und Konzentrationen zulässig sind. Entdimensionalisierung wird hier notwendig, um Konvergenz sicherzustellen.
In der Entdimensionalisierung ersetzt man SI-Einheiten mit charakteristischen Maßen des Systems. Im elektrochemischen Kontext ist z.B. die Debye-Länge eine sinnvolle Wahl, oder auch die Länge der Simulationsdomäne. Die sinnvolle Wahl der Einheiten genießt also Spielraum.
Bei geschickter Wahl der Einheiten werden Parameter in den entdimensionalisierten Gleichungen auf Einheitsgrößen reduziert, sodass eine Simulation über geeignete Skalierung auf unterschiedliche physikalische Situationen bezogen werden kann. Die nun effektiv geringere Anzahl an Parameter macht die Simulation selbst übersichtlicher.
Vorgegebenes Einheitensystem
Nutzen Sie das folgende Einheitensystem:
- Konzentrationseinheit: \(c^\infty \), werden wir als identisch für alle Spezies annehmen.
- Ladungseinheit: \(q_e\), die Ladung eines Elektrons (Absolutwert).
- Einheit für das elektrische Potential: \(k_B T / q_e\), die Temperaturspannung (engl. “thermal voltage”).
- Längeneinheit: \(\lambda = \sqrt {\frac {\varepsilon k_B T }{2 q_e ^2 c^\infty }}\) die Debye-Länge (für zwei Ionenspezies).
- Zeiteinheit: \(\lambda ^2 / D\)
Anmerkung: In diesem und im nächsten Übungsblatt werden wir statische Poisson-Boltzmann-Systeme (PB) und dynamische Poisson-Nernst-Planck-Systeme (PNP) betrachten. In allen Systemen gehen wir von symmetrischen Elektrolyten aus, d.h. Anionen und Kationen sind (im geschlossenen System anfänglich, im offenen System im Bad) in gleicher Menge vorhanden. Wir sprechen daher häufig im Singular von der Referenzkonzentration. In PB-Systemen stellt die in der PB-Gleichung auftretende Konzentration \(c^\infty \) eine natürliche Referenzeinheit dar. Diese Konzentration entspricht im offenen System der Konzentration im Bad. Im (geschlossenen dynamischen) PNP-System werden wir die (initiale) mittlere Konzentration ebenfalls \(c^\infty \) nennen und als Referenzeinheit wählen.
3.1.1 Nichtlineare Poisson-Boltzmann-Gleichung
2 Punkte
Die Poisson-Boltzmann-Gleichung lautet \begin {equation} \nabla ^2\Phi = -\frac {1}{\varepsilon }\sum _{\alpha } q_\alpha \, c^\infty \, \exp \left (\frac {-q_\alpha \Phi }{k_B T}\right )\,\text {.} \end {equation} Schreiben Sie diese Gleichung mit Hilfe der entdimensionalisierten Größen \(\tilde \Phi = \Phi q_e / k_B T\), \(\tilde x = x/\lambda \), \(\tilde y = y/\lambda \) und \(\tilde q_\alpha = q_\alpha / q_e\).
3.1.2 Linearisierte Poisson-Boltzmann-Gleichung
1 Punkte
Bringen Sie die für ein binäres, symmetrisches Elektrolyt einfacher Ladung linearisierte Poisson-Boltzmann-Gleichung in entdimensionalisierte Form.
Poisson-Nernst-Planck-Gleichungen
3 Punkte
Bringen Sie die Poisson-Nernst-Planck-Gleichungen, \begin {equation} \frac {\partial c_\alpha }{\partial t} = \nabla \cdot D_\alpha \left ( \nabla c_\alpha + \frac {q_\alpha c_\alpha }{k_B T} \nabla \Phi \right ) \end {equation} und \begin {equation} \nabla ^2\Phi = -\frac {1}{\varepsilon }\sum _{\alpha } q_\alpha \, c_\alpha , \end {equation} in entdimensionalisierte Form. Wählen Sie \(\tilde c_\alpha = c_\alpha / c^\infty \) und nehmen Sie die Diffusinskonstanten \(D_\alpha \) als konstant im Rechengebiet und equivalent für alle Spezies an, d.h. \(D_\alpha = D\) für alle \(\alpha \).
3.2 Stationäre Lösung der linearisierten Poisson-Boltzmann-Gleichung in 2D - Herleitung des diskretisierten Gleichungssystems
12 Punkte
Lösen Sie die linearisierte Poisson-Boltzmann-Gleichung, die wir in Aufgabe 1 entdimensionalisiert haben, auf dem zweidimensionalen Gebiet \(\Omega \). Auf dem Rand \(\partial \Omega _D\) herrschen Dirichlet-Randbedingungen, und auf dem Rand \(\partial \Omega _N\) gelten Neumann-Randbedingungen.
Der Einfachheit halber verzichten wir darauf, \(\tilde \Phi \) für die entdimensionalisierte Größe zu verwenden, und schreiben direkt \(\Phi \). Bei der Interpretation der Ergebnisse werden wir die Dimensionsbehaftete Größe explizit kennzeichnen: \(\Phi ^D\).
Herleitung des diskretisierten Gleichungssystems (12 Punkte)
3.2.1 Schwache Form
1 Punkte
Bringen sie die lineare Poisson-Boltzmann-Gleichung in die schwache Form und reduzieren Sie die Differenzierbarkeitsanforderungen.
3.2.2 Diskrete Formulierung
1 Punkte
Drücken Sie \(\Phi \) mit Hilfe einer Funktionenbasis \(\{\varphi _i\}\) aus und überführen Sie die Integralgleichung mit Hilfe der Galerkinmethode in ein algebraisches Gleichungssystem.
Im Folgenden werden wir nur noch Dirichlet-Randbedingungen und homogene Neumann-Randbedingungen, \(\nabla \Phi (x,y) = 0\) für \((x,y) \in \partial \Omega _N\), betrachten. Dadurch vereinfacht sich das Gleichungssystem leicht.
3.2.3 Finite-Elemente-Diskretisierung
2 Punkte
Wir betrachten ein rechteckiges Gebiet mit einem regulären Gitter an dreieckigen Elementen. Das folgende Bild zeigt zum Beispiel ein Gitter mit 4 Knoten in x- und y-Richtung.
Eingetragen ist die serielle Indexierung der Knoten, Elemente und Kästen. K Die Knoten und Kästen können auch durch ein Paar kartesischer Indizes identifiziert werden, die wir ebenfalls bei 0 anfangen lassen. Zum Verständnis der kartesischen Indizes ist es nützlich, sich zu überlegen beispielsweise welchen kartesischen Indizes Knoten 5 hat und welchen kartesischen Indizes Kasten 5 hat.
Verwenden Sie stückweise lineare Basisfunktionen und beschränken Sie sich auf den Funktionenraum
\begin {equation} \Phi (x, y) = \sum _i \varphi _i(x, y)\,\Phi _i~\text {.} \end {equation} \(i\) ist der Index des Knotens. Die Basisfunktionen haben die Eigenschaft \(\varphi _i(x_j, y_j) = \delta _{ij}~\text {,}\) wobei es sich bei \(\delta _{ij}\) um das Kronecker-Delta mit \(\delta _{ij} = 1 \) für \(i=j\) und \(\delta _{ij} = 0 \) für \(i\neq j\) und bei \(x_j, y_j\) um die kartesischen Koordinaten des Knoten \(j\) handelt.
Innerhalb eines Elementes \(e\), \( (x, y) \in \Omega _e\), lässt sich \(\varphi _i(x, y)\) mithilfe der linearen Formfunktionen aufschreiben:
\begin {equation} \varphi _{g(e, I)}(x, y) = N^{(e)}_{I}(x, y) = a + b x + c y \end {equation}
Dabei bezeichnet \(g(e, I)\) die Abbildung vom lokalen Knotenindex \(I\in \{0,1,2\}\) auf Element \(e \in \{0, ..., N_e-1\}\) zum globalen Knotenindex \(i \in \{0, ..., N_{n}-1\}\). \(N^{(e)}_{I}(x, y)\) ist eine abgekürzte Schreibweise für \(N_{I}(\xi ^{(e)}(x,y), \eta ^{(e)}(x,y))\).
Geben Sie \(\nabla N^{(e)}_I(x,y)\) an.
3.2.4 Berechnung der Systemmatrix
6 Punkte
Mit der Galerkin-Methode erhält man ein Gleichungssystem der Form \(\sum _j K_{ij} \Phi _j = b_i \), wobei \(K_{ij}\) Integrale der Form \(\int \limits _{\Omega } dx dy \nabla \varphi _i \cdot \nabla \varphi _j\) und \(\int \limits _{\Omega } dx dy \varphi _i \varphi _j\) enthält. Die nicht verschwindenden Elemente dieser Systemmatrix lassen sich elementweise mithilfe der Formfunktionen aufschreiben, z.B.
\begin {equation} \begin {aligned} K_{ij} &=\int \limits _{\Omega } dx dy\ \varphi _{i} \varphi _{j} \\ &= \int \limits _{\Omega } dx dy \left ( \sum _e \sum _I N^{(e)}_I \delta _{g(I,e)i}\right ) \varphi _{j} \\ &= \sum _e \sum _I \delta _{g(e,I)i} \delta _{g(e,J)j} \int \limits _{\Omega _e} dx dy N^{(e)}_I N^{(e)}_J \\ &= \sum _e \sum _I \delta _{g(e,I)i} \delta _{g(e,J)j} K^{(e)}_{IJ} \end {aligned} \end {equation}
Die globale Matrix \(\underline {K}\) kann mit Hilfe von \(3\times 3\) Elementmatrizen \(\underline {K}^{(e)}\) aufgeschrieben werden. Aufgrund des regulären Gitters sind alle Elementmatrizen \(\underline {K}^{(e)}\) identisch. \(g(e,I)\) bezeichnet hier den globalen Knoten, welcher zum lokalen Knoten \(I\) in Element \(e\) gehört.
Die Systemmatrix kann in Beiträge von den unterschiedlichen Elementen, und damit von den einzelnen Formfunktionen, zerlegt werden. Im folgenden sollen Sie diese Beiträge berechnen, die die Elementmarix liefern.
Schritt 1
Berechnen Sie \begin {equation} \int \limits _{\Omega ^{(e)}} dx dy \ \nabla N_I \cdot \nabla N_J \end {equation}
Da in den Gradienten nur Abstände zwischen Gitterpunkten auftreten und die Knoten regelmäßig verteilt sind, werden diese Integrale für jedes Element das gleiche Ergebnis liefern.
Im regelmäßigen Gitter gibt es zwei Arten von Dreiecken. Wenn wir die Ecken \(0, 1, 2\) gegen den Uhrzeigersinn (positiver Umlaufsinn) bennenen, wie unten illustriert, dann sind diese Dreiecke aber identisch bis auf eine 180° Rotation. \(0\) bezeichnet die Ecke mit dem rechten Winkel.
2 | \ | \ dy | \ | \ 0 --- 1 dx oder dx 1 ---0 \ | \ | dy \ | \| 2
Wie lautet die Elementmatrix des Laplace-Operators für den Sonderfall \(\Delta x=\Delta y\)?
\begin {equation} \begin {pmatrix} &&\\ &\int \limits _{\Omega _e} dx dy \ \nabla N^{(e)}_I \cdot \nabla N^{(e)}_J &\\ && \end {pmatrix}_{I=0,1,2; J = 0,1,2} = \frac {1}{2} \begin {pmatrix} 2&-1&-1\\ -1&1&0\\ -1&0&1 \end {pmatrix} \end {equation}
Schritt 2
Implementieren Sie die Elementmatrix für den Laplace-Operator
2 """
3 Construct the element matrix of the partially integrated Laplace operator with
4 homogeneous Neumann boundary conditions.
5
6 The rows and columns of the element matrix correspond to the nodes (corners of
7 the triangle) ordered in geometric positive order, starting with the node with
8 the 90 degree angle.
9
10 Parameters
11 ----------
12 dx : float
13 Grid spacing in x direction
14 dy : float
15 Grid spacing in y direction
16
17 Returns
18 -------
19 element_matrix_ll : numpy.ndarray
20 Array of shape (3,3) containing the element matrix
21 """
22 pass
Schritt 3
Berechnen Sie die Elementmatrix für den Term in \(\Phi \), der bei der Überführung in die schwache Form aus der rechten Seite der linearisierten Poisson-Boltzman Gleichung hervorgegangen ist. Diese Matrix nennen wir Massenmatrix. Schreiben Sie wiederum eine Funktion, die diese Matrix konstruiert:
Die Integrale sind ein aufwendiger zu berechnen, als die des Laplace-Operators.
Berechnen Sie \(\int \limits _{\Omega _e} N^{(e)}_1 N^{(e)}_2 \ dx dy \).
Den Rest der Matrix dürfen Sie aus dem Skript entnehmen.
3.2.5 Zusammenbau der Systemmatrix
1 Punkte
siehe z.B. [Larson, M. G. et al. 10, (Springer Berlin Heidelberg, 2013)] Algo 11
Wir bauen erst die Systemmatrix zusammen, ohne uns Gedanken über die Randbedingungen zu machen. Danach werden die Matrix und die rechte Seite entsprechend der Dirichlet- und Neumann-Randbedingungen (und eventueller Quellterme) justiert.
Sie haben jetzt die einzelnen Integrale der diskretisierten schwachen Formulierung auf Elementebene ausgewertet. Jetzt müssen Sie das gesamte Integral
\begin {equation} \int \limits _{\Omega } \ldots = \sum _{e} \int \limits _{\Omega _e} \ldots \end {equation}
berechnen, also die einzelnen Beiträge der Integrale über die Elemente aufsummieren.
Graphisch kann man sich das folgendermaßen vorstellen:
Die Aufgabe besteht darin, für jedes Element die Einträge der Elementsteifigkeitsmatrix an der entsprechenden Stelle in der Gesamtmatrix aufzuaddieren. Dafür braucht man eine Tabelle oder Funktion, die den Knoten an den Ecken des Elementes einen globalen Index zuweist.
Der Einfachheit halber sollten Sie zuerst die Beziehung zwischen lokalen Elementecken und globalen Knotenindizes für das oben illustrierte 4 x 4 Knoten Gitter von Hand in einer Tabelle aufschreiben:
Damit können Sie im Zweifelsfall auch schon die ersten Testbeispiele weiter unten simulieren. Später sollten die Indizes automatische erstellt werden, um beliebige Gittergrößen zu ermöglichen.
Hier ein Paar Schritte, die bei der Ermittlung der Beziehung zwischen den Indexierungssystemen helfen können:
- Geben Sie die Beziehung zwischen dem seriellen Index \(k\) und dem kartesischen Index-Paar \(i, j\) eines Knotens an.
- Geben Sie die Beziehung zwischen dem seriellen Index \(k\) und den kartesischen Indexen \(i, j\) einer Box an.
Implementieren Sie folgende Funktionen, die zwischen den verschiedenen Indexierungen hin- und herkonvertieren:
2 """
3 Convert the global linear consecutive index of a node to its Cartesian coordinates.
4
5 Parameters
6 ----------
7 g : int
8 global linear index of the node
9 nb_nodes : tuple of ints
10 Number of nodes in the Cartesian directions
11
12 Returns
13 -------
14 i : int
15 x-coordinate (integer) of the node
16 j : int
17 y-coordinate (integer) of the node
18 """
19 pass
20
21
22def node_index(i, j, nb_nodes):
23 """
24 Turn node coordinates (i, j) into their global node index.
25
26 Parameters
27 ----------
28 i : int
29 x-coordinate (integer) of the node
30 j : int
31 y-coordinate (integer) of the node
32 nb_nodes : tuple of ints
33 Number of nodes in the Cartesian directions
34
35 Returns
36 -------
37 g : int
38 Global node index
39 """
40 pass
41
42
43def make_grid(nb_nodes):
44 """
45 Make an array that contains all elements of the grid. The
46 elements are described by the global node indices of
47 their corners. The order of the corners is in order of
48 the local node index.
49
50 They are sorted in geometric positive order and the first
51 is the node with the right angle corner at the bottom
52
53 . Elements within the same box are consecutive.
54
55 This is the first element per box:
56
57 2
58 | \
59 | \
60 dy | \
61 | \
62 0 --- 1
63
64 dx
65
66 This is the second element per box:
67
68 dx
69 1 ---0
70 \ |
71 \ | dy
72 \ |
73 \|
74 2
75
76 Parameters
77 ----------
78 nb_nodes : tuple of ints
79 Number of nodes in the Cartesian directions
80
81 Returns
82 -------
83 triangles_el : numpy.ndarray
84 Array containing the global node indices of the
85 element corners. The first index (suffix _e)
86 identifies the element number and the second index
87 (suffix _l) the local node index of that element.
88 """
89 pass
Für kleine Gitter können Sie folgende Funktion verwenden, um Ihre Indexierung nach Korrektheit zu überprüfen:
2
3
4def show_mesh(nb_nodes, physical_sizes):
5 """
6 Show a simple representation of the grid, showing the order of
7 the linear indices for nodes and elements.
8
9 Parameters
10 ----------
11 nb_nodes : tuple of ints
12 Number of nodes along x and y direction
13 physical_sizes : tuple of floats
14 Length of the simulation domain along x and y direction
15 """
16
17 def nb_elements(nb_nodes):
18 Nx, Ny = nb_nodes
19 return (Nx - 1) * (Ny - 1) * 2
20
21 sx, sy = physical_sizes
22 Nx, Ny = nb_nodes
23 dx = sx / (Nx - 1)
24 dy = sy / (Ny - 1)
25 assert np.prod(nb_nodes) < 40, "grid to big to be reasonably plotted"
26
27 fig, ax = plt.subplots()
28 triangles_el = make_grid(nb_nodes)
29 for k in range(nb_elements(nb_nodes)):
30 nodes_n = triangles_el[k]
31 mean_x = 0
32 mean_y = 0
33 for n in nodes_n:
34 ix, iy = node_coordinates(n, nb_nodes)
35 mean_x = mean_x + ix
36 mean_y = mean_y + iy
37 ax.text(ix * dx, iy * dy, "n{}".format(n))
38 ax.text(mean_x / 3 * dx, mean_y / 3 * dy, "e{}".format(k), color="red")
39
40 ax.set_xlim(- dx, Nx * dx)
41 ax.set_ylim(- dy, Ny * dy)
42
43 ax.set_xlabel("x")
44 ax.set_ylabel("y")
3.2.6 Funktionen für die Darstellung der Ergebnisse
0 Punkte
Farbige 2D “Höhenlinien” Darstellung
2 show_mesh=False, levels=11):
3 """
4 Show a scalar field on a triangulated structured grid using
5 a colormap representation. The field will be linearly interpolated
6 between the node values.
7
8 Parameters
9 ----------
10 field_g : array_like
11 Values of the field in serial form (of length corresponding
12 to number of nodes)
13 nb_nodes : tuple of ints
14 Number of nodes (Nx, Ny) in the Cartesian directions
15 physical_sizes : tuple of floats
16 Length of the simulation domain along x and y direction
17 show_mesh : bool, optional
18 If True, shows nodes and element boundaries in black
19 (Default: False)
20 levels : int, optional
21 Number of levels in the colormap
22 (Default: 11)
23
24 Returns
25 -------
26 figure : matplotlib.Figure
27 Figure object
28 axes : matplotlib.Axes
29 Axes object
30 colorbar : matplotlib.Colorbar
31 Colorbar object
32 """
33 import matplotlib.tri as mtri
34
35 Nx, Ny = nb_nodes
36 sx, sy = physical_sizes
37
38 dx = sx / (Nx - 1)
39 dy = sy / (Ny - 1)
40
41 fig, ax = plt.subplots()
42
43 ax.set_aspect(1)
44
45 x = np.array([node_coordinates(i, nb_nodes)[0] for i in range(np.prod(nb_nodes))])
46 y = np.array([node_coordinates(i, nb_nodes)[1] for i in range(np.prod(nb_nodes))])
47
48 x = x * dx
49 y = y * dy
50
51 triangles_el = make_grid(nb_nodes)
52 triang = mtri.Triangulation(x, y, triangles_el)
53
54 cb = plt.colorbar(ax.tricontourf(triang, field_g, levels=levels))
55 if show_mesh:
56 ax.triplot(triang, ’ko-’)
57
58 ax.set_xlabel(r"x ($\lambda$)")
59 ax.set_ylabel(r"y ($\lambda$)")
60
61 return fig, ax, cb
Zur Interpolation zwischen den Knotenpunkten, z.B. für die Darstellung von Schnitten, können Sie folgende Funktion verwenden
2 """
3 Performs a linear finite element interpolation of the field in x, y
4
5 Parameters:
6 -----------
7 x_i, y_i : array_like
8 Positions where to interpolate the field
9 field_g : array_like
10 Values of the field in serial form (of length corresponding
11 to number of nodes)
12 nb_nodes : tuple of ints
13 Number of nodes (Nx, Ny) in the Cartesian directions
14 physical_sizes : tuple of floats
15 Length of the simulation domain along x and y direction
16
17 Returns:
18 --------
19 interpolated_field_i : numpy.ndarray
20 Values interpolated at position x_i, y_i
21
22 """
23 import matplotlib.tri as mtri
24
25 Nx, Ny = nb_nodes
26 sx, sy = physical_sizes
27
28 dx = sx / (Nx - 1)
29 dy = sy / (Ny - 1)
30
31 _x = np.array([node_coordinates(i, nb_nodes)[0]
32 for i in range(np.prod(nb_nodes))])
33 _y = np.array([node_coordinates(i, nb_nodes)[1]
34 for i in range(np.prod(nb_nodes))])
35
36 _x = _x * dx
37 _y = _y * dy
38
39 triangles_el = make_grid(nb_nodes)
40 triang = mtri.Triangulation(_x, _y, triangles_el)
41
42 interp = mtri.LinearTriInterpolator(triang, field_g)
43
44 return interp(x_i, y_i)
3.2.7 Konstituierendes Gleichungssystem aufstellen
1 Punkte
Sie sollen die einzelnen Schritte beim Erstellen des Gleichungssystems in separaten Funktionen implementieren. Danach sollte sich eine Simulation in einem recht knappen Skript aufsetzen lassen:
2nb_nodes = (4, 4)
3physical_sizes = (3., 3.)
4
5# Construct system matrix ignoring boundary conditions
6# (We call this here the ’bulk’ matrix)
7system_matrix_gg = poisson_boltzmann_bulk_matrix(nb_nodes, physical_sizes)
8
9# Add source term to the right hand side (may not be necessary)
10rhs_g = source(nb_nodes, physical_sizes,
11 source_parameter_1=..., source_parameter_2=...)
12
13# Add boundary condition to matrix and right hand side
14my_interesting_bc(
15 system_matrix_gg, rhs_g, nb_nodes, physical_sizes,
16 bc_parameter_1=..., bc_parameter_2=...)
17
18# Solve linear system
19...
20
21# Visualize solution
22...
Implementieren Sie poisson_boltzmann_bulk_matrix.
3.3 Stationäre Lösung der linearisierten Poisson-Boltzmann-Gleichung in 2D - Modellvalidation gegen das 1D Problem
10 Punkte
Um unser Modell aus der vorherigen Aufgabe zu validieren, betrachten wir wieder den Plattenkondensator aus Übungsblatt 1. Aufgrund der Struktur des Gitters ist es einfacher, wenn die Kondensatorplatten oben (\(y=d\)) und unten (\(y=0\)) liegen statt links und rechts. Bitte setzen Sie auf der unteren Platte ein (entdimensionalisiertes) Potential von \(-1\) an und auf der oberen Platte ein Potential von \(1\) an. Es gibt keinen Quellterm.
3.3.1 Analytische Referenzlösung
0 Punkte
Die analytische Lösung zu \(\frac {\partial ^2\Phi }{\partial x^2} - \Phi = 0\) mit den Randbedingungen \(\Phi (0) = \Phi _0\) und \(\Phi (L) = \Phi _1\) und Abstand \(L\) zwischen den Kondensatorplatten wurde in Übungsblatt 1 aufgestellt. Implementieren Sie diese als
2 """
3 Parameters
4 ----------
5 y : array_like
6 Positions where to evaluate the analytical solution.
7 In unis of the debye length
8 debye^2 = kT*permittivity / (q**2 * concentration_bath)
9 potential_bottom : float
10 Value of the electric potential on the plate at y = 0
11 In units of the thermal_voltage = kb*T / q
12 potental_top : float
13 Value of the electric potential on the plate at y = plate_separation
14 plate_separation : float
15 distance between the bottom and top plate in units of the debye length
16
17 Returns
18 -------
19 potential : np.ndarray
20 Analytical solution of the 1D linearized Poisson-Boltzman equation
21 """
22 pass
3.3.2 Randbedingungen in 2D
2 Punkte
Mit welcher Randbedingungen “links” und “rechts” lässt sich der unendliche Plattenkondensator auch mit der 2D-Diskretisierung modellieren?
Implementieren Sie diese Randbedingungen und vergleichen Sie Ihr Modell mit der analytischen Lösung. Variieren Sie die Größe des Gitters um die Konvergenz der Lösung mit Auflösung der Simulation. Implementieren Sie dafür folgende Funktion, die die Systemmatrix und die rechte Seite (right hand side, rhs) entsprechend modifiziert:
2 """
3 Modifies the global system matrix and the right hand side according to following boundary conditions:
4 bottom and top faces: Dirichlet
5 left and right faces: Neumann homogeneous
6
7 Parameters
8 ----------
9 system_matrix_gg : ndarray of shape (Nx * Ny, Nx * Ny)
10 System matrix with natural boundary conditions
11 rhs_g: ndarray of shape (Nx * Ny)
12 Right hand side with natural boundary conditions
13 nb_nodes : tuple of ints
14 Number of nodes in the Cartesian directions (Nx, Ny)
15 potential_bottom : float
16 Value of the electric potential on the bottom electrode (y = 0)
17 potential_top : float
18 Value of the electric potential on the top electrode (y = plate_separation )
19
20 """
21
22 pass
Anmerkung: Die Dirichlet-Randbedingungen lässt sich am einfachsten implementieren, indem Sie die entsprechenden Randknoten auf den gewünschten Wert zwingen. Das können Sie erreichen, indem Sie in Ihrer Funktion alle Elemente der Systemmatrix \(K_{ii} = 1\) und \(K_{ij} = 0 ~\forall i \neq j\) und den entsprechenden Eintrag im Quellvektor (rhs) \(b_i = \Phi _{0}\) auf dem linken Rand und \(\Phi _{1}\) auf dem rechten Rand setzen, falls \(i\) einen Knoten auf einem dieser Dirichlet-Ränder beschreibt.
3.3.3 Die erste Simulation
4 Punkte
Stellen Sie die Lösung mit der Funktion plot_colormap für \(4\times 4\) Knoten und einer physikalischen Systemgröße von \(5\lambda \times 5 \lambda \) dar.
3.3.4 Gitterkonvergenz und Vergleich mit der analytischen Lösung
4 Punkte
Vergleichen Sie numerische Lösung und analytische Lösung für \(2\times 4\), \(4 \times 4\), \(2 \times 8\), \(2 \times 16\), \(2 \times 32\), \(16 \times 16\) Knoten. Benutzen Sie dabei eine Systemgröße von \(5 \lambda \times 5\lambda \). Stellen Sie die Lösung an den Knotenpunkten bei \(x = 0\) dar.
3.4 Elektrolytkondensator in 2D
8 Punkte
Jetzt widmen wir uns einem Problem, das wir analytisch nicht lösen können. Analog zur Vorlesung schauen wir uns einen Plattenkondensator an, mit dem Unterschied, dass sich jetzt ein Elektrolyt, also bewegliche Ladungen, zwischen den Platten befindet. An der Oberfläche der Platten befinden sich zwei gegenüberliegende Elektroden mit Potiential \(\Phi _u\) und \(\Phi _o\). Die Platten sind unendlich breit und haben den Abstand \(d\). Die Elektrodenhaben eine endliche Länge \(L\).
Links und rechts der Elektroden nehmen wir an, dass das elektrische Feld keine Komponente Senkrecht zur Oberfläche hat. Weit von den Elektroden klingt das elektrische Feld zu Null ab, was wir durch die Bedingung \(\partial \Phi / \partial x = 0\) Auf dem rechten Rand der Simulationsdomäne approximieren. Wir werden später bestimmen, wie groß die “Puffer-Zone” dafür sein muss.
Das System besitzt eine vertikale Symmetrieebene die durch die Mitte beider Elektroden geht. Dies erlaubt uns, nur die Hälfte des Systems zu simulieren, wobei die Symmetriebedingung auch wieder durch eine homogene Neumann Randbedingung erzwungen wird.
Verwenden Sie folgende Parameter
- \(c_{\infty }^D = 100 ~\mathrm {mol}\, \mathrm {m}^{-3}\), wobei \(\mathrm {mol} = N_A\) einfach als Avogadro-Zahl betrachtet wird.
- \(q_e^D = 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}\)
- Elektrodenbreite: \(L^D=10~\text {nm}\)
- Plattenabstand: \(d^D=5~\text {nm}\)
- Potential unten: \(\Phi _u^D = -0.1~\text {V}\)
- Potential oben: \(\Phi _o^D = 0.1~\text {V}\)
3.4.1 Randbedingungen implementieren
4 Punkte
Lösen Sie das obige Problem und stellen Sie die Lösung für \(3\times 3\) Knoten und einer Simulationsdomänengröße von \(10 ~\text {nm}\times 5 ~\text {nm}\). Vergessen Sie nicht, dass nur die Hälfte der Elektroden simuliert werden, i.e. die Simulationsdomänengröße entspricht dem roten Kasten in der Abbildung. Benutzen Sie plot_colormap mit levels=np.linspace(-0.1, 0.1., 11) für die Visualisierung. Stellen Sie die Ergebnisse in SI Einheiten dar.
3.4.2 Gitterkonvergenz
4 Punkte
Stellen Sie den Verlauf des Potentials als Querschnitt bei \(y^D = 0.1~\text {nm}\) dar. Zeigen Sie, wie die Lösung mit
- zunehmender Gitterfeinheit einerseits und
- zunehmendem Abstand des Neumann-Randes zur Elektrode (“Pufferzone”) andererseits
konvergiert.
3.5 Berechnung der Kapazität
10 Punkte
In der Vorlesung wurde der Zusammenhang zwischen der elektrischen Ladung auf der Elektrode und dem elektrischen Feld an dieser Elektrode erklärt. Mithilfe der Ergebnisse aus der vorherigen Aufgabe wollen wir nun untersuchen, wie die Kapizät des Kondensators vom Elektrolyten beeinflusst wird. Wir arbeiten hier wieder mit den dimensionslosen Größen. Die Bedeutung der Ergebnisse erläutern wir auch anhand der dimensionsbehafteten Größen, die mit einem Exponent \(^D\) gekennzeichnet werden. Verwenden Sie die gleichen Parameter wie in der vorhergenden Aufgabe.
3.5.1 Kapazität und Aspektverhältnis
4 Punkte
Stellen Sie die entdimensionalisierte Kapazität \(C = C^D / \varepsilon ^D t^D\) in Abhängigkeit des Aspektverhältnisses \(d / L\) für verschiedene Elektrodenlängen \(L\) auf logarithmischer Skala dar. \(t^D L^D\) ist die Fläche des Kondensators.
Verwenden Sie folgende Parameterkombinationen: - \(\Delta x = \Delta y = 0.0125\), \(N_p=32\), - \(\Delta x = \Delta y = 0.05\), \(N_p=128\), - \(\Delta x = \Delta y = 0.2\), \(N_p=128\),
wobei \(N_p\) die Anzahl an Knoten in der Elektrode ist. Eine Puffer-Zone von einer Länge von etwa \(\text {min}(1, d)\) sollte ausreichend sein, da 1 der Debye Länge entspricht.
Variieren Sie den Abstand \(d\) (\(N_y=8,16,32,64\)) für jede dieser Parameterkombinationen.
Um den Diskretisierungsfehler abzuschätzen können wir noch gröber diskretisierte Simulationen einzeichnen, z.B.:
- \(\Delta x = \Delta y = 0.025\), \(N_p=16\), \(N_y=32\)
- \(\Delta x = \Delta y = 0.1\), \(N_p=256\), \(N_y=16\)
In der Vorlesung wurde gezeigt, dass bei kleinen Aspektverhältnissen \(d/L\), \(L/d\) der Kapazität des Plattenkondensators ohne Elektrolyt entspricht. Zeichnen sie \(C = L/d\) als Anhaltspunkt auch ein.
Bei fixem Aspektverhältnis bedeutet ein größeres \(L\) dass die gesamte Geometrie bei konstanter Debye-Länge vergrößert wird. Das kann aber auch so interpretiert werden, dass die Geometrie fixiert bleibt und sich die Debye-Länge verkürzert, z.B. durch eine höhere Ionenkonzentration.
3.5.2 Analytische Berechnung der Kapazität
2 Punkte
In Übungsblatt 1 haben sie die Poisson-Boltzmann Gleichung im eindimensionalen Fall analytisch gelöst. Berechnen Sie die Kapazität des Kondensators mit Hilfe des analytischen Modells aus. Stellen Sie diese in Zusammen mit Ihren Simulationsergebnissen in 3.5.1 dar.
3.5.3 Praktische Anwendung
0 Punkte
Wie lässt sich dieser Kondensator auch als Sensor verwenden ?