Fallstudie: Stationäre Wärmeleitung - Wissenschaftliches Rechnen (i788)

Lösen ein stationäres Wärmeleitproblem in der Ebene mittels FEM. Dabei bedeutet “stationär”, dass wir nicht die zeitliche Entwicklung simulieren, sondern den nach hinreichend langer Wartezeit vom betrachteten physikalischen System eingenommenen Zustand. Sind Wärmezufluss und -abfluss in einem Körper über einen langen Zeitraum konstant, so wird sich eine stabile Temperaturverteilung im Körper einstellen. Dieser stabile Zustand soll simuliert werden.

12.1Modell¶

Sei $B\subseteq\bbR^2$ die Querschnittsfläche eines Kühlkörpers (z.B. einer CPU). Am unteren Rand

\Gamma_2:=\{(x,0):\;x\in[-1,1]\}

(1)

erfolgt ein zeitlich konstanter Wärmeeintrag

g_2(x,0):=100\,(1-x^2).

(2)

Am restlichen Rand $\Gamma_1$ liegt die Umgebungstemperatur stets bei

g_1(x,y)=20

(3)

(Luftstrom).

import matplotlib.pyplot as plt
import numpy as np

# Liste von Eckpunkten zur Beschreibung des polygonalen Gebiets
outline = [[-1, 0]]
for x in np.linspace(-1 + 0.1, 1 + 0.1, 5, endpoint=False):
    outline.extend([[x, 0.3], [x, 1], [x + 0.2, 1], [x + 0.2, 0.3]])
outline.append([1, 0])

# Randbedingungen
g1 = lambda x, y: 20
g2 = lambda x, y: 100 * (1 - x ** 2)

# Plot der Gebiets
fig, ax = plt.subplots(figsize=(10, 5))
x, y = zip(*(outline + [outline[0]]))
ax.plot(x, y, '-k')
ax.axis('equal')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()

Im Inneren des Kühlkörpers gilt

\Delta u=0

(4)

für die Temperatur $u:B\to\bbR$ .

12.2Schwache Formulierung¶

Die schwache Formulierung erhalten wir als Spezialfall der schwachen Formulierung der Poisson-Gleichung:

Finde $u_0\in V$ , sodass

\int_B \nabla u_0\circ\nabla v\diff b=\int_{\Gamma_2}g_2\,v\diff s-\int_B \nabla u_1\circ\nabla v\diff b

(5)

für alle $v\in V$ gilt. Dabei ist

V:=\{v\in H^1(B):\;v(x)=0\text{ für }x\in\Gamma_1\}

(6)

und $u_1\in H^1(B)$ ist eine Funktion, die auf $\Gamma_1$ gerade $g_1$ entspricht. Die gesuchte Lösung ist dann

u:=u_0+u_1.

(7)

12.3Gebietszerlegung¶

Für die Zerlegung des Gebiets $B$ in Dreieckselemente nutzen wir das Python-Modul pygmsh.

import pygmsh

# Zerlegung in Dreieckselemente
with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(outline, mesh_size=0.05)
    mesh = geom.generate_mesh()

# Knotenkoordinaten als NumPy-Array
#   axis 0: globale Knotennummer (0, 1,...)
#   axis 1: Koordinate (0 (x), 1 (y))
nodes = mesh.points[:, :2]

# globale Knotennummern der Randkanten als NumPy-Array
#   axis 0: Randkantennummer (0, 1,...)
#   axis 1: lokale Knotennummer (0, 1)
boundary_lines = mesh.cells[0].data

# globale Knotennummern der Dreieckselemente als NumPy-Array
# (Elementzusammenhangstabelle)
#   axis 0: Elementnummer (0, 1,...)
#   axis 1: lokale Knotennummer (0, 1, 2)
triangles = mesh.cells[1].data

del mesh  # nicht mehr benötigt

Zur späteren Verwendung trennen wir die von pygmsh generierten Knoten nach inneren Knoten, Knoten auf dem Dirichlet-Rand und Knoten auf dem Neumann-Rand.

