PAMM · Proc. Appl. Math. Mech. 8, 10377 – 10378 (2008) / DOI 10.1002/pamm.200810377

Grundlagen zur Beschreibung der Wärmeleitung mitZellulären Automaten

Michael Müller 1,∗ und Georg-Peter Ostermeyer1, ∗∗1 Technische Universität Braunschweig, Institut für Dynamik und Schwingungen, Schleinitzstraße 20, 38106 Braunschweig

Für die Beschreibung der Temperaturleitung werden oftmals die für die Zeitintegration stabilen impliziten Lösungsverfahren(zum Beispiel in Form von Finite-Elemente-Modellen) verwendet, was für ein Zelluläre-Automaten-Modell jedoch aufgrundder dann extremen Zunahme der Rechenzeit nicht zweckmäßig ist. Stattdessen wird die Dynamik der Wärme im Sinne desZellulären Automaten mit einem expliziten Algorithmus wiedergegeben, der auf dem lokalen Austausch von Energiemengen,also auf Bilanzgleichungen zwischen den Zellen basiert. Dieses Paper zeigt Aussagen zur numerischen Stabilität in Hinblickauf Zeit- wie Ortsdiskretisierung auf.

c© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

1 Ein Algorithmus zur Beschreibung der Wärmeleitung mit Zellulären Automaten

Der hier vorgestellte Algorithmus basiert auf dem Grundgedanken, dass während eines Zeitschrittes Wärmeenergie zwischenbenachbarten Zellen ausgetauscht wird. Zur Anschauung dieses Vorgehens soll Abb.1 links dienen.

Abb. 1 Wärmeaustausch zwischen benachbarten Zellen links: Zelle MZ und ihre Nachbarn Mitte: Schema zum Austausch zwischen Typ-1-Nachbarn rechts: Plot der Funktion aus Gleichung (2)

Die nun folgende Updateregel ist exemplarisch für die Zelle MZ (die Zelle in der Mitte eines 3x3x3-Zellen-Quaders)aufgestellt. Gleichung (1) beschreibt die Temperaturänderung der Zelle MZ zum Zeitpunkt 1 gegenüber dem Zeitpunkt 0 alsFunktion der Temperaturdifferenzen zu ihren 26 Nachbarn i (Abb. 1 links).

TMZ,1 = TMZ,0 −26∑

i=1

ki (TMZ,0 − Ti,0) (1)

Demnach wird eine Angleichung der Temperaturen mit einem zugehörigen Faktor ki vorgenommen. Dieser Faktor steht fürden zeitlichen Gradient dieser Angleichung und hängt im physikalischen und numerischen Sinne von der Wärmeleitfähigkeitsowie Orts- und Zeitdiskretisierung ab. Er muss sorgfältig gewählt werden, um numerische Instabilitäten auszuschließen.

Basis des hier vorliegende Problems ist das bestehende Zelluläre-Automaten-Modell zur Beschreibung der Grenzschicht-dynamik von Bremssystemen. Zur detaillierten Beschreibung des Wachstums harter, dünner und reibintensiver Strukturen,sogenannter ’Patches’ an der Oberfläche des Bremsbelages ist eine Zellengröße von 10 µm festgelegt worden [1]. Die Tem-peraurleitfähigkeit a entspricht der von Grauguss (Bremsscheibenmaterial) und wurde zu a = 14 mm2/s angenommen. Esexistieren bei der implementierten Nachbarschaft nach Moore (alle an MZ angrenzenden Zellen) 3 Typen von Nachbarn.Nachbarn des Typs 1 teilen eine gemeinsame Fläche mit MZ, ihre Gesamtzahl ist 6. Demgegebnüber existieren 12 Zellenvom Typ 2 (eine gemeinsame Kante mit MZ) sowie 8 Zellen vom Typ 3 (ein gemeinsamer Punkt mit MZ). Der ki -Wert vonZellentyp 1 muss größer sein als der von Zellentyp 2, der wiederum größer ist als der von Zellentyp 3, weil mit Zellentyp auchder Abstand zwischen den jeweiligen Zellenmittelpunkten steigt.

Grundlage für eine überschlägige Ermittlung der ki -Faktoren sei zunächst das System eines eindimensionalen Wärme-leiters (Abb.1 Mitte). Zwei benachbarte Zellen haben zu einem Zeitpunkt 0 die unterschiedlichen Temperaturen T1 bzw. T2.

∗ Corresponding author E-mail: mi.mueller@tu-braunschweig.de, Phone: +49 531 391 7005, Fax: +49 531 391 7017∗∗ E-mail: gp.ostermeyer@tu-braunschweig.de, Phone: +49 531 391 7000, Fax: +49 531 391 7017

c© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

10378 Sessions of Short Communications 05: Oscillations

