Modellierung von ultraschnellem Wärmetransport in COMSOL Multiphysics^®

Die Fortschritte in der Technologie ultrakurzer Laserpulse haben dazu geführt, dass es notwendig geworden ist, den Wärmetransport auf extrem kurzen Zeitskalen zu verstehen. Im Subnanosekundenbereich versagt die herkömmliche Wärmeleitungsgleichung, sodass wir verfeinerte Modelle wie hyperbolische Wärmeleitungsgleichungen und Zwei-Temperatur-Modelle einsetzen müssen. In diesem Blog-Beitrag geben wir einen kurzen Überblick über die relevante Theorie für die gängigsten Nicht-Fourier-Modelle des Wärmetransports und zeigen, wie diese in der Software COMSOL Multiphysics^® implementiert werden können.

Ein Blick auf die Wärmeleitungsgleichung

Zunächst geben wir einen kurzen Überblick über eine der in Naturwissenschaften und Ingenieurwesen am häufigsten vorkommenden partiellen Differentialgleichungen: die Wärmeleitungsgleichung für Festkörper. Unser Ausgangspunkt ist die Energiegleichung bzw. die Kontinuitätsgleichung für die Temperatur T:

wobei C=\rho C_p die volumetrische Wärmekapazität, \mathbf{q} der Wärmefluss und S die Volumenwärmequelle ist.

Dies ist genau genommen ein Ausdruck der Energieerhaltung und muss daher stets gelten, unabhängig vom räumlichen und zeitlichen Maßstab. Um eine geschlossene Form der partiellen Differentialgleichung für T zu erhalten, müssen wir einen Ausdruck für \mathbf{q} finden. Die einfachste und intuitivste Wahl ist das Fouriersche Gesetz, das besagt, dass der Wärmefluss proportional zum Temperaturgradienten ist:

Der Proportionalitätsfaktor k wird als Wärmeleitfähigkeit bezeichnet. Setzt man diesen in die Kontinuitätsgleichung ein, erhält man die bekannte Wärmeleitungsgleichung:

Es handelt sich hierbei um eine parabolische partielle Differentialgleichung, daher der Name parabolisches Ein-Schritt-Modell (parabolic one-step model, POS). (Die Bedeutung von „Ein Schritt“ wird später in diesem Blog-Beitrag deutlich werden.)

Wir können einen Einblick in das Verhalten dieser Gleichung gewinnen, indem wir ein halbunendliches Gebiet (x \geq 0) eines Mediums betrachten, wobei die Oberfläche bei x=0 auf der Temperatur T_0 gehalten wird, während die Anfangstemperatur des restlichen Gebiets T_i beträgt. Dieses Problem hat eine analytische Lösung:

wobei \alpha = k/C die thermische Diffusionsleitfähigkeit des Mediums ist und \mathrm{erf} die Fehlerfunktion bezeichnet. Eine Animation dieser Lösung ist unten dargestellt. Aus dem obigen Ausdruck geht klar hervor, dass die Auswirkung der Randbedingung bei x=0 sofort überall spürbar ist, was gegen die spezielle Relativitätstheorie verstößt. Dies ist eine Folge der recht starken Annahme, die bei der Herleitung dieser Gleichung verwendet wurde, nämlich des Fourierschen Gesetzes.

Trotz seines Namens handelt es sich bei dem Fourierschen Gesetz in Wirklichkeit um eine phänomenologische Beziehung, deren strenge Herleitung aus den Grundprinzipien mehrere Näherungen erfordert. Es ist jedoch anzumerken, dass der räumliche Energietransport auf üblichen Zeit- und Längenskalen schnell genug ist, um als augenblicklich angesehen zu werden, sodass dieser Mangel der Wärmeleitungsgleichung für die meisten ingenieurtechnischen Zwecke getrost vernachlässigt werden kann.

Animation der zeitlichen Entwicklung der Lösung der Wärmeleitungsgleichung für ein einfaches halbunendliches Gebiet.