# Knoten nach Lage auswählen (Rand, innen)
b_nodes = set.union(*[{int(n1), int(n2)} for n1, n2 in boundary_lines])
i_nodes = list(set(range(0, nodes.shape[0])) - b_nodes)

# Randknoten in Dirichlet und Neumann trennen
d_nodes = []
n_nodes = []
for n in b_nodes:
    if -1 < nodes[n, 0] < 1 and nodes[n, 1] == 0:
        n_nodes.append(n)
    else:
        d_nodes.append(n)

del b_nodes, boundary_lines  # nicht mehr benötigt

Die Gebietszerlegung sollte immer visuell geprüft werden. Insbesondere sollten die Dreiecke möglichst nah an der Form eines gleichseitigen Dreiecks sein.

# Plot der Gebietszerlegung
fig, ax = plt.subplots(figsize=(10, 5))

# Dreiecke
for n1, n2, n3 in triangles:
    n1231 = [n1, n2, n3, n1]
    ax.plot(nodes[n1231, 0], nodes[n1231, 1], '-k', lw=1)

# Knoten
ax.plot(nodes[i_nodes, 0], nodes[i_nodes, 1], 'og', label='innen')
ax.plot(nodes[d_nodes, 0], nodes[d_nodes, 1], 'or', label='Dirichlet')
ax.plot(nodes[n_nodes, 0], nodes[n_nodes, 1], 'ob', label='Neumann')

# Knotennummern
#for n, (x, y) in enumerate(nodes):
#    ax.text(x + 0.02, y + 0.02, str(n))

ax.axis('equal')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend()

plt.show()

Für den Aufbau der Systemmatrix weiter unten sind die Nachbarschaftsbeziehungen zwischen Knoten sowie zwischen Knoten und Dreieckselementen von Bedeutung.

# alle Knoten außer Dirichlet-Rand
ni_nodes = n_nodes + i_nodes

# Dreiecke für jeden Knoten
#   Listenindex ist globale Knotennummer
#   Listeneinträge sind Mengen von Elementnummern
tris_of_nodes = [set() for _ in nodes]
for tri, (n1, n2, n3) in enumerate(triangles):
    tris_of_nodes[n1].add(tri)
    tris_of_nodes[n2].add(tri)
    tris_of_nodes[n3].add(tri)

# benachbarte Knoten (kleine Knotennummer zuerst)
neighbors = set()
for node, tris_of_node in enumerate(tris_of_nodes):
    for tri in tris_of_node:
        for n in set(triangles[tri, :]) - {node}:
            neighbors.add((node, n) if node < n else (n, node))

12.4Systemmatrix und rechte Seite¶

Für die numerische Berechnung der in der schwachen Formulierung auftretenden Bereichsintegrale verwenden wir stückweise lineare Funktionen auf der Dreieckszerlegung. Für jeden Knoten wird also eine Basisfunktionen $v_k:B\to\bbR$ eingeführt, die am entsprechenden Knoten den Wert 1 annimmt und an allen anderen Knoten den Wert 0 hat.

12.4.1Dirichlet-Randbedingungen¶

Im kontinuierlichen Modell wird die Lösung $u$ zerlegt in die Summe $u=u_0+u_1$ , wobei $u_0$ auf dem Dirichlet-Rand $\Gamma_1$ Null ist und $u_1$ auf dem Dirichlet-Rand die Randbedingung erfüllt. Die Funktion $u_1$ ist abseits von $\Gamma_1$ frei wählbar. Sei $u_1$ also Null im Inneren des Gebiets und auf dem Neumann-Rand.

Im diskretisierten Zustand hat $u_0$ nur Nicht-Null-Koeffizienten for Basisfunktionen $v_k$ , die zu einem inneren Knoten oder zu einem Knoten auf dem Neumann-Rand gehören. Die Systemmatrix wird als weniger Zeilen bzw. Spalten haben als es Knoten (Basisfunktionen) gibt. Umgekehrt genügen für die diskretisierte Darstellung der Funktion $u_1$ die Basisfunktionen die zu einem Knoten auf dem Dirichlet-Rand gehören.

