Diskretisierung des diffusiven Transportterms

Der Anteil des diffusiven Transports zur Temperaturänderung (5.2) wird beschrieben durch
 

Ersetzt man die partiellen Ableitungen durch Differenzenquotienten der Form (5.8) und (5.14), so ergibt sich:

Aufgelöst nach der gesuchten Temperatur zum Zeitpunkt
 
ergibt sich ein explizites Differenzenschema, nach dem sich die Temperatur zu einem Zeitpunkt an einem Punkt (i,j,k) aus den Temperaturen am Punkt selbst und an den sechs nächsten Nachbarpunkten zu einem früheren Zeitpunkt t berechnen lä\3t.

Dies entspricht einer gewichteten Mittelung der Temperaturen einer lokalen Umgebung des gesuchten Punktes. Die Wichtungskoeffizienten sind symmetrisch bezüglich des Punktes (i,j,k). In der Sprache der Digitalen Bildverarbeitung entspricht dies einer Faltung der Temperaturverteilung T mit der Faltungsmaske h der Wichtungskoeffizienten:

Diese Äquivalenz zur Digitalen Bildverarbeitung legt es nahe, die dort entwickelten Konzepte der Optimierung von Filteroperationen zu nutzen. In der Digitalen Bildverarbeitung stellt sich gleicherma\3en das Problem, da\3 mathematische Operatoren auf einem diskreten Gitter möglichst fehlerfrei implementiert werden müssen, um eine optimale Extraktion physikalischer Information aus Bilddaten zu ermöglichen.

Zur Interpretation der Diffusionsgleichung im Sinne der Bildverarbeitung kann Gleichung (5.15) mit den diskreten Koordinaten (5.4) umgeschrieben werden zu
 
wobei nur die zeitliche Ableitung durch den Differenzenquotienten (5.8) ersetzt wurde. Für ein regelmä\3iges Gitter, d. h. , reduziert sich dies auf
 
Dabei ist als diskrete Realisierung des kontinuierlichen Laplaceoperators aufzufassen. Aufgelöst nach der gesuchten Temperatur
 
liefert (5.20) die Operatorgleichung:
 
wobei den Einheitsoperator bezeichnet.

In dieser Darstellung ergibt sich die Temperaturverteilung nach einem Zeitschritt durch die Summe der Temperaturverteilung vor dem Zeitschritt und der Temperaturverteilung nach Anwendung des Laplaceoperators.

Wie bereits in der Herleitung der diskreten partiellen Ableitungen beschrieben wurde, stellen Differenzenquotienten nur Näherungen einer optimalen Lösung dar, da die Taylor-Entwicklung nach dem zweiten bzw. dritten Term abgebrochen wurde. Wie gut ein diskreter Ableitungsoperator die kontinuierliche Lösung approximiert, lä\3t sich an seiner Transferfunktion, d. h. der Fouriertransformierten des Filters, ablesen. Für partielle Ableitungen in Richtung der Koordinate gelten die Zusammenhänge ([Jähne, 93a])
 
zwischen der Darstellung des Operators im Ortsraum und im Fourierraum. Die erste Ableitung einer Funktion f entspricht einer Multiplikation ihrer Fouriertransformierten mit dem imaginären Faktor , wobei die Wellenzahl in Richtung der -Achse bezeichnet. Durch zweimaliges Anwenden der ersten Ableitung ergibt sich die zweite Ableitung einer Funktion f als Multiplikation ihrer Fouriertransformierten mit dem negativen Quadrat der Wellenzahl.

Aus den partiellen, zweiten Ableitungen ergibt sich der dreidimensionale Laplaceoperator

der mit (5.24) im Fourierraum die Form
 
annimmt. Aus Gleichung (5.26) lä\3t sich erkennen, da\3 der ideale, dreidimensionale Laplaceoperator eine isotrope Transferfunktion besitzt. Das komplexe Spektrum einer Funktion wird dabei mit dem negativen Quadrat des Betrages k des Wellenzahlvektors multipliziert.

Um die Forderung der Isotropie zu erfüllen, wird in [Jähne, 93a] vorgeschlagen, die Isotropie von Binomial-Glättungsmasken zur Konstruktion eines dreidimensionalen Laplaceoperators auszunutzen. Die Transferfunktion einer n-dimensionalen Binomialmaske vom Grad 2 kann für kleine Wellenzahlen mit

approximiert werden, wobei k den Betrag des n-dimensionalen Wellenzahlvektors im Fourierraum bezeichnet. Wird das Originalbild von dem mit der Binomialmaske geglätteten Bild subtrahiert (), so ergibt sich für kleine Wellenzahlen eine Transferfunktion

Dabei bezeichnet den Einheitsoperator, der das Bild unverändert lä\3t. Die kombinierte Faltungsmaske
 
ist damit ein n-dimensionaler Ableitungsoperator zweiter Ordnung, d. h. ein n-dimensionaler Laplaceoperator mit isotroper Transferfunktion .