Über das Fouriersche Gesetz hinaus

Wenn wir uns kürzeren Zeitskalen zuwenden, können wir nicht mehr davon ausgehen, dass der räumliche Energietransport augenblicklich erfolgt, wie wir es oben getan haben. Darüber hinaus müssen wir auch die endliche Geschwindigkeit des Energietransfers zwischen den Elektronen und dem Gitter berücksichtigen. Im Folgenden werden wir sehen, wie Ersteres zu hyperbolischem Wärmetransport und Letzteres zu Zwei-Temperatur-Modellen führt (Ref. 1).

Hyperbolischer Wärmetransport

Anstatt anzunehmen, dass der Wärmefluss augenblicklich auf den Temperaturgradienten reagiert, können wir folgendes „verzögertes“ Fouriersches Gesetz vorschlagen:

Die Verzögerungszeit \tau sollte der Elektronenrelaxationszeit entsprechen, d. h. der mittleren freien Weglänge der Elektronen. Sofern diese ausreichend klein ist, lässt sich der obige Ausdruck in eine Taylor-Reihe entwickeln, was zu folgendem Ergebnis führt:

Dieser Ausdruck ist als Cattaneo-Vernotte-Modell bekannt (Ref. 2). Setzt man ihn in die Kontinuitätsgleichung ein (nach Differenzierung dieser Gleichung nach der Zeit), erhält man das hyperbolische Ein-Schritt-Modell (HOS):

Diese hyperbolische partielle Differentialgleichung ähnelt der Wellengleichung, und wie wir im Folgenden sehen werden, führt die endliche Geschwindigkeit des Wärmetransports tatsächlich zum Auftreten von „thermischen Wellen“.

Zwei-Temperatur-Modelle

Die Aufhebung der Annahme, dass sich die Elektronen und das Gitter stets im lokalen thermischen Gleichgewicht befinden, erfordert die Einführung zweier getrennter Temperaturfelder (der Gittertemperatur T_l und der Elektronentemperatur T_e) sowie einer Kopplung zwischen diesen, G. Zunächst wird die Energie des Lasers von den Elektronen absorbiert, was zu einem starken Anstieg von T_e (Schritt 1) führt, der sich dann auf T_l auswirkt (Schritt 2). Die Wärmediffusion durch das Gitter ist im Vergleich zum elektronischen Transport auf Zeitskalen im Subnanosekundenbereich vernachlässigbar, sodass wir schließlich die folgenden Gleichungen erhalten:

Dies ist als parabolisches Zwei-Schritte-Modell (PTS-Modell) bekannt. Wir könnten T_e eliminieren, um eine geschlossene Form der partiellen Differentialgleichung für T_l zu erhalten, doch dies führt zu umständlichen Kreuztermen. Es ist zweckmäßiger, diese Form einer gekoppelten partiellen Differentialgleichung und einer gewöhnlichen Differentialgleichung für das Gebiet beizubehalten. Eine Schätzung für den Elektron-Gitter-Kopplungsparameter G lässt sich durch Integration der Elektron-Photon-Energieaustauschrate über alle möglichen Übergänge erhalten (Ref. 3).

Wir können zusätzlich die endliche Ausbreitungsgeschwindigkeit der Wärmeenergie der Elektronen berücksichtigen, um das hyperbolische Zwei-Schritte-Modell (HTS-Modell) zu erhalten. In Analogie zum obigen Modell lässt sich das Modell, wenn man die Gitterdiffusion vernachlässigt, wie folgt ausdrücken:

Weitere Verallgemeinerungen

Es gibt einige weitere Effekte, die wir möglicherweise berücksichtigen müssen, um genaue Ergebnisse zu erhalten, wie zum Beispiel:

• Berücksichtigung der Temperaturverzögerung zusätzlich zum Wärmefluss, was zum Dual-Phase-Lag-Modell (DPL) führt
• Durch Einbeziehung der Gitterdiffusion in die oben genannten Zwei-Temperatur-Modelle entstehen die Dual-Parabolic-Two-Step- (DPTS) und die Dual-Hyperbolic-Two-Step- (DHTS) Modelle
• Die explizite Modellierung der Ladungsträgerdynamik mittels einer Drift-Diffusions-Formulierung mit zwei Temperaturen, bekannt als das semiklassische Zwei-Schritte-Modell (STS)
• Berücksichtigung des Strahlungswärmetransports

Auch wenn eine detaillierte Besprechung den Rahmen dieses Blog-Beitrags sprengen würde, ist es für viele Anwendungen von entscheidender Bedeutung, komplexere Prozesse zu berücksichtigen, die während der ultraschnellen Erwärmung ablaufen können, wie beispielsweise Ablation und Schmelzen.

Implementierung des Nicht-Fourierschen Wärmetransports in COMSOL Multiphysics^®

Standardmäßig unterstützen die Wärmetransport-Interfaces in COMSOL Multiphysics^® ausschließlich das POS-Modell; die Implementierung hyperbolischer Wärmetransport- und Zwei-Temperatur-Modelle in der Software mithilfe der PDE- und ODE-Interfaces ist jedoch unkompliziert.

Als Beispiel betrachten wir einen Goldfilm mit einer Dicke von 200 nm, der mit einem Laserpuls bestrahlt wird, dessen absorbierte Fluenz \Phi _p = 1\,\mathrm{J}/\mathrm{m}^2 beträgt, wobei diese hier unabhängig von der Elektronen- und Gittertemperatur angenommen wird, und der eine gaußförmige zeitliche Verteilung der Pulsbreite t_p = 1 \, \mathrm{ns}, \; 100 \, \mathrm{ps}, \; 100 \, \mathrm{fs} und 10 \, \mathrm{fs}. \Phi _p enthält Informationen über die Wellenlänge des Lasers, da diese auch von der Reflektivität des Films abhängt. Wir nehmen außerdem an, dass der Strahlradius im Vergleich zur Eindringtiefe d des Strahls groß ist, sodass das Problem eindimensional wird. Somit lässt sich die durch den Puls verursachte Volumenwärmequelle wie folgt schreiben (Ref. 1):

Die Vorfaktoren gewährleisten die korrekte Normalisierung der gesamten absorbierten Pulsenergie. Dieser Puls ist bei t=0 zentriert, daher starten wir die Simulation bei t=-2t_p, um den gesamten Puls zu erfassen. Wir übernehmen die in Ref. 3 angegebenen Materialparameter für Gold und gehen davon aus, dass diese konstant bleiben.

Hyperbolischer Wärmetransport

Um die Zeitableitung zweiter Ordnung des HOS-Modells einzubeziehen, müssen wir das Interface General Form PDE verwenden. Die erforderlichen Einstellungen sind unten dargestellt.

Implementierung des HOS-Modells in COMSOL Multiphysics^®. Um die Zeitableitung zweiter Ordnung einzubeziehen, verwenden wir das Interface General Form PDE.

Zwei-Temperatur-Modelle

Zur Implementierung des PTS-Modells können wir die Elektronentemperatur mit dem vordefinierten Interface Heat Transfer modellieren, und da die Gitterdiffusion in diesem Modell vernachlässigt wird, reicht eine gewöhnliche Differentialgleichung für das Gebiet aus, um die Gittertemperatur zu berechnen. Die erforderlichen Einstellungen sind unten dargestellt.

Implementierung des PTS-Modells in COMSOL Multiphysics^®. Für die Elektronentemperatur (links) können wir das Standardinterface Heat Transfer mit einer zusätzlichen Quelle aufgrund der Kopplung verwenden. Für die Gittertemperatur (rechts) verwenden wir das Interface Distributed ODE, wobei die Kopplung als Quellterm festgelegt wird.