# Abbildung der Knotennummern (ohne Dirichlet-Rand) auf Matrixzeilen/-spalten
node2idx = {n: idx for idx, n in enumerate(ni_nodes)}

Die Systemmatrix wird dünn besetzt sein, sodass wir die Einträge zunächst in Listen sammeln und dann später daraus mit SciPy eine dünn besetzte Matrix generieren.

# leere (dünn besetzte) Matrix, rechte Seite zunächst überall 0
Adata = []
Acols = []
Arows = []
b = np.zeros(len(ni_nodes), dtype=float)

12.4.2Gradient der Basisfunktionen¶

Sei $P_k$ der Knoten zur globale Knotennummer $k$ . Der Gradient der zugehörigen Basisfunktion $v_k$ ist Null auf allen Dreieckselemente, die nicht an den Knoten $P_k$ grenzen. Auf an den Knoten $P_k$ angrenzenden Dreiecken ist der Gradient konstant.

Die Berechnung des Gradients von $v_k$ auf einem Dreieck mit Eckknoten $P_k,P_l,P_m$ kann leicht durch Transformation des Dreiecks auf das Referenzdreieck mit den Eckpunkten $(0,0)$ , $(1,0)$ , $(0,1)$ erfolgen:

Die Transformation der Ecken des Referenzdreiecks auf die Eckknoten $P_k,P_l,P_m$ ist durch

\Phi(\xi,\eta):=P_k+\xi_1\,(P_l-P_k)+\xi_2\,(P_m-P_k)

(8)

gegeben. Die zum betrachteten Dreieck und zur Basisfunktion $v_k$ passende Formfunktion ist $\varphi_1$ (vgl. Section 11.4.3). Mit ihr gilt

(v_k\circ\Phi)(\xi,\eta)=\varphi_1(\xi,\eta).

(9)

Andererseits liefert die Kettenregel

\nabla(v_k\circ\Phi)(\xi,\eta)=J_\Phi(\xi,\eta)^\rmT\,(\nabla v_k)(\Phi(\xi,\eta))

(10)

(IDVID 1210), wobei

J_\Phi(\xi,\eta)={\begin{bmatrix} [P_l-P_k]_1& [P_m-P_k]_1\\ [P_l-P_k]_2& [P_m-P_k]_2\end{bmatrix}}

(11)

die Jacobi-Matrix zu $\Phi$ ist. Somit

{\begin{align} (\nabla v_k)(\Phi(\xi,\eta)) &=J_\Phi(\xi,\eta)^{-\rmT}\,\nabla(v_k\circ\Phi)(\xi,\eta)\\ &=J_\Phi(\xi,\eta)^{-\rmT}\,\nabla\varphi_1(\xi,\eta). \end{align}}

(12)

Analog folgt

(\nabla v_l)(\Phi(\xi,\eta))=J_\Phi(\xi,\eta)^{-\rmT}\,\nabla\varphi_2(\xi,\eta).

(13)

12.4.3Bereichsintegrale¶

Für die Berechnung der Einträge der Systemmatrix und der rechten Seite benötigen wir die Bereichsintegrale

\int_B(\nabla v_k)^\rmT\,\nabla v_l\diff b

(14)

für alle Paare von Basisfunktionen und auch für $k=l$ .

Sind entsprechenden Knoten $P_k$ und $P_l$ nicht benachbart, so ist das Integral Null. Bei benachbarten Knoten ist der Integrand nur auf den beiden Dreiecken von Null verschieden, die sowohl $P_k$ als auch $P_l$ als Eckknoten haben. Bezeichnen wir diese beien Dreiecke mit $T$ und $\tilde{T}$ , so gilt

\int_B(\nabla v_k)^\rmT\,\nabla v_l\diff b=\int_T(\nabla v_k)^\rmT\,\nabla v_l\diff b+\int_{\tilde{T}}(\nabla v_k)^\rmT\,\nabla v_l\diff b.

(15)

