4. Magdeburger Brand- und Explosionsschutztag /3. vfdb Workshop Brandschutzforschung

Magdeburg, 26./27. Marz 2015

Bug oder Feature? Uber numerische Lichtschutzfaktoren imfvDOM-Strahlungsmodell

Michael Plagge1,2, Matthias Munch3,4

1 Anlagentechnik und Anlagensicherheit, FVST, Otto-von-Guericke-Universitat Magdeburg2 CERN, Europaische Organisation fur Kernforschung, Genf, Schweiz3 INURI GmbH, Haderslebener Strasse 9, Berlin4 Numerische Mathematik, FB Mathematik und Informatik, FU-Berlin

Kurzfassung:Bei der Anwendung von CFD-Simulationsprogrammen (CFD: Computational Fluid Dynamics)im Brandingenieurwesen ist eine wesentliche Teilaufgabe die Berechnung der Auswirkung vonVerbrennungsprozessen. Die Simulationsprogramme verfugen daher in der Regel uber ein odermehrere Submodelle fur die Berechnung des Warmetransportes durch die Warmestrahlung. Furdie physikalische Abbildung des Warmeaustausches im Brandfall ist die Warmestrahlung einwichtiger Transportprozess, da der Warmestrahlungstransport auch entgegen der Stromungs-richtung erfolgen und damit den Verlauf und die Ausbreitung entscheidend beeinflussen kann.

Aus betriebswirtschaftlicher Sicht ist hingegen die von den Warmestrahlungsmodellen be-notigte Rechenkapazitat von großem Interesse. Je nach Modellansatz und gewahlter Diskreti-sierung benotigen die Warmestrahlungssubmodelle erhebliche Rechenkapazitaten, die leicht 20Prozent der insgesamt erforderlichen Rechenkapazitat betragen konnen.

Der vorliegende Beitrag berichtet uber Teilergebnisse einer Zusammenarbeit von CERN undINURI, in dem die Einflussfaktoren von Warmestrahlungsmodellen im Rahmen von Brandsimu-lationen untersucht werden. Anhand eines analytischen Testfalls werden im vorliegenden Bei-trag die Ergebnisse einer numerischen Verifikationsstudie des fur finite Volumina entwickeltenDiskrete-Ordinaten-Modells (fvDOM) diskutiert. Dieser Modellansatz kommt in zahlreichen inder Praxis eingesetzten CFD-Programmen zum Einsatz. Unter anderem ist dieser in den Pro-grammcodes ANSYS FLUENT, dem Fire Dynamics Simulator (FDS) des National Institute ofStandards and Technology (NIST) und der Softwarebibliothek OpenFOAM der OpenFOAM-Foundation enthalten. Neben der raumlichen Diskretisierung des Stromungslosers bestimmtsich der Rechenaufwand des fvDOM-Modells maßgeblich uber die Anzahl der Raumwinkel, diezur diskreten Auflosung der Strahlenverfolgung verwendet werden. Hieraus ergibt sich fur denAnwender ein Optimierungsproblem zwischen der diskreten Auflosung und dem noch akzepta-blen Rechenaufwand.

Plagge, Munch

1 Einleitung

In [4] wird ein praxisnahes Simulationsexperiment zur Untersuchung des Einflusses der Git-terweite auf die Modellierung der Brandausbreitung durch die Warmestrahlung vorgestellt. Ineiner Konfiguration zweier planparalleler Platten besitzt eine der Platten eine konstant hoheTemperatur von 1000 Grad Celsius. Die von dieser Strahlerplatte ausgehende Warmestrahlungtrifft auf die gegenuberliegende Platte und entzundet diese schließlich, wenn die in dem Si-mulationsmodell hierfur definierte Oberflachentemperatur auf der jeweiligen Plattenoberflacheerreicht wird. In diesem Simulationsexperiment veranderten sich die ermittelten Entzundungs-zeiten nicht in Analogie mit der variierten Gitterweite, sondern zeigten starke Schwankun-gen. Fur das in diesem Simulationsexperiment verwendete Simulationsprogramm Fire Dyna-mics Simulator (FDS) des National Institute of Standards and Technology (NIST) [5] hat-ten die Programmentwickler einen wesentlichen Modellparameter des verwendeten fvDOM-Strahlungsmodells, die Anzahl der diskreten Strahlungswinkel, mittels eines analytischen Test-falls fest definiert.