Aufgelöst nach liefert (5.29) den Zusammenhang
 

Ein Vergleich von (5.30) mit (5.22) liefert eine modifizierte, explizite Lösung der Diffusionsgleichung:
 

Unter Verwendung des diskreten Laplaceoperators der Form (5.29) ergibt sich die numerische Lösung der dreidimensionalen Diffusionsgleichung (5.15) somit als eine Glättung des Simulationsvolumens mit einer dreidimensionalen Binomialmaske der Ordnung 2. Dies entspricht anschaulich der Tatsache, da\3 bei einer Glättung der Wärmeinhalt eines Gitterpunktes auf die umliegenden Nachbarpunkte verteilt wird, wobei die Koeffizienten der Glättungsmaske die Form der Verteilung, d. h. die Punktantwort, wiedergeben. Die Faltung mit der Punktantwort ist somit die diskrete Realisierung des Konzeptes der Green'schen Funktion zur Lösung von Differentialgleichungen (siehe auch Kapitel 4.5).

Die Green'sche Funktion der Diffusionsgleichung ist eine normierte, gau\3förmige Verteilung, deren Breite mit der Quadratwurzel der Zeit anwächst (4.47). Eine optimale, diskrete Realisierung der Gau\3verteilung wird von den diskreten Binomialmasken dargestellt. Der Grad 2n gibt die Breite der Masken an und repräsentiert somit die Zeitdauer der Diffusion.

Eine bemerkenswerte Eigenschaft von Binomialmasken, die sich zur effizienten Implementierung ausnutzen lä\3t, ist die Tatsache, da\3 sich gro\3e Masken durch wiederholtes Anwenden kleinerer Masken aufbauen lassen. Da die Faltung eine lineare Operation darstellt, kann damit das Simulationsvolumen mehrmals hintereinander mit der elementaren Maske geglättet werden, um effektiv eine Glättung mit der grö\3eren Maske zu erreichen:

Dies entspricht genau der Lösung der Diffusionsgleichung nach Gleichung (5.31). Die Zeit wird in diskrete Schritte eingeteilt. Zu jedem Zeitschritt führt die Glättung mit der Maske zu einer Änderung der Temperatur gemä\3 der Diffusionsgleichung. Die Temperaturverteilung nach einer Zeit ergibt sich durch n-maliges Wiederholen dieses Schrittes.

Die physikalische Grö\3e der realen Zeit , die eine einmalige Glättung repräsentiert, ist durch die Bedingung in Gleichung (5.31) gegeben:

Nur für diese Kombination der diskreten Orts- und Zeitauflösung entspricht die Diffusion während eines Zeitschrittes einer Glättung mit der Binomialmaske .

Für die verwendete Tiefenauflösung von m und die Diffusionskonstante D = 0.0014cms von Wärme in Wasser ergibt sich der Zeitschritt zu s.

Aus der Digitalen Bildverarbeitung ergibt sich eine Methode dreidimensionale Binomialmasken effizient zu implementieren. Dabei wird die Eigenschaft der Separabilität der Masken ausgenutzt. Dies bedeutet, da\3 eine dreidimensionale Glättung mit einer Binomialmaske durch sukzessives Falten der dreidimensionalen Datenstruktur T mit eindimensionalen Binomialmasken entlang der drei Koordinatenachsen ersetzt werden kann:
 
Die eindimensionalen Masken , und haben dabei die Form
 
und laufen entlang der Achsen i, j und k.

Damit ergibt sich die Implementierung des dreidimensionalen Diffusionstermes in Gleichung (4.9) zu drei eindimensionalen Glättungen des dreidimensionalen Simulationsvolumens der Form
 

Diese Möglichkeit der Implementierung approximiert den isotropen, dreidimensionalen Laplaceoperator der Differentialgleichung besser als das explizite Differenzenverfahren der Form (5.17). Beide Implementierungen stellen jedoch diskrete Faltungen der Datenstruktur dar, die in einem Bildverarbeitungsprogramm mit wenig Aufwand realisiert werden können. Im Falle des verwendeten Programmes heurisko konnte die Simulation vollständig in der programmeigenen Hochsprache mit Standard-Bildverarbeitungsoperatoren implementiert werden (siehe Anhang B.1). Dies macht die enge Verwandtschaft zwischen Digitaler Bildverarbeitung und numerischer Mathematik deutlich.

Bemerkungen:

1.

Die Implementierung der dreidimensionalen Diffusion nach (5.36) ergab sich aus dem Wunsch, einen Laplaceoperator zu konstruieren, der die Forderung der Isotropie besser erfüllt als das direkte Differenzenverfahren nach Gleichung (5.17). Für die eindimensionale Diffusionsgleichung

ergibt sich mit (5.8) und (5.14) das explizite Differenzenschema
 
d. h.  eine Faltung mit der Maske . Für stellt dies die eindimensionale Binomialmaske der Form (5.35) dar. Für den eindimensionalen Fall sind somit die beiden Realisierungen (5.36) und (5.17) identisch.

2.