Auf $T$ und $\tilde T$ ist der Integrand jeweils konstant.

Sei $P_m$ der dritte Eckknoten im Dreieck $T$ und sei $\Phi$ die Transformation des Referenzdreiecks auf das Dreieck $T$ (vgl. oben). Dann gilt

\int_T(\nabla v_k)^\rmT\,\nabla v_l\diff b =\int_0^1\int_0^{1-\xi}h(\xi,\eta)\,|\det J_\Phi(\xi,\eta)|\diff\eta\diff\xi

(16)

mit

{\begin{align*} h(\xi,\eta) &=\bigl(\nabla v_k(\Phi(\xi,\eta))\bigr)^\rmT\,\nabla v_l(\Phi(\xi,\eta))\\ &=\bigl(J_\Phi(\xi,\eta)^{-\rmT}\,\nabla\varphi_1(\xi,\eta)\bigr)^\rmT\,J_\Phi(\xi,\eta)^{-\rmT}\,\nabla\varphi_2(\xi,\eta)\\ &=\bigl(\nabla\varphi_1(\xi,\eta)\bigr)^\rmT\,J_\Phi(\xi,\eta)^{-1}\,J_\Phi(\xi,\eta)^{-\rmT}\,\nabla\varphi_2(\xi,\eta)\\ &=\bigl(\nabla\varphi_1(\xi,\eta)\bigr)^\rmT\,\bigl(J_\Phi(\xi,\eta)^\rmT\,J_\Phi(\xi,\eta)\bigr)^{-1}\,\nabla\varphi_2(\xi,\eta) \end{align*}}

(17)

und

\det J_\Phi(\xi,\eta)=\sqrt{\det\bigl(J_\Phi(\xi,\eta)^\rmT\,J_\Phi(\xi,\eta)\bigr)}.

(18)

Die symmetrische und positiv definite Matrix $J_\Phi(\xi,\eta)^\rmT\,J_\Phi(\xi,\eta)$ nimmt für alle $\xi$ und $\eta$ den gleichen Wert

J_\Phi(\xi,\eta)^\rmT\,J_\Phi(\xi,\eta)={\begin{bmatrix}\alpha&\beta\\\beta&\gamma\end{bmatrix}}

(19)

mit Zahlen $\alpha,\beta,\gamma,\delta\in\bbR$ an. Auch die Gradienten der Formfunktionen hängen nicht von $\xi$ und $\eta$ ab:

\nabla\varphi_1(\xi,\eta)={\begin{bmatrix}-1\\-1\end{bmatrix}},\qquad \nabla\varphi_2(\xi,\eta)={\begin{bmatrix}1\\0\end{bmatrix}}.

(20)

Der Integrand ist also konstant und der Integrationsbereich (Referenzdreieck) hat den Flächeninhalt $\frac{1}{2}$ . Somit

{\begin{align} \int_T(\nabla v_k)^\rmT\,\nabla v_l\diff b &=\frac{1}{2}\,\frac{1}{\alpha\,\gamma-\beta^2}\,{\begin{bmatrix}-1&-1\end{bmatrix}}\,{\begin{bmatrix}\gamma&-\beta\\-\beta&\alpha\end{bmatrix}}\,{\begin{bmatrix}1\\0\end{bmatrix}}\,\sqrt{\alpha\,\gamma-\beta^2}\\ &=\frac{\beta-\gamma}{2\,\sqrt{\alpha\,\gamma-\beta^2}}. \end{align}}

(21)

Für $k=l$ erhält man völlig analog (zweimal $\varphi_1$ statt $\varphi1$ mit $\varphi_2$ )

\int_T(\nabla v_k)^\rmT\,\nabla v_k\diff b=\frac{\alpha+\gamma-2\,\beta}{2\,\sqrt{\alpha\,\gamma-\beta^2}}.

(22)

Ist $P_k$ oder $P_l$ ein Knoten auf dem Dirichlet-Rand, so spielt der Wert $a(v_k,v_l)$ keine Rolle in der Systemmatrix, taucht aber bei der Berechnung des Bereichsintegrals $a(u_1,v_k)$ bzw. $a(u_1,v_l)$ in der rechten Seite $b$ auf.