Im Rahmen des Projekts wird untersucht, wie beziehungsweise ob dieser Parameter furBrandsimulationen fest definiert werden kann und welchen Einfluss dies hat. Da das betref-fende fvDOM-Strahlungsmodell auch in anderen Simulationsprogrammen zur Verfugung steht,werden zum Vergleich auch die Simulationsprogramme OpenFOAM [6] und ANSYS Fluent[1] mitbetrachtet.

2 Modellierung der Warmestrahlung

Das Diskrete-Ordinaten-Modell (DOM) wurde 1960 von Chandrasekhar ursprunglich furdie Methode der finiten Differenzen entwickelt. Es folgten mehrere Verallgemeinerungen undOptimierungen, so dass es erst seit rund 30 Jahren verbreitet Anwendung bei der Losung vonWarmetransportproblemen findet. Diese und weitere historische Fakten die Entwicklung be-treffend, finden sich in Modest ”Radiative heat transfer “, Kapitel 17, [3]. In diesem Buchfinden sich auch eine ausfuhrliche Beschreibung des Modells nebst mehreren Beispielen undLosungen. Auf eine Wiederholung der Theorie und der unterschiedlichen Implementierungenin die drei anfangs genannten Programme wird daher an dieser Stelle weitgehend verzichtet.

Die Modellierung der Warmestrahlung im Zusammenhang mit dem fvDOM-Modellansatzlasst sich anhand der Strahlungshalbkugel nach Siegel [7] erklaren. Wie in Abbildung 1 gezeigt,wird um eine ebene strahlende Grundflache dA eine Halbkugel mit dem Radius R angenom-men.

R

ϑ

dω

ϕ

dA

dAk

Abbildung 1: Strahlungshalbkugel nach Siegel [7]

Ausgehend vom Kreismittelpunkt, der gleichzeitig der Mittelpunkt des Flachenstucks dAist, wird ein beliebiger Strahl s in Richtung der Halbkugel betrachtet. Die Lage dieses Strahl

2

Plagge, Munch

s laßt sich durch den Azimutwinkel ϕ und den Polarwinkel ϑ genau definieren. Kleine Aus-lenkungen von Azimut- und Polarwinkel um dϑ und dϕ ergeben auf der Halbkugeloberflacheein Flachenelement senkrecht zum Strahl s. Auf diese Weise lassen sich ausgehend vom Halb-kugelmittelpunkt kegelformige Raumausschnitte bilden. Diese Raumausschnitte werden durchden Raumwinkel dω reprasentiert. Der Raumwinkel dω sei im Folgenden das vom Strahl s be-einflusste Raumvolumen. Je kleiner der Raumwinkel dω, desto feiner die Auflosung der Halb-kugel und umgekehrt.

3 Der Verifikationstestfall

Als Verifikationstestfall wird der Plate-View-Factor Testfall aus dem FDS-Verifikationguideverwendet (vgl. Abbildung 2). Er besteht aus einem quaderformigen Simulationsgebiet, dessenWandrandbedingungen so gestaltet sind, dass die Wandflachen auf 0 Grad Kelvin (blau) gehal-ten werden. Einzig die Deckenflache besitzt eine vom absoluten Nullpunkt abweichende Tem-peratur und wird mit 1273,14 Grad Kelvin (rot) festgelegt. Unter der Annahme schwarzer Strah-lereigenschaften werden die Emissivitatskoeffizienten der Flachen auf 1 gesetzt. Gravitations-und konvektive Warmeubergangseffekte auf den Begrenzungsflachen werden per Konstruktionunterbunden. Diese Konstruktion erzwingt eine Schwarzkorper-Strahlung von der heißen De-cke auf die anderen Berandungsflachen. Da diese auf 0 Grad Kelvin gehalten werden, sendetallein die Decke Warmestrahlung aus.

Abbildung 2: Modifizierter Plate View Factor-Testfall [4].

Der Testfall wurde nach [4] modifiziert. Das Simulationsgebiet ist 2,1m lang, 2,1m breit und1,2m hoch und wird mit einem kartesischen Gitter aquidistanter Gitterweite diskretisiert. Diegroßte Gitterweite betragt 0,3 m, so dass mittig am Boden die in Abbildung 2 dargestellte Refe-renzmessflache (orange) gebildet werden kann. Die Berandung dieser Flache entspricht geradeden Gitterlinien der Diskretisierung. Alle feineren Diskretisierungen sind ein Vielfaches diesesGrundgitters, so dass die Berandung der Referenzflache stets auf den Gitterlinien liegt. Diesermoglicht Diskretisierungen mit aquidistanten Gitterweiten von 1, 2, 6, 10, 15 und 30 cm, diejeweils bundig ineinander passen.