In Anlehnung an die oben genannten Einstellungen können wir das HTS-Modell auch implementieren, indem wir für die Elektronentemperatur anstelle des Interfaces Heat Transfer ein weiteres Interface General Form PDE verwenden. Beachten Sie, dass in diesem Modell die Wärmeleitfähigkeiten als konstante Skalare betrachtet werden; sowohl das Interface Heat Transfer als auch das Interface General Form PDE unterstützen jedoch Koeffizienten, die beispielsweise temperaturabhängig sind und tensorielle Eigenschaften aufweisen, um anisotropen Wärmetransport zu modellieren.

Ergebnisse und Diskussion

Sehen wir uns zunächst die Ergebnisse für t_p = 1 \, \mathrm{ns} an. Die räumlichen Verteilungen der Gittertemperaturen nach der Bestrahlung sind in der folgenden Abbildung dargestellt. Wie zu erwarten war, lassen sich bei dieser Pulsdauer fast identische Ergebnisse mit den verschiedenen Formulierungen feststellen.

Plot der Gittertemperatur als Funktion der Tiefe bei einer Pulsdauer von 1 ns. Alle vier Modelle liefern auf dieser Zeitskala nahezu identische Ergebnisse.

Als Nächstes betrachten wir den gleichen Plot für t_p = 100 \, \mathrm{ps}, wie in der folgenden Abbildung dargestellt. Es zeigt sich eine deutliche Diskrepanz zwischen den Ein-Schritt-Modellen (POS und HOS) und den Zwei-Schritt-Modellen (PTS und HTS): Die Ein-Schritt-Modelle sagen aufgrund der Vernachlässigung des relativ schnelleren Elektronenwärmetransports eine höhere Gittertemperatur in der Nähe der bestrahlten Oberfläche voraus. Zwischen den parabolischen Modellen (POS und PTS) und den hyperbolischen Modellen (HOS und HTS) ist jedoch noch kein wahrnehmbarer Unterschied festzustellen.

Plot der Gittertemperatur in Abhängigkeit von der Tiefe bei einer Pulsdauer von 100 ps. Die Ein-Schritt-Modelle (POS und HOS) sagen höhere Gittertemperaturen in der Nähe der bestrahlten Oberfläche voraus.

Diese Diskrepanz wird umso deutlicher, je weiter wir uns in den Femtosekundenbereich mit t_p = 100 \, \mathrm{fs} bewegen, wie aus dem nachstehenden Plot ersichtlich ist.

Plot der Gittertemperatur als Funktion der Tiefe bei einer Pulsdauer von 100 fs. Der Einfluss der in den hyperbolischen Modellen (HOS und HTS) berücksichtigten Zeitableitung zweiter Ordnung wird allmählich sichtbar.

Nun liegen die Vorhersagen der Ein-Schritt-Modelle um mehrere hundert Grad höher, und die hyperbolischen Modelle weisen aufgrund der endlichen Geschwindigkeit des Elektronenwärmetransports in der Nähe der bestrahlten Oberfläche höhere Temperaturen auf als die entsprechenden parabolischen Modelle.

Bislang haben wir ausschließlich die Gittertemperaturen betrachtet. Wir können uns einen Eindruck vom unterschiedlichen Verhalten der Elektronen- und der Gittertemperatur verschaffen, indem wir die Elektronentemperatur auf demselben Plot, jedoch mit einem anderen Maßstab, darstellen. Das Ergebnis ist im Folgenden als Animation zu sehen.

Animation, die die zeitliche Entwicklung der Elektronen- (links) und Gittertemperaturverteilungen (rechts) für den 100-fs-Puls zeigt.