# Matrix A ohne Diagonale, Dirichlet-Anteil von b
for n1, n2 in neighbors:

    # a(v_n1, v_n2) berechnen
    int_gradgrad = 0
    for tri in tris_of_nodes[n1] & tris_of_nodes[n2]:  # höchstens 2 Läufe
        n3 = next(iter(set(triangles[tri, :]) - {n1, n2}))  # dritter Knoten des Dreiecks
        J = np.vstack((nodes[n2, :] - nodes[n1, :], nodes[n3, :] - nodes[n1, :])).T
        JJ = np.matmul(J.T, J)
        detJJ = JJ[0, 0] * JJ[1, 1] - JJ[0, 1] ** 2
        int_gradgrad += (JJ[0, 1] - JJ[1, 1]) / (2 * np.sqrt(detJJ)) 

    # in Matrix A oder in rechte Seite schreiben?
    n1d = n1 in d_nodes
    n2d = n2 in d_nodes
    if not n1d and not n2d:
        i1 = node2idx[n1]
        i2 = node2idx[n2]
        Acols.append(i1)
        Arows.append(i2)
        Adata.append(int_gradgrad)
        Acols.append(i2)
        Arows.append(i1)
        Adata.append(int_gradgrad)
    elif n1d and not n2d:
        i2 = node2idx[n2]
        b[i2] += -g1(nodes[n1, 0], nodes[n1, 1]) * int_gradgrad
    elif n2d and not n1d:
        i1 = node2idx[n1]
        b[i1] += -g1(nodes[n2, 0], nodes[n2, 1]) * int_gradgrad
    else:  # Dirichlet-Randkante
        pass

# Diagonale von A
for n1 in ni_nodes:
    i1 = node2idx[n1]

    # a(v_n1, v_n1) berechnen
    int_gradgrad = 0
    for tri in tris_of_nodes[n1]:
        n2, n3 = set(triangles[tri, :]) - {n1}
        J = np.vstack((nodes[n2, :] - nodes[n1, :], nodes[n3, :] - nodes[n1, :])).T
        JJ = np.matmul(J.T, J)
        detJJ = JJ[0, 0] * JJ[1, 1] - JJ[0, 1] ** 2
        int_gradgrad += (JJ[0, 0] + JJ[1, 1] - 2 * JJ[0, 1]) / (2 * np.sqrt(detJJ)) 

    # in Matrix schreiben
    Acols.append(i1)
    Arows.append(i1)
    Adata.append(int_gradgrad)

12.4.4Neumann-Randbedingungen¶

In die rechte Seite $b$ müssen noch die Integrale

\int_{\Gamma_2}g_2\,v_k\diff s

(23)

eingebracht werden, wobei $k$ die globalen Knotennummern der Knoten auf dem Neumann-Rand durchläuft.

Das Integrationsgebiet $\Gamma_2$ kann in einzelne Randkanten zerlegt werden, wobei für jede Kante die Endknoten entweder beide auf dem Neumann-Rand liegen oder ein Endknoten auf dem Neumann- und ein Endknoten auf dem Dirichlet-Rand liegt.

In das Integral über eine Randkante gehen nur zwei Basisfunktionen (zwei Neumann-Knoten) oder nur eine Basisfunktion (Dirichlet-Knoten und Neumann-Knoten) ein.

Seien $P_k$ und $P_l$ die Endknoten einer Randkante $K$ . Dann ist durch

\Psi(\zeta):=P_k+\zeta\,(P_l-P_k)

(24)

eine Abbildung des Intervalls $[0,1]$ auf die Randkante gegeben. Damit erhalten wir

\int_K g_2\,v_k\diff s=\int_0^1 g_2(\Psi(\zeta))\,(1-\zeta)\diff\zeta

(25)

und

\int_K g_2\,v_l\diff s=\int_0^1 g_2(\Psi(\zeta))\,\zeta\diff\zeta,

(26)

wobei wir