Fur die Untersuchung der Gitterabhangigkeit des Strahlungslosers wird die einfallende War-mestrahlung RADIATIVE HEAT FLUX (RHF) auf die 0,3m x 0,3m große Referenzflache (oran-ge) betrachtet (vgl. Abbildung 2).

Die Simulationszeit wird 0,01 Sekunden begrenzt und ein Zeitintervall von 0,001 Sekundendefiniert. Auf diese Weise lasst sich uberprufen, ob und wie die numerischen Implementierun-gen der betrachteten Warmestrahlungsmodelle eine stationare Losung erreichen.

3

Plagge, Munch

3.1 Analytische Losung

Die auf die gesamte Bodenflache von der Decke einfallende Warmestrahlungsleistung (Ra-diative Heat Flux (RHF)) ergibt sich fur die beschriebene Testfallkonfiguration zu

q′′RHF = σ 1273,144 = 148,97 kW/m2. (1)

Die analytische Losung fur die auf die kleinere Referenzflache (orange) einfallende Warme-strahlungsleistung q′′RHF,F lache laßt sich auf der Grundlage der optischen Winkelrelationen mitHilfe der Methode der Sichtfaktoren berechnen.

h

a

bA2

A1

Abbildung 3: Parallel quadratische Flachenanordnung zur Bestimmung des Sichtfaktors [4].

Hierbei ergibt sich die von der Decke A1 auf die Referenzbodenflache A2 einfallende Strah-lungsleistung durch Multiplikation eines Sichtfaktors F1−2 mit der von der Decke abgestrahltenWarmestrahlung q′′RHF unter Berucksichtigung der Flachenverhaltnisse.

q′′RHF,F lache = F1−2 q′′RHF

A1

A2

. (2)

Der Sichtfaktor F1−2 fur begrenzte parallele quadratische Flachen berechnet sich nach [2]mit A = a/h und B = b/a zu:

fur A ≥ 0,2 : (3)

F1−2 =1

π A2

{ln

[A2(1 +B2) + 2]2

(Y 2 + 2)(X2 + 2)

+√Y 2 + 4

[Y tan−1 Y√

Y 2 + 4−X tan−1 X√

Y 2 + 4

]+√X2 + 4

[X tan−1 X√

X2 + 4− Y tan−1 Y√

X2 + 4

]}mit X = A (1 +B) und Y = A (1−B)

fur A < 0,2 : (4)

F1−2 =(AB)2

π.

Mit Hilfe der vorgenannten Berechnungsvorschriften ergibt sich die analytische Losung furdie mittlere RHF-Strahlung auf die 0,3m x 0,3m große Referenzflache mit q′′RHF,F lache =72,34 kW/m2.

4

Plagge, Munch

3.2 Simulationsvariationen

Der beschriebene Testfall wird mit der jeweiligen Version des in den Simulationsprogram-men FDS, OpenFOAM und ANSYS FLUENT verfugbaren diskreten Ordinaten Warmestrahl-ungsmodells berechnet. Hierbei wird in der Darstellung der Ergebnisse die Annahme der FDS-Entwickler ubernommen, dass die Anzahl der Strahlungswinkel universell fur alle Anwendun-gen voreingestellt werden kann. Das Simulationsergebnis wird daher in Abhangigkeit von derGitterweite fur jeweils konstante Winkelzahlen dargestellt.

Das beschriebene Simulationsgebiet wird mit aquidistanten Gitterweiten von 1, 2, 6, 10, 15und 30 cm, die jeweils bundig ineinander passen, diskretisiert. Die Anzahl der insgesamt vomWarmestrahlungsmodell verwendeten Raumwinkel wird mit n = 50, 100, 200, 300, 400 und500 variiert, soweit dies die Simulationsprogramme zulassen.

3.3 Ergebnisse mit FDS

Die Simulationen wurden mit dem Fire Dynamics Simulator (FDS) in der vom NIST zurVerfugung gestellten Binarversion FDS 6.1.2 Revision 20564 gerechnet. Soweit nicht anders an-gegeben, wurden die Defaulteinstellungen des Programms verwendet. Die Winkelzahlen konn-ten uber einen entsprechenden Programmparameter vorgegeben werden. Der Vergleich der Si-mulationsergebnisse mit den analytischen Losungen zeigt fur den Vergleich mit der 0,3m x0,3m großen Referenzflache (vgl. Abbildung 4) ein stark gitterabhangiges Losungsverhaltenfur den Warmestrahlungstransport bei festgehaltenem Strahlungswinkel.

