Kristalle wachsen im Rechner

Simulation eines Dendriten, der aus einem kleinen Saatkristall wächst

Das physikalische Problem
Beim Schmelzen oder Auftauen von festen Materialien sind die Wärmeleitung und der Abtransport der entstehenden Schmelzwärme die wesentlichen physikalischen Effekte. Gleiches gilt für langsame Gefriervorgänge. Beim schnellen Gefrieren oder Erstarren von unterkühlten Schmelzen werden zusätzliche Oberflächeneffekte an der Grenze zwischen festem und flüssigem Material wichtig, die zur Entstehung von schnell wachsenden verzweigten Strukturen, sogenannten Dendriten, führen. Phasenübergänge mit dendritischem Wachstum treten in vielen Anwendungen auf, wie z.B. beim Erstarren von Metallschmelzen in der Stahlerzeugung oder beim Gefrieren von Lebensmitteln oder Blutplasmen. Das wesentliche gestaltbildende Prinzip kann schon an reinen, einkomponentigen Materialien studiert werden. Das große Interesse von Anwendern an einem genauen Verständnis der zugrundeliegenden Gesetze wird auch durch die Tatsache belegt, dass mehrmals das Erstarrungsverhalten von Materialien unter Mikrogravitation in Space Shuttle Experimenten beobachtet wurde (http://zeta.lerc.nasa.gov/EXPR2/IDGE.HTM). Materialwissenschaftler sind an der Modellierung des dynamischen Verhaltens der wachsenden Dendriten interessiert. Man beobachtet bei Experimenten dass sich die Spitzen der Dendriten mit einer konstanten Geschwindigkeit bewegen und einen festen Krümmungsradius besitzen. Das Prinzip der Auswahl dieser Wachstumsparameter in Abhängigkeit von Daten wie der Unterkühlung der Schmelze ist ein noch immer nicht ganz gelöstes Problem. Experimente sind teuer und aufwendig; deshalb sind numerische Simulationen von großem Interesse.

Berechenbar machen
Das mathematische Modell für das Dendritenwachstum besteht aus mehreren nichtlinearen partiellen Differentialgleichungen für die Temperatur und für die Position der Phasengrenze. Wenn man diese Gleichungen mit dem Rechner lösen möchte, so müssen sie erst für den Computer in geeigneter Form bereitgestellt werden. Man ersetzt deshalb das kontinuierliche Problem durch ein diskretes. Dies geschieht so, dass man die zu bestimmenden physikalischen Größen nicht überall berechnet, sondern nur noch an Punkten eines Gitters. Dann sind typischerweise 600\,000 Gleichungen mit ebensovielen Unbekannten auf dem Rechner zu lösen. Trotzdem ist diese Zahl noch relativ gering, da der verwendete Algorithmus auf einem adaptiven Verfahren beruht. Würde man keine Adaptivität der Gitter verwenden, so wären 10 Millionen Unbekannte zu berechnen. Eine wichtige algorithmische Idee besteht in der Verwendung sogenannter {\em Finiter Elemente}. Dies ist eine numerische Methode, die dazu führt, dass in jeder der vielen Gleichungen nur wenige Unbekannte vorkommen und die Unbekannten gleichmäßig über die Gleichungen verteilt sind. Eine weitere wichtige Technik ist die schon erwähnte Verwendung {\em adaptiver Gitter}. Das Rechengitter wird dort besonders fein gewählt, wo Strukturen gut aufgelöst werden müssen, z. B. wo die Lösung stark variiert. Dagegen kann man an Stellen, an denen die Lösung fast konstant ist, ein sehr grobes Gitter verwenden. Allerdings hängt die Feinheit des Gitters offensichtlich von der Lösung ab, die man ja gerade ausrechnen will. Mathematische Methoden liefern aber Formeln dafür, an welcher Stelle das Gitter wie fein zu wählen ist, um eine gegebene Fehlertoleranz einzuhalten. Das Besondere an diesen sogenannten Fehlerschätzern ist, dass sie nur Daten verwenden, die auf dem Computer vorhanden sind und nicht etwa Daten, die nur theoretisch bekannt sind. Bei der Berechnung des Wachstums eines Dendriten muß das Gitter in jedem Zeitschritt neu angepaßt werden. Die besondere Schwierigkeit bei der numerischen Simulation ist die Berechnung des freien Randes, das heißt der Phasengrenze.

Adaptiv verfeinertes Tetraedergitter zur Berechnung der Temperatur

Der Freie Rand
Neben der Wärmeleitung in der festen und flüssigen Phase spielen die Verhältnisse an der Phasengrenze zwischen festem und flüssigem Material eine wesentliche Rolle für die Entstehung und das Wachstum von Dendriten. Auf molekularer Ebene betrachtet geschieht das Erstarren, also der übergang vom flüssigen zum festen Zustand, durch Einordnung der im flüssigen Zustand relativ frei beweglichen Moleküle in den festen Verband einer Kristallstruktur. Die Schmelztemperatur, also diejenige Temperatur, bei der der übergang von flüssigem zu festem Zustand eintritt, wird von den kinetischen und geometrischen Eigenschaften der Phasengrenze beinflusst. Je schneller sich die Grenze fortbewegt und je stärker sie gekrümmt ist, desto tiefer ist die Schmelztemperatur. Dieses Verhalten kann durch das scheinbar einfache Gesetz

beschrieben werden. Dabei ist $V$ die Geschwindigkeit der Grenze, $K$ die (mittlere) Krümmung der Grenze, und $T$ die Temperatur beim Phasenübergang. Die Parameter $\alpha$ und $\beta$ hängen vom Material ab und sind darüberhinaus richtungsabängig, was sich aus der geometrischen Struktur des zugrundeliegenden Kristallgitters des Materials begründet. Durch diese Parameter ergeben sich dann auch die zu beobachtenden räumlichen Symmetrien der entstehenden Strukturen, wie man sie auch von Schneeflocken kennt.

Mean Curvature Flow
Dieses Bewegungsgesetz für die Phasengrenze war den letzten zehn Jahren Gegenstand der Forschung in Geometrie und Analysis. Es ist in mathematischer und numerischer Hinsicht ein höchst schwieriges Bewegungsgesetz. Die Forschung hat aber neben vielen theoretischen Resultaten zum Mean Curvature Flow auch effiziente und besonders stabile Algorithmen zur Berechnung der Bewegung der Phasengrenze (des Freien Randes) gemäß diesem Gesetz geliefert. Die Schwierigkeiten mit dem Bewegungsgesetz beginnen damit, dass es gar nicht so klar ist, wie man dem Rechner beibringt, die Krümmung der Phasengrenze zu berechnen. Mit Methoden aus der mathematischen Theorie der Flächen gelingt dies. Außerdem muß ein Algorithmus für das Bewegungsgesetz numerisch stabil sein. Das bedeutet, dass das numerische Verfahren gegenüber den immer auftretenden Rundungsfehlern unempfindlich sein muss. Dies gelang mit einer geeigneten zeitlichen Aufteilung der Berechnungsschritte.

Algorithmik
Das zeitabhängige Problem wird in diskreten Zeitschritten gerechnet. Dabei können die Wärmeleitung und die Bewegung der Phasengrenze getrennt diskretisiert und berechnet werden. In jedem Zeitschritt wird zunächst die Bewegung der diskreten Phasengrenze aus der oben erläuterten Gleichung berechnet. Die Phasengrenze selbst ist durch ein Dreicksnetz dargestellt, für dessen Knoten ein System von gekoppelten Bewegungsgleichungen hergeleitet werden kann, welche die Geschwindigkeit, die Krümmung und die Temperatur an der Phasengrenze beinhalten. Da sich die Phasengrenze im Laufe der Simulation stark vergrößert und verzweigt, ist es notwendig das Dreiecksnetz ständig anzupassen. Dazu werden z.\ B. zu groß gewordene Dreiecke in kleinere zerteilt. Ist die neue Phasengrenze berechnet, so kann die neue Temperaturverteilung und damit die Temperatur an der neuen Phasengrenze berechnet werden. Dabei wird das Temperatur-Rechengitter mit Hilfe der oben erwähnten Fehlerschätzer automatisch der neuen Lösung angepasst. Ein typisches Rechengitter ist im Bild gezeigt. Die Temperatur variiert sehr stark in der Nähe der Phasengrenze, daher muss dort das Gitter auch besonders fein sein. Demgegenüber ist die Temperatur innerhalb der festen Phase fast konstant. In diesem Bereich darf daher ein grobes Gitter verwandt werden. Im Rahmen des geförderten Projekts hat Dr. Andreas Veeser in seiner Promotion Konvergenz und Fehlerabschätzungen für eine Variante des Verfahrens bewiesen.

Dreiecksgitter zur Diskretisierung der Phasengrenze