Die Temperaturdifferenz ∆T12,0 zwischen ihnen beträgt demzufolge ∆T12,0 = T1,0 − T2,0. Die analytische Lösung ihrerTemperaturdifferenz zu einem bestimmten Zeitpunkt t und der zugehörige ki -Faktor beträgt dann:

∆T12 (t) = ∆T12,0 exp(−2a

z2t

)→ k = 1 − exp

(−2a

z2t

)(2)

Der Parameter z entspricht dabei der Zellengröße bzw. dem Abstand zwischen den Zellenmittelpunkten, also 10 µm. Abb.1 rechts zeigt den Plot der linken Funktion aus Gleichung (2). Daraus wird ersichtlich, dass die zu wählenden Zeitschrittefür das Gesamtmodell möglichst klein (also im Sub-Mikrosekundenbereich) zu wählen ist. Bei Zeitschritten von etwa 5 µsbeträgt die entstandene Temperaturdifferenz nur noch 25% der Temperaturdifferenz vor dem Zeitschritt. Da eine Zelle des-weiteren 6 Nachbarn dieses Typs besitzt, wäre eine instabile Numerik beziehungsweise eine Verletzung des 2. Hauptsatzesder Thermodynamik die Folge.

2 Ergebnisse

Die ermittelten Werte für ki werden nun anhand eines Systems mit einem einmaligen Wärmeimpuls Q0 an der Oberflächeeines homogenen und isotropen Halbraumes verifiziert und quantifiziert. Die analytische Lösung dieses Problems (unter derAnfangsbedingung T (t = 0, r) = 0) lautet:

T (r, t) =Q0

4π√

πλt√

t√

aexp

(− r2

4at

)(3)

wobei r der Abstand von der Punktwärmequelle und λ die Wärmeleitfähigkeit des Materials (50 W/Km bei Grauguss)ist. Aus der analytischen Lösung können zeitliche Verläufe der gemittelten Temperaturen in den Zellen berechnet werden.Diese werden denen, die über Gleichung (1) berechnet werden, gegenübergestellt. Dabei zeigt sich, dass die in (2) berechneteBeziehung zwischen Zeit und ki - Wert zu deutlichen Abweichungen zwischen analytischer Lösung und Lösung nach (1)führt, siehe Abb. 2 links. In dieser Abbildung sind exemplarisch die Verläufe für die Zelle AM, in deren Mitte sich die Impuls-Punktwärmequelle befand und für einen Vertreter des Nachbartyps 1 zu Zelle AM geplottet. Die Ursache der Diskrepanz liegtdarin begründet, dass eine Diskretisierung von 10 µm erheblich zu grob ist, um die örtlichen Temperaturgradienten an derWärmequelle hinreichend genau approximieren zu können [2].

Mit Hilfe eines erweiterten Ansatzes, der auf der analytischen Lösung aus Gleichung (3) basiert, kann dieses Problemauch ohne feinere Diskretisierung gelöst werden. Eine Neuberechnung der zugehörigen Faktoren ki zu k1 = 2.227 · 10−3,k2 = 8.733 · 10−5 und k3 = 3.969 · 10−6 für Zeitschritte von 10−8 s führt demnach zum Funktionsverlauf, der in Abb.2 rechts dargestellt ist [2]. Es ist ersichtlich, dass die so berechneten Faktoren ki die Temperaturdynamik erheblich besserbeschreiben. Diese Vorgehensweise kann als eine Art ’Kalibrierung’ des Verfahrens verstanden werden, was eine Ausweitungauf komplexere Systeme (z.B. auch Wechselwirkungen zwischen bewegten Temperaturfeldern, siehe [3]) ermöglicht.

Abb. 2 Vergleich zwischen analytischer Lösung und Lösung nach Gleichung (1) links: über Gleichung (2) berechnetes ki rechts: überverfeinertes Verfahren berechnetes ki, ’an’ entspricht ’analytisch’, ’al’ entspricht ’Algorithmus’, ’N1/N2/N3’ entspricht ’Nachbar 1/2/3’

Literatur

[1] G.-P. Ostermeyer und M. Müller New developments of friction models in brake systems, SAE Paper 2005-01-3942, Proceedings of the23rd SAE Brake Colloquium and Exhibition, Orlando (2005), 187-200, weiterhin in SAE 2005 Transactions Journal of Passenger Cars:Mechanical Systems, 3078-3090

[2] M. Müller, Zur Topographie- und Grenzschichtdynamik in Bremssystemen, Dissertation, Shaker Verlag, Aachen (2007)[3] K. H. Bode und G.-P. Ostermeyer, Wechselwirkung bewegter Temperaturfelder mit Zellulären Automaten, PAMM (2008)c© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.gamm-proceedings.com