60

62

64

66

68

70

72

74

76

0 5 10 15 20 25 30

Rad

iativ

e H

eat F

lux

in k

W/m

2

Gitterweite in cm

Analytische Lösung 50 Winkel100 Winkel200 Winkel300 Winkel400 Winkel500 Winkel

Abbildung 4: Einfallende mittlere RHF-Strahlung auf die 0,3m x 0,3m große Referenzflache in Abhangigkeit vonGitterweite und Strahlungswinkel (FDS-Version 6.1.2).

3.4 Ergebnisse mit OpenFOAM

Im Gegensatz zum Fire Dynamics Simulator muss in der Softwarebibliothek OpenFOAMRevision 2.2.2 die Anzahl der Strahlungswinkel jeweils fur den Azimut- und Polarwinkel ein-zeln vorgegeben werden. Dies fuhrt zu etwas anderen Gesamtwinkelzahlen, die jedoch so ge-wahlt wurden, dass die Großenordnung vergleichbar ist. Die Abbildung 5 zeigt fur den Ver-gleich der Simulationsberechnung des fireFoam-Paketes mit der analytischen Losung fur die0,3m x 0,3m große Referenzflache.

5

Plagge, Munch

71

71.5

72

72.5

73

73.5

74

74.5

75

75.5

76

0 5 10 15 20 25 30

Rad

iativ

e H

eat F

lux

in k

W/m

2

Gitterweite in cm

Analytische Lösung 50 Winkel 98 Winkel200 Winkel288 Winkel392 Winkel512 Winkel

Abbildung 5: Einfallende mittlere RHF-Strahlung auf die 0,3m x 0,3m große Referenzflache in Abhangigkeit vonGitterweite und Strahlungswinkel (OpenFOAM-Version 2.2.2).

Auch hier findet sich das gleiche, schon beim Fire Dynamics Simulator zu beobachtendeSimulationsergebnis.

3.5 Ergebnisse mit ANSYS Fluent

Im kommerziellen Programm ANSYS Fluent, Version 15.0, ist ebenfalls eine Implemen-tierung des Discrete-Ordinate-Model enthalten. Wenn nicht anders angegeben, wurden fur dieModellparameter die Standardeinstellungen verwendet. Alle Rechengitter konnten via Kom-mandozeile von der OpenFOAM-Variante ubernommen werden.