Die Wahl von in Gleichung (5.22) und (5.38) hatte das Ziel, die Binomialmaske als Lösung zu erhalten. Für andere Werte von ergeben sich alternative Realisierungen, die für feste Gitterabständen zu anderen Zeitschritten führen. Der grö\3te, sinnvolle Wert für ergibt sich aus der Bedingung, da\3 nicht mehr Wärme aus der zentralen Box verschwinden kann, als vor dem Zeitschritt in ihr enthalten war, d. h. . Dies führt zu der Bedingung . Genauere Berechnungen zeigen, da\3 dieser Maximalwert bereits zu numerischen Instabilitäten und Oszillationen des Ergebnisses führt ([Cranck, 93]). Der Wert liegt noch im Bereich der stabilen Lösungen. Gleichzeitig ist er gro\3 genug, um nicht zu sehr kleinen Zeitschritten und damit sehr langen Rechenzeiten zu führen. Es zeigt sich hierbei, wie auch bei vielen Anwendungen der Bildverarbeitung, da\3 die Binomialmaske die physikalische Realität in idealer Weise zu approximieren scheint und darüber hinaus sehr effizient zu implementieren ist.

3.

Beim Übergang von Gleichung (5.19) zu Gleichung (5.20) wurde vorausgesetzt, da\3 es sich bei der Diskretisierung des Volumens um ein regelmä\3iges Gitter handelt, d. h. . Nur in diesem Fall ist die Faltungsmaske isotrop bezüglich des diskreten Gitters und führt auf die dreidimensionale Binomialmaske als Realisierung. Im vorliegenden Fall wurde jedoch die vertikale Auflösung um einen Faktor 100 kleiner gewählt als die horizontale Auflösung (Abschnitt 5.3.1). Eine Anwendung der dreidimensionalen Binomialmaske auf diesem Gitter führt zu einer anisotropen Diffusion.

Zum Vergleich der Wirkung der Diffusion in horizontaler und vertikaler Richtung wird die dreidimensionale Binomialmaske nach (5.34) zerlegt und die Ausdehnung der einzelnen Masken bezüglich der kontinuierlichen Koordinaten betrachtet. Eine Binomialmaske vom Grad n, die auf einem Gitter der Gitterkonstante angewendet wird, stellt eine diskrete Realisierung einer Gau\3verteilung der Breite dar, wobei der Zusammenhang

zwischen n und gilt. Für ein unregelmä\3iges Gitter, mit , ergeben sich daher die Ausdehnungen
 
der Verteilung bezüglich der kontinuierlichen Koordinaten, nach einmaligem Anwenden der Binomialmaske in alle Raumrichtungen. Berücksichtigt man noch, da\3 die Breite der Gau\3verteilung bei der Diffusion mit der Quadratwurzel der Zeit anwächst, d. h. , dann entspricht die Glättung mit einer isotropen Gau\3maske auf dem unregelmä\3igen Gitter unterschiedlichen Zeitschritten in horizontaler und vertikaler Richtung. Mit

folgt aus (5.40):

Aufgrund der 100-fach gröberen Auflösung in horizontaler Richtung, entspricht eine Glättung entlang der horizontalen Achsen 10 Glättungen bezüglich eines Gitters der Auflösung in horizontaler Richtung. Um pro Zeitschritt eine isotrope Diffusion zu erhalten, kann die Ausdehnung der Glättungsmasken in horizontaler Richtung eingeschränkt werden. Dazu wird, wie gehabt, der Zeitschritt über die Bedingung (5.38) mit und der vertikalen Auflösung festgelegt. Dieser Zeitschritt liefert, mit der horizontalen Auflösung bzw. , die Bedingung für in horizontaler Richtung (5.38) und damit eine modifizierte Faltungsmaske. Diese stellt jedoch keine Binomialmaske mehr dar. Eine weitere Möglichkeit zur Realisierung einer isotropen Diffusion auf einem unregelmä\3igen Gitter ist es, die horizontale Glättung nur nach jeweils 10 vertikalen Glättungen durchzuführen. Dies reduziert den Rechenaufwand auf etwa 1/3 und führt auf das gleiche Ergebnis.

4.

Die Diskussion der Implementierung der Diffusionsgleichung in diesem Abschnitt zeigt deutlich, wie eine effiziente Implementierung durch die Konzepte der Digitalen Bildverarbeitung erreicht werden kann. Darüber hinaus wird die gewünschte Isotropie der Lösung besser approximiert als durch direkte Differenzenverfahren. Bei einer Implementierung der Form (5.17) sind pro Gitterpunkt 7 Multiplikationen und 6 Additionen nötig. Bei einer Realisierung der Form (5.36) können die Multiplikationen vollständig durch schnelle Shift-Operationen ersetzt werden. Für die drei Faltungen entlang der Koordinatenachsen reduziert sich der Rechenaufwand auf jeweils eine Shift-Operation des zentralen Punktes, zwei Additionen und eine abschlie\3ende Shift-Operation zur Normierung.