VERFAHREN ZUR SIMULATION VON WAREMSTRAHLUNGN ZWISCHEN OBERFLÄCHEN
BESCHREIBUNG
Die Erfindung betrifft ein Verfahren und einer dazugehörigen Algorithmus für die Berechnung der thermischen Oberflächenstrahlung, und ihre Anwendung zur Simulation vom Einfluss der thermisch gekoppelten Oberflächenstrahlung, besonders mit Bezug auf
Gießprozessen .
Weiter wird die Erfindung im Rahmen weitere hilfreiche Verfahren und Algorithmen beschrieben, die hier als Teile des Verfahren und des Algorithmus für die Berechnung der thermischen Oberflächenstrahlung, und ihre Anwendung zur Simulation vom Einfluss der thermisch gekoppelten Oberflächenstrahlung, beschrieben sind, wessen
Anwendungen jedoch nicht nur auf den hiesigen Verfahren und Algorithmus beschränkt sind.
Die Erfindung betrifft deswegen außerdem ein Verfahren zur Diskretisierung eines Raumwinkels zur Anwendung in einem Simulations- oder Berechnungsprozess um Einsparnisse an Computerzeit und verbrauchte
Computermemory zu erzielen.
Weiter betrifft die Erfindung deswegen ein Verfahren des Ray Tracings zur Anwendung in einem Simulations- oder Berechnungsprozess um Beschleunigung und Einsparnisse an Computerzeit und verbrauchte Computermemory zu erzielen. Besonders wird einen unkonventionelle Verfahren mittels
Parallel Computing zur Beschleunigung einer Ray Tracing Berechnung alleine und in seiner Kombination mit Anisotrope Chebychev-Distanz Berechnungen und/oder zusätzlicher Akzeleration durch Kachel-Vereinigung vorgeschlagen.
INHALTSVERZEICHNIS Hintergrund
Beschreibung der Abbildungen
Grundsätze
Energiebilanz an der strahlenden Oberfläche
Berechnung des einfallenden Strahlungsflusses
Strahlungskachel
Diskretisierung des Raumwinkels
Hierarchisches System
Anpassung an die Geometrie und an die Temperaturverteilung
Ray Tracing
Beschleunigung
Ray Tracing
Kachel-Vereinigung
Anisotrope Chebychev-Distanz
Paralleles Ray Tracing
Thermische Berechnung
Anwendungsbeispiele
HINTERGRUND
Viele wichtige Eigenschaften, wie z.B. die Viskosität, die elektrische Leitfähigkeit, das spezifische Volumen oder das chemische Reaktionsvermögen eines Materials oder
Mediums werden von der Temperatur bestimmt. Z.B. im Bereich der Erschmelzung, Läuterung und/oder der formgebende Bearbeitung fester Stoffe, wie z.B. Metalle, Glas, Keramik, kommt der Einhaltung bestimmter Temperaturen, häufig auch zeitabhängig, große Bedeutung zu .
Die Erfindung betrifft ein Verfahren und einer dazugehörigen Algorithmus für die Berechnung der thermischen Oberflächenstrahlung, und ihre Anwendung zur Simulation vom Einfluss der thermisch gekoppelten Oberflächenstrahlung, besonders mit Bezug auf
Gießprozessen . Die Erfindung betrifft weiter ein Verfahren zur Simulation vom Einfluss der thermisch gekoppelten Oberflächenstrahlung durch Berechnen der radiativen Austausch zwischen graue, diffuse Oberflächen dadurch gekennzeichnet dass die zu bestrahlende Oberfläche oder Oberflächen adaptiv, hierarchisch in Strahlungskacheln gleicher oder nahezu gleicher Strahlungsintensität untergeteilt wird und die sich durch Einstrahlung ergebenen Oberflächentemperatur mittels Ray Tracing als die Summe aller Teilbeträge der Strahlungskacheln zu berechnen.
Die Erfindung betrifft weiter ein Verfahren wie oben erwähnt worin der Raumwinkel in seine Teilbereiche durch sphärische Projektion adaptiv und hierarchisch diskretisiert wird. Die Erfindung betrifft weiter einem Verfahren wie oben erwähnt worin die Ray Tracing Prozedur beschleunigt ist.
Weiter wird die Erfindung im Rahmen weitere hilfreiche Verfahren und Algorithmen beschrieben, die hier als Teile des Verfahren und des Algorithmus für die Berechnung der thermischen Oberflächenstrahlung, und ihre Anwendung zur Simulation vom Einfluss der thermisch gekoppelten Oberflächenstrahlung, beschrieben sind, wessen
Anwendungen jedoch nicht nur auf den hiesigen Verfahren und Algorithmus beschränkt sind. Im weiteren betrifft die Erfindung ein Verfahren zur Simulation vom Einfluss thermisch gekoppelter Oberflächenstrahlung an einem Festkörper, welcher Festkörper mindestens eine bestrahlbare Oberfläche aufweist; durch berechnen des radiativen Austauschs zwischen grauen, diffusen Oberflächen; dadurch gekennzeichnet dass die zu bestrahlende Oberfläche oder Oberflächen adaptiv, hierarchisch in Strahlungskacheln gleicher oder nahezu gleicher Strahlungsintensität untergeteilt wird, und die sich durch Einstrahlung ergebene Oberflächentemperatur mittels eines hierarchischen Sichtfaktoren-Verfahrens erfolgt; welche Sichtfaktoren-Verfahren die Auswertung eines
Raumwinkelintegrals unter Anwendung einer primären Raumwinkelunterteilung umfasst; welche primäre Raumwinkelunterteilung eine homogene Sichtfaktoren- Diskretisierung umfasst, wobei jene
Raumwinkelunterteilung in seine Teilbereiche durch sphärische Projektion adaptiv und hierarchisch diskretisiert wird und wobei die Summe aller Teilbeträge jenes Raumwinkelintegrals mittels Ray Tracing zu ermitteln ist.
Die Erfindung betrifft des Weiteren ein Verfahren als Oben angegeben wonach jenes Ray Tracing beschleunigt ist, insbesondere wonach jenes Ray Tracing durch Strahlungskachelvereinigung beschleunigt ist, wonach jenes Ray Tracing durch einen anisotropen Chebychev- Distanz-Verfahren beschleunigt ist, und/oder wonach jenes Ray Tracing mittels Parallel-Computing beschleunigt ist.
Außerdem betrifft die Erfindung des Weiteren einem Computersoftwareprodukt auf einem computerlesbaren
Medium, das einen Softwarecode zum Ausführen eines Verfahren nach dem Oben angegebenen.
Bindick, Ahrenholz und Krafczyk hat in „Heat Transfer - Mathematical Modelling, Numerical Methods and Information Technology", ISBN 978-953-307-550-1 , Aziz Belmiloudi Editor, Gedruckt 2011-02-14 , Kapitel 7 - „Efficient Simulation of Transient Heat Transfer Problems in Civil Engineering" die Stand der Technik innerhalb des technische Gebiet der thermisch gekoppelter Strahlung in drei Dimensionen zusammengefassen (siehe insbesondere Abschnitt 7.3) . Ihre Betrachtungen sind hiermit in ihre Gesamtheit bei Referenz in dieser Arbeit miteinbezogen.
EP 1 667 069 AI beschreibt ein Verfahren zur Bestimmung der Verteilung der lokalen Strahlungsintensität in einem semitransparenten Medium durch Anwendung von Ray Tracing, wobei ein signifikant schnelleres Verfahren mit gleichzeitig geringerem Speicherbedarf zur Bestimmung der Verteilung der lokalen Strahlungsintensität in einem semitransparenten, mindestens eine Grenzfläche aufweisenden Medium hergebracht wurde.
BESCHREIBUNG DER ABBILDUNGEN
Abbildung 1.
Veranschaulichung der Raumwinkel durch eine halbe Einheitssphäre.
Abbildung 2.
Beispiel der fest eingeprägten Diskretisierung des Raumwinkels .
Abbildung 3.
Veranschaulichung vom gewählten Prinzip der festen Unterteilung der Einheitssphäre, die Sicht „von oben" entlang der normalen Richtung der Strahlungskachel. Die Unterteilung ähnelt einer Zielscheibe und besteht aus den azimutal angeordneten Segmenten, die ihrerseits radiale Ringe bilden.
Abbildung 4.
Ein Beispiel der aus der Gleichung (7) ergebenden Systems der Richtungsvektoren für eine Parametrisierung (r=4, n=10) . Die gesamte Anzahl der Richtungen nach (7a) ist gleich N=86. Abbildung 5.
A: Die zweite Verfeinerungsebene, erstellt aus der ersten Ebene mit Unterteilung (r, n) = (4, 10) , siehe Abb. 4. B: die dritte Ebene, erstellt aus der zweiten Ebene, siehe Abb. 4. C: Richtungsvektoren aus allen drei Ebenen. Vektoren aus unterschiedlichen Ebenen sind mit unterschiedlichen Symbolen am Ursprung der Vektoren gekennzeichnet, vergl . mit Abb. 4, 5.A, 5.B. Die Strahlen für alle in Abb. 5.C gezeigten Richtungen werden nach dem Ray Tracing Verfahren verschickt.
Abbildung 6.
Das Ablaufdiagramm für die Auswahl der Strahlungsquellen aus unterschiedlichen Verfeinerungsebenen der Winkel- Diskretisierung .
Abbildung 7.
Ein Beispiel der angepassten Diskretisierung des Raumwinkels. Die sphärischen Vierecke aus unterschiedlichen Ebenen sind in der Abbildung in Richtung der Verfeinerung gekennzeichnet.
Abbildung 8.
Das Schema der Strahlenverfolgung auf dem kartesischen Raster. Sie besteht aus der Detektion einer Sequenz der Schnittpunkte mit dem orthogonalen Raster.
Abbildung 9.
Unterstützte Varianten der Kachelkonstruktion, ausgehend von einzelnen strahlenden Seiten der orthogonalen Gitterzellen (Kreuze) .
Abbildung 10.
Veranschaulichung der Chebychev-Distanz anhand von 2D- Beispiele. 10A: Isotrope Chebychev-Distanz. 10B Anisotrope Chebychev-Distanz.
Abbildung 11.
Ein Beispiel der Strahlenverfolgung in 2D, unterstützt durch anisotrope Chebychev-Distanz.
Abbildung 12.
Beispiel eines mit dem Strahlungsmodell berechneten thermischen Ergebnisses für ein Feingussprojekt. Beispiel einer realen Anwendung.
Abbildung 13.
Elefantenfamilie: Test der Modelleffizienz, Effekte der Verdeckungen und des Abstands zur Wärmequelle. Die Wärme kommt von der Kugel oben links im Bild.
Abbildung 14.
Modellierung der gleichzeitigen Wärmebehandlung von 4 Zahnrad-Aggregaten in einem Ofen mit halb geschlossenen Hohlraum. A: das gesamte Modell. B: temperierte Aggregate. Die Wärmequellen sind an der linken und rechten Wand angebracht.
Abbildung 15.
Berechnete Szene mit drei Wägen. Die Wärmequelle ist vorne oben, im Bild nicht gezeigt.
Abbildung 16.
Wärmebehandlung von 5 gegossenen Pumpengehäusen in einem Ofen.
GRUNDSÄTZE
Energiebilanz an der strahlenden Oberfläche
Die Berücksichtigung der Oberflächenstrahlung in einem thermischen Modell erfolgt durch die Anpassung der Gleichung für die Bedingung der Energiebilanz an der strahlenden Oberfläche.
Ohne Strahlung wird diese Bedingung allein durch die Wärmeleitungsflüsse ausgedruckt, (l). ÄL T-ÄR T)N = 0
Die Größen mit den Indizes L und R in (1) kennzeichnen die entsprechend links und rechts von der Oberfläche evaluierten Wärmeleitfähigkeiten und Gradienten der Temperatur. Der Vektor N ist der Normalenvektor der Oberfläche .
Ein zusätzlicher Term, die Netto-Strahlungsflussdichte, kommt durch Strahlung in die Bilanzgleichung (1) :
( 2 ) . (AL LT - λΚΤκΤ )iV + qnet = 0
Die Gleichung (2) muss in einem thermisch gekoppelten Modell mit Oberflächenstrahlung gelöst werden.
Der Netto-Strahlungsfluss besteht aus der Differenz der absorbierten und der emittierten Strahlungswärme
(3) . qnet = e(Rin ~ «"Γ4) e in (3) ist die Emissivität der Oberfläche. Sie wird im Falle der sogenannten grau und diffus strahlenden Oberfläche über das gesamte elektromagnetische Spektrum und über den Raumwinkel hemisphärisch gemittelt. Sie wird auch gleich dem ebenfalls gemittelten Absorptions- koeffizienten angenommen. qin in (3) ist die einfallende Strahlungsflussdichte (incident heat flux density) . Die Aufgabe eines
numerischen Modells für Oberflächenstrahlung ist, die Werte der einfallenden Strahlung an jedem Stück der strahlenden Oberfläche zu ermitteln. Berechnung des einfallenden Strahlungsflusses
Der einfallende Strahlungsfluss wird durch ein Integral an den ausgehenden Strahlungsfluss gekoppelt: (4). qin = fi >0(ÜN)qoutd[l
Es wird über den Raumwinkel 2π integriert. Das Ergebnis gilt für das Zentrum des Stücks der strahlenden Oberfläche .
Das Verfahren der Sichtfaktoren ersetzt das Integral in (4) durch eine Summe über die Beiträge der einzelnen Unterteilungen der strahlenden Oberfläche, für die die direkte Sichtbarkeit durch lineare Optik gegeben ist. Als einzelne Elemente der Oberfläche gelten in einem normalen Fall die Flächen des numerischen Gitters, wobei auch die Vereinigung von mehreren benachbarten Flächen des Gitters in solche strahlenden Elemente nach dem hierarchischen Sichtfaktoren-Verfahren bekannt ist.
Der Raumwinkel wird oft durch eine halbe Einheitssphäre veranschaulicht wie in Abbildung 1 gezeigt. Eine Einheitssphäre wird um das Zentrum der strahlenden Oberfläche gelegt. Das umliegende sichtbare strahlende Oberflächengitter wird auf die Einheitssphäre zentral projiziert. Nach dem Nusselt-Prinzip ist die zur äquatorialen Ebene parallele Komponente der Projektion auf die Einheitssphäre gleich dem jeweiligen Sichtfaktor.
Die Berechnung der Sichtfaktoren ist in Abbildung 1 geometrisch veranschaulicht. Um den Wert eines Sichtfaktors zwischen der Oberfläche j draußen und der Oberfläche i (grau schattiert im Bild) zu ermitteln, soll zuerst der sichtbare Teil von Fläche j auf die Einheitssphäre um das Zentrum der Fläche i projiziert werden. Die Berechnung des gesamten Integrals (4) bedeutet die Projektion des gesamten strahlenden Gitters auf die Einheitssphäre, wie es in Abb. (1) verdeutlicht ist .
Die Komplexität der Quadratur (4) ist generell proportional das Quadrat der Unterteilungen der strahlenden Oberfläche durch das numerische Gitter. Deswegen kann die Berechnung nach (4) zu einem übermäßigen Speicher- und Rechenaufwand besonders für große, komplexe Geometrien führen. Das geometrische Schema der Sichtfaktoren-Verfahren nach Abb. (1) kann umgedreht werden, was im Prinzip zu einem kleineren Speicher- und Rechenaufwand führen kann.
Im Falle der Sichtfaktoren ist die Unterteilung (=Diskretisierung) des Raumwinkels durch das zuvor erzeugte numerische Gitter fest vorgegeben. Die Feinheit der Unterteilung ist komplett durch das numerische Gitter gegeben .
Es ist aber umgekehrt möglich, eine beliebige Unterteilung des Raumwinkels zunächst unabhängig von dem numerischen Gitter am Rande des transparenten Hohlraums im ersten Schritt fest vor ;ugeben. Im zweiten Schritt wird diese Unterteilung der Einheitssphäre dann auf das
umliegende Gitter hinaus projiziert. Das Prinzip ist in Abb. 2 erklärt.
Dieses Schema führt zu einer anderen Logik bei der Integration des Strahlungsflusses nach Gleichung (4). Während die Aufgabe in der Sichtfaktoren-Verfahren aus der Ermittlung der einzelnen Unterteilungen des Raumwinkels besteht, wobei die Unterteilungen ein Abbild der sichtbaren Unterteilungen der Oberfläche sind, sind die einzelnen Unterteilungen des Raumwinkels in der alternativen Verfahren fest vorgegeben und deswegen im Voraus bekannt. Die eigentliche Aufgabe ist jetzt die Suche nach einem repräsentativen Element des strahlenden Gitters, das die Rolle der Strahlungsquelle für das gegebene Stück der Unterteilung des Raumwinkels übernimmt. Es wird also angenommen, dass die von der gefundenen Strahlungsquelle ausgehende Strahlungs¬ intensität in der gesamten Unterteilung des Raumwinkels homogen vorliegt (Siehe Abb. 2) .
Bei der fest eingeprägten Diskretisierung des Raumwinkels wird im folgendem Vorgegangen. Im ersten Schritt wird eine bestimmte Unterteilung der Einheitssphäre um das Zentrum der ausgewählten strahlenden Fläche i unabhängig vom numerischen Gitter vorgenommen. Im zweiten Schritt wird eine Zentralprojektion dieser Unterteilung auf das umliegende Gitter durchgeführt. Die Richtung der Projektion ist durch Pfeile nach oben verdeutlicht, vergleiche mit der umgekehrten Projektionsrichtung in der Verfahren der Sichtfaktoren in Abb. 1. Die Projektion des Mittelpunktes jedes Teils des Raumwinkels wird einem Element des strahlenden Gitters zugeordnet, so dass eine Abbildung (Mapping) der gesamten Unterteilung des
Raumwinkels auf die strahlenden Flächenelemente zustande kommt .
Ein Fehler bei der Diskretisierung der Verteilung der Strahlungsflusses über den Raumwinkel entsteht, falls der betrachtete Teil des Raumwinkels, projiziert auf das strahlende Gitter, Flächenelemente mit stark unterschiedlichen ausgehenden Strahlungsflussdichten hat, was zu einer großen Variation der Strahlungsintensität innerhalb des Raumwinkels führt.
Der Vorteil dieses Verfahren ist, das die komplizierte geometrische Analyse der Verdeckungen entfällt. Für jede Unterteilung des Raumwinkels ausgehend vom Zentrum der Fläche wird ein einziger Probestrahl verschickt. Die Strahlenrichtung entspricht dem Mittelpunkt des gegebenen Teils des Raumwinkels.
Auf das Verschicken der Strahlen, oder Ray Tracing, ausführlicher im Abschnitt „Ray Tracing" eingegangen.
Strahlungskachel
Das Verfahren der Unterteilung der strahlenden Flächen ist unten exemplarisch dargestellt mittels des aktuell verwendeten numerischen Gitters im Programm MAGMAsoft (MAGMA Gießereitechnologie GmbH) . Das Verfahren ist aber im Prinzip für beliebige Gittertypen ohne Begrenzung anwendbar .
MAGMAsoft benutzt ein Tensorprodukt-Gitter. Das 3D-Gitter besteht aus drei Systemen der Gitterlinien in 3 kartesischen Raumrichtungen X, Y, und Z, die sich durch das gesamte Modell durchziehen und es in einen Quader
einbettet. Das Gitter besteht also aus den orthogonalen quaderförmigen Zellen.
Das numerische Gitter besteht also aus einem orthogonalen Raster und wird daher allein durch die Materialverteilung über die Gitterzellen in diesem Raster und drei Reihen der Koordinaten entlang der drei kartesischen Richtungen komplett vorgegeben. Die einzelnen strahlenden Oberflächen im Strahlungsmodell werden bei ihrer Initialisierung durch Materialnachbarschaften im Gitter gefunden. Eine strahlende Oberfläche wird als eine rechtwinklige Facette einer Gitterzelle mit dem Normalenvektor in einer der 6 Raumrichtungen +X, -X, +Y, -Y, +Z, -Z definiert, falls eine der 2 Bedingungen erfüllt ist:
1. Die Facette trennt zwei Gitterzellen, wovon eine mit einem opaken Material und eine andere mit einem transparenten Material belegt sind. Der Normalenvektor zeigt in Richtung der transparenten
Zelle. Das transparente Material in MAGMAsoft ist normalerweise Luft, wobei in diesem Fall allein die Wärmeleitung berechnet wird.
2. Die Facette trennt zwei Gitterzellen, wovon eine mit einem opaken Material und eine andere mit dem
Material ID „boundary" belegt sind. Die Zelle mit dem Material „boundary" liegt außerhalb des Rechengebiets von MAGMAsoft. Eine solche Facette wird als eine strahlende Fläche definiert, falls sie nicht an der Grenze der Bounding Box des
Gitters liegt, und zwar deshalb, weil dort von der Oberfläche aus gesehen keine andere Facette sichtbar ist. In diesem Fall ist die Bestimmung des Netto-Strahlungsflusses trivial.
Die auf solche Weise definierten strahlenden Flächen tragen zur Energiebilanz nach Gleichung (2) bei. Sie werden im Weiteren als Strahlungskachel oder Kachel bezeichnet.
DISKRETISIERUNG DES RAUMWINKELS
Der komplette Raumwinkel wird in diesem Verfahren so unterteilt, dass jede Unterteilung dem gleichen Sichtfaktor VF entspricht, nämlich VFj = 1/N = const, wobei N die gesamte Anzahl der Unterteilungen ist. Eine solche im Bezug auf Sichtfaktoren homogene Unterteilung der Einheitssphäre ist keine homogene Unterteilung ihrer Fläche wegen des Terms (HN) im entsprechenden Integral:
(5) . /Ω (HV)an = / /θ sin(G) cos(6)dcpdQ = ^ = const
Der Vorteil einer solchen Unterteilung ergibt sich bei der Berechnung des einfallenden Strahlungsflusses nach Gleichung (4) . Der stets derselbe Sichtfaktor kann ausgeklammert werden, nachdem die Ermittlung des einfallenden Flusses in eine Mittelwertbildung des ausgehenden Flüsse überführt wird.
Eine nach den Sichtfaktoren homogene Unterteilung der Einheitssphäre ist nicht eindeutig und kann auf unzähligen Wegen erfolgen. Das hier gewählte Verfahren bittet zusätzlich eine gewisse Symmetrie in der Unterteilung im Bezug auf den Normalenvektor der Strahlungskachel und ist einfach in der Handhabung.
Die Einheitssphäre wird zunächst axialsymmetrisch, von einem Kreis am Nordpol angefangen, in eine Reihe von radial aufeinander folgenden Ringe unterteilt. Anschließend wird jeder Ring in eine unterschiedliche Anzahl der Ringsegmente in azimutaler Richtung unterteilt. Jedes Segment ist ein sphärisches Viereck, begrenzt durch 2 azimutalen und 2 meridionalen Koordinatenlinien (Kreisbögen) des sphärischen Koordinatensystems. Die Anzahl der Unterteilungen in einzelnen Ringen bilden eine arithmetische Progression (Siehe Abb. 3) .
Die Unterteilung kann komplett durch die Anzahl der Ringe n in meridionaler Richtung und die Anzahl der azimutalen Segmente des ersten Rings am Nordpol r parametrisiert werden. Für jede mögliche Parametrisierung (n,r) gibt es eine einzige Lösung für die meridionalen Koordinaten der Ringe, so dass jedem sphärischen Segment immer derselbe Sichtfaktor zugeordnet werden kann:
,^, 1 ,Λ (n-£+l)(2r+n+£-2\
7 . Θ; = -2arccosfvl -- - -)
N in (7) ist die gesamte Anzahl der Unterteilungen nach der arithmetischen Progression,
(7a). N= 2^ln + 1
2
Die Anzahl der Unterteilungen N und damit auch die Anzahl der Strahlen pro eine Kachel wächst nach (7a) quadratisch mit der Anzahl der vorgegebenen meridionalen Ringen n (Siehe Abb. 4) .
Die Segmente in einem Ring aus der Abb. 3 können in der azimutalen Richtung um einen freien Winkel gedreht werden, damit ein größerer Winkelabstand zwischen den Segmenten der benachbarten Ringe entsteht.
Hierarchisches System
Eine zu grobe Unterteilung der Einheitssphäre führt zu den numerischen Fehlern in der Integration des Strahlungsflusses, wie es im Abschnitt „Energiebilanz an der strahlenden Oberfläche" erklärt wurde. Um eine höhere Winkelauflösung im Bezug auf die Geometrie des strahlenden Gitters zu erreichen, wird in diesem Verfahren mit einem hierarchischen System der Diskretisierungsebenen gearbeitet. Dieses Verfahren ist den unterschiedlichen Verfeinerungsstufen des numerischen Gitters in einem expliziten Mehrgitterverfahren ähnlich.
Die erste Ebene stellt das nach der Gleichung (7) erzeugte System der Raumrichtungen, siehe Abb. 4. Die nächste Ebene wird durch Unterteilung von jedem Segment der ersten Ebene durch die doppelte Halbierung in azimutaler und in der meridionaler Richtung, also durch eine Viertelung erstellt.
Die Ausnahme gibt es bei der ersten Verfeinerung des kreisförmigen Gebiets am Nordpol, die mit 4 azimutalen Unterteilungen in 4 sphärischen Dreiecken geteilt wird. Bei weiteren Unterteilungen gelten keine Ausnahmeregeln, alle sphärischen Gebiete werden als sphärische Vierecke, wie oben beschrieben verfeinert. Die 4 sphärischen Dreiecke rund um den Nordpol werden als entartetem sphärischem Vierecke ab der dritten Verfeinerungsebene behandelt .
Danach wird die nächstfeinere Ebene weiter rekursiv nach gleichen Regeln unterteilt. So kommt es zu einer geometrischen Progression von der Anzahl der Strahlen auf jeder nächsten Ebene. Für k Ebenen der Verfeinerung ergeben sich insgesamt
, ^ „ N»rtot = N»r—4k-l = ,2r+n-l +, „.—k-l
(7b).
Strahlen auf allen Verfeinerungsebenen.
Die Richtungsvektoren nach der Verfeinerung sind in Abbildung 5 gezeigt. Drei Verfeinerungsebenen werden in der aktuellen Implementierung benutzt.
In dem hierarchischen Verfahren von Ray Tracing werden zunächst Strahlen von allen Verfeinerungsebenen nacheinander verschickt und alle dadurch erhaltenen Ergebnisse gespeichert. Die Ergebnisse werden dann ausgewertet .
Anpassung an die Geometrie und an die Temperatur¬ verteilung
Es wäre durchaus möglich, die im vorangehenden Abschnitt beschriebene hierarchische Unterteilung des Raumwinkels virtuell vorzunehmen, und ausschließlich die Strahlen für die Vektoren der feinsten erreichten Ebene zu versenden. Die Konsequenz wäre, alle nach dem Ray Tracing gefundenen Strahlungsquellen müssen im thermischen Modell prozessiert werden, was bei einer großen Anzahl der Strahlen zur einem unnötig hohen Rechenaufwand geführt hätte. Wenn zum Beispiel die eingestellte Winkelauflösung
durch die Anzahl der Strahlen in einer Raumrichtung zu hoch ausfällt, würde die gleiche durch eine Strahlungskachel definierte Strahlungsquelle durch unterschiedliche Strahlen mehrfach getroffen. An dieser Stelle wäre also die Winkelauflösung aus den gröberen Verfeinerungsebenen ausreichend.
Wenn allerdings die Strahlen von allen definierten Verfeinerungsebenen versendet werden, besteht die Möglichkeit für die lokale Anpassung der Winkelauflösung. Das Ergebnis von Ray Tracing für die thermische Berechnung ohne eine Anpassung wäre eine komplette Liste der Strahlungsquellen, die sich alle entweder aus der primären Unterteilung (Abb. 4) oder aus einer der Verfeinerungsebenen (Abb. 5) ergeben.
Die lokale Anpassung der Winkeldiskretisierung ersetzt eine solche Liste durch Unterlisten aus allen Diskretisierungsebenen . Das Vorgehen ist dabei das Folgende.
Die Unterteilungen von allen Verfeinerungsebenen werden als „nicht verfeinert" vormarkiert. Zunächst wird angenommen, dass die Diskretisierung nur Vektoren der feinsten Ebene enthält. Dann geht man die Unterteilungen der nächst gröberer Ebene durch, bis die gröbste primäre Ebene erreicht ist.
Jedes Mal werden die Strahlungsquellen der vorangegangenen feineren Ebene geprüft, die den vier Vierteln des unterteilten sphärischen Vierecks der feineren Ebene entsprechen. Die zuvor definierten aktiven Strahlungsquellen von den 4 Vierteln werden entweder
durch eine Strahlungsquelle aus dem aktuellen sphärischen Segment ersetzt, oder beibehalten.
1. Wenn mindestens eine der 4 Unterteilungen bereits als „verfeinert" markiert ist, werden alle bereits definierten Strahlungsquellen auf den feineren Ebenen innerhalb des aktuellen sphärischen Segments beibehalten, die Quelle aus der aktuellen Ebene wird verworfen .
2. Wenn alle 4 Strahlen der feineren Ebene das gleiche opake Material getroffen haben, oder alle 4 Strahlen „ins Leere" gingen, d. h. den halbgeöffneten transparenten Hohlraum verlassen haben, werden die bereits definierten Strahlungsquellen der feineren Ebene verworfen und durch die Strahlungsquelle der aktuellen Ebene ersetzt.
3. Wenn 4 Strahlen sowohl den Hohlraum sowohl verlassen als auch ein opakes Material treffen, werden die bereits definierten Strahlungsquellen aus der feineren Ebenen beibehalten.
4. Wenn 4 Strahlen unterschiedliche opake Materialien treffen, werden die Temperaturen der getroffenen Strahlungskacheln miteinander verglichen. Wenn dabei das hierarchische System der Richtungen zum ersten Mal vor der thermischen Simulation angepasst wird, kommen Temperaturen der Quellen aus der ersten Initialisierung der Temperatur für unterschiedliche Materialien. Ansonsten, wenn die Anpassung im Laufe der thermischen Berechnung dynamisch erfolgt, sind das die aktuell berechneten Temperaturen. Falls die absolute Differenz zwischen der maximalen und der minimalen Temperatur, bezogen auf die maximale Temperatur kleiner als die fest vorgegebene Grenze ist, werden die bereits definierten
Strahlungsquellen der feineren Ebene ebenfalls verworfen und durch eine Quelle aus der aktuellen Ebene ersetzt. Ansonsten werden sie beibehalten.
5. Nachdem eine Unterteilung in der aktuellen Ebene bearbeitet ist, wird sie je nach dem Ergebnis entweder als „verfeinert" oder „nicht verfeinert" markiert. Wenn dabei zur Ersetzung der Quellen der feineren Ebenen durch eine Quelle der aktuellen Ebene kommt, wird diese in die Liste der berücksichtigten Quellen der aktuellen
Verfeinerungsebene eingetragen. Wenn eine der Quellen als „nicht verfeinert" markiert ist, wird sie Quelle in die Liste der berücksichtigten Quellen für die feinere Ebene eingetragen. Die Eintragung erfolgt, nur falls die neue Quelle einer Oberfläche des opaken Materials gehört. Ansonsten wird der Sichtfaktor des äußeren nicht verdeckten Raums um einen konstanten Wert der aktuellen
1
Verfeinerungsebene ι gleich—— erhöht.
Der oben beschriebene Ablauf ist in Abbildung 6 dargestellt. Der Algorithmus liefert k Listen der Strahlungsquellen für k Ebenen der homogenen
Diskretisierung . Die angepasste Winkeldiskretisierung wird dabei so gewählt, dass eine Verdichtung der repräsentativen Strahlungsquellen in den Winkelbereichen entsteht, wo eine unstetige Abhängigkeit der Strahlungsintensität von der Raumrichtung potenziell zu erwarten ist.
Der Sprung in der Strahlungsintensität kommt in den halbgeöffneten Hohlräumen an einer Grenze zwischen dem heißen opaken Material und dem offenen Raum vor. Ein
solcher Sprung wird auch durch den thermischen Kontrast von zwei im Winkelraum benachbarten opaken Materialien verursacht. Die mit unterschiedlichen IDxs gekennzeichneten Materialien können entlang ihrer Grenze im Winkelraum entweder sich im direkten thermischen Kontakt befinden, oder ein anderen zum Teil verdecken.
Die angepasste Diskretisierung des Raumwinkels für die betrachtete Strahlungskachel erfolgt auf diese Weise durch die Auswahl der Strahlungsquellen. Die unterteilte Oberfläche der Einheitssphäre kann für die Darstellung der Winkeldiskretisierung verwendet werden. Sie wird nach der Anpassung lückenlos und ohne Überschneidungen durch sphärische Segmente aus unterschiedlichen Verfeinerungs- ebenen bedeckt.
Ein Beispiel der angepassten Diskretisierung des Raumwinkels ist in Abbildung 7 gezeigt. Ein Sichtfaktor von einer H-förmigen Fläche, in Abbildung links oben gezeigt, soll durch Ray Tracing angenähert werden. Dafür wurden 3 Diskretisierungsebenen mit der ersten Ebene parametrisiert mit (n=15,r=4) nach Gleichung (7) benützt. Die Anzahl der Strahlen ist gleich 166, 664 und 2656 für die entsprechend 1-e, 2-e und 3-e Diskretisierungsebene, so dass insgesamt 3486 Strahlen für das Abtasten der Hemisphäre verschickt wird. Nach der Durchführung der geometrischen Anpassung sind 218 Strahlen ausgewählt worden, was um einen Faktor 16 kleiner ist als an der feinsten Ebene. Die sphärischen Vierecke aus unterschiedlichen Ebenen sind in der Abbildung sortiert in Richtung der Verfeinerung durch ihre Größe gekennzeichnet .
Die ausgehenden Strahlungsflüsse, vorgegeben durch Quellen aus unterschiedlichen Ebenen, müssen in der Berechnung des einfallenden Strahlungsflusses unterschiedlich gewichtet werden. Die Gleichung (6) für den einfallenden Strahlungsfluss wird durch Vorfaktoren von unterschiedlichen Ebenen ergänzt:
Für drei Ebenen ist k=3. Drei Summen mit den Gewichtsfaktoren 1, H und 1/16 ergeben nach (8) den gesamten einfallenden Strahlungsfluss . Die Terme q
o ] ut in der inneren Summe sind die Beiträge aus der Liste der Strahlungsquellen der jeweiligen Diskretisierungsebenen „level".
RAY TRACING Ein Voxel-gestütztes Ray Tracing Verfahren mit der Rückwärtsverfolgung der Strahlen wird eingesetzt.
Ein Voxel unterstützt den Test, ob der Strahl im gegebenen Voxel ein Objekt trifft. Ein Voxel ist ein würfelförmiges Volumen, wo die Informationen über geometrische Objekte in ihm enthalten sind. Bei der Anwendung mit MAGMAsoft bietet sich ein Tensorprodukt- Gitter als aus den einzelnen Voxeln bestehender Raster selbst an. Jede Gitterzelle wird zu einem Voxel. Das Verfahren gemäß der Erfindung ist jedoch nicht an das Programm MAGMAsoft gebunden sondern ist für beliebige Voxel-gestütztes Ray Tracing Verfahren einsetzbar.
Die Rückwärtsverfolgung bedeutet, dass zunächst alle Strahlen von der Kachel, die Strahlung empfängt, hinaus verschickt werden. Die von Strahlen getroffenen Strahlungsquellen werden dabei ermittelt. Die Strahlungsquellen senden aber Energie physikalisch zum Empfänger. Die Energie erreicht die Kachel auf dem „geraden" Weg, was einer im Strahlengang von Ray Tracing umgekehrten Richtung entspricht. Deswegen wird das als Ray Tracing mit Rückwärtsverfolgung gekennzeichnet.
Die geometrischen Objekte, die die Strahlen absorbieren und reflektieren, sind die definierten Strahlungskacheln, die mithilfe der Materialnachbarschaft in den Gitterzellen definiert werden.
Es ist wegen des strukturierten Charakters des numerischen Gitters ausreichend, 3 IDxs der Kachel je Gitterzelle zu speichern. Wenn man 6 Seiten einer Gitterzelle als front, back, west, east, nord, south bezeichnet, verweisen die 3 IDxs auf drei möglichen Kacheln an der Back-, East- und North-Seite der Zelle. Die Information über die Kachel an drei verbleibenden Seiten der Zelle Top, West und South können aus der Nachbarzelle entnommen werden. Die positiven Werte werden für die reell existierenden Kacheln an einer Seite der Zelle zugewiesen. Die Zuweisung erfolgt während der Definition der Kachel.
Vor dem Beginn der Strahlenverfolgung wird der gesamte einmal definierte Bündel der Vektoren, die Richtungen einzelner Strahlen vorgeben, um das Normalenvektor der jeweiligen Kachel zentriert. Das wird durch die Multiplikation der Vektoren mit einer Drehmatrix erreicht, die den zentralen Vektor des Bündels in
kartesischer Richtung +Z in den Normalenvektor der Kachel überführt .
Danach werden die Strahlen aus allen Verfeinerungsebenen die Reihe nach verschickt. Die Strahlenverfolgung durch das numerische Gitter besteht aus der Suche nach Schnittpunkten zwischen der Fortsetzung des Strahls und den einzelnen Seiten der Gitterzelle, in der sich der Strahl gerade befindet (Siehe Abb. 8) .
Dafür werden drei möglichen Seiten der Gitterzelle für die nächste Überschneidung geprüft. Die möglichen Seiten sind durch die Vorzeichen der drei Komponenten der
Strahlenrichtung Ω gegeben. Der Schnittpunkt mit dem minimalen Abstand von der aktuellen Position des Strahls gilt als der nächste Punkt auf dem Strahl. Die Seite der Zelle mit der minimalen Länge bis zum Schnittpunkt wird dabei gesucht, und der Strahl bis zum gefundenen
Schnittpunkt in Richtung Ω um die gefundene Länge Xray fortgesetzt:
(9a) . XMIN = M/NL=(ij-fc)-^—
(9b). l^X™* = xMINnt
Nachdem der nächste Schnittpunkt gefunden ist, wird ID der entsprechenden Kachel an der Seite der Gitterzelle abgefragt . · Falls ID einem reellen Kachel entspricht, d. h. die
Seite der Zelle trennt ein opakes und ein transparentes Material, wird die Verfolgung des
Strahls gestoppt und das gefundene globale ID der Strahlungsquelle zurückgeliefert .
• Falls ID einer Symmetrieebene mit einem kartesischen Normalenvektor entspricht, wird die zur Symmetrieebene normale Richtungskomponente invertiert und der reflektierte Strahl weiterverfolgt .
• Ansonsten wird die Prozedur der Strahlenverfolgung in der nächsten Gitterzelle wiederholt, bis der Strahl eine Kachel trifft oder die Grenzen des numerischen Gitters verlässt. In letzten Fall wird ein festes ID des äußeren Raums zurückgeliefert.
BESCHLEUNIGUNG Ray Tracing Kachel-Vereinigung
Die Ebenen fein vernetzten strahlenden Oberflächen (Platten) kommen oft bei den modellierten Geometrien vor. Sie bringen sehr viele strahlende Elemente hervor. Der Rechenaufwand bei der Strahlungsmodellierung ist mit dem vorgestellten Modell linear proportional der Anzahl der Strahlungskacheln. Rechenzeit kann gespart werden, falls die Feinheit der Vernetzung die nötige Auflösung in der Verteilung des integrierten Strahlungsflusses
(=irradiance) entlang der Oberfläche übertrifft. In diesem Fall wäre es ausreichend, weniger gröbere Kachel zu haben, als durch das numerische Gitter eingeprägt.
In veröffentlichten Verfahren wie zum Beispiel das DTRM- Modell von FLUENT wird auf clustering von
Strahlungskacheln hingewiesen. Dort wird die Vereinigung von mehreren benachbarten Strahlungskacheln vorgenommen, falls sie von einer Ebene nur um einen minimalen Winkel abweichen (Planarität) . Die Strahlengänge und die Strahlungsflüsse werden nicht für jeden Bestandteil des Clusters, sondern für die gesamte Vereinigung am Ort ihres geometrischen Zentrums berechnet. Eine ähnliche Technik ist für eine „billige" Beschleunigung in aktuellem Verfahren implementiert.
Alle Strahlungskacheln haben eine kartesische Orientierung im aktuellen Modell. Die Vereinigung der Nachbarn wird zugelassen, wenn diesen sowohl in einer Ebene liegen als auch ihre Gitterzellen das gleiche opake Material haben.
Es wird bis zu 3 vereinigten Kacheln in jeder der zwei lateralen Richtungen entlang der Oberfläche zugelassen. So können von 1 bis zu 9 Seiten der Gitterzellen in eine Kachel vereinigt werden. Die zugelassenen Varianten sind In Abb. 9 gezeigt.
Clustering der Strahlungskacheln kommt meistens nur an flachen Oberflächen massiv zustande. Bei gekrümmten Geometrien ist das Oberflächengitter stufig, wodurch sich die Anzahl der komplanaren Seiten der Gitterzellen verringert .
Eine Beschleunigung um Faktor 2 oder leicht darüber wird bei komplexen Geometrien durchschnittlich durch das implementierte Verfahren der Kachelvereinigung erreicht. Programm-technisch besteht die Kachelvereinigung in Cluster aus folgenden Schritten:
• Die Liste der zu einer cluster-Strahlungskachel gehörenden Gitterzellen wird in der Datenstruktur einer Strahlungskachel abgelegt.
• Die Strahlenverfolgung startet für jeden Cluster aus seinem geometrischen Zentrum.
• Informationen über die gefundenen und angepassten Strahlungsquellen werden ebenfalls für das gesamte Cluster anstatt für jeden kleinsten Bestandteil gespeichert .
· Der emittierte Strahlungsfluss des Clusters wird mit
Berücksichtigung von Temperaturverteilung berechnet:
(10). qem = e^- mit S Fläche der Seite i und T absolute Temperatur an der Oberfläche der Gitterzelle.
· Der berechnete einfallende Strahlungsfluss ist für jedes Mitglied der Clusters derselbe, was zu einer lokalen Verschmierung, Glätterung des
Strahlungsflusses führt. Dies ist aber der Preis für die Beschleunigung.
· Der gemittelte einfallende Strahlungsfluss wird auf gleiche Weise für jede Seite der Gitterzelle aus dem Cluster in die Energiebilanz (2) einbezogen.
Anisotrope Chebychev-Distanz
Die Strahlenverfolgung nach dem Voxel-basierten Ray Tracing Verfahren setzt den Besuch jeder auf dem Strahlengang liegenden transparenten Gitterzelle voraus, siehe Abb 8. Die Bearbeitung von vielen leeren (d. h. transparenten) Zellen kann einen wesentlichen Teil der Rechenzeit nehmen, besonders wenn das Modell große mit Luft gefüllten Zwischenräume bzw. die feine Vernetzung in den Zwischenräumen hat. Ein Verfahren würde das Ray
Tracing beschleunigen, dass einen Strahl statt den Besuch jeder sich unterwegs auftreffenden Gitterzelle sofort über einen größeren Block der leeren Zellen „tunnelt". Die Information über die Größe von solchen leeren Zellenblöcken soll dafür vorhanden sein.
Solche leeren Blöcke können explizit im Gitter definiert werden. Ein besseres Verfahren ist so genannte Chebychev- Distanz, auch als Schachbrett-Distanz genannt. Es geht dabei um ein Distanzmaß für diskrete Objekte wie Gitterzellen, das in ganzen Zahlen gemessen wird.
Eine Veranschaulichung der Chebychev-Distanz anhand von 2D-Beispielen ist in Abb. 10 gezeigt. Dabei wird wie folgt vorgegangen.
Die Distanz von den gestrichenen Gitterzellen (Nullzellen) wird gemessen. Abb. 10A zeigt die isotrope klassische Chebychev-Distanz. Die Zellen mit gleicher Distanz um die gestrichenen Nullzellen mit Distanz Null sind in einer Schicht am Rande eines Quadrats angeordnet. Abb. 10B zeigt die anisotrope Chebychev-Distanz für den linken oberen Quadrant gezeigt. Die anisotrope Distanz ist entsprechend ausschließlich nur für die Zellen definiert, die in diesem Quadrant in Bezug auf eine der Nullzellen liegen. Ansonsten wird sie genauso wie die isotrope Distanz in Abb. 10A gemessen. Die betrachteten Richtungen (Quadrant) sind in beiden Fällen mit roten Pfeilen in den Nullzellen gekennzeichnet. Es fällt auf, dass die anisotrope Distanz größere Werte als die isotrope erreicht.
Die Benutzung der Chebychev-Distanz lässt Schritte bei der Strahlenverfolgung adaptiv, abhängig vom aktuellen
Abstand zu einer Oberfläche gestalten. Wenn die aktuelle Position des verfolgten Strahls weit von der nächsten Oberfläche entfernt ist, wird der Sprung im Gitter auch entsprechend groß sein. Wenn der Strahl sich einer Oberfläche nähert, nimmt die Chebychev-Distanz ab, und mit ihr auch die Sprunggröße. Je näher ein Strahl der Oberfläche kommt, desto langsamer wird seine Fortbewegung . Die isotrope klassische Chebychev-Distanz wird im Bezug zu den markierten Gitterzellen mit der Distanz Null gemessen. Solche Zellen werden im Weiteren als Null- Zellen genannt. Die Zellen mit gleicher Chebychev-Distanz von einer ausgewählten Nullzelle sind in quadratischen Schichten um die Null-Zelle angeordnet.
Die Distanz erhält man durch Suche nach den noch nicht markierten Gitterzellen. Wenn am Ende des ersten Schritts des Algorithmus alle nächsten Zellen aus der ersten Schicht mit Distanz gleich 1 gefunden und markiert sind, wird die nächste Schicht mit Distanz 2 belegt. Dafür werden alle nicht markierten Zellen mit direkten Nachbarn oder Nachbarn „over edge" mit Distanz 1 gesucht. Man geht für die weiteren Schichten ähnlich vor, bis keine Zellen ohne eingetragene Distanz mehr gefunden werden.
Der Aufwand für eine direkte Berechnung der Chebychev- Distanz verhält sich entsprechend als wo N Anzahl der Gitterzellen ist.
Diese Abhängigkeit ist aber nicht zwingend und lässt sich optimieren, wenn bei der Belegung jeder nächsten Schicht der Zellen mit Distanz nur die zuvor gefundenen und zwischengespeicherten Gitterzellen mit Distanz i
prozessiert werden. Es werden nur ihre Nachbarzellen untersucht. So wird der Front der zuletzt gefundenen „aktiven" Zellen propagiert. Das Verfahren setzt aber voraus, dass der Schreibzugriff in die äußeren Zellen des Stencils möglich ist. Deswegen setzt es Beschränkungen für die Parallelisierung durch die Gebietszerlegung.
Das Vorgehen bei der Berechnung der anisotropen Chebychev-Distanz ist sehr ähnlich. Eine Distanz für eine Zelle im isotropen Fall wird jetzt aber durch 22=4 Distanzen in einem 2D Fall und 23=8 Distanzen in einem 3D Fall ersetzt. Jede der 4 Distanzen wird im 2D Fall nur für Gitterzellen berechnet, die im entsprechenden Quadrant, von Null-Zellen aus gesehen, liegt, siehe Abbildung 10B. Die Belegung der Distanzen dort wird für den Quadrant zwischen den kartesischen Richtungen -X und +Y berechnet. Im 3D wird zwischen den 8 Oktanten unterschieden. Es findet in beiden Fällen im Gegenzug zur isotropen Chebychev-Distanz keine Belegung der Distanzen in den Zellen statt, die im gegebenen Quadrant (2D) oder Oktant (3D) keine Null-Zellen haben, von denen die Distanz gemessen wird.
Die Bestimmung eines Oktanten ist der Benutzung der anisotropen Distanz in der Strahlenverfolgung vorausgesetzt. Es wird zunächst geprüft, in welchem
Oktant die Strahlenrichtung Ω liegt. Dann wird auf eine der 8 anisotropen Distanzen im 3D Fall zugegriffen, die einem gegenüberliegenden Oktant entspricht. Wenn z. B. ein Strahl aus der Richtung (0.5, 0.2, -0.1) kommt, liegt er im Oktant (+X, +Y, -Z) . Die entsprechende Distanz, von den Null-Zellen ausgehend wurde in einem Oktant entgegen der Strahlenrichtung gemessen. Die nötige Distanz ist also im Oktant (-X, -Y, +Z) gespeichert.
Der Vorteil einer anisotropen Distanz über einer isotropen ist eine genauere Vorgabe für einen möglichen Sprung über einen Block der leeren Zellen in der Strahlenverfolgung. Der Wert der anisotropen Distanz in einem Oktant fällt oft größer aus, als bei der isotropen, vergleiche Abb. 10A und 10B. Die Wahrscheinlichkeit, dass eine Nullzelle in einem Oktant-Würfel des Gitters liegt, ist kleiner, als im isotropen Falle, wo der entsprechende Würfel aus 8 solchen Oktanten zusammengesetzt ist. Daher kann die Strahlenverfolgung größere Sprünge im Gitter machen .
Ein Nachteil der anisotropen Distanz ist der dafür zusätzlich erforderliche Speicher bei der
Strahlenverfolgung, was extra 96 Bytes je eine Gitterzelle bei dem Typ „int" der Distanz ausmacht. Um Speicher zu sparen, werden Distanzen mit dem Typ „unsigned char" in diesem Verfahren belegt. Dieser Typ lässt eine maximale Größe der Distanz von 128 Gitterzellen zu. Theoretisch können größere Werte in einer Anwendung vorkommen. Deswegen wird die Belegung der Distanzen bei ihrer Berechnung nach dem Erreichen der Schicht mit Wert 128 unterbrochen, und die verbliebenen nicht belegten transparenten Zellen mit demselben Wert von 128 belegt. Dadurch wird die maximale mögliche Länge eines Sprungs im Gitter bei der Strahlenverfolgung beschränkt . Als Nullzellen werden im Verfahren die opaken Zellen mit definierten Strahlungskacheln und die Zellen am äußeren Rand des numerischen Gitters im Kontakt mit einer transparenten Zelle vorgegeben. Durch das letzte wird sichergestellt, dass alle transparenten Zellen lückenlos
von einer Schicht der Nullzellen umschlossen sind, und für jede transparente Zelle mindestens eine Referenz- Nullzelle auch im anisotropen Fall vorliegt. Ein Beispiel der Strahlenverfolgung in 2D, unterstützt durch anisotrope Chebychev-Distanz , ist in Abbildung 11 gezeigt. Das Gitter und die berechnete Distanz sind aus der Abb. 10 rechts übernommen. Der Strahl im Oktant (+X,- Y) wird mit Distanz im Oktant (-X,+Y) bearbeitet. Die Startposition des Strahls im Schnittpunkt mit einer Zellenseite oben links liegt in der Zelle mit Distanz 3. Deswegen wird der nächste Schnittpunkt mit der Seite des 3x3 Würfels, gesucht. Nach dem Verschieben der Position zu dem Schnittpunkt wird die Chebychev-Distanz aus nächster Zelle in Richtung des Strahlenganges genommen. Sie ist wieder gleich 3, und der nächste Sprung führt den Strahl direkt zum Ziel unten rechts. In diesem Beispiel waren 2 Schnittpunkte gesucht; ohne Anwendung der Chebychev-Distanz wären das ganze 9.
Die Mehrkosten bei der Strahlenverfolgung mit Chebychev- Distanz sind in der Bestimmung der diskreten Position des Strahls nach dem Sprung über einen Block der Gitterzellen. Die Koordinaten des Schnittpunktes lassen sich ähnlich nach den Gleichungen (9) berechnen. Anstelle der kartesischen Koordinaten am Rande der Gitterzelle XL kommt jedoch die Koordinate am Rande des Oktanten- Würfels, die der Koordinate einer weiter entfernten Gitterzelle gleich ist.
Die Indizes der nach dem Sprung getroffenen Zelle in 3 kartesischen Richtungen sind unbekannt. Sie werden nach dem binären Suchalgorithmus bestimmt. Der Index im Gitter am Schnittpunkt wird ausgehend aus einem äquidistanten
Gitter im Würfel in jeder Richtung abgeschätzt. Dann wird durch den Koordinatenvergleich geprüft, ob der gefundene Schnittpunkt in der Zelle mit dem geschätzten Index liegt. Wenn nein, wird die entsprechende Hälfte des Intervalls mit Länge der Chebychev-Distanz weiter halbiert. Die Prozedur wird auf diese Weise rekursiv wiederholt, bis der gefundene Index der Zelle dem Schnittpunkt gehört. Der Aufwand der binären Suche in einem festen Intervall hängt logarithmisch von der Intervalllänge ab, d. h. für eine maximale Sprunglänge von 128 Zellen würde ~ln(128) = 81n(2) ~ 5 Vergleichsoperationen je Richtung nötig.
Trotz dem Mehraufwand für einen Gittersprung liefert die Chebychev-Distanz eine beträchtliche Beschleunigung des Verfahrens, besonders für Modelle, wo viel Raum mit transparentem Material belegt ist.
Paralleles Ray Tracing
Lastenverteilung bei Ray Tracing:
Eine gute Skalierbarkeit im parallelen Betrieb ist die wichtige Voraussetzung der Effizienz bei der Initialisierung des Modells (Chebychev-Distanz, Ray Tracing) und ihrer Anwendung in thermischer Berechnung. Entscheidend ist dabei die Rechenzeit für Ray Tracing.
Die gängige Praxis für die Parallelisierung mit verteiltem Speicher wie MPI bei der Lösung von partiellen Differentialgleichungen ist die Zerlegung des Rechengebiets. Ihre Anwendung für Ray Tracing hat aber eine schlechte Skalierbarkeit in dieser Entwicklung gezeigt. Der Grund dafür sind hohe Kosten der
Kommunikation. Daten der Millionen von Strahlen müssen immer wieder zwischen den einzelnen Partitionen übertragen werden. Das spezifische Gewicht der Kommunikation wächst mit steigender Anzahl der CPUs, was schließlich zu einer schlechten Skalierbarkeit führt. Deswegen ist in diesem Modell ein anderes Verfahren der Parallelisierung entwickelt.
Ray Tracing und die Berechnung der Chebychev-Distanz werden auf jedem CPU für das komplette Gitter gleichzeitig durchgeführt. Die Strahlungskacheln bleiben aber dabei auf jedem CPU lokal definiert.
Das globale Gittermodell wird im ersten Schnitt erstellt. Die Gitterdaten von dem gesamten Gitter werden zu allen CPUs kommuniziert. Einerseits, werden drei Reihen der Koordinaten in drei kartesischen Richtungen mitgeteilt. Andererseits wird ein Flag für die Transparenz jeder Gitterzelle von Root ermittelt und in einem numerischen Feld vom Typ „unsigned char" an alle CPUs kommuniziert. Diese Daten reichen aus, um die Durchführung von Ray Tracing auf dem Tensorprodukt-Gitter unabhängig voneinander auf jedem CPU durchzuführen. Die Symmetriegrenzen im thermischen Modell werden an dieser Stelle auch berücksichtigt.
Die durch Ray Tracing lokal ermittelten Daten über die Strahlungsquellen der Strahlungskacheln sollen aber schließlich zu dem CPU gelangen, wo die entsprechende Kachel lokal definiert ist. An dieser Stelle entstehen Kosten für die Kommunikation. Um sie zu minimieren, wird wie folgt vorgegangen.
1. Die Anzahl der vorhandenen auf einem CPU lokal definierten Strahlungskacheln N wird ermittelt.
2. Dann wird der arithmetische Mittelwert Nav über die CPUs berechnet.
3. Von einem CPU-Donator mit dem Überschuss der Kacheln
> Nav wird eine Portion AN der Kacheln von dem Donator virtuell abgezogen und an das nächste gefundene CPU-Akzeptor mit N2 < Nav zugeschrieben, so dass nach dieser Operation entweder die Bedingung N1 = Nav oder N2 = Nav erfüllt wird. Die übertragene
Portion ist dabei gleich AN = M1N(N1— Nav, N2— Nav) . Die Indizes der zu übertragenden Kacheln und die IDxs der CPUs-Akzeptoren werden jedes Mal gespeichert. 4. Der Schritt 3 wird wiederholt, bis kein Ausgleichen der Kacheln zwischen der CPUs mehr möglich wird.
Nachdem die Parameter für den Lastenausgleich wie oben beschrieben jedem CPU bekannt werden, werden die lokalen Daten zur Lage der Strahlungskacheln im Gitter von den CPUs-Donatoren an die CPUs-Akzeptoren kommuniziert. Die CPUs-Akzeptoren erhalten auf diese Weise eine Liste der zu bearbeitenden importierten Kacheln zusätzlich zur Liste der eigenen Kacheln. Die Listen der zu bearbeitenden eigenen Kacheln auf den CPUs-Donatoren werden entsprechend gekürzt.
Paralleles Ray Tracing:
1. Vor dem Ray Tracing erfolgt die Berechnung der anisotropen Chebychev-Distanz auf jedem CPU unabhängig voneinander. Keine Kommunikation ist hier erforderlich, Berechnungen und Ergebnisse auf jedem CPU sind an dieser Stelle identisch. Der Aufwand für
diese Berechnung beträgt wenige Sekunden, so dass die Organisation einer Arbeitsverteilung zwischen den einzelnen CPUs und anschließende Kommunikation hier kaum lohnen würde.
Während des Ray Tracing werden auf jedem CPU unabhängig voneinander alle Strahlen ausgehend von zuerst eigenen zu bearbeitenden Kacheln verschickt, die Strahlungsquellen gefunden und nach ihrer geometrischer und thermischer Anpassung die ausgewählten Quellen gespeichert. Die Speicherung der Strahlungsquellen erfolgt in diesem Schritt direkt in die Datenstruktur der eigenen Strahlungskacheln .
Nachdem alle eigenen Kacheln bearbeitet sind, sind die importierten Kacheln auf den CPUs-Akzeptoren an der Reihe. Die gleiche Prozedur der Strahlenverfolgung und der Ermittlung der Strahlenquellen wird bei der Bearbeitung von diesen Kacheln verwendet. Die gefundenen und durch die Prozedur der Anpassung komprimierten Strahlungsquellen werden in einem Buffer zwischengespeichert .
Wenn alle Kacheln bearbeitet sind, erfolgt die Kommunikation der Ergebnisse. Die CPUs-Akzeptoren senden die gefundenen Strahlungsquellen der importierten Kacheln an die CPUs-Donatoren zurück. Die MPI-Kommunikation erfolgt asynchron, wobei bei einer unregelmäßigen Verteilung der Kacheln sowohl jeder Donator mehrere Akzeptoren, als auch jeder Akzeptor mehrere Donatoren haben darf. Die vom Donator empfangenen Strahlungsquellen werden in diesem Schritt gleich nach dem Empfang von einem CPU-Akzeptor in die Datenstruktur der eigenen Kacheln geschrieben.
5. Am Ende von Ray Tracing werden nicht mehr benötigten Daten wie Chebychev-Distanzen und Zwischenspeicher für die Kommunikation gelöscht., THERMISCHE BERECHNUNG
Um die Summierung der ausgehenden Strahlungsflüsse nach dem Mehrlevel-Verfahren nach Gleichung (8) zu lösen wird ein passender Algorithmus zur näherungsweisen Lösung von linearen Gleichungssystemen eingesetzt. Die Summierung ist unten anhand der Gauß-Seidel-Iteration exemplifiziert, könne aber auch mittels anderen Algorithmen, beispielsweise das Jacobi-Verfahren oder das SOR-Verfahren, gelöst werden.
In einer Ausführungsform wird die Gauß-Seidel-Iteration wird für die Berechnung des absorbierten
Strahlungsflusses an jeder Strahlungskachel benutzt. Dieser Fluss stellt das Endergebnis für das thermische Modell von der Strahlungsseite dar. Das Wissen der ausgehenden Strahlungsflüsse wird an dieser Stelle CPU- übergreifend erforderlich, da der berücksichtigte Strahlengang für jeden Strahl ohne Rücksicht auf die Gebietszerlegung berechnet wurde.
Es ist vom Vorteil dass die Gauß-Seidel-Iteration auf jedem CPU nur für die lokalen Strahlungskacheln erfolgt. Die Summierung der ausgehenden Strahlungsflüsse nach dem Mehrlevel-Verfahren nach Gleichung (8) beinhaltet die Beiträge von allen durch Ray Tracing ermittelten Kacheln, auch von anderen CPUs .
Dadurch können die ausgehenden Strahlungsflüsse im Algorithmus global sichtbar gemacht werden. Jeder
Strahlungskachel wird ein durchgehender globaler Index zugewiesen, wobei vorteilhaft dieselben Indizes anstelle der Strahlungsquellen für jede Kachel während Ray Tracing gespeichert werden können.
Die globalen ausgehenden Strahlungsflüsse können in einem Array nach dem globalen Index indiziert werden, wobei vor jedem Durchlauf der Gauß-Seidel-Schleife die ausgehenden Flüsse von allen lokalen Kacheln in das globale Array kommuniziert werden.
POTENZIELLE ANWENDUNGSGEBIETE, BEISPIELE
Das oben vorgestellte Modell der thermisch gekoppelten Oberflächenstrahlung für diffuse graue Strahler kann die Anwendung in vielen technischen Verfahren finden, wo Hohlräume oder halbtransparente Materialien und hohe Temperaturen vorliegen. Die attraktive Seite des Verfahrens liegt in der Leichtigkeit und sehr kurzer Rechenzeit, in der die anspruchsvollsten Aufgaben in der Modellierung der Fernwirkung der thermischen Strahlung bearbeitet werden.
Die interessantesten Anwendungen in der Simulation der Gießprozesse sind der Feinguss (Investment casting) , Wärmebehandlung, Blockguss bei den angereihten Anordnungen der Anlagen.
Weitergehende Anwendungen hängen mit der Produktion und Betrieb der Industrieöfen zusammen, Beispielsweise bei der Glasfertigung oder Keramikproduktion, in die Chemie, der Kristallzüchtung, der Metallurgie oder auch in die Lebensmittelindustrie, z.B. Backöfen.
In den Abbildungen 12 bis 16 sind verschiedene Beispiele der Anwendung des Verfahrens gezeigt. Wie gezeigt ist das hiesige Verfahren nicht auf nur eine Wärmequelle beschränkt sondern lässt sich auch für eine Vielzahl, z.B. zwei (Abb. 14) und eine höhere Anzahl von Wärmequellen verwenden.
Die hier angegebene Verfahren zur Diskretisierung eines Raumwinkels zur Anwendung in einen Simulationsprozess ist in seiner Anwendung nicht nur in Zusammenhang mit thermisch gekoppelter Oberflächenstrahlung zu denken, sondern kann auch in andere Berechnungs- und Simulationsprozesse eingesetzt werden wo Beschleunigung und Einsparnisse an Computerzeit und verbrauchte Computermemory, wie z.B. Computerdarstellungen (rendering Software) in Allgemeinen und in Computerspiele ins Besonders oder wissenschaftliche Darstellungswerkzeuge, von Bedeutung ist. Die hier angegebene Verfahren des Ray Tracings zur Anwendung in einem Simulationsprozess ist gleichweise in seiner Anwendung nicht nur in Zusammenhang mit thermisch gekoppelter Oberflächenstrahlung zu denken, sondern kann auch in andere Berechnungs- und Simulationsprozesse eingesetzt werden wo Beschleunigung und Einsparnisse an Computerzeit und verbrauchte Computermemory, wie z.B. Computerdarstellungen (rendering Software) in Allgemeinen und in Computerspiele ins Besonders oder wissenschaftliche Darstellungswerkzeuge, von Bedeutung ist. Besonders die hier vorgeschlagene unkonventionelle Verfahren des Parallel Computing zur Beschleunigung einer Ray Tracing Berechnung alleine und in seiner Kombination mit Anisotrope Chebychev-Distanz Berechnungen und/oder
zusätzlicher Akzeleration durch Kachel-Vereinigung ist dem Erfinder aus dem Stand der Technik nicht bekannt.