Abbildung 6: Beispielrechengitter dargestellt in ANSYS Fluent (hier Zellen mit mit je 30 cm Kantenlange in x-,y- und z-Richtung.

Die Simulationen wurden seriell instationar (transient) unter Verwendung des Druck-basiertenLosers (Pressure-Based) mit absoluten Geschwindigkeiten (Absolute Velocity Formulation) ge-

6

Plagge, Munch

rechnet. Aus der Liste der verfugbaren Modelle wurde nur das Warmestrahlungsmodell (Radia-tion) aktiviert. Hierbei wurde das Discrete-Ordinate Submodell ausgewahlt. Die Diskretisierungder Raumwinkel in ANSYS Fluent erfolgt gemaß Gl. 5

nRaumwinkel = 8 nΦ nΘ . (5)

Hier unterscheidet sich ANSYS Fluent von OpenFOAM. Im Ergebnis sind auch hier die imTestfall vorgesehenen Winkelauflosungen nur naherungsweise erreichbar. Tatsachlich wurdendaher 48, 96, 200, 288, 400 und 504 Raumwinkel verwendet.

ANSYS Fluent definiert automatisch die feste Phase Aluminium (Solid - aluminum) und dieGasphase Luft (Fluid - air). Bei beiden wurden die Standardwerte ubernommen. Insbesonderebei der Luft wurde der Absorptionskoeffizient gepruft. Dieser ist standardmaßig Null, was imSinne des Vergleichs mit der analytischen Losung ist.

Das zu berechnende Volumen wurde mit der definierten Luft gefullt (Cell Zone Conditions -Type fluid). Die Randbedingungen fur die sechs Aussenflachen sowie die Zielflache im Bodenwurden analog zu den beiden anderen Programmen definiert. Fur die Geschwindigkeit wurdenalle Flachen als stationare Wand (Stationary Wall) ohne Wanddurchstromung (No slip) definiert.Die Temperatur fur die Decke wurde auf 1273,14 Grad Kelvin gesetzt, die der anderen Flachenauf 0 Grad Kelvin. Der Emissionsgrad wurde durch das Programm auf eins gesetzt. Fur dasWarmestrahlungsmodell setzt ANSYS Fluent alle Flachen automatisch als lichtdicht (opaque)mit einem diffusivem Anteil von eins.

Die Erfassung der einfallenden Warmestrahlung wurde mit Hilfe virtueller Messgroßen rea-lisiert (Area-Weighted Average - Wall Fluxes - Radiation Heat Flux). Fur die Berechnung wurdedie konstante Zeitschrittweite von 0,001 s festgelegt. Die Berechnungen erfolgten fur jeweilszehn Zeitschritte, i.e. uber zehn Iterationen. Die Abbildung 7 zeigt das Ergebnis der 36 Simula-tionen in der Ubersicht.

66

68

70

72

74

76

78

80

0 5 10 15 20 25 30

Rad

iativ

e H

eat F

lux

in k

W/m

2

Gitterweite in cm

Analytische Lösung 48 Winkel 96 Winkel200 Winkel288 Winkel400 Winkel504 Winkel

Abbildung 7: Einfallende mittlere RHF-Strahlung auf die 0,3m x 0,3m große Referenzflache in Abhangigkeit vonGitterweite und Strahlungswinkel (ANSYS Fluent Version 15.0).

Abbildung 8 zeigt exemplarisch den auf die Decke und Bodenflache der betrachteten Geo-metrie abgegebenen bzw. einfallenden Warmestrom.

7

Plagge, Munch

Abbildung 8: Beispielergebnis bei 15 cm Gitterweite und 200 Raumwinkeln nach 10 Iterationen.

4 Zusammenfassung

Die vorgestellten Resultate sind ”Work-in-progress“. Alle betrachteten diskrete Ordinaten-modelle der verschiedenen Simulationsprogramme zeigen unabhangig voneinander das gleicheVerhalten. Bis auf wenige Ausnahmen nahert sich die auf die Referenzflache einfallende mittle-re Warmestrahlung mit zunehmend kleiner werdender Gitterweite dem analytisch berechnetenWert an und durchbricht diesen dann.

Interessanterweise zeigt das in ANSYS Fluent verwendete Modell fur eine Gitterweite von10 cm fur eine Gesamtwinkelzahl von 100 Winkeln eine sehr gute Wiedergabe des analytischberechneten Wertes. Dieses Verhalten ist auch beim Fire Dynamics Simulator zu beobachtenund veranlasste die FDS-Entwickler zur Festlegung einer universellen Gesamtwinkelanzahl von100 Winkeln.

Die Abbildungen 4, 5 und 7 zeigen bereits, dass die universelle Festlegung einer fixen Win-kelzahl unabhangig von der Gitterweite keine so gute Idee zu sein scheint. Hier ist vielmehr einabgestimmtes Konzept erforderlich, dass die Winkelanzahl in Relation zur Gitterweite betrach-tet.

LITERATUR

[1] ANSYS, Inc.: ANSYS Fluent, Version 15.0.

[2] John R. Howell. A catalog of radiation heat transfer configuration factors, 2010. Thecatalog is a resource for Thermal Radiation Heat Transfer, 5th ed., John R. Howell, RobertSiegel and M. Pinar Menguc, Taylor and Francis/CRC, New York, 2010.

[3] Michael F. Modest. Radiative Heat Transfer, 3. Auflage. Academic Press, 2013.

[4] Matthias Munch. Konzept zur Absicherung von CFD-Simulationen im Brandschutz undin der Gefahrenabwehr. INURI Verlag, zugleich Dissertation der Otto-von-Guericke-Universitat Magdeburg, Juni 2013.

[5] NIST: Fire Dynamic Simulator (FDS) and Smokeview, Version 6.1.2 Rev. 20564.

8

Plagge, Munch

[6] OpenFOAM Foundation: Open Field Operation and Manipulation (OpenFOAM), Version2.2.2-9739c53ec43f.

[7] Robert Siegel, John R. Howell, and Joachim Lohrengel. Warmeubertragung durch Strah-lung, Grundlagen und Materialeigenschaften, Volume 1. Springer-Verlag, 1988.

9