v_k(\Psi(\theta))=1-\zeta\qquad\text{und}\qquad v_l(\Psi(\theta))=\zeta

(27)

verwendet haben.

Diese Integrale können leicht mittels numerischer Integration berechnet werden.

from scipy.integrate import quad

# Neumann-Anteil von b
for n1, n2 in neighbors:

    # mindestens ein Neumann-Randknoten
    n1n = n1 in n_nodes
    n2n = n2 in n_nodes
    n1d = n1 in d_nodes
    n2d = n2 in d_nodes
    if not((n1n and n2n) or (n1n and n2d) or (n1d and n2n)):
        continue  # keine Neumann-Randkante

    # Integral von g_2 * v_n1 und von g_2 * v_n2
    P1 = nodes[n1, :]
    P2 = nodes[n2, :]
    distP1P2 = np.sqrt((P1[0] - P2[0]) ** 2 + (P1[1] - P2[1]) ** 2)
    if n1n:  # Basisfunktion bei D-Knoten trägt nicht zum Integral bei
        i1 = node2idx[n1]
        func = lambda t: g2(*(P1 - t * (P2 - P1))) * (1 - t)
        b[i1] += distP1P2 * quad(func, 0, 1)[0]
    if n2n:
        i2 = node2idx[n2]
        func = lambda t: g2(*(P1 - t * (P2 - P1))) * t
        b[i2] += distP1P2 * quad(func, 0, 1)[0]

12.5Lösen des Gleichungssystems¶

Das Gleichungssystem kann mit Standardverfahren für dünn besetzte Matrizen gelöst werden.

import scipy.sparse as ss

# Matrix A als SciPy-Objekt
A = ss.csr_array(
    (Adata, (Arows, Acols)), 
    shape=(len(ni_nodes), len(ni_nodes))
)

# Gleichungssystem lösen
u0 = ss.linalg.spsolve(A, b)

12.6Zusammensetzen der Lösung¶

Das Array u0 enthält nun die Funktionswerte der Lösungsfunktion an den inneren Knoten und an den Randknotes des Neumann-Randes. Auf dem Dirichlet-Rand ist die Lösung bekannt.

# Lösung zusammensetzen
u = np.zeros(len(nodes), dtype=float)
u[d_nodes] = g1(nodes[d_nodes, 0], nodes[d_nodes, 1])
u[ni_nodes] = u0

Für die Darstellung der Lösung kann der Bereich zwischen den Knoten linear interpoliert werden. Dies entspricht der kontinuierlichen Darstellung der Lösung als Summe stückweise linearer Ansatzfunktionen.

# Plot der Lösung
fig, ax = plt.subplots(figsize=(10, 5))

# Lösung
contour_plot = ax.tricontourf(
    nodes[:, 0], nodes[:, 1], triangles, u, cmap='jet', levels=100
)
fig.colorbar(contour_plot, ax=ax, orientation='vertical')

# Dreiecke
for n1, n2, n3 in triangles:
    n1231 = [n1, n2, n3, n1]
    ax.plot(nodes[n1231, 0], nodes[n1231, 1], '-w', lw=0.5)

ax.axis('equal')
ax.set_xlabel('x')
ax.set_ylabel('y')

plt.show()

12.7Realitätsnähe¶

Die Modellierung mit Dirichlet-Randbedingungen ist wenig realitätnah. Um die Oberflächentemperatur konstant bei 20°C zu halten, müsste das Bauteil beispielsweise von einer 20°C warmen Luftströmung umströmt werden, die beliebig viel Wärme aufnehmen kann. Realistischer sind Robin-Randbedingungen, die die Wärmeabgabe proportional zur Differenz zwischen Oberflächen- und Umgebungstemperatur modellieren (vgl. Diffusionsgleichung).

Bei der schwachen Formulierung mit Robin-Randbedingungen erhält man sowohl in der Bilinearform $a$ als auch in der rechten seite $b$ zusätzlich ein Kurvenintegral über den Robinrand.

12 Fallstudie: Stationäre Wärmeleitung