Numerische
Simulation von Explosionsereignissen
Numerical
simulation of explosions
1
Inhaltsverzeichnis
2..... Einleitung
3..... Einfache Methoden zur Beschreibung von
Explosionsszenarien
3.1 Spitzenüberdruck der einfallenden Druckwelle
3.2 Spitzenüberdruck mit Reflexion
3.3 Positive Druckphase - Überdruckphase
3.4 Impuls
3.5 Sogphase
4..... Berechnung mit AT-Blast
4.1 Maximaler Überdruck
4.2 Impuls
5..... Explizite Zeitintegration
5.1 Bewegungsgleichung
5.2 Berechnungsverfahren
5.3 Stabilität und Genauigkeit
5.4 Kritischer Zeitschritt
5.5 Vergleich zur impliziten Zeitintegration
6..... Numerische Simulation
6.1 Materialmodelle
6.2 Vorgehen
6.3 1-D Modell – sphärische Ausbreitung der
Druckwelle
6.4 2-D Modell – sphärische Ausbreitung
6.5 2-D Modell – hemisphärische Ausbreitung, ohne
Reflexion
6.6 3-D Modell – hemisphärische Ausbreitung
6.7 Vergleich der Ergebnisse
7..... Fazit
8..... Ausblick
2
Einleitung
Vor dem Hintergrund der zunehmenden Bedrohung
durch terroristische Anschläge, ist es die Aufgabe von Ingenieuren, insbesondere
von Bauingenieuren, unsere Bauwerke und somit unsere Umwelt unter Berücksichtigung
dieser außergewöhnlichen Einwirkung zu bemessen. Aber nicht nur terroristische
oder militärische Angriffe üben die in dieser Arbeit untersuchten dynamischen
Beanspruchungen auf unsere Bauwerke aus. So sind auch zum Beispiel Unfälle in
technischen Anlagen oder Gasexplosionen in Wohnhäusern als Szenarien denkbar.
Bei einer Explosion wird in einem Bruchteil
einer Sekunde eine enorme Energie freigesetzt. Hierdurch entsteht eine Volumenzunahme,
die wiederum eine Druckwelle auslöst. Diese Druckwelle breitet sich mit bis zu
36.000 km/h aus. Durch die schlagartige Volumenzunahme wird die Luft stark
verdichtet, wodurch sich die Umgebungsluft drastisch erhitzt. In der eigentlich
Definition beinhaltet eine Explosion eine Detonation (Ausbreitung über
Schallgeschwindigkeit) sowie in eine Deflagration (Ausbreitung unter
Schallgeschwindigkeit). Die Begriffe werden in dieser Arbeit jedoch synonym
Verwendet.
Neben dem maximalen Druck interessieren auch
die Dauer der positiven Druckphase und der mit ihr einhergehende Impuls. Die
Bedeutungen dieser Größen sollen in dieser Arbeit beschrieben werden.
Um unsere Bauwerke ein Stück weit sicherer zu
gestalten, soll diese Arbeit dabei helfen, die Einwirkung durch eine Explosion
mit einfachen Mitteln und Methoden zu beschreiben. Hierzu erfolgen die
Auswertung und der Vergleich von in der Literatur angegebenen Formeln zur
Berechnung einer Druckwelle, Rechnungen mit einem einfachen Simulationsprogramm
AT-Blast sowie die numerische Simulation mit ANSYS Autodyn.
Es sei jedoch angemerkt, dass die hier
dargestellten Verfahren zur Ermittlung der Einwirkung einer Explosion nicht
bedingungslos exakte Ergebnisse liefern werden. Jedes Explosionsszenario hat
eine Vielzahl von Randbedingungen, die nur schwer bzw. niemals vollständig zu
erfassen sind.
3 Einfache Methoden zur Beschreibung von Explosionsszenarien
In diesem Kapitel sollen verschiedene
Verfahren zur Beschreibung von Explosionsszenarien vorgestellt werden. Hierzu
gehören die Berechnung des Spitzenüberdrucks, die Ermittlung der positiven
Druckdauer sowie die Berechnung des maximalen Impulses.
Die nachfolgenden Berechnungen zur Ermittlung
des Spitzenüberdrucks einer Druckwelle sollen alle unter Verwendung eines
skalierten Faktors, welcher den Abstand zur Detonationsquelle sowie die Masse
der Sprengladung berücksichtigt, erfolgen. Dieser Faktor ist wie folgt
definiert:
R:
Abstand zum Detonationsursprung [m]
W: Masse
der Sprengladung [kg TNT-Äquivalent]
Mit Hilfe des Faktors Z lassen sich
Explosionen zudem klassifizieren. Diese Klassifizierungen werden häufig bei der
Beschreibung von Explosionen verwendet. So gilt nach Mills [1]
Kontaktexplosion
|
Nahexplosion
|
Fernexplosion
|
Um die Berechnung für die Vielzahl von
verschiedenen Sprengstoffen, und der damit verbundenen Eigenschaften, allgemein
gültig zu machen, erfolgt in der Regel eine Umrechnung der Sprengstoffmasse in
ein TNT-Äquivalent. Hierdurch wird die unterschiedliche massenbezogene
spezifische Energie berücksichtigt. In Tabelle 1 ist beispielhaft ein Auszug der
Umrechnungsfaktoren gegeben.
Tabelle 1:
Umrechnungsfaktoren für explosives Material in TNT-Äquivalent nach [2]
Explosives
Material
|
Massenbezogene spezifische Energie
Qx [kJ/kg]
|
TNT-
Äquivalent (Qx/QTNT)
|
Chemische Verbindung B (60 % Hexagon (RDX) und
40 % TNT)
|
5190
|
1,148
|
RDX (Cyclonit)
|
5360
|
1,185
|
HMX (Oktogen)
|
5680
|
1,256
|
Nitroglycerin (flüssig)
|
6700
|
1,481
|
TNT
|
4520
|
1,000
|
60 % Nitroglycerin Dynamit
|
2710
|
0,600
|
Semtex (Plastiksprengstoff)
|
5660
|
1,250
|
3.1 Spitzenüberdruck
der einfallenden Druckwelle
In diesem Abschnitt sollen verschiedene in der
Literatur angegebene Formeln zur Berechnung des Spitzenüberdrucks bei freier
Ausbreitung betrachtet werden. Unter freier Ausbreitung ist hier das Fehlen von
Hindernissen (Reflexionsflächen) zu verstehen. Dieses Szenario tritt in der
Realität nur sehr selten auf und ist für die Bemessung von Bauteilen i.d.R. irrelevant.
Die Berechnungen liefern jedoch erste, wichtige Grundlagen, um das Verhalten
einer Druckwelle zu beschreiben. Das Ergebnis aller in diesem Kapitel angegeben
Formeln ist der Spitzenüberdruck, d.h. um den tatsächlichen herrschenden
Luftdruck zu ermitteln, ist der Umgebungsluftdruck P0 zu addieren.
3.1.1
Berechnung nach Brode
Nachdem Neumann und Bethe [3] analytische Lösungen zur Berechnung von Druckwellen veröffentlichten,
stellte Brode [4] im Jahr 1954 weitere, aber diesmal numerische Untersuchungen zur
Berechnung des Spitzenüberdrucks einer Explosion an.
Brode gibt die Formel zur Berechnung des
Spitzenüberdrucks bei freier, sphärischer Ausbreitung wie folgt an:
[atmos]
bzw.
[atmos]
für
In der Formel von Brode wird nicht der unter
Abschnitt 3 angegebene
Skalierungsfaktor Z verwendet. Daher soll die Formel nach Brode zunächst so
umgeformt werden, dass der allgemeine Skalierungsfaktor verwendet werden kann. Der
ursprüngliche Skalierungsfaktor nach Brode lautet:
r : Radialer Abstand zum Detonationsursprung [m]
ε = (Etot
: Energiegehalt [J] ; P0: Umgebungsluftdruck [Pa])
Der Energiegehalt pro kg TNT beträgt wie in Tabelle 1 angegeben 4,520 J, der
Luftdruck auf Meereshöhe beträgt durchschnittlich 101325 Pa.
Bildet man den Quotienten aus diesen Werten und zieht die dritte Wurzel, so
erhält man den Faktor f = 3,547. Die Formel nach Brode ist somit mit dem Faktor
f zu modifizieren um mit dem Skalierungsfaktor Z arbeiten zu können.
Die modifizierte Formel nach Brode ergibt sich
somit zu:
[bar]
bzw.
[bar]
für
Bei der Berechnung nach Brode ist zu beachten,
dass diese laut Rutner [2] keine adäquate
Lösung im Nahbereich liefert. Hierzu sollte man andere Verfahren, die im Folgenden
noch vorgestellt werden, heranziehen (z.B. nach Henrych Kap. 0) .
Es ist zu erwähnen, dass die Modifizierung der
Formeln nach Brode bereits in anderen Quellen ([2],[5]) erfolgt ist, und
leicht abweichende Ergebnisse liefert. Die dort angegebenen Ergebnisse lauten
wie folgt:
[bar]
bzw.
[bar]
für
Die genaue Herleitung der Formeln ist nicht im
Detail nachzuvollziehen. Unterschiede können durch unterschiedliche
Bezugsgrößen für Etot (Energiegehalt) und P0
(Umgebungsluftdruck) hervorgerufen werden. Betrachtet man die umgeformten
Gleichungen nach [2] bzw. [5] so stellt man jedoch fest, dass nicht alle Terme der Gleichung mit
einem konstanten Faktor umgeformt wurden, d.h. die Unterschiede können nicht nur
durch einen unterschiedlichen Energiegehalt oder Umgebungsluftdruck hervorgerufen
werden. Die genaue Umformung kann daher nicht nachvollzogen werden und bleibt an
dieser Stelle ungeklärt.
Der Spitzenüberdruck, mit der Berechnung nach
Brode, wird in Abbildung 1, in Abhängigkeit von der Entfernung zum
Detonationsursprung, verdeutlicht.
Es werden die Verläufe beider hier erwähnten
Umformungen dargestellt. Die Unterschiede in den Ergebnissen beider Umformungen
sind, wie man sieht, sehr gering, lediglich im Nahbereich, in dem die Formeln
nach Brode ohnehin nur unzureichende Ergebnisse liefern, kommt es zu
unwesentlichen Abweichungen. In weiteren Betrachtungen werden die Ergebnisse
nach Henrych [5] bzw. Rutner [2] verwendet.
Die Berechnung erfolgte für eine Sprengladung
von 100 kg TNT-Äquivalent.
Abbildung 1: Spitzenüberdruck nach
Brode, hergeleitet und nach [2],[3]
3.1.2
Berechnung nach Henrych
Henrych [5] baut auf den
Ergebnisse von Brode auf und modifiziert dessen Formeln, indem er die von Brode
aufgestellten Differentialgleichungen verändert. Als Ergebnis stellt er drei verschiedene
Gleichungen, abhängig vom skalierten Faktor Z, zur Berechnung des
Spitzenüberdrucks bei freier sphärischer Ausbreitung. Sie lauten wie folgt:
[bar]
für
[bar]
für
[bar]
für
In der folgenden Grafik sollen nun sowohl die
Ergebnisse von Henrych als auch zum Vergleich, die von Brode veranschaulicht
werden. Die Berechnung erfolgte ebenfalls für eine Sprengladung von 100 kg TNT-Äquivalent.
Abbildung 2: Spitzenüberdruck nach
Brode und Henrych
Es zeigt sich, dass die beiden
Berechnungsverfahren sehr ähnliche bis nahezu identische Ergebnisse liefern.
3.1.3
Weitere sphärische Berechnungen
Weitere sphärische Betrachtungen nach Kinney
& Graham [6], Mills [1] oder Naumyenko,
Petrovski & Sadovski [7] liefern ebenfalls sehr ähnliche Ergebnisse. Nachstehend sind die verschiedenen
Formeln sowie deren Ergebnisse in Abbildung 3 dargestellt. Es zeigt sich, dass alle Verfahren im Fernbereich ähnlich
gute Ergebnisse liefern. Im Nahbereich kommt es vereinzelt zu größeren
Abweichungen.
Kinney & Graham:
Die Berechnungen nach Kinney & Graham
beruhen auf semi-empirischen Untersuchungen. Sie geben den Spitzenüberdruck wie
folgt an:
[bar]
mit [bar]
Mills:
[bar] für
Naumyenko, Petrovski und Sadovski:
Die Koeffizienten in den Formeln nach Naumyenko
& Petrovski wurden experimentell ermittelt. Ihre Formeln lauten, in
Abhängigkeit vom skalierten Faktor Z, wie folgt:
[bar] für
bzw.
[bar] für
Abbildung 3: Weitere sphärische
Berechnungen im Vergleich
3.1.4
Kontaktexplosion – Hemisphärische Ausbreitung
Die in den Abschnitten zuvor dargestellten
Formeln gingen von einer sphärischen, d.h. kugelförmigen Druckwellenausbreitung
aus. In der Realität, wird dieses idealisierte Szenario nur sehr selten
eintreten. Wahrscheinlicher ist es, dass der Sprengkörper auf oder nahe einer annähernd
starren Unterlage liegt (Beispiel Autobombe). Hier spricht man von einer
hemisphärischen Ausbreitung der Druckwelle. Eine hemisphärische Ausbreitung
soll hier als quasi freie Ausbreitung bezeichnet werden, da sie lediglich an der
Symmetrieachse der freien Ausbreitung reflektiert wird.
In der Veröffentlichung von Mays & Smith [8] wird ein Faktor zur
Erhöhung des Spitzenüberdrucks bei hemisphärischer Ausbreitung eingeführt. Geht
man von einer ideal ebenen, starren Oberfläche aus, so bildet diese eine Reflexionsebene
und der Faktor wäre somit mit 2 anzusetzen. Da diese Idealisierung jedoch nur
äußerst selten zutrifft, empfehlen Mays und Smith den Faktor 1,8 zu verwenden.
Abbildung 4: Hemisphärische
Druckwellenausbreitung
Der Spitzenüberdruck für eine hemisphärische
Druckwellenausbreitung errechnet sich also zu:
In Kap. 4.1.3 wird ein weiterer Ansatz zu Berechnung des
Spitzenüberdrucks bei hemisphärischer Ausbreitung vorgestellt.
3.2
Spitzenüberdruck mit Reflexion
In dieser Arbeit wurden bisher nur freie (sphärisch)
bzw. quasi freie (hemisphärisch) Ausbreitungen betrachtet. Bei der Ermittlung
von Bemessungslasten können jedoch nicht die unter 3.1 angegebenen Formeln zur
freien Ausbreitung der Druckwelle, verwendet werden. Gesucht sind nun die Größen,
die auf ein Bauteil wirken. Daher muss nun zusätzlich die Reflexion der
Druckwelle durch das Bauteil selbst berücksichtigt werden. Hierzu soll
vereinfacht von einer ebenen Reflexionsfläche ausgegangen werden.
3.2.1
Senkrechte Reflexion an starren Oberflächen
Geht man von einer unendlich großen, unendlich
starren, ideal ebenen und senkrecht angeströmten Oberfläche aus, liegt zunächst
der Gedanke nahe, den einfallenden Druck durch die Reflexion zu verdoppeln.
Diese Näherung trifft wie schon unter 3.1.4 beschrieben für eine Sprengladung
direkt auf bzw. nahe dieser Oberfläche ausreichend genau zu.
Zusätzlich zu dem beschriebenen
Verdopplungseffekt, muss jedoch noch der Staudruck berücksichtigt werden. Der
Staudruck entsteht dadurch, dass auf die reflektierten Luftteilchen
nachströmende Luftteilchen treffen und die reflektierten Teilchen in ihrer
Bewegung behindern. Dies ist schematisch in Abbildung 5 dargestellt.
Abbildung 5: Entstehung des Staudruck,
schematisch
Es ist somit mit einem resultierenden
Überdruck zu rechnen, der mehr als das Zweifache des einfallenden Drucks beträgt.
Rankine & Hugonoit [9] haben hierzu erste Näherungen
zur Druckwellenausbreitung mit Reflexion aufgestellt. Sie gehen bei ihren
Berechnungen von einer unendlich großen, starren und ideal ebenen Oberfläche
aus. Ihre Gleichung lautet wie folgt:
mit
Die Gleichung lässt sich jedoch nur verwenden,
wenn, wie bereits erwähnt, die Druckwelle in Normalenrichtung auf die
Reflektionsfläche trifft. Abbildung 6 zeigt, wie sich die Form der Stoßfront mit
dem Abstand verändert. Man erkennt, dass die Druckwelle mit zunehmendem Abstand
vertikal verläuft. Daher können Flächen, die weit vom Detonationsursprung
liegen, näherungsweise als in normalenrichtung angeströmt betrachtet werden.
Abbildung 6: Entwicklung einer Stoßfront
nach Cooper [10]
Die Ergebnisse der Berechnung nach Rankine
& Hugonoit sollen im Folgenden veranschaulicht werden. Als Eingangsgröße
für den einfallenden Überdruck werden die Berechnungen nach Brode (vgl. 3.1.1) verwendet.
Abbildung 7: Reflektierter
Spitzenüberdruck nach Rankine & Hugonoit
Henrych [5] stellte
ebenfalls eine Gleichung zur Berechnung des reflektierten Spitzenüberdrucks
auf.
Die Gleichung liefert nahezu identische
Ergebnisse zu denen nach Rankine & Hugonoit. Auf eine weitere grafische
Darstellung der Ergebnisse wird daher verzichtet.
Bildet man den Quotienten () aus
reflektierter Druckwelle und einfallender Druckwelle, stellt man fest, dass man
keine Werte kleiner 2 und keine Werte größer 8 erhält (siehe Abbildung 8). Dies gilt nach Gebbeken und Döge [11] nur unter der Annahme, dass sich Luft wie ein ideales Gas verhält (Isotropenexponent). Diese
Annahme trifft jedoch nur bedingt zu. Bei einem hohen einfallenden Druck kann
sich der Isotropenexponent auch verringern, was einen höheren Quotienten zur Folge
hat.
Nach [11] sind diese Abweichungen bis zu einem einfallenden Überdruck von 2 MPa (entspricht bspw. einer Explosion von 2700 kg TNT
in 10 m Entfernung) jedoch vernachlässigbar.
Abbildung 8: Reflexionsfaktor Cr
Dieser Abschnitt sollte einen ersten Eindruck
über die Reflexion von Druckwellen liefern. Für viele Bemessungsfälle ist
jedoch nicht mit einer senkrechten Anströmung zu rechnen. Daher ist für weitere
Betrachtungen der Anströmwinkel von großer Bedeutung. Im Folgenden sollen
dessen Auswirkungen näher betrachtet werden.
3.2.2
Nichtsenkrechte Reflexion
Um den unter 3.2.1 eingeführten
Reflexionsfaktor nun für
allgemeinere Fälle anzuwenden entwickelte Schindler [12] das in Abbildung 9 dargestellte Diagramm. In diesem Diagramm wird wie
schon in den vorherigen Kapiteln der Isotropenexponent zu angesetzt
(Annahme ideales Gas). Es sei noch einmal angemerkt, dass dies in vielen Fällen
eine gute Näherung liefert.
Abbildung 9: Reflexionsfaktor nach Schindler [12]
Es gehen der Auftreffwinkel sowie der
einfallende Spitzenüberdruck in atü, hier als bezeichnet,
ein. Bestimmt man den Reflexionsfaktor nach dem hier dargestellten Diagramm für
einen Auftreffwinkel von 90° (senkrechte Reflexion), kann man die Ergebnisse
mit denen nach Rankine und Hugonoit vergleichen.
Beispiel:
Entfernung zum Detonationsursprung: 4,5 m
Masse TNT-Äquivalent: 100 kg
Nach Brode (vgl. Kap 3.1.1) ergibt sich somit
ein einfallender Druck von:
Nach Gleichung (3.17) ergibt sich der reflektierte Überdruck zu:
Und somit der Reflexionsfaktor zu:
Nach Abbildung 9 erhält man für den einfallenden Druck und den
Faktor
In einem weiteren Diagramm nach TM 5-855-1 [13] wird der
veränderliche Isotropenexponent, abhängig vom einfallenden Druck berücksichtigt
und liefert dementsprechend leicht veränderte Ergebnisse. Es sind die
unterschiedlichen Eingangsgrößen atü (Atmosphären Überdruck – veraltete
Einheit) und MPa () zu
unterscheiden. Außerdem ist der Einfallwinkel unterschiedlich
definiert.
Abbildung 10: Reflexionsfaktor nach TM 5-855-1 [13]
Ermittelt man auch hier für einen einfallenden
Druck und dem
Einfallswinkel (senkrechte
Reflexion) den Reflexionsfaktor so erhält man ebenfalls
3.3
Positive Druckphase - Überdruckphase
Neben der bisher vorgestellten Größe Spitzenüberdruck
gibt es noch weitere für die Bemessung relevante Größen. Hierzu soll zunächst
der Verlauf einer Druckwelle betrachtet werden.
Abbildung 11: Überdruck-Zeit-Verlauf, schematisch
Die Dauer der positiven Druckphase ist eine
wichtige Größe für die Bemessung. Aus ihr geht z.B. der Impuls, welcher
grafisch der Fläche der positiven Druckphase nach Abbildung 11entspricht, hervor. Insbesondere für
Glasbauteile (und andere spröde Werkstoffe) kann auch die Sogphase, d.h.
Unterdruck, von Bedeutung sein.
Zur Berechnung der Dauer der positiven
Druckphase haben
Kinney & Graham [6] eine Formel,
abhängig vom skalierten Faktor Z und der Masse des TNT-Äquivalents, aufgestellt.
mit: mTNT
in [kg]
z in [-]
Aus der Dauer der positiven Druckphase kann
nun der Überdruck-Zeit-Verlauf ermittelt werden. Dieser wird nach Kinney &
Graham [6] idealisiert
durch eine Exponentialfunktion, der so genannten Friedlander-Gleichung,
beschrieben.
Hierbei ist der ermittelte Spitzenüberdruck [hPa] sowie
die Dauer der positiven Druckphase [ms] sowie
ein Formfaktor einzusetzen.
Dieser Formfaktor ist folgender Tabelle zu entnehmen (hier ist nur beispielhaft
ein Auszug angegeben, weitere Werte sind Tabelle XI in [6] zu
entnehmen).
Tabelle 2: Formfaktor α nach [6]
Z [m]
|
α
|
2,4
|
1,04
|
2,5
|
0,99
|
2,6
|
0,94
|
2,7
|
0,90
|
2,8
|
0,86
|
2,9
|
0,82
|
3,0
|
0,79
|
3,5
|
0,67
|
4,0
|
0,60
|
4,5
|
0,54
|
5,0
|
0,50
|
Nach Gebbeken [11] kann bei der Bemessung jedoch auf die Ermittlung dieses Faktors
verzichtet werden, und vereinfacht mit angesetzt
werden. Hierdurch wird ein linearer Druck-Zeit-Verlauf unterstellt. Diese
Vereinfachung wird in Normen angesetzt und soll auch im Folgenden verwendet
werden.
Beispiel:
Sprengung von 100 kg TNT bei sphärischer Ausbreitung im Abstand von 15 m
Nach [6] ergibt
sich:
Nach Gleichung (3.23) ergibt sich:
Der nach Gleichung (3.24) resultierende Druck-Zeit-Verlauf mit ist in Abbildung 12 graphisch dargestellt.
Abbildung 12: Linearisierter Druck-Zeitverlauf,
100 kg TNT – 15 m Entfernung
3.4 Impuls
Anschaulich kann man sich den Impuls als die
„Wucht“, mit der zwei Körper aufeinander treffen vorstellen. Als physikalische
Einheit trägt der Impuls die Größe Masse mal Geschwindigkeit.
Des Weiteren ist der Impuls grafisch gesehen
die Fläche unter der Druck-Zeit-Kurve. Das heißt durch Integration dieser Kurve
erhält man den Impuls. Der maximale Impuls sei mit bezeichnet.
Unter Beachtung der unter 3.3 genannten
Vereinfachung ergibt
sich der vereinfachte maximale Impuls nach Kinney & Graham [6] zu:
Es handelt sich somit um eine Dreiecksfläche
und der Impuls kann einfach anhand des einfallenden Spitzenüberdrucks und der
Dauer der positiven Druckphase bestimmt werden.
Obwohl diese Annahme sehr drastisch ist, wird
auch in technischen Regelwerken wie der DIN EN 13123-1 [14] zur
Klassifizierung von Sprengwirkungshemmungen der Formfaktor als ein Idealfall
vereinfacht zu null angenommen.
3.5
Sogphase
Es soll nun noch ein einfaches Verfahren zur
Berechnung der Sogphase dargestellt werden. Brode [4] hat folgende Gleichung
zur Berechnung des maximalen Unterdrucks aufgestellt:
für Z> 1,6
Ferner wird der negative Impuls wie folgt
angegeben:
Alle hier gezeigten Methoden liefern nur unter
idealisierten Randbedingungen hinreichend genaue Ergebnisse. Zur
Berücksichtigung von z.B. Mehrfachreflexion, wie sie z.B. im städtischen
Bereich vorkommt, müssen komplexere numerische Lösungsverfahren angewandt
werden. Im folgenden werden zunächst verschiedene Berechnungen mit dem Programm
AT-Blast durchgeführt und vorgestellt. Im Anschluss daran werden einige
Grundlagen zur expliziten Zeitintegration erläutert sowie die Unterschiede
zwischen expliziter und impliziter Zeitintegration aufgezeigt. Letztlich
erfolgt eine numerische Simulation von einfachen Explosionsereignissen mit dem
Programm ANSYS Autodyn sowie ein Vergleich der Ergebnisse der verschiedenen
Berechnungen.
4
Berechnung mit AT-Blast
AT-Blast wurde von der Firma “Applied Research
Associates, Inc.“ entwickelt und ermöglicht einfache Berechnungen der
Druckwelle in Folge einer Explosion. Als Eingangsgrößen werden die Masse der
Sprengladung, die minimale und maximal Entfernung zur Sprengladung, die
Schrittweite der Auswertungspunkte sowie der Reflexionswinkel benötigt. Als
Ergebnis erhält man die Geschwindigkeit, den Druck, den Impuls sowie die Dauer
der positiven Druckphase. Alle Ergebnisse beruhen auf einer hemisphärischen
Ausbreitung der Druckwelle. Die Grundlage der Berechnungen in AT-Blast ist
nicht veröffentlicht, daher können nur Vermutungen darüber angestellt werden,
wie AT-Blast rechnet.
In diesem Kapitel sollen die Ergebnisse der
Berechnungen aus AT-Blast zum einen mit denen aus Kapitel 3 und zum anderen mit Werten aus technischen
Regelwerken (ISO 16934) verglichen werden. Hierdurch lassen sich, wie in Kap. 4.1.3 gezeigt, weitere Erkenntnisse über das
Verhalten bei hemisphärischer Ausbreitung gewinnen.
4.1
Maximaler Überdruck
Zunächst soll der maximale Überdruck
betrachtet werden. Hierbei werden die Berechnungen aus AT-Blast mit den
Berechnungen nach Brode bzw. der ISO 16934 verglichen.
4.1.1
Vergleich zu Brode – ohne Reflexion
Die unter 3.1.1 aufgeführten Werte wurden unter Annahme
einer sphärischen Druckwellenausbreitung ermittelt. Wie bereits erwähnt,
rechnet AT-Blast mit einer hemisphärischen Ausbreitung. In Tabelle 3 sind die Ergebnisse nach Brode [4],
deren Umrechnung in hemisphärische Werte nach G. Mays [8]
(vgl. Kap. 3.1.4), sowie die
Ergebnisse aus AT-Blast dargestellt.
Tabelle 3: Vergleich Brode - AT-Blast, hemisphärisch, ohne
Reflexion, 100 kg TNT-Äquivalent
Entfernung
[m]
|
Brode-sphärisch
[kPa]
|
Brode-hemisphärisch mit Mays
[kPa]
|
AT Blast
[kPa]
|
5,00
|
682,00
|
1227,60
|
1157,22
|
7,50
|
252,84
|
455,10
|
464,02
|
10,00
|
133,20
|
239,76
|
239,45
|
12,50
|
84,32
|
151,77
|
145,69
|
15,00
|
59,54
|
107,16
|
99,08
|
20,00
|
35,88
|
64,58
|
56,47
|
25,00
|
24,96
|
44,93
|
37,99
|
30,00
|
18,83
|
33,90
|
28,2
|
35,00
|
14,95
|
26,92
|
22,2
|
40,00
|
12,29
|
22,12
|
18,27
|
45,00
|
10,35
|
18,62
|
15,44
|
50,00
|
8,87
|
15,97
|
13,38
|
Es zeigt sich, dass durch die unter 3.1.4 erwähnte Umrechnung (Faktor 1,8) Werte in der
Größenordnung 6-20 % größer als nach AT-Blast
entstehen. Die Abweichungen können an der Ungenauigkeit der sphärischen
Eingangsgrößen sowie an der Umrechnung in hemisphärische Größen liegen. Die
Ergebnisse sind zudem in Abbildung 13 grafisch dargestellt.
Abbildung 13: Vergleich Brode - AT-Blast, hemisphärisch, ohne
Reflexion, 100 kg TNT-Äquivalent
4.1.2
Vergleich zur ISO 16934 sowie Brode – mit Reflexion
In der internationalen Norm ISO 16934 [15] oder der Deutschen Norm DIN EN 13123-1 [14] werden zur Klassifizierung von Glasscheiben verschiedene
Explosionsszenarien dargestellt. Alle Szenarien beruhen auf einer
hemisphärischen Druckwellenausbreitung sowie einer senkrecht angeströmten
Reflexionsfläche und wurden nach [13] berechnet.
Vergleicht man die in der Norm angegebenen
Werte mit den Ergebnissen aus AT-Blast, so stellt man fest, dass die Werte sich
exakt gleichen. Nun liegt bereits die Vermutung nah, dass die Berechnungen in
AT-Blast ebenfalls nach [13] erfolgen. Zum Vergleich wurden außerdem die mit dem Faktor 1,8
umgerechneten Ergebnisse nach Brode, d.h. hemisphärisch ohne Reflexion, durch
die unter 3.2.2 dargestellten
Diagramme in Werte für hemisphärische Ausbreitung mit Reflexion umgerechnet. In
Tabelle 4 sind der Ergebnisse
dargestellt.
Tabelle 4: Vergleich Spitzenüberdruck ISO, Brode und
AT-Blast, hemisphärisch, mit Reflexion
Ent-fernung [m]
|
TNT-Äquivalent
[kg]
|
ISO
16934
[kPa]
|
AT Blast
[kPa]
|
Brode ohne Ref.
[kPa]
|
Brode mit Ref. nach [12]
[kPa]
|
Brode mit Ref. nach [13]
[kPa]
|
RL nach [12]
[-]
|
RL nach [13]
[-]
|
33,00
|
30,00
|
29,00
|
28,96
|
16,31
|
35,55
|
34,24
|
2,18
|
2,10
|
34,00
|
100,00
|
51,00
|
50,75
|
28,10
|
63,22
|
60,41
|
2,25
|
2,15
|
33,00
|
160,00
|
70,00
|
69,71
|
37,23
|
86,00
|
81,90
|
2,31
|
2,20
|
39,00
|
500,00
|
104,00
|
104,25
|
52,02
|
125,89
|
117,56
|
2,42
|
2,26
|
41,00
|
1000,00
|
153,00
|
153,34
|
70,24
|
178,42
|
175,61
|
2,54
|
2,50
|
46,00
|
2000,00
|
201,00
|
200,91
|
85,93
|
230,30
|
218,27
|
2,68
|
2,54
|
49,00
|
2500,00
|
206,00
|
206,43
|
87,66
|
236,68
|
223,53
|
2,70
|
2,55
|
Auch hier fallen die gegenüber den Werten aus
AT-Blast bzw. ISO 16934 höheren Werte auf. Die Verwendung anderer
Eingangsgrößen als Brode z.B. nach Kinney würde zu noch höheren Werten führen. Die
Ergebnisse sind noch einmal in Abbildung 14 grafisch dargestellt. Man stellt einen sehr ähnlichen Verlauf fest.
Abbildung 14: Vergleich ISO, Brode und AT-Blast,
hemisphärisch, mit Reflexion
Allgemein wird es als unzureichend betrachtet,
eine pauschale Umrechnung unabhängig von der Entfernung durch einen konstanten
Faktor vorzunehmen, wie es in Kap. 3.1.4 erfolgte. Weitere Überlegungen zu dem Thema
erfolgen in Kap. 4.1.3.
Ermittelt man nun die Ergebnisse aus AT-Blast
einmal mit und einmal ohne Reflexion und errechnet das Verhältnis, so fällt
auf, dass diese Werte sehr stark den abgelesenen Werten nach Abbildung 10 ähneln. Auch dies ist ein weiteres
Indiz dafür, dass die Berechnungen in AT-Blast nach [13] erfolgen.
Tabelle 5: Reflexionsfaktor aus AT-Blast
AT Blast mit
Reflexion [kPa]
|
AT Blast ohne
Reflexion [kPa]
|
RL
errechnet
|
RL
abgelesen
|
28,96
|
13,65
|
2,12
|
2,08
|
50,75
|
23,24
|
2,18
|
2,15
|
69,71
|
31,10
|
2,24
|
2,25
|
104,25
|
44,54
|
2,34
|
2,39
|
153,34
|
61,91
|
2,48
|
2,44
|
200,91
|
77,50
|
2,59
|
2,59
|
206,43
|
79,22
|
2,61
|
2,60
|
4.1.3
Neuer Ansatz zur hemisphärischen Ausbreitung einer
Druckwelle
In diesem Abschnitt wird ein eigener Ansatz
zur Berechnung von hemisphärischen Ausbreitungen vorgestellt. Im Gegensatz zu
den Überlegungen nach Mays soll hier jedoch nicht der Druck für sphärische
Ausbreitung mit einem Faktor versehen werden, sondern die Masse der
Sprengladung verändert werden. Die Sprengladung wurde unter Annahme von
Idealbedingungen mit dem Faktor 2 erhöht. In Tabelle 6 sind die Ergebnisse
einer Berechnung für 200 kg TNT-Äquivalent, bei
sphärischer Ausbreitung nach Brode (Annahme: entspricht 100 kg TNT-Äquivalent bei hemisphärischer Ausbreitung), 100 kg TNT-Äquivalent bei sphärischer Ausbreitung nach Brode modifiziert nach Mays sowie 100 kg TNT-Äquivalent nach AT-Blast dargestellt.
Tabelle 6: Hemisphärische
Ausbreitung Brode 200 kg
Entfernung [m]
|
Brode (200kg)
[kPa]
|
Brode (100kg) – mit Mays (1,8)
[kPa]
|
AT-Blast
[kPa]
|
5,00
|
1247,18
|
1227,60
|
1157,22
|
7,50
|
439,92
|
455,10
|
464,02
|
10,00
|
221,88
|
239,76
|
239,45
|
12,50
|
135,47
|
151,77
|
145,69
|
15,00
|
92,89
|
107,16
|
99,08
|
20,00
|
53,67
|
64,58
|
56,47
|
25,00
|
36,36
|
44,93
|
37,99
|
30,00
|
26,97
|
33,90
|
28,20
|
35,00
|
21,18
|
26,92
|
22,20
|
40,00
|
17,29
|
22,12
|
18,27
|
45,00
|
14,28
|
18,62
|
15,24
|
50,00
|
14,51
|
15,97
|
13,38
|
Man stellt fest, dass die Ergebnisse der
Berechnung für 200 kg TNT-Äquivalent im Schnitt
nur um 4% von den Ergebnissen aus AT-Blast abweichen. Die Ergebnisse für 100 kg unter Berücksichtigung des Faktors nach Mays weichen
durchschnittlich um 11% ab. Die Vereinfachung, die Masse der Sprengladung für
eine hemisphärische Ausbreitung gegenüber einer sphärischen Ausbreitung zu
verdoppeln, scheint unter Verwendung der Berechnungen nach Brode eine gute
Näherung zu sein.
Abbildung 15: Hemisphärische Ausbreitung Brode 200 kg
4.2
Impuls
Zur Berechnung des maximalen Impulses benötigt
man wie bereits unter Kap. 3.4
erwähnt zunächst die Dauer der positiven Druckphase. Hierzu wurde unter
Verwendung der unter Kap. 3.3
angegebenen Formel (Kinney & Graham) für die Dauer der positiven Druckdauer
sowie nach Gl. (3.28) der maximale Impuls ermittelt.
Die Ergebnisse sollen nun mit den Berechnungen aus AT-Blast verglichen werden.
Abbildung 16: Impuls nach Kinney & Graham und AT-Blast
Es zeigt sich, dass die Vereinfachung nach
Kinney und Graham in größerer Entfernung sehr ähnliche Ergebnisse wie AT-Blast
liefert. Kommt man näher an den Detonationsursprung, so sind die Abweichungen
nach Kinney und Graham zum Teil erheblich (bis zu 70 %), für die Bemessung von
Bauteilen, liegt der vereinfachte Ansatz auf der sicheren Seite.
AT-Blast bietet die Möglichkeit, mit
einfachsten Mitteln wichtige Kenngrößen einer Explosion zu ermitteln. Die
AT-Blast zugrunde liegende Rechenvorschrift ist zwar nicht veröffentlicht, aber
die Vermutung, dass [13] als Grundlage dient, ist äußerst wahrscheinlich.
5
Explizite Zeitintegration
Wie im vorherigen Abschnitt beschrieben, ist
das Lösen von komplexeren Problemen in der Regel nur noch durch numerische
Verfahren möglich. Grundlage für die numerischen Berechnungen sind
mathematische Modelle, welche in der Regel aus Differentialgleichungen
bestehen. Zum besseren Verständnis soll an dieser Stelle die allgemeine
Bewegungsgleichung vorgestellt werden. Durch schrittweises Lösen der Bewegungsgleichung
erhält man sowohl die Verschiebung, die Geschwindigkeit als auch die
Beschleunigung eines Körpers.
5.1
Bewegungsgleichung
Der Übersicht halber wird ein Ein-Freiheitsgrad-Schwinger
mit Dämpfung betrachtet. Gegeben sei folgendes System:
Abbildung 17: Ein-Freiheitsgrad-Schwinger
Aus der Mechanik ist bekannt, dass sich die
Kräfte wie folgt ergeben:
Mit dem Gleichgewicht nach dem Prinzip von
d’Alambert folgt:
Zum numerischen Lösen der Bewegungsgleichung
ist eine Diskretisierung der Zeit erforderlich. Hierbei kann sowohl eine
implizite als auch eine explizite Zeitintegration verwendet werden. Liegt nun
ein stark nichtlinearer Fall vor, d.h. z.B. große Deformationen oder sehr kurze
Belastungszeit, so stößt die implizite Zeitintegration an ihre Grenzen (vgl. 5.5). Im Weiteren sollen nun verschiedene explizite
Integrationsverfahren zur Berechnung solcher Probleme vorgestellt werden.
5.2
Berechnungsverfahren
Die explizite Zeitintegration fällt unter die Integrationsverfahren
für Anfangswertprobleme. Explizite Verfahren werten die Bewegungsgleichung am alten,
bekannten Zeitschritt aus, implizite Verfahren dagegen verwenden den neuen,
unbekannten Zeitschritt. Die Unterschiede beider Verfahren sollen später gesondert
betrachtet werden.
Bei den expliziten Verfahren gibt es
verschiedene Ansätze von denen einige im Folgenden betrachtet werden sollen.
5.2.1
Euler‘sches Polygonzugverfahren (Euler-vorwärts)
Dieses Verfahren ist ein sehr einfaches
Verfahren zur Lösung von Anfangswertproblemen. Wir betrachten zur
Verdeutlichung folgendes Anfangswertproblem 1. Ordnung.
Entwickelt man die Funktion in eine
Taylor-Reihe so ergibt:
Durch Vernachlässigung des Fehlerterms erhält
man eine erste Näherungslösung. Es liegt nah, dass mit zunehmender Schrittweite
die
Genauigkeit abnimmt.
Jede Differentialgleichung höherer Ordnung
lässt sich in ein System von Differentialgleichungen 1. Ordnung
überführen. Daher lassen sich auch Systeme höherer Ordnung durch dieses
Verfahren berechnen.
Das Euler’sche Polygonzugverfahren ist ein
sehr einfaches Verfahren mit dem man schnell zu Ergebnissen kommt. Es ist
jedoch im Gegensatz zu anderen Verfahren mit einem größeren Fehler behaftet.
5.2.2
Mehrschrittverfahen (Zentrale Differenzenmethode)
Das Mehrschrittverfahren greift im Gegensatz
zum Euler-vorwärts Verfahren auf mehrere zurückliegende Ergebnisse zurück.
Dadurch kann eine höhere Genauigkeit erzielt werden. Das Problem dabei ist
jedoch, dass zu Beginn häufig nicht ausreichend Ergebnisse vorliegen. Hierzu
müssen weitere Berechnungen erfolgen.
Wie bereits geschildert, werden bei diesem
Berechnungsverfahren zur Berechnung des Zeitschritts die
Ergebnisse aus dem Zeitschritten sowie benötigt.
Abbildung 18: Integration mit der zentralen
Differenzenmethode
Setzt man einen linearen Verschiebungsverlauf voraus, so erhält man
nach Blaut [16]
Stellt man
diese Gleichung um, so erhält man:
Dieses Verfahren bietet im Vergleich zum
Euler-vorwärts Verfahren etwas genauere Ergebnisse, ist jedoch auch etwas
aufwändiger in der Anwendung. Trotzdem gilt das Verfahren als ein sehr
einfaches und wird z.B. in Programmen wie ANSYS Autodyn oder LS DYNA
eingesetzt.
5.2.3
Prädiktor-Korrektor & Runge-Kutta-Verfahren
Das Prädiktor-Korrektor-Verfahren bzw. das
Runge-Kutta-Verfahren sind Mischformen aus impliziten und expliziten Integrationsverfahren.
5.2.3.1
Prädiktor-Korrektor Verfahren (Heunsches Verfahren)
Hierbei wird ein expliziter Startwert
(Prädiktor) sowie
eine implizite Integrationsregel (Korrektor) verwendet. Der Prädiktor ist im
Prinzip ein Schätzwert, der dann im impliziten Teil des Verfahrens verwendet
wird und sich nach Müller [17]
wie folgt ergibt:
Anschließend erhält man über die Trapezregel
die eigentliche Lösung zum Zeitpunkt t+Δt durch
Abbildung 19: Integration mit dem Heunschen Verfahren
5.2.3.2
Runge-Kutta-Verfahren
Das
Runge-Kutte-Verfahren verwendet im Gegensatz zum Prädiktor-Korrektor-Verfahren
nicht nur einen, sondern mehrere Schätzwerte. Hierdurch kann die Genauigkeit
erhöht werden. Zum Runge-Kutta-Verfahren gibt es diverse Varianten. An dieser
Stelle soll das Verfahren 4.Ordnung oder auch Standard-Runge-Kutta-Verfahren
dargestellt werden.
Nach Müller [17]
gilt
Euler-Prädiktor
(Halbschritt)
Euler
Korrektor (Halbschritt)
Mittelpunkts-Prädiktor
(Ganzschritt)
Simpson-Regel-Korrektor
(Ganzschritt)
5.3 Stabilität
und Genauigkeit
Je höher die Ordnung eines Verfahrens desto
mehr Therme der Taylorreihe werden mitgenommen. Daher hängt die Genauigkeit der
Berechnung mit der Ordnung des Verfahrens zusammen. So ist beispielsweise, das
Euler-Vorwärts-Verfahren ein Verfahren erster Ordnung,
das Mehrschrittverfahren eines zweiter Ordnung und das Standard-Runge-Kutta-Verfahren
eines vierter Ordnung.
Wenn sich bereits kleine Ungenauigkeiten im
Laufe des Integrationsprozesses aufsummieren und zu untragbaren Fehlern führen
spricht man von einem instabilen System. Allgemein lässt sich sagen, dass ein
explizites Verfahren nie als bedingungslos stabil angesehen werden kann. Dies
kann nur durch implizite Verfahren erreicht werden. Es ist zu beachten, dass diese
Eigenschaft jedoch keine direkte Aussage über die Genauigkeit zulässt.
Ein wichtiges Kriterium für die Stabilität
expliziter Berechnungen ist der gewählte Zeitschritt pro Integrationsschritt.
Dieser Zeitschritt wird bei expliziten Verfahren durch den kritischen Zeitschritt
begrenzt.
5.4 Kritischer
Zeitschritt
Der kritische Zeitschritt zur Berechnung eines
Systems hängt von dessen Eigenfrequenz ab und wird mit bezeichnet.
Der verwendete Zeitschritt muss kleiner als der kritische Zeitschritt gewählt
werden.
Moderne Berechnungsprogramme wie z.B. LS-DYNA
oder ANSYS-Autodyn berechnen die kritischen Zeitschritte eines Systems automatisch.
Es gibt dennoch Möglichkeiten, die Größe des Zeitschritts zu steuern. Möchte
man einen Zeitschritt wählen, der größer als die kritischen Zeitschritte
einzelner Elemente des Systems ist, so werden diese Elemente mittels
Massenskalierung bearbeitet.
Über die Massenskalierung wird die Dichte der
Elemente so verändert, dass der errechnete kritische Zeitschritt dem gewählten
Zeitschritt entspricht.
5.5 Vergleich
zur impliziten Zeitintegration
Wie schon einleitend in diesem Kapitel
beschrieben, wird bei expliziten Verfahren die dynamische Bewegungsgleichung
zum Zeitpunkt und bei
impliziten zum Zeitpunkt ausgewertet.
Ein Vorteil der expliziten Rechnung gegenüber
der impliziten ist der geringere Rechenaufwand und damit einhergehend eine
höhere Rechengeschwindigkeit. Dabei sollte jedoch beachtet werden, dass die
explizite Rechnung weniger stabil ist als die implizite.
Die unter Kap. 5.4 genannte Einschränkung der Schrittweite gilt
für implizite Rechnungen nicht. Hierdurch lassen sich größere Schrittweiten
generieren, was durch den erhöhten Rechenaufwand pro Zeitschritt vorteilhaft
ist. Explizite Verfahren finden bei starken quasistatischen Verformungen,
transienten Analysen mit hoher nichtlinearer Dynamik (Crashtests) sowie
Schockwellenberechnung ihren Einsatz. Implizite dagegen eher bei statischen
Simulationen oder transienten Berechnungen mit linearer Dynamik.
6
Numerische Simulation
Neben den bereits in Kapitel 3 und 4 betrachteten Berechnungen erfolgt nun eine
Simulation mit Hilfe des Programms ANSYS Autodyn. Dies ist ein explizites
FE-Programm zur Beschreibung nichtlinearer dynamischer Vorgänge von
Festkörpern, Flüssigkeiten, Gasen sowie deren Interaktion.
In diesem Kapitel sollen verschiedene Modelle
sowie verschiedene Parameter innerhalb der Modelle untersucht und verglichen
werden. Die Ergebnisse werden anschließend in Kapitel 6.7 mit den Ergebnissen der Berechnungen der
anderen Kapitel gegenübergestellt.
6.1 Materialmodelle
In Autodyn liegt eine Vielzahl von Materialien
in einer Bibliothek als Modelle bereits vor. Für die im Rahmen dieser Arbeit zu
untersuchenden Szenarien sind die Materialen „Air“ sowie „TNT“ bzw. „TNT-2“
ausreichend. Durch Einbringen von geeigneten Randbedingungen lassen sich
sphärische und hemisphärische Ausbreitungen bzw. Reflexionen simulieren.
Die Materialmodelle von „TNT“ sowie „TNT-2“ beinhalten
verschiedene Parameter, die als Eingangsgrößen für die der Rechnung zu Grunde
liegenden Zustandsgleichung von Jones-Wilkins-Lee (JWL) (Gl. (6.1) nach Baudin [18]) benötigt
werden. Dazu zählen die Parameter . sowie die
Dichte des Sprengstoffs
und die
Dichte der Detonationsprodukte .
Weiterhin gilt: .
Die für „TNT“ angegebenen Parameter
entsprechen den üblichen Literaturwerten nach Dobratz [19].
Der Unterschied der Materialmodelle „TNT“ bzw.
„TNT-2“ wird innerhalb des 1-D Modells untersucht. Desweiteren sind die Materialdaten in Abbildung 20 gegeben.
Das Materialmodell „Air“ beruht auf der
Annahme eines idealen Gases. Daher ist der Isotropenexponent zu angegeben
(vgl. Kap. 3.2.1). Neben ihm sind die
Dichte, die spezifische Wärmekapazität, sowie eine Bezugstemperatur gegeben.
Abbildung 20: Materialdaten TNT und TNT-2
6.2
Vorgehen
In Autodyn lassen sich 1-D, 2-D als auch 3-D
Modelle erstellen. Je nach Dimension und Größe des Systems variieren die
Rechenzeiten enorm. Um den Aufwand insbesondere durch die lange Rechenzeiten in
einem angemessenen Maß zu halten, wurde entschieden, die ausführliche Parameterstudie
an einem 1-D Modell durchzuführen.
Bei den 2-D Modellen werden zwei verschiedene
Grundtypen, „Box“ und „Circle“ näher betrachtet. Diese sollen nun noch unter
Berücksichtigung der Erfahrungen aus den 1-D Modellen in ihrem Netztyp variiert
werden. Weiterhin soll die Funktion des Remappings untersucht werden, bei dem
die Ergebnisse eines 1-D Modells in ein 2-D überführt werden. In Abbildung 21 ist das Vorgehen grafisch dargestellt.
Abbildung 21:
Vorgehensweise numerische Simulation der sphärischen
Druckwellenausbreitung
Im Anschluss an diese Untersuchungen zur sphärischen
Druckwellenausbreitung erfolgt noch eine Simulation für den hemisphärischen
Fall an dem 2-D Modell „Circle“. Abschließend wird beispielhaft ein 3-D
Szenario in einer realen Beispielumgebung simuliert.
6.3 1-D
Modell – sphärische Ausbreitung der Druckwelle
Der Grundaufbau dieses
Modells, sowie aller anderen Modelle in dieser Arbeit beruhen auf folgendem
Vorgehen:
Zunächst werden die Materialien aus der in
Autodyn hinterlegten Bibliothek geladen. Danach wird die Grundgeometrie als einer
neuer „Part“, welcher den Betrachtungsraum darstellt, erstellt. Als Berechnungsansatz
oder auch „Solver“ wird „Euler, 2D Multi-Material“ gewählt. Nun muss die
Geometrie definiert werden, in diesem Fall ein „Wedge“ (Keil).
Der „Wedge“ kann auf Grund seiner eindimensionalen Netzausbreitung als ein 1-D
Modell angesehen werden. In den 2-D Modellen werden die Formen „Box“ und
„Circle“ gewählt. Nach Eingabe der Abmessungen wird das Netz definiert. Die
dort möglichen Einstellungen sind von der gewählten Geometrie abhängig. Bei einem
„Wedge“ wird lediglich die Anzahl der Elemente über den Radius definiert.
Abschließend wird die Geometrie mit dem Material „Air“ gefüllt. Nachdem nun die
Umgebung fertiggestellt ist wird der Sprengstoff platziert. Hierzu wird an dem
Ort, an dem die Sprengung stattfinden soll ein Teil der Geometrie mit „TNT“
gefüllt. Über das Volumen bzw. in 1-D und 2-D Betrachtungen die Fläche, wird
die Menge und somit die Masse des Sprengstoffs definiert (vgl. Gl. (6.2)). In allen Modellen wird stets eine kugel- bzw. kreisförmige
Sprengladung platziert.
Anschließend wird der Detonationsursprung im
Zentrum der Sprengladung definiert. Abschließend werden noch, neben weiteren
Randbedingungen wie zum Beispiel Symmetrieebenen oder Reflexionsflächen,
„Gauges“ (Messpunkte) zu Erfassung der entstehenden Druckwelle platziert.
In Abbildung 22 ist der Modellaufbau der 1-D-Simulation
schematisch dargestellt. Es handelt sich um einen keilförmigen Körper dessen
untere Kante eine Rotationsachse bildet. Der Keil hat eine Länge von 20 m und entlang der Symmetrieachse sind in 5 m bzw.
15 m Abstand vom Detonationsursprung Messpunkte positioniert.
Am rechten Rand des Keils wurde als
Randbedingung ein so genannter „Flow-Out“ gewählt. Dies bedeutet, dass die
Druckwelle nicht gestaut wird und somit ungehindert entweichen kann. Hierdurch
wird eine freie Ausbreitung der Druckwelle simuliert. In Abbildung 23 ist außerdem der Rotationskörper
dargestellt.
Abbildung 22: 1-D Modell schematisch
Abbildung 23: 1-D Modell rotiert
Einer der zu untersuchenden Parameter und
dessen Auswirkung ist die Netzdichte. In dem vorliegenden Modell handelt es
sich um ein sehr einfaches Netz, dass, wie bereits erwähnt, nur über die Anzahl
der Zellen über den Radius definiert wird.
Ein weiterer Parameter ist der Zeitschritt.
Wie bereits in Kapitel 5.4 beschrieben, erfolgt die Berechnung des kritischen
Zeitschritts in Autodyn automatisch. Zur Untersuchung des Einflusses des
Zeitschritts auf die Ergebnisse soll dieser jedoch auch manuell eingegeben bzw.
nach oben begrenzt werden.
6.3.1
Variation der Netzdichte
Bei der Wahl des Netzes ist darauf zu achten,
dass ausreichend Zellen mit der Sprengladung gefüllt werden. Nach [20] sollten mindesten 10 Zellen gefüllt sein.
Zur Untersuchung des Einflusses der Netzdichte
auf die Ergebnisse wurden folgende Netze untersucht:
Netz A: 500
Zellen über Radius (20 m) à Zellgröße 4,00 cm
Netz B:
750 Zellen über Radius (20 m) à Zellgröße 2,67 cm
Netz C:
1000 Zellen über Radius (20 m) à Zellgröße 2,00 cm
Netz D:
2000 Zellen über Radius (20 m) à Zellgröße 1,00 cm
Netz E: 4000
Zellen über Radius (20 m) à Zellgröße 0,50 cm
Netz F:
8000 Zellen über Radius (20 m) à Zellgröße 0,25 cm
Die Netze A und B erfüllen zwar nicht die eingangs
des Abschnittes genannte Bedingung zur Zahl der mit TNT gefüllten Zellen,
wurden dennoch zur Verdeutlichung einer Konvergenz mit untersucht.
In Tabelle 7 sind die Ergebnisse der Berechnungen
aus Autodyn angegeben. Es zeigt sich, dass mit Zunahme der Netzdichte der
Spitzenüberdruck größer wird. Wie schon in Kapitel 3 festgestellt werden konnte, nehmen die
Schwankungen von mit
zunehmenden Abstand zum Detonationsursprung ab.
Tabelle 7: Ergebnisse
Variation Netzdichte, 1-D
Netz
|
(5 m)
[kPa]
|
(15 m)
[kPa]
|
Netz A
|
627,39
|
61,31
|
Netz B
|
655,80
|
62,84
|
Netz C
|
676,06
|
63,29
|
Netz D
|
718,30
|
64,95
|
Netz E
|
748,32
|
66,01
|
Netz F
|
759,31
|
66,28
|
In Abbildung 24 ist einmal der Druck Zeitverlauf für Netz F
beispielhaft dargestellt. Der unter Kapitel 3.3 vorgestellte typische Verlauf (Friedlander-Kurve)
ist gut zu erkennen. Es ist zu beachten, dass die in Tabelle 7 angegebenen Werte den Spitzenüberdruck
darstellen, d.h. der Umgebungsluftdruck (101,33 kPa) wurde abgezogen, in Abbildung 24 ist er enthalten.
Abbildung 24: 1-D Modell, Druck-Zeitverlauf aus Autodyn, Netz
F, 5 m
In Abbildung 25 sind Plots zu verschiedenen Zeitpunkten des 1-D Modells gegeben.
Hierbei ist der Rotationskörper des Keils, ein Kegel, dargestellt. Man erkennt
eine klare Spitze der Werte an der Druckwellenfront, die danach schnell
abfällt. Dies spiegelt den in Abbildung 24 angegebenen Verlauf gut wieder. Außerdem ist die negative Druckphase
gut zu erkennen.
Abbildung 25: Plot Druckverlauf, 1-D Modell, Netz C
Zur Feststellung einer möglichen Konvergenz
der Ergebnisse, sind diese in Abhängigkeit der Netzdichte in Abbildung 26 für
eine Entfernung von 5 m sowie in Abbildung 27 für 15 m aufgetragen.
Abbildung 26: 1-D Modell, Spitzenüberdruck abhängig von
Netzdichte, 5 m
Abbildung 27: 1-D Modell, Spitzenüberdruck abhängig von
Netzdichte, 15 m
Es zeigt sich, dass die Ergebnisse mit
zunehmender Netzdichte konvergieren. Eine weitere
Verdichtung des Netzes wurde auf Grund der hohen Rechenzeiten und der erzielten
Konvergenz nicht weiter untersucht.
6.3.2
Variation des Zeitschritts
In Kapitel 5.4
wurde bereits der Begriff „kritischer Zeitschritt“ erläutert. Autodyn berechnet
standardmäßig den kritischen Zeitschritt automatisch und versieht diesen noch
mit einem Sicherheitsfaktor von 0,66. Zur Steigerung der Genauigkeit der
Ergebnisse soll der kritische Zeitschritt nun manuell begrenzt werden. Hierzu
wird eine obere Grenze festgelegt. Die untersuchten Zeitschrittweiten sowie die
Ergebnisse sind in Tabelle 8
dargestellt.
Tabelle
8: Ergebnisse Variation
Zeitschritt, 1-D
Δtmax
[ms]
|
(5 m)
[kPa]
|
(15 m)
[kPa]
|
Auto Timest.
|
655,80
|
62,84
|
0,1
|
655,80
|
62,84
|
0,05
|
655,80
|
62,84
|
0,01
|
655,44
|
62,84
|
0,005
|
643,04
|
62,29
|
0,001
|
619,73
|
60,74
|
0,0005
|
618,22
|
60,61
|
0,0001
|
616,79
|
60,47
|
0,00005
|
616,71
|
60,45
|
Man sieht, dass die manuelle Begrenzung erst
ab einem maximalen Zeitschritt von 0,01 zu Abweichungen in den Ergebnissen
führt. Davor liegt der automatisch errechnete Zeitschritt von Autodyn unter dem
manuell definierten Zeitschritt.
Eine grafische Übersicht der Ergebnisse für 5 m bzw. 15 m in logarithmischer Darstellung liefern Abbildung 28 und Abbildung 29.
Abbildung 28: 1-D Modell,
Spitzenüberdruck abhängig von Zeitschritt, 5 m
Abbildung 29: 1-D Modell, Spitzenüberdruck abhängig von Zeitschritt,
15 m
Im Gegensatz zur Konvergenz bei Verfeinerung
der Netzdichte entsteht bei der Begrenzung des Zeitschritts eine Konvergenz die
scheinbar eine untere Grenze besitzt.
6.3.3
Variation des Sprengstoffes
Wie schon unter 6.1 beschrieben liegen in der Materialbibliothek
von Autodyn zwei Modelle für TNT vor. Zur Untersuchung der Auswirkungen des
TNT-Modells auf die Ergebnisse wurden verschiedenen Berechnungen einmal für
„TNT“ und einmal für „TNT-2“ durchgeführt. Die Ergebnisse für den maximalen
Druck sowie deren Verhältnis sind in Tabelle 9 gegeben.
Tabelle 9: Ergebnisse Variation TNT, 1-D
|
A
|
B
|
C
|
D
|
|
|
Netzdichte
|
(5 m,TNT)
[kPa]
|
(5 m,TNT-2)
[kPa]
|
(15 m,TNT)
[kPa]
|
(15 m,TNT-2)
[kPa]
|
B / A
|
D / C
|
500
|
627,39
|
681,80
|
61,31
|
64,98
|
1,087
|
1,060
|
750
|
655,70
|
718,34
|
62,84
|
66,65
|
1,095
|
1,061
|
1000
|
676,06
|
742,30
|
63,29
|
67,13
|
1,098
|
1,061
|
2000
|
718,10
|
792,10
|
64,84
|
68,58
|
1,103
|
1,058
|
Es zeigt sich, dass die Ergebnisse für „TNT-2“
in einer Entfernung von 5 m um ca. 10 % und in einer Entfernung von 15 m um ca. 6 % größer sind als für „TNT“. In
Abbildung 30 und Abbildung 31 sind die Ergebnisse noch einmal grafisch
dargestellt.
Abbildung 30: 1-D Modell, Vergleich TNT mit TNT-2, 5 m
Abbildung 31: 1-D Modell, Vergleich TNT mit TNT-2, 15 m
In den folgenden Berechnungen wird „TNT“
verwendet.
6.4
2-D Modell – sphärische Ausbreitung
6.4.1
Modell A „Box“:
Zur Untersuchung einer sphärischen Ausbreitung
der Druckwelle mit Erfassung der Werte in 5 m und
15 m Entfernung muss eine „Box“ mit den Mindestmaßen 30 m x 30 m (gewählt 40 m x 40 m) und einer mittig platzierten Sprengladung erstellt werden. Durch
Verwendung einer Symmetrieachse kann das Modell in zwei Teile geteilt und es
verbleibt eine „Box“ mit den Maßen 40 m x 20 m. Rein optisch lässt dieses Modell eine hemisphärische
Druckwellenausbreitung vermuten. Durch die Symmetrieachse entspricht dieses
reduzierte Modell jedoch dem großen Ausgangsmodell und simuliert somit eine
sphärische Ausbreitung. Durch Rotation entlang der Symmetrieachse entsteht ein
Zylinder welcher eine anschauliche Darstellung (Abbildung 33) der Druckwellenausbreitung zulässt.
Neben Messpunkten („Gauges“) im Abstand von 5 m bzw. 15 m entlang der horizontalen Achse wurden auch Messpunkte entlang einer
um 45° geneigten Achse sowie entlang der vertikalen Achse im jeweils gleichen Abstand
(5 m und 15 m) installiert. In Abbildung 32 ist der beschriebene Modellaufbau
schematisch gezeigt.
Abbildung 32: 2-D Modell A, schematisch
Abbildung 33: 2-D Modell A, rotiert 270°
6.4.1.1
Variation der Netzdichte
Zur Untersuchung des Einflusses der Netzdichte
auf die Ergebnisse wurden folgende Netze untersucht:
Netz A: Zellgröße = 10,0 cm
Netz B: Zellgröße = 5,00 cm
Netz C: Zellgröße = 2,50 cm
Eine weitere Verdichtung des Netzes führt zu unverhältnismäßig
großen Rechenzeiten. Dennoch werden weitere Netzverfeinerungen an einem
kleineren, und somit schneller zu berechnenden, Modell untersucht. Hierzu wurde
ein Modell mit den Kantenlängen 20 m und 10 m erstellt. Die Messpunkte in 15 m Entfernung entfallen
daher.
Hierbei wurden folgende Netze weiterhin
betrachtet:
Netz D: Zellgröße = 2,50 cm
(Vergleichsnetz zu Netz C)
Netz E: Zellgröße = 1,25 cm
Die Ergebnisse der Berechnungen sind in Tabelle 10 aufgetragen.
Tabelle 10: Ergebnisse Variation Netzdichte, Modell A, 2-D
sphärisch
Netz
|
(5 m)
horizontal
[kPa]
|
(15 m)
horizontal
[kPa]
|
(5 m)
45°
[kPa]
|
(15 m)
45°
[kPa]
|
(15 m)
vertikal
[kPa]
|
(15 m)
vertikal
[kPa]
|
Netz A
|
561,60
|
67,39
|
660,85
|
56,32
|
607,36
|
65,85
|
Netz B
|
600,29
|
69,33
|
696,08
|
59,70
|
640,27
|
67,97
|
Netz C
|
643,18
|
70,13
|
692,28
|
62,46
|
669,47
|
68,88
|
Netz D
|
643,18
|
-
|
692,28
|
-
|
669,47
|
-
|
Netz E
|
656,35
|
-
|
704,88
|
-
|
680,06
|
-
|
In Abbildung 34 bzw. Abbildung 35 sind die
Ergebnisse der Berechnungen aus den verschiedenen Netzen grafisch aufgetragen.
Man erkennt, dass sich die Werte entlang der unterschiedlichen Achsen immer
weiter annähern. Das bedeutet, dass mit zunehmender Netzdichte die Annahme
einer kreisförmigen Druckwellenausbreitung immer besser erfüllt wird.
Abbildung 34: Spitzenüberdruck abhängig von Netzdichte, 5 m
Abbildung 35: Spitzenüberdruck abhängig von Netzdichte, 15 m
Der in Abbildung 36 dargestellte Verlauf der Druckwelle
unter Verwendung von Netz B zeigt die besonders zu Beginn der Explosion nicht
ideal kreisförmige Ausbreitung. Mit Voranschreiten der Explosion nimmt die
Ausbreitung eine immer bessere Kreisform an. Auch hier ist die Sogphase
deutlich zu erkennen.
Abbildung 36: Plot Druckverlauf, 2-D Modell A, Netz B
6.4.1.2
Variation des Netzaufbaus
(Zoning)
Für die im Modell A gewählte Geometrie „Box“
bietet Autodyn neben der Wahl der Anzahl der Elemente über eine Kantenlänge die
Möglichkeit, diese in bestimmten Bereichen weiter zu verdichten. In Abbildung 37 ist das zugehörige Eingabefeld
dargestellt.
Abbildung 37: Eingabefeld
Netz "Box"
Man erkennt, dass durch diese Einstellungen eine
höhere Netzdichte im Bereich des Detonationsursprungs geschaffen werden kann. Man
muss allerdings beachten, dass eine Verdichtung auf der einen Seite eine
gegensätzliche Wirkung auf der anderen Seite nach sich zieht. Daher ist
gegebenenfalls die totale Anzahl der Zellen anzupassen. Auch hier wurde zur
Verringerung der Rechenzeit auf die Messreihe in 15 m
Entfernung verzichtet und mit einem Modell mit den Maßen 20 m mal 10 m gerechnet.
Bei der Wahl der Netze wurde zunächst ein
Bereich von 1 m um die Explosion ()
verfeinert. Die Anzahl der Zellen wurde entsprechend angepasst. Zur
Untersuchung der Auswirkungen wurden folgende Netze erstellt:
Netz A:
Cells I = 810, Cells J = 410, dx = dy = 20 mm, nI = nJ =
50
Netz B:
Cells I = 860, Cells J = 460, dx = dy = 10 mm, nI = nJ =
100
Netz C:
Cells I = 960, Cells J = 560, dx = dy = 5 mm, nI = nJ =
200
In einem weiteren Schritt wurde der zu
verfeinernde Bereich um die Explosion vergrößert. Hierzu wurde eine Rechnung
durchgeführt.
Netz D: Cells I = 1120, Cells J = 720, dx = dy = 5 mm, nI
= nJ = 400
Die Ergebnisse der Berechnungen sind in Tabelle 11 abgetragen.
Tabelle 11: Ergebnisse Variation Netzaufbau (Zoning),
Modell A, 2-D sphärisch
Netz
|
(5 m)
horizontal
[kPa]
|
(5 m)
45°
[kPa]
|
(5 m)
vertikal
[kPa]
|
Netz A
|
646,36
|
699,80
|
670,59
|
Netz B
|
671,00
|
689,54
|
695,33
|
Netz C
|
627,73
|
674,89
|
685,79
|
Netz D
|
654,19
|
690,97
|
699,57
|
Betrachtet man in Abbildung 38 den Spitzenüberdruck abhängig von der
Netzverfeinerung dx bzw. dy, so lassen sich aus den Ergebnissen nur
eingeschränkt Schlüsse über die Auswirkung der Netzverfeinerung ziehen.
Abbildung 38:
Spitzenüberdruck abhängig vom Netzaufbau (Zoning), 5 m
Die Untersuchungen an Modell A lassen vermuten,
dass das für eine „Box“ zur Verfügung stehende Netz ob mit oder ohne Zoning
nicht ideal ist, um kreisförmige Druckwellenausbreitungen zu simulieren. Daher
wurde entschieden, ein weiteres, kreisförmiges Modell zu untersuchen.
6.4.2
Modell B „Circle“
Da sich die Druckwelle einer freien Explosion im
Idealfall kugelförmig ausbreitet, wurde bei Modell B ein kreisförmiges Modell
erstellt. Hierdurch entfallen die für den Druckverlauf prinzipiell
irrrelevanten Elemente in den oberen Eckbereichen. Das Ziel ist es hierbei, die
Rechengeschwindigkeit zu steigern. Des Weiteren wurde wie bei Modell A, unter
Ausnutzung der Symmetrie, nur das halbe System erstellt. Der Aufbau in einer
kreis- bzw. halbkreisförmigen Struktur bietet außerdem neue Möglichkeiten in
der Netzmodellierung. Abbildung 39 zeigt
den Modellaufbau schematisch, Abbildung 40 den Rotationskörper.
Abbildung 39: 2-D Modell B, schematisch
Abbildung 40: 2-D Modell B, rotiert
6.4.2.1
Variation des Netztyps
Für einen „Circle“ gibt es in Autodyn zwei
verschiedene Netztypen. Diese sind in Abbildung 41 schematisch dargestellt.
Abbildung 41: 2-D Modell B, Netztypen
Zur Untersuchung der Unterschiede beider Typen
wurden für beide folgende Netze untersucht:
Netz A: 100
Zellen über Radius (10 m) à Zellgröße 10,00 cm
Netz B: 250
Zellen über Radius (10 m) à Zellgröße 4,00 cm
Netz C:
500 Zellen über Radius (10 m) à Zellgröße 2,00 cm
Netz D: 800
Zellen über Radius (10 m) à Zellgröße 1,25 cm
Die Ergebnisse sind in Tabelle 12 angegeben sowie in Abbildung 42 grafisch dargestellt.
Tabelle
12: Ergebnisse Variation Netztyp,
Modell B, 2-D sphärisch
Netz
|
(5 m)
horizontal
[kPa]
|
(5 m)
45°
[kPa]
|
(5 m)
vertikal
[kPa]
|
|
Typ 1
|
Typ 2
|
Typ 1
|
Typ 2
|
Typ 1
|
Typ 2
|
Netz A
|
548,11
|
562,24
|
668,75
|
693,41
|
571,67
|
578,19
|
Netz B
|
610,23
|
611,53
|
688,97
|
709,13
|
633,86
|
647,74
|
Netz C
|
655,74
|
656,31
|
711,49
|
706,73
|
675,63
|
680,87
|
Netz D
|
687,23
|
686,90
|
729,25
|
724,75
|
699,52
|
703,36
|
Abbildung 42: Spitzenüberdruck abhängig von Netztyp, 5 m
Es zeigt sich, dass die Unterschiede beider
Netze jeweils sehr gering sind und mit zunehmender Netzdichte immer ähnlichere
Ergebnisse liefern. Die Werte sind betragsmäßig größer als die aus Modell A. Da
somit eine weitere Annäherung an die Werte aus dem 1-D Modell erfolgt ist, wird
Modell B als besser geeignet beurteilt. Außerdem erkennt man wie schon in Kap.
6.4.1.1, dass die Ergebnisse entlang
der verschiedenen Achsen mit zunehmender Netzdichte immer weniger voneinander
abweichen. Dennoch bleibt auch unter Verwendung des „Circle“ die beim 1-D Modell festgestellte Konvergenz der Ergebnisse aus. Es wird Vermutet,
dass hierzu eine größere Zahl an Elementen benötigt wird. Die Untersuchungen
hierzu konnten jedoch im Rahmen dieser Studienarbeit, auf Grund der hohen
Rechenzeiten, nicht durchgeführt werden. In Abbildung 43 ist ebenfalls die Annäherung der
Kreisform mit voranschreitender Zeit gut zu erkennen.
Abbildung 43: Plot Druckverlauf, 2-D Modell B, Netz B
6.4.3
Remapping
Sowohl in Kap. 6.4.1 als auch in Kap. 6.4.2 wurden nicht die erwünschten Ergebnisse,
d.h. eine ideal kugelförmige Ausbreitung der Druckwelle sowie konvergierende
Ergebnisse mit Zunahme der Netzdichte erreicht. In diesem Abschnitt wird nun
die Funktion des „Remappings“ untersucht. „Remapping“ ist die Überführung eines
1-D Modells in ein 2-D Modell bzw. ein 2-D Modell in ein 3-D Modell. Konkret
wird in diesem Abschnitt die sphärische Ausbreitung des 1-D Modells in ein 2-D Modell
implementiert. Hierzu wurde im 1-D Modell (Netz F, vgl. Kap. 6.3.1) nach 1 ms ein so genanntes „Datafile“ geschrieben. In diesem „Datafile“ sind
alle Ergebnisse der Simulation zum Zeitpunkt t = 1 ms gespeichert.
Der Vorteil der Verwendung des „Remappings“
ist, dass man die hoch dynamischen Vorgänge der ersten Millisekunde in einem
sehr feinen 1-D Modell berechnet und dann erst in ein weniger stark
verfeinertes 2-D Modell überführt. Dabei ist jedoch darauf zu achten, dass der
Bereich in den das 1-D Modell implementiert wird auch den Randbedingungen
entspricht, unter denen das 1-D Modell erstellt wurde. Die Funktion des „Remappings“
wurde an Modell A sowie Modell B untersucht. Dabei wurden folgende Netze betrachtet:
Modell A:
Netz A-(A): Zellgröße = 10,0 cm
Netz B-(A): Zellgröße = 2,50 cm
Netz C-(A): Zellgröße = 1,60 cm
Netz D-(A): Zellgröße = 1,25 cm
Modell B:
Netz A-(B): 200 Zellen über Radius (10 m) à Zellgröße 5,00 cm
Netz B-(B): 400 Zellen über Radius (10 m) à Zellgröße 2,50 cm
Netz C-(B): 800 Zellen über Radius (10 m) à Zellgröße 1,25 cm
Netz D-(B): 1250 Zellen über Radius (10 m) à Zellgröße 0,80 cm
Der resultierende Spitzenüberdruck in 5 m Entfernung entlang der verschiedenen Messachsen ist in Tabelle 13 angegeben sowie in Abbildung 44 und Abbildung 45 grafisch dargestellt.
Tabelle
13: Ergebnisse Remapping, 2-D
sphärisch
Netz
|
(5 m)
horizontal
[kPa]
|
(5 m)
45°
[kPa]
|
(5 m)
vertikal
[kPa]
|
Netz A-(A)
|
536,08
|
648,77
|
585,26
|
Netz B-(A)
|
688,81
|
700,33
|
693,96
|
Netz C-(A)
|
708,75
|
724,07
|
711,42
|
Netz D-(A)
|
722,29
|
735,37
|
723,94
|
Netz
A-(B)
|
644,68
|
707,15
|
652,46
|
Netz
B-(B)
|
688,96
|
705,71
|
685,52
|
Netz
C-(B)
|
722,27
|
735,77
|
720,74
|
Netz
D-(B)
|
737,75
|
747,23
|
739,53
|
Abbildung 44: Spitzenüberdruck Modell A, abhängig von Netz
mit Remapping , 5 m
Abbildung 45: Spitzenüberdruck Modell B, abhängig von Netz
mit Remapping , 5 m
In beiden Modellen wird für ein feines Netz eine
quasi ideal kugelförmige Ausbreitung der Druckwelle erreicht. So beläuft sich
die Abweichung des maximalen Überdrucks der verschiedenen Messpunkte bei Modell
A auf maximal 1,8 % und bei Modell B auf 1,3 %.
Des Weiteren lässt sich beobachten, dass die
maximalen Überdrücke sich immer stärker den Berechnungen des 1-D Modells angleichen. Eine Konvergenz gegen einen oberen Grenzwert bleibt
dennoch aus. Auch hier wird, wie in Kapitel 6.4.2.1, der Grund in der zu groben Netzeinteilung
gesehen.
6.5 2-D
Modell – hemisphärische Ausbreitung, ohne Reflexion
Nachdem in den vorherigen Abschnitten dieses
Kapitels nur sphärische Ausbreitungen betrachtet wurden, soll an dieser Stelle
die hemisphärische Ausbreitung der Druckwelle untersucht werden. Bei den hier
durchgeführten Simulationen soll eine Sprengladung 1 m über
der Erdoberfläche betrachtet werden. Dieser Aufbau entspricht beispielsweise
dem Szenario Bombe in Kleintransporter. Die Erdoberfläche wird als ideal
starre, ideal ebene und unendlich große Fläche angenommen.
Das in Abbildung 46 schematisch dargestellte Modell besteht
aus einem Viertelkreis dessen waagrechte Kante die Symmetrieachse bildet und
dessen senkrechte Kante die Reflexionsebene darstellt.
Abbildung 46: 2-D Modell, hemisphärisch, schematisch
Das System muss man sich als an der unteren
Kante gespiegelt und um 90° gegen den Uhrzeigersinn gedreht vorstellen, um den
hemisphärischen Fall vorliegen zu haben.
Da nun der maximale Druck in Oberflächennähe
zu erwarten ist, wurden die Messpunkte in 5 m und
15 m Entfernung vom Detonationsursprung und in einem Abstand zur
Reflexionsebene von 1 m angebracht.
In dieser Simulation soll ebenfalls die
Funktion des Remappings verwendet werden. Hierzu wurde aus dem 1-D Modell (Netz
F, 100 kg) eine „Datafile“ nach 0,175 ms geschrieben. Zu diesem Zeitpunkt
beläuft sich die Druckwellenausbreitung noch unter 1 m und
kann somit im 2-D Modell verwendet werden.
Es wurde auf Grund annehmbarer Rechenzeiten folgendes
Netz verwendet:
Netz A: 800 Zellen über Radius (20 m) à Zellgröße 2,50 cm
Tabelle 14: Ergebnisse, 2-D hemisphärisch
Netz
|
(5 m)
[kPa]
|
(15 m)
[kPa]
|
Netz A
|
1348,20
|
112,17
|
Betrachtet man nun den in Abbildung 47
gegebenen Druck-Zeit Verlauf nach N. Gebbeken [21], so fällt auf, dass dort ein Spitzenüberdruck von ca. 110 bis 120 kPa
für die freie hemisphärische Ausbreitung auftritt. Der mit Autodyn simulierte
Wert beträgt 112,17 kPa. Daraus lässt sich schließen, dass die hier gewählten Einstellungen
gut für die hemisphärische Ausbreitung funktionieren.
Abbildung 47: Druck-Zeit Verlauf, 100 kg, 15 m, nach [21]
6.6
3-D Modell – hemisphärische Ausbreitung
Nachdem eine ausführliche Betrachtung der 1-D
und 2-D Modelle erfolgte, wird nun beispielhaft eine 3-D Simulation
durchgeführt. Die Sprengladung ist in einem Fahrzeug und befindet sich somit 1 m über der Erdoberfläche. Das Fahrzeug wird jedoch nicht weiter
berücksichtigt. Als Umgebung wurde der Campus Lichtwiese Abschnitt L5 der TU
Darmstadt gewählt. Abbildung 48 zeigt
ein Luftbild der Umgebung. Als Sprengladung wurden 100 kg TNT
gewählt. Die Auswertung der Druckwelle erfolgt beispielhaft an einem Büro im
sechsten Stockwerk des Gebäudes L5|06 (MP1) sowie an einer Gebäudekante des Gebäudes
L5|07 (MP2).
Abbildung 48: Luftbild,
Campus Lichtwiese L5 Bauingenieurwesen [22]
Bei der 2-D Betrachtung wurde bereits die
Funktion des Remappings eingeführt. Da eine genaue Rechnung eine sehr hohe
Netzdichte und somit sehr hohe Rechenzeiten erfordert, soll das Remapping auch
hier angewendet werden. In Abbildung 49 ist der Vorgang des Remappings schematisch dargestellt. Zunächst wird
das 1-D Modell (Netz F, vgl. Kap. 6.3.1) in ein 2-D Modell (hemisphärische
Ausbreitung) überführt. Hierzu kann jedoch nicht das Modell aus Kap. 6.5 verwendet werden, da dort ein „Circle“ zum
Einsatz kam. Bei dem Versuch eine solche Geometrie in ein 3-D Modell zu führen,
entsteht eine Fehlermeldung. Daher wurde eine weitere Simulation der
hemisphärischen Ausbreitung mittels einer „Box“ durchgeführt. Die Ergebnisse
wurden bei einer Zeit von 1,5 ms in eine „Datafile“ geschrieben. Anschließend wurde das 3-D Modell
erstellt. Hierbei wurde zunächst die gesamte Umgebung als eine „Box“ mit den
Maßen 160 m x 100 m x 40 m als Luft gefüllter Raum erstellt. Anschließend wurden die Bereiche,
in denen sich die Gebäude befinden, als „unused Part“ definiert. Hierdurch
wurden die Geometrien der Gebäude stark vereinfacht und als ideal starr angenommen.
Abbildung 49: Remapping schematisch
Bei dem Remapping eines 2-D Modells in ein 3-D
Modell ist darauf zu achten, dass die Rotationsachse des 2-D Modells in den
Koordinatenursprung des 3-D Modells fällt. Gegebenenfalls ist eine Translation
des Koordinatensystems notwendig. In Abbildung 50 ist der Modellaufbau dargestellt.
Abbildung 50: Modellaufbau, 3-D Modell
An allen Grenzflächen, mit Ausnahme der
Bodenfläche, die als starr definiert wurde, wurde als Randbedingung ein „Flow-Out“
definiert. Die Zellgröße wurde auf Grund der Größe des Modells relativ grob mit
1 m eingestellt. Die maximalen Überdrücke an den Messpunkten sind in Tabelle
15 angegeben.
Messpunkt 1 befindet sich in einer Entfernung von ca. 29 m und Messpunkt 2 in ca. 14 m zum
Detonationsursprung.
Tabelle
15: Ergebnisse, 3-D
Messpunkt
|
[kPa]
|
1
|
44,62
|
2
|
113,96
|
Der Druck-Zeit Verlauf an beiden Messpunkten
ist in Abbildung 51 angegeben. Auch hier ist der typische Verlauf, wie er in
Kap. 3.3 angegeben wurde, zu
erkennen.
Abbildung 51: 3-D Modell, Druck-Zeitverlauf
Da der Druckverlauf nur in 1-D bzw. 2-D Plots
anschaulich darstellbar ist, sind in Abbildung 53 die Geschwindigkeits-Vektoren
zu verschiedenen Zeitpunkten t dargestellt. Anhand dieser Darstellung lässt
sich zwar keine quantitative Aussage über den Druckverlauf machen, aber wie man
in Abbildung 52 erkennt, sind die
Verläufe für Druck und Geschwindigkeit qualitativ sehr ähnlich, daher
beschreiben die Plots qualitativ den Druckwellenverlauf sehr gut. In den Plots
für t = 145,10 ms sowie t = 314,60 ms erkennt man zudem, dass zwei Schockwellen entstanden sind. Die erste
wurde durch die eigentliche Detonation und die zweite durch die Reflexion der
ersten an dem Gebäude L5|06 verursacht.
Abbildung 52: 3-D Modell, Druck-Zeit- und
Geschwindigkeit-Zeit-Verlauf
Abbildung 53: 3-D Modell, Plot Geschwindigkeitsvektoren zum
Zeitpunkt t
Die beispielhafte Berechnung am 3-D Modell
ermöglicht nur erste Eindrücke und Erfahrungen mit der Simulation von 3-D
Szenarien. Weiterführend sind auch hier die verschiedenen Parameter wie z.B.
die Zellgröße zu untersuchen. Des Weiteren besteht auch die Möglichkeit ein 1-D
Modell direkt in ein 3-D Modell zu überführen, diese Funktion sollte noch
weiter betrachtet werden. Ziel hierbei ist es konvergente, und damit numerische
genaue Ergebnisse bei möglichst geringem Rechenaufwand zu erhalten.
6.7 Vergleich
der Ergebnisse
In diesem Kapitel sollen die Ergebnisse der numerischen
Simulationen mit den sonstigen in dieser Arbeit behandelten Berechnungsansätzen
verglichen werden.
Betrachtet man die sphärische Ausbreitung der
Druckwelle für 100 kg TNT-Äquivalent, so wurde dieses Szenario zum einen über verschiedene,
einfache Berechnungsverfahren sowie numerisch in einem 1-D sowie 2-D Modell
berechnet. Die Ergebnisse sind auszughaft in Tabelle 16 aufgetragen.
Tabelle
16: Vergleich der Ergebnisse
sphärische Ausbreitung
Berechnungsansatz
|
In Kapitel
|
(5 m)
[kPa]
|
(15 m)
[kPa]
|
Brode
|
3.1.1
|
682,00
|
59,54
|
Henrych
|
3.1.3
|
673,51
|
69,01
|
1-D (Netz F)
|
6.3.1
|
759,31
|
66,28
|
2-D (Netz E)
|
6.4.1.1
|
680,43
|
-
|
2-D (Netz D)
|
6.4.2.1
|
705,33
|
-
|
Bei der Betrachtung der Ergebnisse stellt man
fest, dass der kleinste und größte ermittelte Wert für den Spitzenüberdruck in
5 m Entfernung nur um rund 11 %
voneinander abweichen. Diese Abweichung wird hinsichtlich der Vielzahl an
unkalkulierbar oder schwer kalkulierbaren Randbedingungen als plausible
Streuung betrachtet. Bei den vorliegenden Ergebnissen in 15 m Entfernung beläuft sich die Abweichung auf maximal
ca. 15 %.
Bei dem Vergleich der Ergebnisse einer
hemisphärischen Druckwellenausbreitung ohne Reflexion liegen sowohl Ergebnisse
aus den einfachen Berechnungen, aus AT-Blast sowie der numerischen Simulation
vor. Diese sind in Tabelle 17 angegeben.
Tabelle
17: Vergleich der Ergebnisse hemisphärische
Ausbreitung
Berechnungsansatz
|
In Kapitel
|
(5 m)
[kPa]
|
(15 m)
[kPa]
|
Brode
|
3.1.4
|
1227,60
|
107,16
|
AT-Blast
|
4.1.1
|
1157,22
|
99,08
|
2-D Modell
|
6.5
|
1348,20
|
112,17
|
Auch bei diesen Berechnungen betragen die Abweichungen maximal 16 %. Bei
der hemisphärischen Ausbreitung kommt zur Problematik der Erfassung der
Randbedingungen noch das Problem des unklaren Versuchsaufbaus hinzu. Es ist
z.B. bei AT-Blast nicht klar nachzuvollziehen, in welcher Höhe sich die
Sprengladung über der Erdoberfläche befindet. Darüber hinaus ist die Anordnung
der Messpunkte (in Achse des Detonationsursprungs oder an der Erdoberfläche)
unklar. Schließlich wird wie in Kapitel 4.1.2 erwähnt die pauschale Umrechnung nach Mays
als unzureichend erachtet.
7
Fazit
Schon vor über 50 Jahren beschäftigte man sich
eingehend mit dem Thema der Berechnung von Explosionsereignissen. In dieser
Zeit entstand bereits eine Vielzahl an Berechnungsansätzen, deren Ergebnisse zum
einen der Spitzenüberdruck der Explosionsdruckwelle, bei sphärischer und
hemisphärischer Ausbreitung, und zum anderen der maximale positive Impuls sind.
Auf Grund der vielfältigen Ansätze, die alle
einen skalierten Faktor (Z), welcher die Entfernung zur Sprengladung sowie
deren Masse berücksichtigt, galt es diese miteinander zu vergleichen. Durch diesen
Vergleich wurde festgestellt, dass gerade im Nahbereich unterschiedlichen
Ergebnissen für den Spitzenüberdruck bei sphärischer Ausbreitung erzielt werden.
Mit zunehmender Entfernung gleichen sich die Ergebnisse an. Das betrachtete Verfahren
zur Berechnung hemisphärischer Ausbreitungen vollzieht lediglich eine
Umrechnung der Ergebnisse der sphärischen in die hemisphärische Ausbreitung
mittels eines konstanten Faktors. Diese Pauschalisierung wird jedoch kritisch
gesehen und sollte lediglich als erste Näherung betrachtet werden. Der maximale
positive Impuls lässt sich ebenfalls unter Verwendung des Faktors Z berechnen.
Mit der vereinfachten Annahme eines linearisierten Druck-Zeit-Verlaufs ist der
Impuls leicht zu ermitteln. Diese wird auch in technischen Regelwerken angewendet.
Für die Berücksichtigung von Reflexionen an
starren Flächen wurden zum einen Formeln zur senkrechten Reflexion sowie zum
anderen zwei Diagramme zur nicht senkrechten Reflexion untersucht. Diese
Verfahren liefern sehr gute Ergebnisse.
Durch die ersten Betrachtungen konnte ein
Eindruck von der Entstehung sowie der Ausbreitung einer Druckwelle in Folge
einer Explosion erlangt werden. In weiteren Schritten wurden dann die
vereinfachten Ergebnisse unter Verwendung von den Programmen AT-Blast sowie
ANSYS Autodyn verglichen.
Da die bisher erwähnten Berechnungen
prinzipiell alle von einer sphärischen Druckwellenausbreitung ausgingen, mit
Ausnahme der genannten konstanten Umrechnung, war es zudem interessant, die
Ergebnisse für den hemisphärischen Fall aus AT-Blast zu erhalten. AT-Blast
benötigt als Eingangsgrößen, ähnlich wie die oben beschriebenen Verfahren,
lediglich die Masse der Sprengladung sowie den Abstand zur Sprengladung.
Außerdem lässt sich der reflektierte Druck für einen beliebigen Reflexions- bzw.
Anströmwinkel ermitteln. Der genaue Rechenalgorithmus des Programms AT-Blast
scheint ebenfalls auf analytischen Berechnungsverfahren, wie sie auch in dieser
Arbeit beschrieben werden, zu beruhen.
Auf Grund seiner Einfachheit ist AT-Blast auf
eine kleine Zahl von Szenarien, nämlich eine hemisphärische Ausbreitung mit
Einfachreflexion, beschränkt. Dennoch liefert es für diese Sachverhalte
hinreichend genaue Ergebnisse. Für eine Vordimensionierung sind die Ergebnisse
ausreichend und sollten dort Verwendung finden. Dies sieht man allein daran,
dass die in AT-Blast errechneten Werte für den Spitzenüberdruck sowie den
Impuls den in den Normen angegeben Werten entsprechen.
Als Grundlage zu den folgenden numerischen
Berechnungen wurden verschiedene Verfahren der expliziten Zeitintegration sowie
die Unterschiede zur impliziten Zeitintegration aufgezeigt. Der Grundlegende
Unterschied liegt darin, dass in expliziten Verfahren das Gleichgewicht zum
bekannten Zeitpunkt und bei den impliziten zum unbekannten, gesuchten Zeitpunkt
aufgestellt wird. Ein weiteres Merkmal der expliziten Zeitintegration ist die
Tatsache, dass die maximale Schrittweite pro Integrationsschritt beschränkt ist,
bei der impliziten ist dies nicht der Fall. Um bei der Berechnung von stark
nichtlinearen dynamischen Vorgängen zu einem Ergebnis zu kommen, müssen
explizite Verfahren, wie das z.B. in Autodyn eingesetzte Mehrschrittverfahren,
verwendet werden.
In der numerischen Simulation konnten nun für
alle vorangegangenen Berechnungen Vergleichswerte ermittelt werden. Hierzu
wurde zunächst eine eingehende Parameterstudie anhand eines 1-D Modells
durchgeführt. Es wurden die Parameter Netzdichte, maximaler
Zeitschritt sowie das Materialmodell des Sprengstoffs untersucht.
Die Berechnungen am 1-D Modell können nur zur Simulation von sphärischen Druckwellenausbreitungen
verwendet werden. Es konnte eine Konvergenz der Ergebnisse bei zunehmender
Netzdichte sowie bei der Begrenzung des Zeitschritts festgestellt werden. Ebenfalls
sphärische Untersuchungen wurden in zwei unterschiedlichen 2-D Modellen durchgeführt.
Hierbei lag das Interesse darin, den Einfluss der verschiedenen Netzgeometrien
der Modelle zu erfassen. Eine weitere, zunächst am 2-D Modell betrachtete, Funktion
ist die des Remappings. Hierbei werden Ergebnisse aus Modellen niederer
Dimension in ein Modell höherer Dimension überführt. Hierdurch konnte eine
Steigerung der Genauigkeit sowie eine Verringerung der Rechenzeit erreicht
werden.
Die hemisphärische Druckwellenausbreitung
erfolgte ebenfalls am 2-D Modell. Durch Einbringen einer Randbedingung wurde
eine ideal starre Reflexionsfläche geschaffen. Der anschließende Vergleich der
Ergebnisse mit anderen Simulationen bzw. mit in Versuchen ermittelten Werten
erlaubt es, die Simulation als sehr gut einzustufen.
In einem weiteren Schritt wurde nun
beispielhaft in einer realen Umgebungsgeometrie eine Simulation an einem 3-D
Modell durchgeführt. Auch an dieser Stelle wurde das Remapping verwendet.
Hierdurch konnte das Netz im 3-D Modell gröber eingestellt und somit Rechenzeit
eingespart werden. In einem 3-D Modell können komplexe Sachverhalte wie z.B.
die Mehrfachreflexion bei enger Bebauung berücksichtigt werden.
Im Allgemeinen hat sich gezeigt, dass die
Ergebnisse stark von der Wahl des Modells sowie der darin verwendete
Netzfeinheit abhängen. Das Modell sollte immer nur soweit „verkompliziert“
werden, wie es nötig ist. So ist zur Berechnung einer freien, sphärischen
Druckwellenausbreitung das 1-D Modell ausreichend und unter der Annahme von
Idealbedingungen mit der größten Genauigkeit behaftet. Möchte man jedoch den
hemisphärischen Fall berücksichtigen, muss ein 2-D Modell verwendet werden. Bei
komplexerer Fragestellung ist in der Regel ein 3-D Modell notwendig. ANSYS
Autodyn bietet einfache Möglichkeiten zur Simulation von solchen
Explosionsereignissen. Dabei sollte zur Steigerung der Genauigkeit die
vorgestellte Funktion des Remappings stets Verwendung finden. Bei den durchgeführten
Berechnungen konnte schnell festgestellt werden, dass die benötigte Rechenzeit stark
mit der Feinheit der Modelle zusammenhängt.
Abschließend lässt sich sagen, dass alle in
dieser Arbeit behandelten Berechnungsverfahren für einfache Explosionsszenarien,
d.h. eine überschaubare Anzahl von Randbedingungen, ähnliche Ergebnisse liefern. Möchte man jedoch komplexere
Szenarien nachbilden, so sind die einfachen Berechnungsverfahren nach
beispielsweise Brode [4]
oder Henrych [5]
nicht mehr ausreichend. Für solche Szenarien wird eine Simulation mit einem
numerischen Berechnungsprogramm zwingend erforderlich sein. Bei der Berechnung
von großen, aufwendigen Modellen sollte zudem über den Einsatz eines
Hochleistung-Rechners nachgedacht werden.
8
Ausblick
Aufbauend auf diese Arbeit sollten noch
weitere, an dieser Stelle noch ungeklärte Sachverhalte betrachtet werden.
Da bei der manuellen Begrenzung des
Zeitschritts zwar eine Konvergenz der Ergebnisse, jedoch nicht dieselbe wie bei
der Netzverfeinerung festgestellt wurde, sollten hierzu noch weitere
Untersuchungen angestellt werden.
Die Untersuchungen am 3-D Modell wurden im
Rahmen dieser Arbeit nur anhand eines Modells beispielhaft durchgeführt. Hier
sollten auch die Auswirkungen durch eine Veränderung der Netzdichte genauer betrachtet
werden.
Es wäre außerdem interessant, sich eingehender
mit dem Technical Manual 5-855-1 [13] zu befassen. Dies ist eine der meist zitierten Quellen im Zusammenhang
mit Explosionsereignissen. Hierdurch wären möglicherweise die genauen
Rechenoperationen in AT-Blast nachzuvollziehen.
In Autodyn besteht neben den hier betrachteten
Ansätzen zur Simulation einer Explosionsdruckwelle noch die Möglichkeit, die
Beanspruchung auf eine Struktur durch eine Explosion automatisch zu berechnen.
Hierbei wird eine Randbedingung generiert, in der die Koordinaten einer
Sprengladung sowie deren Masse hinterlegt sind. Nach dem Aufbringen dieser Randbedingung
auf die Struktur wird diese dem fiktiv entstandenen Druckverlauf ausgesetzt. Die
Grundlage dieser Berechnung ist ebenfalls TM 5-855-1. Dieses Verfahren sollte
getestet und mit den bereits ermittelten Werten verglichen werden.
Nachdem in dieser Arbeit ausschließlich die Einwirkungsseite
betrachtet wurde, ist nun der nächste Schritt die Widerstände zu betrachten. Hierzu
gibt es zwei Möglichkeiten:
1. Man ermittelt die Einwirkungen anhand der hier vorgestellten Verfahren
und überträgt diese als Belastung auf ein beliebiges System bzw. Material.
Da die Druckwelle an einer
ideal starren Fläche ermittelt wurde entfällt in diesem Fall jedoch die
Interaktion zwischen Bauteil und Druckwelle.
2. Man ersetzt die ideal starren Flächen durch die zu untersuchende
Struktur.
Hierdurch wird die
Bauteil-Druckwellen-Interaktion berücksichtigt.
Abbildungsverzeichnis
Abbildung 1:
Spitzenüberdruck nach Brode, hergeleitet und nach [2],[3]. 5
Abbildung 2:
Spitzenüberdruck nach Brode und Henrych. 6
Abbildung 3: Weitere
sphärische Berechnungen im Vergleich. 7
Abbildung 4:
Hemisphärische Druckwellenausbreitung. 8
Abbildung 5: Entstehung
des Staudruck, schematisch. 9
Abbildung 6: Entwicklung
einer Stoßfront nach Cooper [10]. 9
Abbildung 7:
Reflektierter Spitzenüberdruck nach Rankine & Hugonoit 10
Abbildung 8:
Reflexionsfaktor Cr. 11
Abbildung 9:
Reflexionsfaktor nach Schindler [12]. 11
Abbildung 10:
Reflexionsfaktor nach TM 5-855-1 [13]. 13
Abbildung 11:
Überdruck-Zeit-Verlauf, schematisch. 13
Abbildung 12:
Linearisierter Druck-Zeitverlauf, 100 kg TNT – 15 m Entfernung. 15
Abbildung 13: Vergleich
Brode - AT-Blast, hemisphärisch, ohne Reflexion, 100 kg TNT-Äquivalent 18
Abbildung 14: Vergleich
ISO, Brode und AT-Blast, hemisphärisch, mit Reflexion. 20
Abbildung 15:
Hemisphärische Ausbreitung Brode 200 kg. 21
Abbildung 16: Impuls
nach Kinney & Graham und AT-Blast 22
Abbildung 17:
Ein-Freiheitsgrad-Schwinger. 23
Abbildung 18:
Integration mit der zentralen Differenzenmethode. 25
Abbildung 19:
Integration mit dem Heunschen Verfahren. 27
Abbildung 20:
Materialdaten TNT und TNT-2. 31
Abbildung 21:
Vorgehensweise numerische Simulation der sphärischen Druckwellenausbreitung. 32
Abbildung 22: 1-D Modell
schematisch. 33
Abbildung 23: 1-D Modell
rotiert 33
Abbildung 24: 1-D
Modell, Druck-Zeitverlauf aus Autodyn, Netz F, 5 m... 35
Abbildung 25: Plot
Druckverlauf, 1-D Modell, Netz C.. 36
Abbildung 26: 1-D
Modell, Spitzenüberdruck abhängig von Netzdichte, 5 m... 36
Abbildung 27: 1-D
Modell, Spitzenüberdruck abhängig von Netzdichte, 15 m... 37
Abbildung 28: 1-D
Modell, Spitzenüberdruck abhängig von Zeitschritt, 5 m... 38
Abbildung 29: 1-D
Modell, Spitzenüberdruck abhängig von Zeitschritt, 15 m... 38
Abbildung 30: 1-D
Modell, Vergleich TNT mit TNT-2, 5 m... 39
Abbildung 31: 1-D
Modell, Vergleich TNT mit TNT-2, 15 m... 39
Abbildung 32: 2-D Modell
A, schematisch. 40
Abbildung 33: 2-D Modell
A, rotiert 270°. 41
Abbildung 34:
Spitzenüberdruck abhängig von Netzdichte, 5 m... 42
Abbildung 35:
Spitzenüberdruck abhängig von Netzdichte, 15 m... 42
Abbildung 36: Plot
Druckverlauf, 2-D Modell A, Netz B.. 43
Abbildung 37:
Eingabefeld Netz "Box". 44
Abbildung 38:
Spitzenüberdruck abhängig vom Netzaufbau (Zoning), 5 m... 45
Abbildung 39: 2-D Modell
B, schematisch. 46
Abbildung 40: 2-D Modell
B, rotiert 46
Abbildung 41: 2-D Modell
B, Netztypen. 47
Abbildung 42:
Spitzenüberdruck abhängig von Netztyp, 5 m... 48
Abbildung 43: Plot
Druckverlauf, 2-D Modell B, Netz B.. 49
Abbildung 44:
Spitzenüberdruck Modell A, abhängig von Netz mit Remapping , 5 m... 50
Abbildung 45:
Spitzenüberdruck Modell B, abhängig von Netz mit Remapping , 5 m... 51
Abbildung 46: 2-D
Modell, hemisphärisch, schematisch. 52
Abbildung 47: Druck-Zeit
Verlauf, 100 kg, 15 m, nach [21]. 53
Abbildung 48: Luftbild,
Campus Lichtwiese L5 Bauingenieurwesen [22]. 54
Abbildung 49: Remapping
schematisch. 55
Abbildung 50:
Modellaufbau, 3-D Modell 56
Abbildung 51: 3-D
Modell, Druck-Zeitverlauf 57
Abbildung 52: 3-D
Modell, Druck-Zeit- und Geschwindigkeit-Zeit-Verlauf 57
Abbildung 53: 3-D
Modell, Plot Geschwindigkeitsvektoren zum Zeitpunkt t 58
Tabellenverzeichnis
Tabelle 1:
Umrechnungsfaktoren für explosives Material in TNT-Äquivalent nach [2]. 2
Tabelle 2: Formfaktor α nach [6]. 14
Tabelle 3: Vergleich
Brode - AT-Blast, hemisphärisch, ohne Reflexion, 100 kg TNT-Äquivalent 18
Tabelle 4: Vergleich
Spitzenüberdruck ISO, Brode und AT-Blast, hemisphärisch, mit Reflexion. 19
Tabelle 5:
Reflexionsfaktor aus AT-Blast 20
Tabelle 6:
Hemisphärische Ausbreitung Brode 200 kg. 21
Tabelle 7: Ergebnisse
Variation Netzdichte, 1-D.. 34
Tabelle 8: Ergebnisse
Variation Zeitschritt, 1-D.. 37
Tabelle 9: Ergebnisse
Variation TNT, 1-D.. 39
Tabelle 10: Ergebnisse
Variation Netzdichte, Modell A, 2-D sphärisch. 41
Tabelle 11: Ergebnisse
Variation Netzaufbau (Zoning), Modell A, 2-D sphärisch. 44
Tabelle 12: Ergebnisse
Variation Netztyp, Modell B, 2-D sphärisch. 47
Tabelle 13: Ergebnisse
Remapping, 2-D sphärisch. 50
Tabelle 14: Ergebnisse,
2-D hemisphärisch. 52
Tabelle 15: Ergebnisse,
3-D.. 56
Tabelle 16: Vergleich
der Ergebnisse sphärische Ausbreitung. 59
Tabelle 17: Vergleich
der Ergebnisse hemisphärische Ausbreitung. 59
Literaturverzeichnis
[1] C. Mills, The design of
concrete structure to resist explosions and weapon effects. Edinburgh, 1987.
[2] Gebbeken Rutner, Mangerig
"Stahlbaukonstruktionen unter Explosionsbeanspruchung," in Stahlbau
Kalender 2008, ed, 2008.
[3] H. A. Bethe J. Neumann, Shock
hydrodynamics and blast waves: Los Alamos National Laboratory, 1944.
[4] H.L. Brode,
"Numerical solution of spherical blast waves.," Journal of Applied
Physics, 1955.
[5] Henrych, The Dynamics
of Explosion and Its Use. Prag, 1979.
[6] K. Graham G. Kinney, Explosive
shocks in air. New York, 1985.
[7] I.G. Petrovski I.A.
Naumyenko, M.A. Sadovski, The shock wave of a nuclear explosion. Moskau,
1956.
[8] P. Smith G. Mays, Blast
effects on buildings vol. 2. London, 1995.
[9] W.J.M. Rankine, On the
thermodynamic theory of waves of finite longitudinal disturbance.
Philossophical Transactions of the Royal Society, 1870.
[10] P.W. Cooper,
"Explosives Engineering," Wiley-VCH, 2000.
[11] T. Döge N. Gebbeken,
"Vom Explosionszenario zur Bemessungslast," Der Prüfingenieur 2006.
[12] G. Schindler, "Handbuch
der Waffenwirkung für die Bemessung von Schutzbauten," Bundesamt für
Zivilschutz, 1964.
[13] Defense Special Weapons
Agency and Joint Department of the Army, "Design and Analysis of Hardened
Structures to Convential Weapons Effects," TM 5-855-1, 1997.
[14] DIN EN 13123-1,
"Sprengwirkungshemmung - Anfroderung und Klassifizierung," ed, 2001.
[15] ISO-16934, "Glass in
building - Explosion-resistant security glazing," ed, 2007.
[16] M. Blaut. (
04.12.2011). Algorithmen zur Dynamik im Zeitbereich.
[17] R. Müller, "Numerische
Methoden der Mechanik," Skript zur Vorlesung, 2007.
[18] Serradeill Baudin,
"Review of Jones-Wilkins-Lee equation of state," 2010.
[19] B. M. Dobratz, LLNL
Explosive Handbook, Properties of Chemical Explosives and Explosive Simultants.
Livermore CA, 1981.
[20] Century Dynamics,
"Autodyn Remapping Tutorial," ed, 2005.
[21] T. Döge N. Gebbeken, I.
Mangerig, S. Beucher, "Von der Explosion zur Bemessungslast und zur
Brandeinwirkung durch Folgebrand," Universität der Bundeswehr München.
[22] Google Maps. (2012, 30.01.).