Um den Effekt der Einbeziehung der Zeitableitung zweiter Ordnung weiter zu verdeutlichen, vergleichen wir abschließend die von den Zwei-Schritt-Modellen vorhergesagte Elektronentemperaturverteilung für t_p = 10 \, \mathrm{fs} . Nun lässt sich deutlich eine sich in den Film ausbreitende thermische Welle erkennen, die in der folgenden Abbildung dargestellt ist. Es sei darauf hingewiesen, dass die Elektronentemperatur bei dieser Pulslänge Werte erreicht, die hoch genug sind, um die Einbeziehung des Strahlungswärmetransports zu rechtfertigen, der im vorliegenden Modell jedoch vernachlässigt wird.

Plot der Elektronentemperatur in Abhängigkeit von der Tiefe bei einer Pulsdauer von 10 fs. Wir beobachten die Ausbreitung einer vom HTS-Modell vorhergesagten Wärmewelle.

Um schließlich einen Eindruck von der Dynamik des Wärmetransports durch den Dünnfilm zu gewinnen, können wir die zeitliche Entwicklung der Temperatur an der Vorderseite (bestrahlte Oberfläche) und an der Rückseite als Plot darstellen. Aufgrund der großen Diskrepanz im Maßstab müssen wir jedoch jede Datenreihe anhand ihrer jeweiligen Maximaltemperatur normieren.

Das Ergebnis ist in der folgenden Abbildung für eine Pulsdauer von 100 fs dargestellt. Auf dieser Zeitskala ist der Unterschied zwischen dem parabolischen und dem hyperbolischen Modell moderat, doch die Ein-Schritt-Modelle können das qualitative Verhalten der Temperatur an der Rückseite nicht vorhersagen. Die Ein-Schritt-Modelle zeigen eine monoton ansteigende Temperatur, während die Zwei-Schritt-Modelle vorhersagen, dass die Temperatur an der Rückseite ein Maximum erreicht, da die Elektron-Gitter-Relaxation langsamer abläuft als der Elektronenwärmetransport. Diese Ergebnisse stimmen gut mit den Abb. 3 und 4 aus Ref. 3 überein.

Zeitliche Entwicklung der normierten Elektronentemperatur (bei Zwei-Schritt-Modellen) und der Gittertemperatur (bei Ein-Schritt-Modellen) bei einer Pulsdauer von 100 fs an der Vorder- und Rückseite des Films.

Fazit

In diesem Blog-Beitrag haben wir gesehen, wie COMSOL Multiphysics^® zur Simulation ultraschneller Wärmetransportprozesse verwendet werden kann, indem vier einfache und wichtige Nicht-Fourier-Modelle für den Wärmetransport implementiert werden, nämlich das parabolische Ein-Schritt-Modell (POS), das hyperbolische Ein-Schritt-Modell (HOS), das parabolische Zwei-Schritt-Modell (PTS) und das hyperbolische Zwei-Schritt-Modell (HTS).

Nächste Schritte

Laden Sie die Modelldatei für das hier vorgestellte Beispiel herunter:

Referenzen

1. V.E. Alexopoulou and A.P. Markopoulos, “A Critical Assessment Regarding Two-Temperature Models: An Investigation of the Different Forms of Two-Temperature Models, the Various Ultrashort Pulsed Laser Models and Computational Methods,” Arch Computat Methods Eng, Vol. 31, S. 93–123, 2024. https://doi.org/10.1007/s11831-023-09974-1
2. Y. Zhang, D.Y. Tzou, and J.K. Chen, “Micro- and Nanoscale Heat Transfer in Femtosecond Laser Processing of Metals”, High-Power and Femtosecond Lasers: Properties, Materials and Applications, eds. P.H. Barret and M. Palmerm, Nova Science Publishers, Inc., Hauppauge, NY, K. 5, S. 159–206, 2009. https://arxiv.org/abs/1511.03566
3. T.Q. Qiu and C.L. Tien, “Heat Transfer Mechanisms During Short-Pulse Laser Heating of Metals.” ASME J. Heat Transfer, Vol. 115, A. 4, S. 835–841, 1993. https://doi.org/10.1115/1.2911377