UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 EVALUATION OF 3 D RADIANT HEAT TRANS FER IN STREET CANYON S By TUCKER FLAGG STACHITAS A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2009 PAGE 2 2 2009 Tucker Flagg Stachitas PAGE 3 3 To the American people, m ay your light always shine brightly and freely PAGE 4 4 ACKNOWLEDGMENTS First and foremost I thank my Lord and Savior Jesus Christ for the wisdom and patience to complete this research. Furthermore, I thank my advisor Dr. Glenn Sjoden for his tireless guidance and unwavering dedication to my success. I am also deeply indebted to Todd Mock, Kevin Manalo, and the other students at the Florida Institute for Nuclear Detection and Security (FINDS) for their valuable conceptual advice and programming assistance Lastly, I thank the U nited States Navy for the opportunity to pursue my graduate degree. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ...................................................................................................... 4 LIST OF TABLES ................................................................................................................ 8 LIST OF FIGURES ............................................................................................................ 10 ABSTRACT ........................................................................................................................ 12 CHAPTER 1 INTRODUCTION ........................................................................................................ 13 2 THERMAL RADIATION FUNDAMENTALS ............................................................... 17 Thermal Partition ........................................................................................................ 17 Thermal Pulse Timing ................................................................................................. 17 Fireball Size ................................................................................................................ 19 Radiant Heat Trans fer ................................................................................................ 20 Surface Property Definitions ................................................................................ 21 Black Body Radiator Approximation .................................................................... 22 View Factor Definition .......................................................................................... 22 Radiant Ex change between Black Bodies .......................................................... 24 Gray Diffuse Emitters, Reflectors, and Absorbers .............................................. 24 Semi Infinite Solids .............................................................................................. 26 Unit Conventions: Total Power vs. Emissive Power ........................................... 28 3 MODEL CONCEPTUALIZATION ............................................................................... 31 Geometry .................................................................................................................... 31 Material Properties ...................................................................................................... 32 Source Term Considerations ...................................................................................... 33 Statio nary Plane S ource Approximation ............................................................. 34 Direct Temperature Method ................................................................................. 35 Black Body Temperature Method ........................................................................ 35 Shadowing and Ground Effects ........................................................................... 36 Atmo spheric Attenuation ............................................................................................ 37 Initial Modeling ............................................................................................................ 39 Open Field Model ................................................................................................. 40 Sparrow Passage Model ...................................................................................... 41 Initial Model Comparison ..................................................................................... 42 4 CANYON CODE DEVELOPMENT ............................................................................ 51 Reference System ...................................................................................................... 52 PAGE 6 6 View Factors ............................................................................................................... 52 Whole Surface View F actors ............................................................................... 53 Differential View Factors ...................................................................................... 54 Numerical Issues ........................................................................................................ 57 Surface Meshing Constraints ..................................................................................... 58 Temperature Initialization and Source Definition ....................................................... 58 Time Steps .................................................................................................................. 59 Flux Profile Resolution ................................................................................................ 60 Direct Inversion .................................................................................................... 60 Gauss Seidel Iteration ......................................................................................... 61 Temperatur e Updating ................................................................................................ 63 Non Reflecting Boundary Conditions ......................................................................... 63 Benchmarking ............................................................................................................. 64 Steady State Benchmark ..................................................................................... 64 Steady Sta te Integration Benchmark ................................................................... 65 Transient Benchmark ........................................................................................... 65 5 CANYON ANALYSIS & RESULTS ............................................................................ 73 Baseline Parameters Summary ................................................................................. 73 Conv ergence Studies ................................................................................................. 73 Mesh ..................................................................................................................... 74 Time Step ............................................................................................................. 74 Tolerance.............................................................................................................. 75 Transient Wall Temperatures ..................................................................................... 75 Parameter Studies Outline ......................................................................................... 76 Material Property Studies ........................................................................................... 76 Canyon Dimension Studies ........................................................................................ 77 Source T emperature Selection Study ........................................................................ 77 Canyon Effect and Model Comparisons .................................................................... 78 Yield Studies ............................................................................................................... 79 Visibility Studies .......................................................................................................... 79 6 CONCLUSIONS .......................................................................................................... 89 7 FUTURE WORK ......................................................................................................... 91 CANYON Improvements ............................................................................................ 91 Further Anal yses ......................................................................................................... 92 APPENDIX A TRANSMITTANCE CURVE FIT COEFFICIENTS .................................................... 94 B CANYON SOURCE CODE ........................................................................................ 95 C CANYON SAM PLE INPUT FILE .............................................................................. 143 PAGE 7 7 D CANYON SAMPLE OUTPUT FILE (ABRIDGED) ................................................... 145 REFERENCES ................................................................................................................ 149 BIOGRAPHICAL SKETCH .............................................................................................. 151 PAGE 8 8 LIST OF TABLES Table page 3 1 Baseline thermal properties for concrete ............................................................... 43 3 2 Fireball diameters for 2 kT and 20 kT surface bursts ........................................... 44 3 3 Temperature method comparison at peak power ................................................. 45 3 4 International Visibility Code .................................................................................... 47 3 5 Time independent OFM parameters and result .................................................... 47 3 6 Time dependent OFM parameters and results ..................................................... 48 3 7 Sparrow Passage Model parameters and results ................................................. 49 4 1 Surface definitions for CANYON ............................................................................ 67 4 2 Whole surface view factor solution method summary .......................................... 68 4 3 Whole surface view factors for debugging dimensions ......................................... 68 4 4 Differential view factors from surface 5 to 3 for debugging dimensions ............... 70 4 5 View factor method comparison for approximation of zero as small number ...... 70 4 6 Number of time steps at various yields for t of 1x105 s. .................................... 71 4 7 Steady state benchmark results for CANYON ...................................................... 71 4 8 Transient benchmark parameters .......................................................................... 72 4 9 Transient benchmark results .................................................................................. 72 5 1 Baseline parameters summary for CANYON ........................................................ 80 5 2 Results summary for mesh convergence study in CANYON ............................... 80 5 3 Results summary for time step convergence study in CANYON ......................... 81 5 4 Results summary for tolerance convergence study in CANYON ......................... 83 5 5 Selected output temperatures for the baseline case in CANYON ........................ 83 5 6 Material property sensitivity study results summary for CANYON ....................... 83 5 7 Canyon dimension studies results for CANYON ................................................... 84 PAGE 9 9 5 8 Total fluence for 20 kT, 10 kT, & 5 kT at 1610m ................................................... 86 5 9 Canyon effect for 20 kT, 10 kT, & 5 kT at 1610m ................................................. 86 5 10 Total fluence for sea level visibilities of 5, 10, 15, & 23 km compared to baseline. ................................................................................................................. 87 5 11 Description of each parameter study ..................................................................... 87 5 12 Summary of results for parameter studies and model comparisons .................... 88 A1 Transmittance curve fit coefficients ....................................................................... 94 PAGE 10 10 LIST OF FIGURES Figure page 1 1 Flash burns on wooden poles 1.17 mi from ground zero at Nagasaki. ................ 15 1 2 Thermal effects and the street canyon. ................................................................. 16 2 1 Variation of apparent fireball surface temperature with time in a 20 kT air burst. ....................................................................................................................... 29 2 2 Normalized fireball power and fraction of thermal ener gy emitted versus normalized time in the thermal pulse of an air burst below 100,000 ft. ................ 29 2 3 Notation for view factor between two finite surfaces. ............................................ 30 2 4 Radiative exchange in a gray diffuse enclosure. .................................................. 30 3 1 The cont inuous passage simplification. ................................................................. 43 3 2 Plane source street canyon schematic. ................................................................. 44 3 3 Fireball surface flatness example. ......................................................................... 44 3 4 Reproduced observed surface temperature history for a 20 kT air burst. ............ 45 3 5 Direct and black body temperature method comparison. ..................................... 46 3 6 Schematic of a fireball destroying surrounding buildings. .................................... 46 3 7 Atmosp heric attenuation for urban aerosols and zero altitude at various sea level visibilities. ....................................................................................................... 47 3 8 Flux profile from the time dependent OFM for a 20 kT burst at 1mi. .................... 48 3 9 Passage diagram for SPM. .................................................................................... 48 3 10 SPM energy flow normalized to passage width. ................................................... 49 3 11 Flux profile from the SPM for a 20 kT burst at 1 mi. ............................................. 50 3 12 OFM and SPM flux profile comparison for 20 kT burst at 1 mi. ............................ 50 4 1 Coordinate system for CANYON. .......................................................................... 66 4 2 View fa ctor notation for parallel plates of equal size. ............................................ 67 4 3 View factor notation for perpendicular plates with a common edge of equal length. ..................................................................................................................... 67 PAGE 11 11 4 4 Differential surface meshing schematic. ................................................................ 68 4 5 View factor notation for rectangles in parallel planes. .......................................... 69 4 6 View factor notation for rectangles in perpendicular planes. ................................ 69 4 7 Equivalent whole surface and differential view factors example. ......................... 70 4 8 Surface meshing example for mxsfc = 4. .............................................................. 71 4 9 A 5 m cube of black surfaces. ................................................................................ 71 5 1 Baseline geometry for CANYON ........................................................................... 80 5 2 Average flux on the exit plane for baseline case mesh convergence in CANYON. ............................................................................................................... 81 5 3 Average flux on the exit plane for a mesh parameter of 4 and t of 1x105 s compared to mesh parameters of 4, 6, & 7 and t of 1x106 s. ............................ 82 5 4 Average flux on the exit plane for a mesh parameter of 7 and t of 1x106 s compared to mesh parameters of 8 & 10 and t of 5x107 s ................................ 82 5 5 Cany on entrance dimension changes. .................................................................. 84 5 6 Comparison of average exit plane flux vs. time for black body and direct source temperatures. ............................................................................................. 85 5 7 Comparison of average flux on exit plane vs. time for SPM, OFM, and enclosed & baseline CANYON models. ................................................................ 85 5 8 Average exit plane f lux vs. time at 1610 m for 20 kT, 10 kT, & 5 kT .................... 86 PAGE 12 12 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science EVALUATION OF 3 D RADIANT HEAT TRANS FER IN STREET CANYON S By Tucker Flagg Stachitas December 2009 Chair: Glenn Sjoden Major: Nuclear Engineering Sciences In the aftermath of a domestic urban nuclear detonation, investigators will seek to determine all technical details related to the explosion. In doing so, it is of interest to evaluate the propagation of thermal energy down an urban street canyon. Hence, this investigation ex amined the threedimensional transmission of thermal energy originating from an assumed nuclear explosion in an urban setting, where results were compared to existing public domain open field nuclear test data. Under the assumption of gray diffuse surfaces and black body source emission, the effects of various street canyon parameters were examined through the developm ent and validation of a Fortran 90 program called CANYON. It was shown that the total energy delivered was insensitive to material propertie s, but highly dependent on canyon geometry and source strength. The most important unknown was how the nuclear driven fireball interacts with immediately surrounding buildings as it releases energy into the street canyon. Assuming a well characterized source and average street canyon parameters the street canyon amplified the delivered energy by roughly 25% compare d to the same source in an open field PAGE 13 13 CHAPTER 1 INTRODUCTION In the event of a domestic urban nuclear detonation, one of the many issues investigators will seek to resolve is the origin of the device. A reliable data point in that post event assessment is the estimated yield of the device. As part of a National Technical Nuclear Forensics (NTNF) research effort funded through the Air Force Institute of Technology, a new estimation method is under development to correlate explosive yield with thermal damage to painted surfaces as a function of distance from the hypocenter in an urban street canyon environment. An evaluation of how the urban environment will affect the transmission of thermal radiation is c entral to obtaining an accurate es t imate of the yield. By simplifying the highly dynamic and multivariate phenomenon of thermal radiation from a nuclear detonation this investigation quantified the street canyon effect and compared it to existing open field data from testing carried out prior to the Limited Test Ban Treaty (LTBT) of 1963 ( Medalia, 2008) The correlation of explosive yield to thermal damage requires an accurate understanding of three general areas: the nuclear driven fireball source, the transmission down the street canyon and the receiver. This investigation focused on the transmission of thermal energy based on an approximation of the fireball source; r eceiver effects are the subject of a different investigation that focuses on thermal damage to p ainted surfaces Publically available nuclear fireball data was limited and approximate and was derived primarily from The Effects of Nuclear Weapons (Glasstone and Dolan, 1977) N uclear weapons were first developed and tested when sensi ng equipment was crude compared to todays standards which lends uncertainty PAGE 14 14 to the available data Unclassified weapons effects reports including this one, generally revert to the public source cited above for most of their background information. Existing thermal ef fects data are largely based on testing done in flat open areas. Trends based on yield, height of burst, and range have been characterized for rural scenarios through the combination of experimental results and theoretical calculations (Glasstone and Dolan, 1977). Empirical data for the urban environment are more limited because of modern treaty restrictions (e.g. 1963 LTBT) and the impracticality of testing on full scale cities, yet it is clear that the urban environment has the potential to significantl y change the propagation of weapon effects. For example, Figure 31 shows wooden utility poles from Nagasaki with uncharred portions that were protected from the thermal radiation by a fence. The fence was subsequently knocked down by the blast wave ( Gla sstone and Dolan, 1977). The poles were 1.17 miles from ground zero and the charred areas received 5 to 6 cal/cm2. Like the walls of a canyon funneling water, one can imagine long blocks of highrise buildings channeling thermal energy down a city street hence the term street canyon The burning car in Figure 1 2A illustrates the potential for thermal flash damage ( Glasstone and Dolan, 1977). Were this car found in an environment like that of Figure 1 2B, the effects of the thermal flash could have been very different in both range and magnitude Therefore, the major question this thesis addresses is how the transmission of thermal radiation will change in the presence of a street canyon. Chapter 2 covers thermal radiation from a nuclear fireball and the basics of radiant heat transfer. Chapter 3 describes how those concepts were merged into a general model of the source and PAGE 15 15 the street canyon. Chapter 4 describes the development of a computer program named CANYON that calculates the transmis sion of thermal energy based on that model. Chapter 5 shows the results from a variety of tests in CANYON and discusses how they compare to each other and to other models. Chapter 6 discusses conclusions about the governing assumptions and the overall st reet canyon effect. Chapter 7 discuss es future work. Figure 11. Flash burns on wooden poles 1.17 mi from ground zero at Nagasaki. PAGE 16 16 Figure 12. Thermal effects and the street canyon. A) Thermal radiation ignites a car. B) Sunset a ligned with a New York City street canyon showing conceptual similarity to a fireball source. (Source: http://en.wikipedia.org/wiki/File:Manhattanhenge2.jpg. Last accessed September, 2009) A. B. PAGE 17 17 CHAPTER 2 THERMAL RADIATION FUNDAMENTALS The enormous energy released by a nuclear explosion manifests itself in many forms which include the thermal flash, shock wave, and nuclear radiation. The fraction of energy released in each modality depends on a variety of factors such as the height of burst (HOB), the explosive energy yield ( W ), and the weapon design. The explosive energy yield, or simply yield, of a nuclear weapon is typically expressed in terms of the equivalent mass TNT that would release the same amount of explosive energy (i.e. a 1 kT nuclear detonation equals the explosive force of 1000 tons of TNT) (Glasstone and Dolan, 1977) Thermal Partition The fraction of the total yield released as thermal energy is known as the th ermal partition, or f For low altitude air bursts (below 15000 ft ), f was experimentally shown to be 0.35 for all yields ranging from 1 kT to 10 MT (Glasstone and Dolan, 1977) The thermal partition can be as high as 0.45 for air bursts above 15000 ft, depending on yield. According to Glasstone and Dolan, s urface bursts, which are defined as having a fireball that makes direct contact with the ground, generally deliver less thermal exposure than air bursts of the same overall yield. Dust and water vapor kicked up by the explosion can also attenuate the radiant energy of surface bursts resulting in an effective thermal partition of roughly 0.18. Chapter 3 contains a more complete discussion of ground effect s and other source term considerations. Thermal Pulse Timing Based on a complex combination of hydrodynamic forces, x ray interactions with air, and the principles of heat transfer, a nuclear fireball from a low altitude air burst PAGE 18 18 shows an apparent temperature history with two distinct peaks before a final cooling phase For a 20 kT air burst, the first and second peak s occur at roughly 0.44 ms and 150 ms respectively, as shown in the observed curve of Figure 21 (Glasstone and Dolan, 1977) The temperature m inimum after the first peak occurs at about 15 ms, which is very short compared to total length of the pulse. Though the temperatures in the first pulse are hi gh, it represents roughly 1% of the total thermal energy emitted and is typically ignored in calculations considering the e ntire pulse (Marrs et al., 2007) Further references to the thermal maximum are thus only referring to the second peak even though the first peak is technically at a higher temperature. The temperatures of Figure 21 are i ndependent of the yield, but the pulse timing can be scale d to the yield (Glasstone and Dolan, 1977 ). When power is considered instead of temperature, the thermal output can be scaled to yield in both magnitude and time For the second peak, maximum ther mal power ( maxP ) and the time at which it occurs ( maxt ) are governed by Equations 21 and 2 2 for bursts below 15,000 ft (Glasstone and Dolan, 1977). Equation 23 is a fit for the normalized power curve in Figure 22 (Bridgman, 2001). The curve is plotted together with the percent of thermal energy emitted versus normalized time, which shows that by 10 maxt 80% of thermal energy is emitted and the curve has mostly leveled off Both curves are accurate to within about 25% ( Glasstone and Dolan, 1977). 56 0 12 max10 18 3 W P (2 1) 44 0 max0417 0 W t (2 2) 4 2 max1 2 ) ( P t P (2 3) PAGE 19 19 Where; W = yield [ kT ] maxP = maximum thermal power [ cal/s ] maxt = time of thermal maximum [s] t = time after detonation [s] maxt t / = normalized time. Fireball Size The intensity of the thermal pulse is closely related to the size of the fireball which changes rapidly with time Though the true fireball shape can be asymmetric based on HOB, terrain, and weapon design, it is generally assumed to be spherical. Even so, i t is difficult to write a single expression for the radius versus time b ecause fireball expansion depends on both yield and HOB (Marrs et al., 2007 ) Instead, empirical fits have been derived for the radius based on the yield for three specific stages of fireball development : thermal minimum, breakaway, and thermal maximum Breakaway is when the hot opaque shock front becomes transparent enough as it cools to reveal the still hot fireball behind it. A t thermal minimum the radius is similar for air bursts and surface bursts (Equation 23). At breakaway, hydrodynamic forces dominate fireball development so that the radius is smaller for air bursts (Equation 24) than for surface bursts (Equation 2 5) ( Glasstone and Dolan, 1977 ). Calculations with the RADFLO computer code indicate a radi us at the thermal maximum expressed by Equation 2 6 (Marrs et al., 2007 ). In light of the difficulty of accurately describing the fireball size versus time, this investigation considered the fireball to be a steady sphere of max thermalR as was d one in a similar report from the Lawrence Livermore National Laboratory (Marrs et al., 2007). This assumption skews the effective fireball temperature. For a monotonically increasing fireball size, temperature will be artificially low before the thermal maximum and artificially high afterward This effect is discussed in more detail PAGE 20 20 in Chapter 3. Unless otherwise noted, max thermalR is referred to as only R for simplicity of notation. 4 0 min thermal27 R W (2 3) 4 0 breakaway air34 R W (2 4) 4 0 breakaway surface44 R W (2 5) 4 0 max thermal50 R W (2 6) Where; W = yield [ kT ] R = radius [ m ] Radiant Heat Transfer Two parameters allow the thermal flash to be characterized by radiant heat transfer alone.1 First the fireball develops temperatures of thousands of degrees as shown in Figure 21 Conduction and convection are driven by linearly proportional temperature d ifferences whereas radiation depends on 4th power temperature differences. Therefore, radiation tends to dominate in high temperature applications like emission from nuclear fireball s (Modest 2003 ). Secon d the short nature of the thermal flash precludes conduction or convection from having a significant effect on the deposition of thermal energy at long distances. Under the proposed methodology, paint samples would likely be collected at ranges of a mile or more from the detonation point Photons from radiation travel at the speed of light whereas the molecular collisions of conduction and convection are much slower Radiation should therefore dominate conduction and convection d uring the roughly 3 s pulse of a 20 kT air burst (Figure 22). 1 Because this thesis does not address ionizing radiation (i.e. X rays and Gamma rays), the term radiation is used only in reference to the non ionizing type associated with thermal heat transfer. PAGE 21 21 Considering the 4th power temperature dependence of radiation, the long ranges involved, and the short pulse duration, fireball thermal emission is reasonably characterized b y radiant heat transfer alone. Surface Property Definitions In radiant heat transfer, s urfaces are generally characterized by four fundament al properties that are defined as follows (Modest 2003 ): Reflectance, Absorptance, Transmittance Emittance e to black body emission at the same temperature The first three definitions cover all possibilities for incoming radiation, whi ch are summarized in Equation 27 where 0 Emittance ( ) is also a fraction between 0 and 1. 1 (2 7) For black surfaces (perfect absorbers and emitters) it follows that = = 1 and = = 0. Opaque surfaces are defined as surfaces that do not transmit (i.e. = 0 ), which reduces Equation 27 to Equation 28. Classifying a surface as opaque does not imply anything further about its emittance. 1 (2 8) For a given material, the aforementioned properties can vary with surface temperature, wavelength of incoming radiation, and direction of incoming or outgoing radiation. Chapter 3 addresses these dependencies in more detail PAGE 22 22 Black Body Radiator Approximation The t hermal emission spectrum from a nuclear fireball closely approximates the Planckian distribution of a black body radiator at the same temperature, meaning that it is also a near perfect absorber and emitter (Glasstone and Dolan, 1977). The emissive power of a black body is governed by the St efan Boltzmann law (Equation 29 ). Therefore, knowing the fireball temperature was a good approximation for knowing its emissive power. 4 bE T (2 9 ) Where; bE = black body emissive power [W/m2] = 5.6704x108 [W/(m2K4)] (Stefan Boltzmann constant) T = temperature [K] A black body emits energy diffusely, or equally in all directions. In the absence of other interactions or reflecting surfaces diffusely emitted energy decreases in intensity away from the source by spherical spreading which can be accounted for by the diffuse view factor. View Factor Definition Radiant exchange between two surfa ces is described in part by a geometric relationship known as the diffuse view factor. The diffuse view factor Fij is defined as the fraction of diffuse radiation leaving surface i that strikes surface j directly ( engel, 2007). The diffuse view factor is referred to as such to differentiate it from the specular view factor which considers the diffuse radiation received by surface j from surface i directly or after any number of specular reflections (Modest 2003 ). Because the ass umption of diffuse surfaces is common in many engineering applications, the diffuse PAGE 23 23 view factor is generally and henceforth referred to as simply the view factor with references to the specular view factor clarified as necessary G eneral formulae for view factor s are commonly derived on a differential basis One result is an integral expression for the view factor between t wo finite surfaces (Equation 2 10 ) based on the notation of Figure 23 (Modest 2003 ). Similar expressions are derived for arrangements involving differential surfaces or infinite surfaces. Di rect integration of Equation 210 can be exceedingly difficult for all but the simplest geometries. Considerable effort has therefore been put into tabula ting expressions for the view factor in a wide variety of geometries. A thorough catalogue of view factors is available online (Howell 2001) i jA A i j j i i ijdA dA S A F2cos cos 1 (2 10 ) The definition of the view factor leads to two important relationships. The first is known as the summation relation (Equation 211) which states that all the view factors leaving a particular surface must sum to 1 (i.e. 100% of the radiation leaving the surface). The second is the law of reciprocity (Equation 2 12), which states that the view factors ijF and jiF are proportionally related by the areas iA and jA Both these relations are useful when determining view factors in an enclosure. An enclosure of N surfaces will have N2 v iew factors. N j ijF11 (2 11) ji j ij iF A F A (2 12) PAGE 24 24 Radiant Exchange between Black Bodies The net rate of radiant heat transfer between two finite black bodies is governed by Equation 2 13, which was derived from the view factor (Equation 2 10), the emissive power (Equation 2 9), and a simple application of reciprocity (Equation 212). The complete derivation can be found in a typical heat transfer textbook ( engel, 2007). ) (4 2 4 1 12 1 12T T F A Q (2 13) Where; 12 Q = Net radiant transfer from surface 1 to 2 [W] 1A = Area of surface 1 [m] 12F = View factor from surface 1 to 2 = Stefan Boltzmann Constant, 5.6704x108 [W/(m2K4)] nT = Surface temperature [K] Gray Diffuse Emitters Reflectors and Absorbers While the black body approximation is appropriate for the fireball, it would not have been appropriate for the other surfaces of the street canyon. R eal sur faces tend to have large enough reflectances to invalidate simplifying them as black But, another class of idealized surfaces is reasonably applicable to the street canyon walls. Gray diffuse sur faces are hypothetical surfaces that emit equally at all w avelengths and emit, reflect and absorb diffusely (Modest 2003 ). D iffuse emission is a good assumption for most real surfaces while d iffuse reflection is a good assumption only for rough surfaces like sandpaper Very smooth surfaces like a mirror tend to reflect specularly and would not qualify as diffuse although the weathering of real surfaces may tend to make them more diffuse. Gray surfaces are those for which emittance is equal over the entire spectrum (Modest 2003) Local thermodynamic equilibr ium dictates that gray surfaces also emit and absorb equally ( ). Though uniform spectral emission is never truly PAGE 25 25 the case, many surfaces emit equally over a particular wavelength band of interest which allows the assumption of gray surfaces to be valid. By assuming the canyon walls to be gray diffuse surfaces, the radiant analysis was greatly simplified compared to what would have been required for real surfaces This assumption ignore d any specular reflections off smooth surfaces like windows Since windows and concrete building construction materials are spaced intermittently in a typical building, specular reflections will likely be followed by diffuse reflections, making the gray diffuse assumption reasonable in most cases. In an enclosure of gray diffuse surfaces, Equation 214 relates temperature T and heat q flux in integral form bas ed on the notation in Figure 24 T he emissive power of a gray body is related to bE by (Equation 2 1 5 ) A dA dA b b o dA dA AdF E E H dF q q' ') ( ) ( ) ( ) ( 1 ) ( 1 ) ( ) ( r' r r r' r' r r (2 1 4 ) 4T E Eb (2 15 ) As is common practice in many engineering applications, Equation 214 is recast in discrete summation form for computer calculations By introducing the Kronecker delta ( ij ) and breaking the hypothetical enclosure in Figure 23 into N subsurfaces Equ ation 21 4 is recast as Equation 21 6 (Modest 2003) N j N j oi bj ij ij j ij j j ijH E F q F1 11 1 (2 1 6 ) Where; j i = surface indices j i j iij, 0 1 PAGE 26 26 Equation 216 is expressed in matrix fo rm as Equation 2 1 7 with the elem ents of the coefficient matrices defined in Equations 2 1 8 and 21 9 For a gray diffuse enclosure at steady state, Equation 2 16 can be solved for heat flux given temperature or for temperature given heat flux using a variety of matrix solution methods that are discussed in Chapter 4 o bh e A q C (2 1 7 ) ij j j ij ijF C 1 1 (2 18 ) ij ij ijF A (2 1 9 ) Semi Infinite Solids For a transient analysis of radiant heat transfer, it was necessary to account for change s in street canyon wall temperatures with time based on the flux incident from the fireball T emperature in a solid is governed by th e heat equat ion (Equation 2 20 ), which is typically solved in 3D using numerical methods like finite differencing. Given the computational effort required to resolve the 3 D steady state flux profiles alone, a thorough numerical treatment of transient conduction int o every canyon surface would have rendered simulation times unfeasibly long. Therefore, a simpler treatment of conduction was highly desirable. t T k q T 1 '2 (2 20 ) Where in Cartesian coordinates ; 2 2 2 2 2 2 2, z y x q = volumetric heat generation [ W/m3] k = thermal conductivity [W/(mK)] = thermal diffusivity [m2/s] t = time [s]. PAGE 27 27 Equation 220 was greatly simplified by assuming the canyon walls to be semi infinite solids with no internal heat generation. Semi infinite solids are hypothetical bodies that are infinite in all directions with a single plane surface They are useful when the region of interest is close to one surface of a thick solid an d the observation time is short enough that the other boundaries of the solid are inconsequential ( engel, 2007) These conditions held for a roughly 3 s thermal pulse (Figure 21) on city buildings that are many meters thick especially because only the surface temperatures were of interest. Under the semi infinite condition expressions are derived for the transient temperature distribution into a body for different boundary conditions. For the case of a known constant surface heat flux ( sq ) the temperature profi le is described by Equation 2 21 which calls for the complementary error function (Equation 222). Equation 222 cannot be analytically integrated. Fortunately, only needing the surface temperatures ( 0 x ) reduced Equation 2 21 to a simple analytic expression (Equation 2 23). t x x t x t k q T t x Ts i 2 erfc 4 exp 4 ) (2 (2 21 ) 022 1 erfc du eu (2 22 ) t k q T t Ts i4 ) 0 ( (2 23) Where; x = distance into the solid [m] iT = initial temperature [K] sq = constant surface flux [W/m2] PAGE 28 28 Unit Conventions: Total Power vs. Emissive Power Equations 21 and 29 show an apparent unit inconsistency between maximum thermal power ( maxP ) and black body emissive power ( bE ). In physics, power is defined as a unit of energy per unit time, which agrees with the definition of maxP in cal/s. In radiant heat transfer, emissive power is defined as a unit of energy per unit time per unit area, which agrees with the definition of bE in W/m2 (where W = J/s). A more descriptive term for bE would be the black body emissive areal power T he conventional terminology for bE as simply emissive power is maintained throughout this document. Confusion between further references to power ( P ) and emissive power ( E ) is avoided by the inclusion of units and by referring to power ( P ) as total power when necessary. The thermal characteristics of nuclear fireballs and the principles of radiant heat transfer presented in this chapter were combined and applied to develop a model for transmission through the street canyon. The development of that model is addressed in the next chapter. PAGE 29 29 Figure 21. Variation of apparent fireball surface temperature with time in a 20 kT air burst. Figure 22. Normalized fireball power and fraction of thermal energy emitted versus normalized time in the thermal pulse of an air burst below 100,000 ft. PAGE 30 30 Figure 23. Notation for view factor between two finite surfaces.(Source: Modest, M.F., 2003. R adiative Heat Transfer. Academic Press, San Diego, CA., p. 135). Figure 24. Radiative exchange in a gray diffuse enclosure. (Source: Modest, M.F., 2003. Radiative Heat Transfer. Academic Press, San Diego, CA., p. 169). PAGE 31 31 CHAPTER 3 MODEL CONCE PTUALIZATION The urban environment is a complex mixture of geometries and materials that must be considered for nuclear weapons effects modeling. The potential variety from city to city and even block to block necessitated a range of simplifying assumptio ns whose significance could be investigated once a workable model was established. The fireball source is also highly dynamic and difficult to describe in full detail. The following sections address how these aspects of the model were treated with additi onal discussion of atmospheric effects and initial modeling results. Geometry The first simplification was considering a street canyon of multiple blocks as a continuous passage with planar walls (Figure 3 1 ). The rectangular street grids present in many cities in the United States support the planar wall assumption While cross streets represent potential leakage paths for radiation being channeled by the street canyon, coding and analysis w ere greatly simplified by assuming a continuous passage. Cross streets are also generally narrow compared to the width of the blocks between them which will help minimize the leakage in a real scenario. For example, the buildings in Figure 31 A are roughly 5 times wider than the cross streets along the direction of the main avenue. The continuous passage assumption would, of course, be less valid for wider cross streets. Standard street canyon dimensions were selected to facilitate model comparisons. The standard height considered in this work was 100 m, which is about the height of a 30 story building. Many city blocks are dominated by shorter buildings, but this parameter can easily be changed in any given CANYON model The standard width PAGE 32 32 considered was 40 m, which is about the width of a large city street inc luding the sidewalks. The baseline length was 1610 m, or 1 mi. The peak overpressure from a 20 kT detonation will be about 5 psi at 0.75 mi, which is enough to crumple small cars (Glasstone and Dolan, 1977).1 Furthermore, c loser than a mile, the thermal radiation becomes so intense that any paint would be charred beyond measurement. A crumpled and c ompletely c harred car would likely not provide useful forensic samples for the method under investigation, but it would provide an indicator of t he minimum energy delivered at that distance. Material Properties The next simplification was assuming uniform material properties for each surface along the length of the street canyon. The walls and floor of a real street canyon are a diverse mixture of different stone s glass, and asphalt whose relative proportions can vary widely from one building faade to the next. The heat transfer properties of each material dictate how they would interact with the thermal ra diation from a nuclear fireball G la ss can be highly specularly reflective but it is also transmissive, whereas concrete is a diffus e reflector and is completely opaque. A sophisticated representation of spatially varying material properties was not appropriate for the early stages of modeling given the unknown importance of small property changes on the total energy transmission. A rough estimate of average material properties was sufficient to characterize the street canyon with the understanding that results would dictate the necessity of more detailed property definitions. The walls and floor were assumed to be 1 Larger modern cars with plastic composite exterior paneling will likely crumple more easily than the cars originally subject to nuclear testing. This will tend to extend the range at which structurally intact samples are available. PAGE 33 33 entirely concrete. This assumption ignored potentially large property differences, but was an acceptable uncertainty for initial modeling; Chapter 5 shows that 30% property changes have minimal effect on the result Concrete itself is available in many different forms whose thermal properties vary widely based on constituent materials and construction methods. Fortunately, a large database of build ing material properties revealed some typical values that are summarized in Table 31 and were used as the baseline material properties (Clarke et al., 1991) The database contained 3 entries for the emittance ( ) of concrete that ranged from 0.56 to 0.95. Among 98 different entries for various cements, clays, & concretes, 82 types were listed with a specific heat ( Cp) of 840 or 837 J/ (k gK ) with sporadic values up to 2000 J/ (k gK ). The density ( ) ranged from 400 to 2500 k g/m3. The thermal conductivity ( k ) ranged from 0.1 to 2.3 W/ (mK ). Material properties were also assumed to be constant for all surface temperatures. While real surface properties would change as they are heated by the fireball, it was more important to establish a reasonable starting point from which a sensitivity analysis could show if more thorough material specifications were necessary Source Term Considerations By combining the background information in Chapter 2 with the assumed geometry, the street canyon was initially framed as a plane walled enclosure of gray diffuse surfaces. A rising and expanding fireball does not readily fit this framework, so i t was necessary to develop a reliable method for applying the fireball source term The fireball was modeled as one of the enclosures surfaces so that its temperature defined its strength. The details and implications of this approach are discussed below. PAGE 34 34 Stationary Plane Source A pproximation The spherical fireball surface was modeled as a stationary plane source filling the entr ance of the canyon, represented by the red s urface on the left of Figure 32 An analysis of Equations 23, 2 5, and 26 confirmed that a fireball will almost always fill a canyon entrance of the baseline dimensions for the range of yield considered. Table 32 shows the fireball diameter at each phase of development for 2 kT and 20 kT surface bursts All the noted diameters exceed the 100 m canyon entrance height except for the 2 k T diameter at thermal minimum. The plane source assumption is justified because over the entrance area of the canyon the fireball surface will appear reasonably flat. For example, a 20 kT burst will have a radius of 135.7 m (Equation 26). As shown in Figure 33 and Equation 3 1 i f a circle of that radius is drawn so that it just intersects the sides of the canyon entrance and a chord of 40 m is drawn between the intersections, the distance from the center of the circle to the chord ( ) will be 134.2 m, which is nearly the entire radius Therefore, the surface of the spher e is nearly planar from the perspective of the street canyon. m x R 2 134 2 40 7 .1352 2 2 2 (3 1) The stationary source assumption was supported by a comparison of the fireball radius to the length of the street canyon. A radius growth from 0 to 135.7 m represented only an 8.4% shortening of a 1610 m street canyon. The assumption appeared reasonable but its effect on the result was analyzed in Chapter 5 This change becomes less important for longer canyons. PAGE 35 35 Direct Temperature Method As shown in Figu re 21, a specific temperature history was available for a 20 kT air burst One method for determining the time varying source strength was directly using this temperature history along with Equation 29 to calculate the emissive power of the fireball. T he data was extracted from Figure 21 with a Java shareware program that works by setting reference points on a scanned image of a plot and taking user inputs for values of the origin and axes limits (Tummers, 2006) It then uses an algorithm to trace the curve of interest and output the data to a text file which can easily be imported into Excel Figure 34 shows that program faithfully reproduced the observed temperatures from Figure 21. Black Body Temperature Method Another prospective method for impl ementing the source term was to infer a surface temperature from the normalized power curve (Figure 2 2) and the black body radiator assumption. The assumed fireball radius was very important to this calculation I t required converting the total power ( P ) to emissive power ( bE ) through dividing by the surface area of the fireball so a n oversized or undersized radius had the potential to skew the effective temperature dramatically. The conversion is shown in Equat ion 3 2 where Equations 21 through 2 3 were used to convert normalized power to total power The resulting black body temperature ( bbT ) was compared with the direct temperature (dirT) at peak power (Table 33 ) and over the length of the pulse (Figure 3 7) The uncertainty on bbT comes from Glasstone and Dolans 25% error estimate for normalized power (Figure 22 ). Table 3 3 shows that dir T was within the uncertainty for bbT at peak power, but Figure 34 shows that this was not the case elsewhere. If the black body PAGE 36 36 assumption is valid for fireball emission, dir T should be equal to bbT at all times. The discrepancy can be ac counted for mostly by the steady state radius assumption, but it is difficult to thoroughly explain the discrepancy without knowledge of the original data collection methods and uncertainties involved. These specifics were not addressed in the source and are not publically available (Glasstone and Dolan, 1977). Because other related publications have used the normalized power curve, this investigation used bbT as the baseline source t emperature (Bridgman, 2001; Marrs et al., 2007) Chapter 5 addresses the influence of th is choice on the final result. 4 24 ) ( ) ( R t P t Tbb (3 2) Shadowing and Ground Effects In a flat open field, an entire hemisphere of the fireball is visible from the receiver and there are no terrain or building features to distort fireball development. In an urban environment, large buildings surrounding the detonation could ether absorb o r reflect blast energy and would block some of the view of the fireball from the other end of the street canyon. If e nergy is spent destroying the immediately surrounding buildings it would not be available f or release as thermal radiation (Figure 33) On the other hand, if the buildings remain intact the fireball itself will be channeled by the street canyon thereby moving the thermal source closer to potential receivers and possibly intensifying the radiation by reducing the effective surface area. In either case, the buildings will be channeling the fireball across an area defined by the street canyon. The hydrodynamics and materials analysis required to accurately model the interaction of the fireball with the immediately surrounding buildings w as beyond the PAGE 37 37 scope of this investigation. Instead, the focus was on the interaction of the thermal radiation with the street canyon walls assuming a well characterized thermal source. Nonetheless, a primitive example was demonstrative of the importance of the source term effects. A circle with the baseline 20 kT fireball radius of 135.7 m would have an area of 57870 m2. The assumed 40 m by 100 m canyon entrance has an area of 4000 m2, or roughly 6.9% of the area of the circle. If the fireball is perfe ctly shadowed b y the entrance, only this small fraction of it is visible down the street canyon. This effect is made clear by some of the results in Chapter 5. Even in an open field, a surface burst will have different thermal characteristics than an air burst because of reflection off the ground. To account for this effect, surface bursts are commonly modeled as air bursts with twice the yield. This doubling is reasoned from the 4 steradians of a sphere being reflected into the 2 hemis phere. On the other hand, the ground debris and buildings surrounding an urban surface burst will attenuate the fireball as noted. Ground reflection and debris attenuation compensate for each other so that it was reasonable to assume no multiplication on the yield ( Marrs et al., 2007) Atmos pheric Attenuation If the interior of the street canyon was a vacuum, then spherical spreading via view factors would be sufficient to describe the degradation of radiant intensity away from the source. Radiant heat t ransfer problems of this sort (where the space between surfaces is considered a vacuum) are known as nonparticipating media problems. I n reality, thermal radiation will interact with air, dust, and other atmospheric particles as it travels. If they interact, the photons of thermal radiation will either be scattered into new directions or absorbed by the atoms they hit. Both phenomena tend to lessen the PAGE 38 38 amount of energy ultimately transmitted through the atmosphere.2 The probabilities of these interactions are known as cross sections. They can vary by orders of magnitude depending on the energy of the incoming photon, the temperature of the target, and the molecular, atomic, and nuclear properties of the target. Problems that consider these p henomena are known as participating media problems, which are considerably more difficult and computationally intensive to solve, particularly in 3D Doing so directly was beyond the scope of this investigation especially considering the variation of at mospheric properties with season, weather, and altitude Fortunately, the pertinent cross sections for atmospheric transmittance have been compiled into a computer program called LOWTRAN that was developed by US Air Force Research Laboratory, Geophysical Research Directorate (Bridgman, 2001). LOWTRAN allows the user to specify atmospheric properties and visibilities. It also includes a s eries of standard atmospheres that can be selected. A 1985 thesis from the Air Force Institute of Technology (AFIT) u sed LOWTRAN to calculate transmittances specifically for nuclear thermal sources (Leonard, 1985). Its results were directly applicable to this investigation. Leonards thesis calculated transmittance vs. distance from ground zero at a variety of altitud es, aerosol conditions, and sealevel visibilities. It then found 5th order polynomial curve fits for the results of the form in Equation 3 3 and reported the appropriate constants 5 6 4 5 3 4 2 3 2 1) (ln ) (ln ) (ln ) (ln ) (ln x C x C x C x C x C C (3 3 ) Where; 2 Only a black body so urce spectrum was considered in this work. Potential spectral variations are left for future work. PAGE 39 39 = atmospheric transmittance x = distance from ground zero [km] Four particular cases from the thesis were applicable to this investigation: the set calculated for urban aerosols and zero altitude at visibilities of 5 km, 10 km, 15 km, and 23 km, which roughly correspond to haze, light haze, clear, and very clear conditions in the International Visibility Code (Table 34) (Glasstone and Dolan, 1977) The fits were only valid for 1 km x 25 km and are plotted together in Figure 37 The appropriate constants are available in Appendix A. For the purposes of initial model comparisons, transmittance was assumed to be 100% Transmittance was included in the final model as a multiplier on the view factor. Only the 4 distinct visibilities could be selected. The effect of variable transmittance is addressed in Chapter 5. Initial Modeling A literature review revealed two prospective methods to use as baselines for comparison to the final model. The first w as derived from a combination of Glasstone and Dolans The Effects of Nuclear Weapons (1977) and Chapter 6 excerpted from Bridgmans Introduction to the Physics of Nuclear Weapon Effects (2001) The model developed a radiant exposure relationship based on spherical spreading from a point source fireball It is herein referred to as the Open Field Model (OFM) because it does not directly consider urban environmental effects. A paper from the 1964 Symposium on Thermal Radiation of Solids by E. M. Sparrow t itled Radiant Emission, Absorption, and Transmission Characteris tics of Cavities and Passages addresses the transport of thermal radiation through passages. It is referred to as the Sparrow Passage Model (SPM). These two models served as the lower and upper bounds of expected results for the final street canyon model. PAGE 40 40 Open Field Model The OFM assumed the fireball to be an isotropically radiating point source. A simple equation provided a time independent estimate for the total thermal fluence received by a target based on yield, thermal partition, atmospheric transm ittance, and range (Equation 3 4 ). For the input values in Table 3 5 the predicted total thermal fluence was 21.49 cal/cm2. 2 124 10 D W f q (3 4 ) Where; q = thermal fluence [cal/cm2] = atmospheric transmittance D = slant range (target to receiver distance) [cm] A time depe ndent OFM was constructed in TK Solver 5.0 using Equations 21 through 2 3 The f lux was calculated by inserting the time dependent power into Equation 34 to get Equation 35 which was then integrated to find the fluence. This time dependent model might seem superfluous given Equation 3 4 but the time independent fluence was not suitable for two reasons. First it provided no indication of the total time for emission, which prevented comparison to later models based on integrated fluxes. Second integrating until 100% of the thermal energy was emitted would not have been realistic for th e street canyon scenario because at longer time scales other weapons effects come into play that were beyond the scope of this investigation. For example, before the thermal pulse is over the blast wave could kick up enough debris to effectively end the l ong range transmission of thermal energy from increased atmospheric attenuation The integration time limit was therefore set at the length of the direct temperature history (Figure 2 1), which ends at about 3.2 s or PAGE 41 41 roughly 20.5tmax. The same time was used for integrations based on the black body temperatures. 24 ) ( ) ( D t P t q (3 5 ) For the baseline case of 20 kT at 1 m i the time dependent OFM predicted a total thermal fluence of 17.29 cal/cm2 and a peak flux of 52.24 cal/( cm2s) As expect ed, the flux profile (Figure 33 ) reflects the shape of the normalized power curve (Figure 2 2). The model parameters and res ults are summarized in Table 36 The total fluence was calculated using 7point Gaussian Quadrature integration with a tolerance of 1x105. The time dependent OFM fluence was about 80.4% of the time independent OFM fluence, which agree s with Figure 2 2. Sparrow Passage Model The SPM calculates the energy flow in a passage between opposing isothermal black body radiators given the passage dimensions and the temperatures of the environments at either end. The model assumes adiabatic plane walls that are gray di ffuse ref lectors and emitters (Figure 39 ). The parameters used in the model are: = Half angle of tap er T1 = Temperature of environment 1 [K] T2 = Temperature of environment 2 [K] x = Distance along passage dA(x) = Differential wall area at x. h = Half width of passage L = Length of passage The net rate of energy flow through the passage was given in int egral form by Equation 36 With Equation 36 normalized to passage width, the effect of passage length at varying angles of taper is shown in F igure 3 10. Equation 3 6 was simplified into Equation 37 for long parallel walled passages (which correlated well to street PAGE 42 42 canyons), where L/h was large and = 0. Large was defined as the ratio where the curve in Figure 39 levels off for a particular For = 0, p assages were considered large at L/h ) ( ) ( 1 1 ) (1) ( 0 4 2 4 1 4 2 4 1 4 2 4 1 1x dA F T T T x T A T T A QA x dA L (3 6 ) ) / ( ] 1 ) / ln( 2 [ ) (4 2 4 1 1h L h L T T A Q (3 7 ) Using the 20 kT black body fireball temperatures, the SPM predicted a total fluence of 234. 76 cal/cm2 and a peak flux of 709.50 cal/(cm2s) for the parameters in Table 3 7 The total fluence was once again calculated using 7point Gaussian integr ation with a tolerance of 1x105. The SPM flux profile (Figure 311) also resembles the normalized power curve because it was based on the black body temperatures of the fireball as a function of time. The exit temperature was maintained at 300 K to represent an unchanging ambient environment. This condition is addressed in more detail in Chapter 4. Initial Model Comparison The SPM predicted a fluence 13.59 times higher than the time dependent OFM. The same ratio appeared in a comparison of the peak fluxes. T his result agreed with qualitative expectations. The walls of the SPM prevent the spherical spreading of source radiation. Because they are adiabatic gray diffuse surfaces, all the incident en ergy is returned to the passage through reflection or absorption and reemission. A direct comparison of the flux profiles emphasized the predicted e ffect of the canyon (Figure 312 ). PAGE 43 43 The results of the SPM were useful but the model has an important limit ation. The wall, exit, and entrance boundary conditions were well defined, but there was no direct discussion of the vertical dimension (Sparrow, 1965). It was not clear whether the passage was considered infinitely tall or if the ceiling and floor were supposed to be adiabatic gray diffuse surfaces like the walls. Figure 3 9 implies that the model was 2D but Equation 36 calls for the entrance area, implying a third dimension. It is possible that Equation 36 required an area per unit height and was i ntended to be a 2 D formula, but this was not specifically clarified. In light of this significant ambiguity, the results of the SPM were considered cautiously. It was then necessary to apply a more detailed method for calculating the 3 D radiant heat tr ansfer, which is the subject of the next chapter. A. B. Figure 31. The continuous passage simplification. A) An aerial view of the Manhattan street grid near the Westside Rail Yards (Source: http://www.nyc.gov/html/dcp/gif/hyards/aerial detail.jpg Last accessed September, 2009). B) A schematic overlay of the continuous planar passage. Table 3 1. Baseline thermal properties for concrete Property Symbol Value Units Emittance 0.70 1 Specific Heat C p 840.0 J/(kgK) Density 1800.0 kg/m 3 Thermal Conductivity k 1.50 W/(mK) PAGE 44 44 Figure 32. Plane source street canyon schematic. Table 3 2. Fireball diameters for 2 kT and 20 kT surface bursts Fireball phase Diameter at 2 kT [m] Diameter at 20 kT [m] Thermal minimum 71 179 Surface breakaway 116 292 Thermal maximum 126 271 Figure 33. Fireball surface flatness example. Source Receiver 40 m 100 m 1610 m Street Canyon Fireball PAGE 45 45 Figure 34. Reproduced observed surface temperature history for a 20 kT air burst.3 Table 3 3. Temperature method comparison at peak power Parameter Symbol Value Units Yield W 20 kT Radius R 135.7 m Max power maxP 7.122E13 W Time to max power maxt 0.1558 s Black body Temperature bbT 8583 + 492 596 K Direct Temperature dirT 8140 K The 25% error on maxP produces a lopsided uncertainty on bbT because of the 4T relationship. 3 The temperature is plotted in C to maintain similarity with Figure 2 1 but was converted to K for all calculations. 1000 10000 1.E 04 1.E 03 1.E 02 1.E 01 1.E+00 1.E+01Temperature [C]Time After Explosion [s] PAGE 46 46 Figure 35. Direct and black body temperature method comparison. Figure 36. Schematic of a fireball destroying surrounding buildings. 1000 10000 1.E 04 1.E 03 1.E 02 1.E 01 1.E+00 1.E+01Temperature [K]Time After Explosion [s] Direct Black Body Black Body 25% uncertainty PAGE 47 47 Table 3 4. International Visibility Code Atmospheric Condition Upper limit of visibility [km] Light to thick fog 1 Thin fog 2 Haze 4 Light haze 10 Clear 20 Very clear 50 Exceptionally clear 280 Figure 37. Atmospheric attenuation for urban aerosols and zero altitude at various sealevel visibilities. Table 3 5. Time independent OFM parameters and result Variable Value Units f 0.35 1 W 20 kT 1 1 D 161000 cm q 21.49 cal/cm 2 0 2.5 5 7.5 10 12.5 15 17.5 20 22.5 25 .1 0 .1 .2 .3 .4 .5 .6 .7 .8 Distance from Ground Zero [km] 5km 10km 15km 23km PAGE 48 48 Table 3 6. Time dependent OFM parameters and results Parameter Symbol Value Units Yield W 20 kT Range D 1610 m Transmittance 1 1 Total Fluence q 17.29 cal/cm 2 Peak Flux maxq 52.24 cal/(cm 2 s) Figure 38. Flux profile from the time dependent OFM for a 20 kT burst at 1mi. Figure 39 Passage diagram for SPM. (Source: Sparrow, E.M., 1965. Radiant emission, absorption, and transmission characteristics of cavities and passages. Symposium on Thermal Radiation of Solids San Francisco 1964. Scientific and Technical Information Division, NASA, Washington D.C., p. 111. ) 0 .5 1 1.5 2 2.5 3 3.5 0 5 10 15 20 25 30 35 40 45 50 55 Time [s]Thermal Flux [cal/(cm^2*s)] PAGE 49 49 Figure 310. SPM energy flow normalized to passage width. (Source: Sparrow, E.M., 1965. Radiant emission, absorption, and transmission characteristics of cavities and passages. Symposium on Thermal Radiation of Solids San Francisco 1964. Scientific and Technical Information Division, NASA, Washington D.C., p. 111. Table 3 7. Sparrow Passage Model parameters and results Parameter Symbol Value Units Half wi dth of passage h 20 m Passage length L 1610 m Height of passage H 100 m Total fluence q 234.76 cal/cm 2 Peak flux maxq 709.50 cal/(cm 2 s) PAGE 50 50 Figure 311. Flux profile from the SPM for a 20 kT burst at 1 mi. Figure 312. OFM and SPM flux profile comparison for 20 kT burst at 1 mi. 0 .5 1 1.5 2 2.5 3 3.5 0 100 200 300 400 500 600 700 800 Time [s]Thermal Flux [cal/(cm^2*s)] 0 .5 1 1.5 2 2.5 3 3.5 0 100 200 300 400 500 600 700 800 Time [s] Thermal Flux [cal/(cm^2*s)] SPM Flux OFM Flux PAGE 51 51 CHAPTER 4 CANYON CODE DEVELOPMENT As a general topic, the calculation of 3D radiant heat transfer is not unique or new. Sophisticated commercial packages are available that provide user friendly simulations of multimode heat transfer and computational fluid dynamics.1 T hese programs were built for use in typical engineering applications like product design and performance analysis As this chapter makes clear, analyzing the thermal flash from a nuclear fireball introduced numerical issues not common to such applications. The short time scales, long ranges, and large energy magnitudes involved could have easily pushed commercial programs beyond their computational limits. The abili ty to customize the inputs and outputs of commercial programs was also a concern because of the desire to make computational results easily available for incorporation into further research efforts related to this investigation. Therefore, the CANYON program was written from scratch because of the unique application and the likelihood that commercial codes would have been inadequate for the purposes of this work CANYON was written in Fortran 90 primarily using Compaq Visual Fortran Version 6.6 and Microso ft Developer Studio. It used the geometry, surface properties, and source emission characteristics discussed in Chapters 2 and 3 to calculate time dependent heat flux es fluence s and temperature s along the street canyon walls and through the exit plane. By careful validation during development, CANYON delivered results that accurately captured the physics included in the model. This chapter describes the logic and assumptions implemented in CANYON. The sourc e code is 1 Two examples are SolidWorks COSMOSWorks and FLUENT. PAGE 52 52 included in Appendix B For the baseline case developed in this and the next chapter, s ample input and output file s are included in Appendices C and D respectively Reference System The first step in code development was to select a coordinate system and establish a reference conventi on. As implied by Figures 3 1 and 32, the street canyon was considered a right rectangular prism Because every surface was modeled as planar, Cartesian coordinates were most appropriate. The origin was placed at the lower left corner of the source plane with axes obeying the right hand rule so that the entire canyon was in the first octant (Figure 4 1 ). Opposing surfaces were numbered in consecutive pairs with odd surfaces intersecting the origin. The X, Y, & Z dimensions were defined as dx, dy, & dz and given the respective values of 40 m, 200 m, and 80 m for debugging purposes. The surface numbers and dimension labels are summarized in Table 4 1. View Factors The next phase of code development was to begin calculating the appropriate view factors which is typically the most difficult aspect of radiant heat transfer computations because of the double area integral required by Equation 210. A variety of approaches were available for calculating view factors including direct integration, statistical determination, and special methods based on particular geometries (Modest, 2003). The geometry of the street canyon model was simplified enough to avoid more computationally intensive methods like numerical integration and statistical simulation. Every required view factor was calculated by reciprocity, summation, and previously determined analytic relationships (Howell, 2001; Ehlert and Smith, 1992). PAGE 53 53 Whole Surface View Factors Each of the 6 surfaces was first considered in its entirety, allowing for relatively simple view factor calculations. The 6 primary surfaces were later meshed into subsurfaces, so these relationships are referred to as whole surface view factors and the latter are referred to as differential surface view factors Only two basic relationships exist in Figure 41: parallel plates of equal size ( Figure 4 2 & Equation 4 1) and perpendicular plates with a common edge o f equal length (Figure 4 3 & Equation 42) After employing Equation 4 1 once and Equation 42 twice, summation, reciprocity and symmetry were then sufficient to calculate the remaining view factors Y Y X X X Y X Y Y X Y X Y X Y X XY F1 1 2 1 2 2 1 2 2 1 2 2 2 2 2 1tan tan 1 tan 1 1 tan 1 1 1 1 ln 2 (4 1) Where; c a X / c b Y / ) (2 2 2 2 1 2 2 1 2 2 2 2 2 1 2 2 1 2 2 2 1 2 1 2 1 ln 4 1 2 2 1 1 tan 2 2 1 1 tan 1 1 tan 1 2 1 H W H H W H H W H W W H W W H W H W W H W H H H W W W F (4 2) Where; l h H / l w W / The whole surface view factors 2 1 F 3 1 F & 4 3 F were calculated directly. From there it was possible to determine the 33 other whole surface view factors. Many specific paths were available to complete the set by combining symmetry, summation, PAGE 54 54 and reciprocity in different sequences. The particular solution method used is briefly described here, with the full details av ailable in Appendix B. The whole surface view factors were stored in a 6x6 matrix indexed by the row number as the originating surface and the column number as the receiving surface. After the three direct calculations, the self seeing view factors were set to zero by inspection because all the surfaces were flat. Starting with surface 1, symmetry and summation were then used to calculate all the view factors leaving a surface. Lastly, reciprocity was applied over the entire matrix so that if a view fac tor value had been determined above the matrix diagonal, its reciprocal could be calculated. The program then returned to the second row of the matrix to calculate the remaining view factors leaving surface 2 by symmetry and summation before cycling back through the reciprocity loop. This process was repeated until the entire matrix was filled and all the whole surface view factors were calculated. Table 42 is an expression of the whole surface view factor matrix that indicates how each specific entry w as calculated. Table 4 3 shows the whole surface view factor matrix for the 40 m x 80 m x 200 m debugging dimensions. It reflects the relationships noted in Table 42 and the sum of each row is 1.00000 as expected from the summation rule Differential View Factors Whole surface view factors would have been insufficient to describe an accurate temperature profile because CANYON could only calculate a set of temperatures as large as the set of view factors. By meshing the whole surfaces into differential surfaces, finer view factors were calculated to provide the basis for better temperature profiles. Whole surfaces were meshed in each direction as shown in Figure 44. As with the whole surface view factors, two basic relationships were required to calc ulate PAGE 55 55 all the differential view factors. The first relationship was for one rectangle to another in parallel planes (Figure 4 5 Equation 4 3) (Howell, 2001) The second was for one rectangle to another in perpendicular planes with all boundaries paralle l to one of the axes (Figure 4 6 Equation 44) (Ehlert and Smith, 1992) 2 1 2 1 2 1 2 1 ) ( 1 2 1 2 2 1) ( 1 ) )( ( 1l k j l k j i i l k j iy x G y y x x F 2 2 2 2 2 / 1 2 2 1 2 / 1 2 2 2 / 1 2 2 1 2 / 1 2 2ln2 tan tan 2 1 z y x z z y x z y x z x y z x y G (4 3) 2 1 2 1 2 1 2 1 1 2 1 2 2 1) ( 1 1l k j i l k j i l k j iy x G y y x x F 2 2 2 2 2 2 2 / 1 2 2 1 2 / 1 2 2ln 4 1 tan 2 1 y x y x x y x y G (4 4) Once again, there were a variety of solution paths available to solve for t he set of differential view factors. A complex procedure of cycling through symmetry, summation, and reciprocity could have been used to decrease the number of direct calculations, but that approach would have significantly increased the difficulty of the p rogramming. Instead, every differential view factor was calculated directly using Equation 43 or 44 except for coplanar differential view factors, which were zero by inspection. Although this approach was computationally intensive, it only had to be performed once at the beginning of a simulation. It also left CANYON open to future development of more PAGE 56 56 complicated geometries where an intricate cycle would have tied it very closely to sp ecific indexing relationships T he indexing of differential surfaces for CANYON was still fairly complex. To describe a particular differential surface pair required keeping track of the originating whole surface, a 2 part index on the originating whole surface, the receiving whole surface, and a 2part index on the receiving whole surface. The differential view factors were stored in a 6dimensional array called VFdiff with indices for element VFdiff (iFROM, n1, n2; iTO, m1, m2) defined as follows: iFROM: Originating whole surface n1: Leading alph abetical dimension index of differential surface on iFROM (i.e. n1= y index if iFROM is in y z plane) n2: Trailing alphabetical dimension index of differential surface on iFROM (i.e. n2 = z if iFROM is in y z plane) iTO: R eceiving whole surface m1: Leading alphabetical dimension index of differential surface on iTO m2: Trailing alphabetical dimension index of differential surface on iTO To validate the differential view factor calculations the relationships between surfac e 5 (the street) and surface 3 (the left wall) were closely examined for the debugging dimensions with two meshes in each direction. A view factor is based on relative dimensions and distances so it should be the same for a surface pair that has been scaled equally in all directions from an equivalent arrangement. For example, with an equal number of meshes in each direction, the whole surface view factor F5 3 should be the same as the differential view factor VFdiff(5,1,1;3,1,1) (Figure 47 ). As expect ed, CANYON reported F5 3 equal to VFdiff(5,1,1;3,1,1) as well as the other differential view PAGE 57 57 factor directly along the y axis (Table 44) Other expected relationships like symmetry also held true All the differential view factors for Figure 4 7 are reported in Table 44. The verification of these relationships indicated that CANYON calculated the differential view factors correctly. The values are reported in double precision with 10 significant figures because other aspects of CANYON required very precise values. These numerical issues are addressed in the next section. Numerical Issues Implementing Equations 43 and 44 in CANYON revealed some numerical issues that had to be overcome. Calculating ) (l k j iy x G for both cases required adding and subtracting nearly equal large numbers As a consequence, t he 8 significant digits of single precision would probably not have been sufficient to avoid significant truncation error. Refining the mesh and lengthening the street canyon only exacerbated the problem. The issue was resolved by casting every variable, array, and function as double precision, which provided 16 significant digits plus an exponent Despite the strict use of double precision, CANYON often report ed infinity and overflow errors during initial differential view factor calculations. A closer inspection of Equations 43 and 44 revealed two sources of the problems. For both relationships the equation is undefined i f the denominator of the arctangent argument is z ero. Similarly, the natural log goes to negative infinity as its argument goes to zero. These instances tended to occur when an edge of one of the differential surfaces coincided with one of the axes or when the surfaces overlapped in one or more dimens ions. These problems were circumvented by approximating zero as a very small number so that the expressions would always return finite values. This approach was PAGE 58 58 recommended by the team that derived Equation 4 4 (Ehlert and Smith, 1992) and was validated in CANYON by confirming agreement between the whol e surface view factors and the differential surface view factors calculated at a single mesh (Table 4 5) Surface Meshing Constraints To expedite the programming, CANYON was written with a single mesh param eter that appli ed to all three dimensions. The surface mesh parameter is called mxsfc in Appendix B and can be any integer greater than zero It defines the number of differential elements in each direction on each whole surface (as opposed to defining the number of times each direction is sliced). For example, if mxsfc is 4 there will be 16 differential surfaces on each whole surface for a total of 96 differential surfaces (Figure 4 8 Even with a single mesh parameter, CANYON required many sets o f nested DO loops with levels that terminated at mxfc If the X Y , and Z direction meshes had been different, each one of those loop terminations would have depended on which whole surface was involved Nearly every one of the nested do loops would have then been repeated through a series of IF statements to account for the particular mesh parameter. While a single mesh parameter was not the most flexible numerical approach, it simplified the programm ing enough to complete CANYON in a timely manner. Temperature Initialization and Source Definition After calculating and verifying all the differential view factors, the next step in developing CANYON was to initialize the wall temperatures and implement the source term Initial t emperatures were first read from the input file (Appendix C ) for each whole PAGE 59 59 surface and then applied uniformly to the appropriate set of differential surfaces .2 For simplicity, the baseline initial temperature f or surfaces 2 t hrough 6 was 300K, but any initial temperature could have been specified. Results are not significantly affected by m inor changes in the initial wall temperature s due to the high temperatures of the fireball source. The source temperature (surface 1) wa s defined in a similar manner. If the black body source temperature method was selected, Equation 32 was used to define the temperature for any particular time. If the direct temperature method was selected, the corresponding fireball temperature history data was stored in a text file to be read by CANYON The data was stored as normalized time in the first column and the corresponding tem perature in the second column.3 At the beginning of each time step CANYON performed a linear interpolation on the two nearest data points in the temperature history and then applied the interpolated temperature to every differential source surface. With a temperature defined for every differential surface, black body emissive powers ( bE ) were cal culated using Equation 29 as the next step towards determining the flux profile. Time Step s The desired time step ( t ) was read from the CANYON input file and converted to a normalized time step ( tnorm) using tmax. With the total pulse length defined as previously discussed, the number of steps was determined by dividing the total normalized time (20.5) by the normalized time step ( tnorm). For example, at the 2 The other surface and material properties were also defined in this manner (i.e Cp and k ) 3 By storing the source temperatures vs. normalize time, tmax was easily used to scale the profile to yield. PAGE 60 60 standard pulse length of 20.5tmax and a t of 1x105 s the number of time steps at various yields are shown in Table 46.4 Calculations for smaller yields were expected to have shorter run times because fewer steps were required to simulate the same normalized pulse length. Flux Profile Resolution Given the view factors, emissive powers, and emittances for every differential surface, the gray diffuse enclosure assumption was then implemented to determine the steady state flux profile at each time step. A variety of matrix methods were available to solve Equation 2 17 The first approach was direct inversion of the C matrix by Gauss Jordan elimination followed by matrix multiplication. When that method proved insufficient for finely meshed problems, Gauss Seidel iteration was used t o solve for the flux profiles. Direct Inversion Direct inversion via Gauss Jordan elimination is a method whereby the identity matrix (I ) is appended to the right side of a square matrix (C ) that is to be inverted. Elementary r ow operations are then used to switch the identity matrix to the left side of the augmented matrix, leaving the inverse of the original matrix on the right (C1). The inverse matrix can then be used to solve the system of equations via matrix multiplication. This process is summarized below in Equations 45 and 46 CANYON did not consider any external sources of radiation, so oh was dropped from Equation 2 17. 1 1 1C C C C C C ] [] [ ] [ I I I (4 5) 4 A time step of 1x105 s was used during code development. An appropriate t was verified through the convergence studies in Chapter 5. PAGE 61 61 b 1 b 1 1 be A C q e A C q C C e A q C (4 6) This method was implemented in CANYON through a nested loop of row operations. After the identity matrix was appended to C rows of the augmented matrix were sorted in descending order by the value of the first column. Each row was then divided by its le ading value if it was nonzero so that the first column became entirely 1s or 0s. The top row was then subtracted from each of the rows below with non zero leading values leaving the top row with a leading value of 1 and lower rows with leading value s of zero. This process was repeated down the diagonal of the matrix until the diagonal was all 1s with 0s below it. A similar process was used to resolve the upper half of the matrix. The inverse matrix was then detached and used as shown in Equation 46. This algorithm was extensively validated on simple test matrices. Early testing of CANYON show ed that direct inversion was only useful for very coarse meshes. For mesh parameters greater than 2, the code could not resolve the flux profiles after mo re than a few time steps reporting overflow errors instead. Because the direct inversion algorithm was independently validated, these errors suggested that CANYON was trying t o solve an ill conditioned system of equations which is where a small change in the coefficient matrix produces a large change in the solution vector ( Kaw and Kalu, 2009). An iterative solution method was then investigated so that finer meshes could be used. GaussSeidel Iteration Gauss Seidel iteration is a method where a system of equations is solved progressively starting from an initial guess of the solution vector. As each row of the system is solved, it is used in the solution of the next row so that the algorithm is always PAGE 62 62 o perating with most up to date information. A new solution vector is thus calculated and then stored as the next guess The process is repeated until the answer no longer changes within a certain tolerance. Equation 4 7 shows t he general expression for a particular element of the solution vector (Duderstadt and Hamilton, 1976) Equation 48 shows how that expression was applied to Equation 216. 1 1 1 ) ( ) 1 ( ) 1 (1i j N i j m j ij m j ij i ii m ia a S a (4 7) Where; ) 1 ( m i = Element being solved in current (m+1) iteration ) 1 ( m j = Previously solved element in current (m+1) iteration ) ( m j =Element from previous (m) iteration. N i j m j ij i j m j ij N j bj ij ii m iq C q C E A C q1 ) ( 1 1 ) 1 ( 1 ) 1 (1 (4 8) The Gauss Seidel method was implemented in CANYON by calculating the coefficient matr ices and guessing a flat flux profile of unity before beginning the time steps The converged flux profile from one time step was then stored as the initial guess for the next time step, which helped reduce the required number of iterations as the simulat ion progressed. Preliminary test cases showed consistent solutions when the Gauss Seidel tolerance (Equation 4 9) was tightened, indicating that the algorithm properly converged. Chapter 5 addresses this convergence in more detail. GS m i m i m i qtol q q qi 1 (4 9) Where; iq = relative error for element qi in current (m+1) iteration tolGS = Gauss Seidel tolerance (e.g. 1x105) PAGE 63 63 Temperature Updating Once the flux profile was known, a new temperature was calculated for each differential surface using the semi infinite solid method. The heat flux within each time step was assumed to be constant and t was replaced by t in Equation 221. The new temperatures were substituted for the old temperatures and then the code calculated the flux profile for the next time step. Fluence was calculated for each time step as the flux multiplied by t It was then stored for summation with the fluence of the next time step. A more sophisticated integration routine was not required because this summation was done at every time step (i.e. hundreds of thousands of times) Non Reflecting Boundary Conditions Even under the continuous passage assumption, a street canyon is not actually completely enclosed as it has been described so far. In reality a canyon is open to the sky and does not just end with an opaque surface at the chosen exit plane. Energy will stream out of the canyon at these surfaces instead of being reflected or emitted back into the problem. It would be more realistic to cast the sky and exit as nonreflecting boundar ies .5 This could have easily been accomplished by skipping surfaces 2 and 6 in the flux profile calculation, but that approach would not have worked because the ultimate result depended on tr acking the flux at sur face 2 Instead, a workaround was designed to allow CANYON to still take advantage of the gray diffuse enclosure equation without modeling the problem as completely enclosed. Another way to describe a nonreflecting surface is to say that it absorbs every thing and emits nothing. Recalling that = for gray diffuse surfaces, this perfect 5 A partially reflective and partially transparent surface would be best, but such a condition could not have been adapted as easily to the assumpt ion of gray diffuse surfaces. PAGE 64 64 absorption with no emission concept was used to model the nonreflecting condition. A flag was added to the input file to designate each surface as on or off. If a surface was on, it was treated normally with the designated from the input file. If a surface was off, its was automatically set to 1 so that it absorbed all incident energy. The only way for an off surface to then return energy to the problem was b y emission at its current temperature. Without further treatment, off surfaces would ha ve been black body radiators that increased in temperature by the semi infinite solid assumption. Instead, the temperatures of off surfaces were held constant at their initial values In this way, absorbed energy was not reemitted but the off surfaces still functioned as ambient radiant boundaries. Benchmarking Before CANYON was used to analyze t he fireball scenario, a series of simpler benchmark problems were run to validate the programming logic and accuracy. CANYONs results for the fireball scenario could be believed with more confidence if the code was shown to match progressively more invol ved cases with known or independently verified solutions. Ideally, CANYON would be compared to experimental results, but such a benchmark is impractical for the reasons discussed in Chapter 1. Steady State Benchmark A textbook example was used to verify the steady state flux calculation ( engel, 2007). The example was for a 5 m cube of black surfaces where the base, top, and sides were at uniform temperatures of 800 K, 1500 K, and 500 K respectively (Figure 4 9 ). According to the benchmark, the net rat e of heat transfer to surface 5 ( 5Q ) was 925 kW. CANYON reported a value of 924.3 kW which was lower only because PAGE 65 65 CANYON used a more precise Stefan Boltzmann constant and more precise view factors. The CANYON result was the same for both flux methods at 1, 2, and 4 meshes in each direction (Table 4 7). This benchmark indicated that CANYON correctly calculated the flux profile for a simple steady state case for a variety of meshes and either flux method .6 Steady State Integration Benchmark To test the integration method in CANYON, the steady state benchmark temperatures were held constant over the length of the 3.2 s pulse. The total energy transfer should have been 925 kW x 3.2 s, or 2960 kJ. CANYON consistently reported a total en ergy transfer of 2958 kJ. This benchmark indicated that CANYON correctly integrated a steady state flux profile. Transient Benchmark The final benchmark was a comparison against a transient model independently constructed in Tk Solver 5.0 The whole surface view factor for parallel plates (Equation 4 1) was used with the 20 kT black body fireball temperatures and the equation for radiant exchange between two black bodies (Equation 213) The plates were given the baseline canyon entrance dimensions. W ith no other view factors in the model, it was like the street canyon with no walls, street or sky The exit plane tempe rature was held constant at 300K. Steady state fluxes were calculated for each source temperature from the black body profile. The f luxes were then put in a list function with cubic mapping. The fluence was calculated via 7 point Gaussian integration like in previous Tk Solver models. 6 Direct inversion worked for meshes greater than 2 because the cubic benchmark geometry did not introduce the same numerical difficulties that came from the high aspect ratio of the street canyon. PAGE 66 66 To achieve parity from CANYON, the nonreflecting boundary condition was applied to every surface but the fireball. The walls, street, and sky were held at 0 K and the exit plane was held at 300K. The scenario was run for a single mesh on each surface as additional validation of the view factor calculation. The parameters and results for both models are summarized in Table 4 8 The results in Table 4 9 are quoted to extra significant digits to emphasize the numerical agreement between the models. A flux plot is not included here because the two methods were visually indistinguishable. This benchmar k indicated that CANYON correctly integrated a transient flux profile It also indicated that the whole surface view factor for equal parallel plates matches the single mesh differential view factor fo r rectangles in parallel planes Collectively, t he be nchmarks showed that CANYON delivered reliable results for simple test models. It was then necessary to examine the results at finer m eshes with the walls in place. Figure 41. Coordinate system for CANYON. Z X Y 1 2 3 4 5 6 dx dy d z PAGE 67 6 7 Table 4 1. Surface definitions for CANYON Surface number Symbolic dimensions Description 1 dx, dz Plane source 2 dx, dz Exit plane 3 dy, dz Left wall 4 dy, dz Right wall 5 dx, dy Street 6 dx, dy Sky Figure 42. View factor notation for parallel plates of equal size. (Source: Howell, J.R., 2001. A Catalog of Radiation Heat Transfer Configuration Factors. 2nd Ed. University of Texas at Austin, Austin, TX. (http://www.me.utexas.edu/~howell/sectionc/C 11.html Last accessed September 2009)). Figure 43. View factor notation for perpendicular plates with a common edge of equal length. (Source: Howell, J.R., 2001. A Catalog of Radiation Heat Transfer Configuration Factors. 2nd Ed. University of Texas at Austin, Austin, TX. (http://www.me.utexas.edu/~howell/sectionc/C 14.html Last accessed September 2009)). PAGE 68 68 Table 4 2. Whole surface vi ew factor solution method summary F n m m n 1 2 3 4 5 6 1 0 # # SY SM/SY SM/SY 2 R 0 SY SY SY SY 3 R R 0 # SM/SY SM/SY 4 R R R 0 SM/SY SM/SY 5 R R R R 0 SM 6 R R R R R 0 # : direct calculation; SY: symmetry; SM: summation, R: reciprocity Table 4 3. Whole surface view factors for debugging dimensions F n m m n 1 2 3 4 5 6 1 0 0.02392 0.31780 0.31780 0.17024 0.17024 2 0.02392 0 0.31780 0.31780 0.17024 0.17024 3 0.06356 R 0 0.52993 0.17147 0.17147 4 0.06356 R 0.52993 0 0.17147 0.17147 5 0.06809 0.06809 0.34295 0.34295 0 0.17792 6 0.06809 0.06809 0.34295 0.34295 0.17792 0 Figure 44. Differential surface meshing schematic. Z X Y 3 dx dy d z 5 PAGE 69 69 Figure 45. View factor notation for rectangles in parallel planes. (Source: Howell, J.R., 2001. A Catalog of Radiation Heat Transfer Configuration Factors. 2nd Ed. University of Texas at Austin, Austin, TX. (http://www.me.utexas.edu/~howel l/sectionc/C 13.html Last accessed September 2009)). Figure 46. View factor notation for rectangles in perpendicular planes. (Source: Howell, J.R., 2001. A Catalog of Radiation Heat Transfer Configuration Factors. 2nd Ed. University of Texas at Austin, Austin, TX. (http://www.me.utexas.edu/~howell/sectionc/C 15.html Last accessed September 2009)). PAGE 70 70 Figure 47. Equivalent whole surface and differential view factors example. Table 4 4. Differential view factors from surface 5 to 3 for debugging dimensions VFdiff(iFROM, n1, n2; iTO, m1, m2) Value VFdiff( 5, 1, 1; 3, 1, 1) 0.3429471859D+00 VFdiff( 5, 1, 1; 3, 1, 2) 0.3783171193D 01 VFd iff(5, 1, 1; 3, 2, 1) 0.1921299123D 01 VFd iff(5, 1, 1; 3, 2, 2) 0.8594659693D 02 VFd iff(5, 1, 2; 3, 1, 1) 0.1921299123D 01 VFd iff(5, 1, 2; 3, 1, 2) 0.8594659693D 02 VFd iff(5, 1, 2; 3, 2, 1) 0.3429471859D+00 VFd iff(5, 1, 2; 3, 2, 2) 0.3783171193D 01 VFd iff(5, 2, 1; 3, 1, 1) 0.1576919242D+00 VFd iff(5, 2, 1; 3, 1, 2) 0.7731022698D 01 VFd iff(5, 2, 1; 3, 2, 1) 0.2222177005D 01 VFd iff(5, 2, 1; 3, 2, 2) 0.2008390183D 01 VFd iff(5, 2, 2; 3, 1, 1) 0.2222177005D 01 VFd iff(5, 2, 2; 3, 1, 2) 0.2008390183D 01 VFd iff(5, 2, 2; 3, 2, 1) 0.1576919242D+00 VFd iff(5, 2, 2; 3, 2, 2) 0.7731022698D 01 Table 4 5. View factor method comparison for approximation of zero as small number View Factor Whole surface method Differential surface method F 1 2 0.00049 0.4904689700D 03 F 1 3 0.34523 0.3452283617D+00 F 1 5 0.15453 0.1545264038D+00 By symmetry, F1 3 = F1 4 and F1 5 = F1 6. Z X Y 3 d x =40 m d y =200 m dz=80 m 5 mz=1 nx=1 my=1 ny=1 PAGE 71 71 Figure 48. Surface meshing example for mxsfc = 4. Table 4 6. Number of time steps at various yields for t of 1x105 s. Yield [kT] t max [s] t norm Number of steps 20 0.1558 6.418E 05 319406 10 0.1149 8.707E 05 235445 5 0.0847 1.181E 04 173554 Figure 49. A 5 m cube of black surfaces. Surface numbers here reflect the conventions of the source publication but were changed within CANYON to match the aforementioned reference system. Table 4 7. Steady state benchmark results for CANYON Heat transfer to surface 5 [kW] Mesh parameter Direct Inversion Gauss Seidel Iteration 1 924.3 924.3 2 924.3 924.3 4 924.3 924.3 Z X Y 5 2 3 T 6 = 1500 K T1,2,3,4 = 500 K T 5 = 800 K 5 m 5 m 5 m PAGE 72 72 Table 4 8. Transient benchmark parameters Parameter Tk Solver CANYON Units Canyon length 1610 1610 m Canyon width 40 40 m Canyon height 100 100 m Exit temperature 300 300 K 7 pt Gauss tolerance 1.0E 05 1 Number of subintervals 13 1 Time step 1E 05 s Gauss Seidel tolerance 1E 05 1 Table 4 9. Transient benchmark results Result Tk Solver CANYON Units View Factor 4.90468970E 04 4.90468971E 04 1 Fluence 1.1932 1.1929 cal/cm 2 PAGE 73 73 CHAPTER 5 CANYON ANALYSIS & RESULTS Baseline Parameters Summary In order to evaluate parameter changes in CANYON, it was desirable to have a baseline set of input variables. These base line parameters were developed and justified throughout Chapters 2, 3 & 4. They are collectively summarized in Table 51. Unless otherwise noted, all cases in this chapter used the open sky and exit plane geometry (Figure 51). An initial mesh parameter is not specified in Table 51 because it was determined from subsequent convergence studies, along with refined values for the time st ep and Gauss Seidel tolerance. Convergence Studies Before meaningful results could be extracted from CANYON, the simulation parameters had to be refined until the solution converged. The differential view factor approach assumed a uniform temperature acr oss each differential area. In reality, surface temperatures would change continually as a function of distance from the source, so the accuracy of the steady state flux profile depend ed on the size of the differential areas relative to the slope of the a ctual temperature profile along the walls. Solving the flux profiles with Gauss Seidel iteration meant that the tolerance also factored into steady state convergence. Convergence in the transient solution depended additionally on the size of the time step. The source temperature changed rapidly by thousands of degrees. If the time step were too large, the temperature history would not have been sampled enough to accurately capture its steep slope s PAGE 74 74 Mesh The first convergence test was to increase the mesh parameter and observe the effects on the average flux history and total fluence on the exit plane CANYON was run with the baseline inputs at each mesh parameter from 1 through 8. The resulting fluences and runtimes are summarized in Table 5 2 along with percent changes from the previous run. The decreasing differences between each run show ed a converging trend. The same trend was observed in plots of average flux incident on the exit plane versus time (Figure 52) As expected, the flux curves all had the same basic shape as the normalized p ower curve (Figure 22). But, w hen the mesh parameter was set to 4 an unphysical flattening of the peak appeared. This anomaly indicated numerical instability in the run and c ast doubt on the runs at finer meshes. The time step was then modified in an attempt to eliminate the irregularity. Time Step The next convergence technique was refining the time step to more closely capture the temperature history. T he time step was reduced to 1x106 s and CANYON was rerun for mesh parameters of 4, 6, and 7. The resulting runtimes and fluences are shown in Table 5 3. The flux profile for a mesh parameter of 4 no longer showed the numerical anomaly at the peak and the finer meshes s till appeared to be converging (Figure 53) CANYON was then run at a time step of 5x107 s for mesh parameters of 8 and 10 to investigate how tightly the result s converged for progressively finer meshes and time steps The resulting flux profiles were a lmost indistinguishable (Figure 54) A mesh parameter of 8 was selected for use in further runs because of the long runtime for the 10 mesh case and the small change in total fluence (3.36%) from 8 to 10 meshes. PAGE 75 75 Tolerance The final convergence test was a check on the Gauss Seidel tolerance. It was first reduced to 1x106 to see if tighter tolerances would be required, but t here was no appreciable change in the total fluence (Table 5 4). The tolerance was then loosened to 1x104 to see if the same result could be achieved in shorter runtimes. The total fluence under the looser tolerance only changed by 0.002% from the 1x105 case and saved about 30 minutes of runtime so 1x104 was used as the tolerance for all subsequent runs As a result of these convergence studies, the baseline parameters (Table 51) were modified to include a mesh parameter of 8, a time step of 5x107 s, and a Gauss Seidel tolerance of 1x104. The results of subsequent parameter studies were compared to a baseline total fluence of 1.5353 cal/cm2. Transient Wall Temperatures As implied by the benchmark problems and convergence studies, the primary output of CANYON was the total fluence averaged over the exit plane. This value is what would probably be determined in a real forensic scenario. By default, temperature profiles in CANYON are not written to the output file (Appendix D). In the last block of the input file (Appendix C), the user may designate locations on the left wall, right wall, or floor for which he would like to see temperatures.1 The temperatures for the differential surfaces containing each requested location are then printed over the length of the pulse to a separate output file called temps.put For example, temperatures were printed from the baseline case for the middle of the floor and for 7 m high on the 1 When the sky and exit plane are turned off under the non reflecting approximation, their temperatures are not meaningful because they are false surfaces. PAGE 76 76 left wall at 300 m, 900 m, and 1500 m down the canyon at times of 0.5tmax, tmax, and 20.5 tmax (Table 55 ). Table 5 5 emphasizes the large temperature gradient s present. Parameter Studies Outline With a benchmarked and converged model in place, each input parameter was examined for its effect on the total fluence. The first goal was to determine the relative importance of each parameter, which would help assess the validity of the corresponding assumptions. The second goal was to determine the expected fluence under a variety of conditions that might be encountered in a real scenario. Fluences were calculated for changes from baseline for each material property The canyon dimensions were then changed to examine the stationary source assumption, the entrance dimensions, and the decrease of fluence with range. The importance of source temperature selection was then evaluated by comparing r esults for b oth temperature methods. Next, the overall canyon effect was determined through model comparisons. The final studies examined changes in yield and visibility. Material Property Studies One of the potential oversimplifications in CANYON was the uniform material properties on each surface. To test the sensitivity to these properties, each one was altered by 30% and resulting fluence s were compared to the baseline value (Table 5 6 ) The largest change in fluence was less than 0.10% of the baseline value, showing that the CANYON model was largely insensitive to material property changes. This result showed that it is probably not necessary to add a detailed material map or temperature dependent properties to CANYON because doing so would not significantly affect the result. PAGE 77 77 Canyon Dimension Studies The effect of the stationary source assumption was also of interest along with the other canyon dimension assumptions Model r uns were calculated at canyon lengths plus and minus the radius of the fireball T he fluence decreased by 16.37% and increased by 21.93% respectively (Table 5 7 ). The canyon length was then extended to 2415 m (1.5 mi) and 3220 m (2 mi). The total fluences showed large decreases as expected. The most dramatic change in fluence came from doubling the canyon entrance height (a 279.27% increase). Increasing the entrance dimensions directly increases the size of the source, so a larger fluence makes sense. One migh t also expect the doubling of the entrance width to show the same relative increase because the entrance area was doubled in either case. But, the wider entrance only showed an 89.33% increase in fluence because doubling the width also doubled the leakage area through the sky whereas the taller canyon had the same sky dimensions as the baseline case. These changes are emphasized in Figure 55. Source Temperature Selection Study The source temperature selection was examined next for its influence on the re sult. As discussed in Chapter 3, the fireball could have been defined by the direct temperature history ( Tdir) or by the black body temperatures derived from the normalized power curve (Tbb). When the direct temperatures were used instead of the baseline black body temperatures the fluence dropped by 36.62 % to 0.9730 cal/cm2. Except for the first peak (not modeled by the black body temperatures), the direct temperature flux was lower than the black body flux for the entire length of the pulse (F igure 56 ). This result showed that the total fluence was very se nsitive to source temperature, which PAGE 78 78 was expected because temperature was a direct measure of the energy emitted by the source. Canyon Effect and Model Comparisons After examining the sensit ivity of the CANYON result to the baseline parameters, it was compared to the time dependent Open Field Model ( OFM ) and Sparrow Passage Model ( SPM ) results. For the baseline case, the time dependent OFM predicted a fluence of 17.29 cal/cm2 and CANYON predicted a fluence of 1.5353 cal/cm2, or 8.9% of the OFM result. This reduction is similar to the entrance area compared to the cross sectional fireball area (6.9%) covered in Chapter 3 This fluence comparison emphasizes the importance of the source treatm ent. It was initially expected that the street canyon would increase the total fluence delivered to the exit plane relative to the open field scenario, but that expectation did not account for shadowing of the source. Even so, there was still a signific ant canyon effect. When all the walls were set as non reflecting boundaries at 300K (i.e. open field with in CANYON) the total fluence was 1.2282 cal/cm2. The 1.5353 cal/cm2 for the baseline case was a 25 % increase from this value, which showed that the presence of the canyon walls and street should increase the fluence for a given source plane, even if the buildings have shadowed a significant portion of the actual source. As an example, the black body t emperatures were linearly scaled for the open field CANYON case until the fluence matched the time dependent OFM result. The results matched for a temperature factor of about 1.95. When the walls were turned back on, the fluence showed the same 25% increase. In other words, CANYON predicted a 25% canyon effect for the baseline parameters regardless of the source strength. PAGE 79 79 For comparison to the SPM, all the walls (including the sky & exit) were then turned on to see the result of modeling the canyon as a complete enclosure. This scenario produced dramatically higher fluxes and fluence. The total fluence was 2410.7 cal/cm2, or roughly 10 times higher than the SPM fluence and 1570 times higher than the baseline CANYON fluence. The flux profiles showed similar differences (Figure 5 7 ) The large values for the enclosed CANYON model were expected because there was nowhere for the radiated heat to go except into the walls via conduction. Therefore, much more was channeled down the canyon without spreadi ng Yield Studies CANYON would ultimately only be useful if it could predict measurable differences in fluence for different yields. To test this sensitivity, the model was also run at 10 kT and 5 kT which should have scaled the power output and pulse time length according to the normalized power curve (Figure 22) The resulting fluences did show marked reductions from the 20 kT baseline (Table 58 ). The flux profiles also showed shorter pulses for the smaller yields, indicating that CANYON properly normalized time (Figure 5 8 ) The overall canyon effect was also examined for 10 kT and 5 kT cases, which again showed about a 25% increase when the walls were in place (Table 59 ). Visibility Studies A final set of runs was conducted for the visibilities discussed in Chapter 3. The resulting fluences and relative changes were noteworthy but less than one might expect (Table 510 ). Because the transmittance fits were not valid for distances less than or equal to 1 km, CANYON set the short range transmittance to 100% regardless of the specified visibility. Therefore, any exchange between differential surfaces separated by 1 km or less was not subject to atmospheric attenuation. All the preceding param eter PAGE 80 80 studies are described in Table 5 1 1 and t he resul ts are summarized in Table 51 2 Overall, CANYON enabled a variety of tests of the street canyon model. General conclusions from these tests are presented in the next chapter. Table 5 1. Baseline parameters summary for CANYON Parameter Symbol Value Units Yield W 20 kT Canyon dimensions Length dy 1610 m Width dx 40 m Height dz 100 m Material properties Initial wall temperature T 300 K Emittance 0.7 1 Specific heat Cp 840.0 J/(kgK) Density 1800.0 kg/m 3 Conductivity k 1.500 W/(mk) Simulation parameters Time step 1.00E 05 s Gauss Seidel tolerance tol 1.00E 05 1 Figure 51. Baseline geometry for CANYON Table 5 2. Results summary for mesh convergence study in CANYON Mesh parameter Total fluence [cal/cm 2 ] Change from previous run Runtime [s] 1 100.2034 1.89 2 16.8398 83.19% 3.14 3 5.1640 69.33% 8.34 4 2.6056 49.54% 35.62 5 1.8847 27.67% 42.27 6 1.6376 13.11% 110.39 7 1.5371 6.14% 507.31 8 1.4891 3.12% 932.17 40 m 100 m 1610 m PAGE 81 81 Figure 52. Average flux on the exit plane for baseline case mesh convergence in CANYON. Table 5 3. Results summary for time step convergence study in CANYON Mesh parameter Time step [s] Total fluence [cal/cm 2 ] Runtime [hrs] 4 1.00E 05 2.6056 0.01 4 1.00E 06 2.6845 0.07 6 1.00E 06 1.6962 0.37 7 1.00E 06 1.5847 0.87 8 5.00E 07 1.5352 2.86 10 5.00E 07 1.4854 11.23 0.01 0.1 1 10 100 1000 0 0.5 1 1.5 2 2.5 3 3.5Thermal Flux [cal/(cm^2*s)]Time [s] 1 2 3 4 5 6 7 8 PAGE 82 82 Figure 53. Average flux on the exit plane for a mesh parameter of 4 and t of 1x105 s compared to mesh parameters of 4, 6, & 7 and t of 1x106 s. Figure 54. Average flux on the exit plane for a mesh parameter of 7 and t of 1x106 s compared to mesh parameters of 8 & 10 and t of 5x107 s 0.01 0.1 1 10 0 0.5 1 1.5 2 2.5 3 3.5Thermal Flux [cal/(cm^2*s)]Time [s] 4(old) 4 6 7 0.01 0.1 1 10 0 0.5 1 1.5 2 2.5 3 3.5Thermal Flux [cal/(cm^2*s)]Time [s] 7 8 10 PAGE 83 83 Table 5 4. Results summary for tolerance convergence study in CANYON Mesh parameter Time step Tolerance Total fluence [cal/cm 2 ] Runtime [hrs] 8 5.00E 07 1.00E 05 1.5352 2.86 8 5.00E 07 1.00E 06 1.5352 4.07 8 5.00E 07 1.00E 04 1.5353 2.45 Table 5 5. Selected output temperatures for the baseline case in CANYON Normalized time Time [s] Fireball temp [K] Surface T [K] (y=300 m) T [K] (y=900 m) T [K] (y=1500 m) 0.5 0.077 7057.21 Left wall 3737.78 924.68 338.15 Floor 3628.33 949.54 359.95 1.0 0.156 8581.70 Left wall 4562.07 1722.00 575.98 Floor 4428.11 1675.28 623.32 20.5 3.194 2254.04 Left wall 1202.81 494.32 339.96 Floor 1167.56 482.15 336.51 Left wall temperatures were taken 7 m above the floor. Floor temperatures were taken in the middle of the canyon. Table 5 6. Material property sensitivity study results summary for CANYON Cp [J/(kgK)] [kg/m 3 ] k [W/(mk)] Fluence [cal/cm 2 ] Change from baseline 0.70 840.0 1800.0 1.50 1.5353 0.91 1.5351 0.007% 0.49 1.5351 0.012% 1092.0 1.5343 0.065% 588.0 1.5363 0.068% 2340.0 1.5343 0.065% 1260.0 1.5363 0.068% 1.95 1.5343 0.065% 1.05 1.5363 0.068% A indicates no change from the baseline value in the first row. Other properties are changed by 30% from baseline. PAGE 84 84 Figure 55. Canyon entrance dimension changes. A) Baseline geometry. B) Double height geometry. C) Double width geometry. Width and height are to scale, but length is not. Table 5 7. Canyon dimension studies results for CANYON Lengt h (dy) [m] Width (dx) [m] Height (dz) [m] Fluence [cal/cm 2 ] Change from baseline 1610 40 100 1.5353 1474 1.8719 21.93% 1746 1.2840 16.37% 2415 0.6464 57.89% 3220 0.3661 76.15% 80 2.9067 89.33% 20 0.9129 40.54% 200 5.8828 279.27% 50 0.0207 98.65% A indicates no change from the baseline value in the first row. 40 m 40 m 8 0 m 10 0 m 20 0 m 100 m 161 0 m 1610 m 161 0 m A sky 2 A sky A entrance 2 Aentrance A B C PAGE 85 85 Figure 56. Comparison of average exit plane flux vs. time for black body and direct source temperatures. Figure 57. Comparison of average flux on exit plane vs. time for SPM, OFM, and enclosed & baseline CANYON models. 0.01 0.1 1 10 0 0.5 1 1.5 2 2.5 3 3.5Thermal Flux [cal/(cm^2*s)]Time [s] Black Body Direct 0.01 0.1 1 10 100 1000 10000 0 0.5 1 1.5 2 2.5 3 3.5Thermal Flux [cal/(cm^2*s)]Time [s] CANYON enclosed SPM OFM CANYONbaseline PAGE 86 86 Table 5 8. Total fluence for 20 kT, 10 kT, & 5 kT at 1610m Yield [kT] Fluence [cal/cm 2 ] Change from baseline 20 1.5353 10 1.1297 26.42% 5 0.8302 45.92% Table 5 9. Canyon effect for 20 kT, 10 kT, & 5 kT at 1610m Yield [kT] No wall fluence [cal/cm 2 ] Baseline geometry fluence [cal/cm 2 ] Relative increase with walls in place 20 1.2282 1.5353 25.00% 10 0.9053 1.1297 24.78% 5 0.6674 0.8302 24.41% Figure 58. Average exit plane flux vs. time at 1610 m for 20 kT, 10 kT, & 5 kT 0.01 0.1 1 10 0 0.5 1 1.5 2 2.5 3 3.5Thermal Flux [cal/(cm^2*s)]Time [s] 20 kT 10 kT 5 kT PAGE 87 87 Table 5 10. Total fluence for sea level visibilities of 5, 10, 15, & 23 km compared to baseline. Visibility [km] Total fluence [cal/cm 2 ] Change from baseline 1.5353 23 1.3861 9.71% 15 1.3698 10.78% 10 1.3545 11.77% 5 1.3340 13.11% Table 5 11. Description of each parameter study Test Number Type of Study Description A1 Material property Emittance increased by 30% A2 Emittance decrease by 30% A3 Specific heat increased by 30% A4 Specific heat decreased by 30% A5 Density increased by 30% A6 Density decreased by 30% A7 Conductivity increased by 30% A8 Conductivity decreased by 30% B1 Canyon dimension Length decreased by fireball radius (135.7 m) B2 Length increased by fireball radius (135.7 m) B3 Length increased to 1.5 mi B4 Length increased to 2.0 mi B5 Width doubled B6 Width halved B7 Height doubled B8 Height halved C1 Temperature method Direct Temperatures ( T dir ) used instead of T bb D1 Canyon Effect All walls off at 300K D2 All walls on D3 Transient Open Field Model D4 Sparrow Passage Model E1 Yield Yield reduced to 10 kT E2 Yield reduced to 5 kT E3 Yield reduced to 10 kT with all walls off at 300K E4 Yield reduced to 5 kT with all walls off at 300K F1 Visibility Sea level visibility set to 23 km F2 Sea level visibility set to 15 km F3 Sea level visibility set to 10 km F4 Sea level visibility set to 5 km All relative parameter changes are in reference to the baseline values in Table 51. PAGE 88 88 Table 5 12. Summary of results for parameter studies and model comparisons Test Number Type of study Parameter change from baseline Total fluence [cal/cm 2 ] Change from baseline Baseline 1.5353 A1 Material property 1.5354 0.007% A2 1.5351 0.012% A3 Cp = 1092 J/(kgK) 1.5343 0.065% A4 Cp = 588 J/(kgK) 1.5363 0.068% A5 3 1.5343 0.065% A6 3 1.5363 0.068% A7 k = 1.95 W/(mk) 1.5343 0.065% A8 k = 1.05 W/(mk) 1.5363 0.068% B1 Canyon dimension dy = 1474 m 1.8719 21.93% B2 dy = 1746 m 1.2840 16.37% B3 dy = 2415 m 0.6464 57.89% B4 dy = 3220 m 0.3661 76.15% B5 dx = 80 m 2.9067 89.33% B6 dx = 20 m 0.9129 40.54% B7 dz = 200 m 5.8228 279.27% B8 dz = 50 m 0.0207 98.65% C1 Temp method Direct temperatures 0.9730 36.62% D1 Canyon effect and model comparison All walls off (at 300K) 1.2282 20.00% D2 All walls on 2410.7 157018.17% D3 Transient OFM 17.29 1126.16% D4 SPM 234.76 15290.82% E1 Yield W = 10 kT 1.1297 26.42% E2 W = 5 kT 0.8302 45.92% E3 W = 10 kT, Walls OFF 0.9053 N/A E4 W = 5 kT, Walls OFF 0.6674 N/A F1 Visibility Vis = 23 km 1.3861 9.71% F2 Vis = 15 km 1.3698 10.78% F3 Vis = 10 km 1.3545 11.77% F4 Vis = 5 km 1.3340 13.11% N/A: These two test are not directly comparable to the baseline case. PAGE 89 89 CHAPTER 6 CONCLUSIONS Through the theoretical development and numerical examination of a wide range of scientific assumptions this investigation created a program called CANYON for modeling the transmission of thermal energy from nuclear detonations in urban environments. Results showed that while most of the assumptions turned out to be reasonable, some of them require further development before consistent fluence estimates can be made. In a real forensic scenario, material properties would be difficult to measure and highly variable throughout the street canyon. Fortunately, the street canyon model was in sens itive to material property changes as large as 30% On the other hand, the model was more sensitive to canyon dimensions and atmospheric conditions U nlike material properties, actual canyon dimensions could easily be measured in a real scenario. The gr owth of the fireball could also be modeled with a quasi static approach, where the length of the canyon is changed at distinct stages. Atmospheric cond itions could be accounted for using visibility estimates from the local weather authority. T he total t hermal fluence at the far end of the canyon was ultimately driven by the source strength regardless of street canyon properties. The fluence dropped by about 26% for each halving of the yield from 20 kT to 10 kT to 5 kT. Recalling that regardless of yield the fireball will nearly always fill the canyon entrance, simple 50% reductions were not seen because smaller yields produce a tradeoff between a shorter pulse that decreases the integration time and a smaller fireball that increases black body temperat ure. For the stationary plane source fireball approximation, the street canyon increased the fluence by about 25% compared to the same source in an open field over PAGE 90 90 the range of yields considered. An accurate forensic measurement would depend on knowing how the fireball interacts with the immediately surrounding buildings to either degrade or enhance the intensity of thermal energy released into the street canyon. The CANYON program could be modified to model the transmission of that energy. Therefore, t he utility of this investigation is contingent upon an accurate determination of the source term via some other research effort. PAGE 91 91 CHAPTER 7 FUTURE WORK This investigation produced meaningful results but also established a framework for many additional developments. Future work should focus in two basic areas: improving the efficiency and flexibility of CANYON and using it to perform more detailed analyses. Pending the accomplishment of some or all of these suggestions, CANYON will be a powerful tool for post event forensic assessments. CANYON Improvements The mesh parameter was limited in practice by run times CANYON was written in a straightforward manner without direct consideration of memory usage or calculation optimization. A variety of changes could help reduce runtimes to allow finer meshes and facilitate more detailed parameter studies. The matrices developed from finite differencing schemes are frequently sparse with simple diagonal dominances ( e g. tridiagonal or pentadiagonal matrices) that allow condensed operations instead of storing and working with large arrays The coefficient matrices in CANYON do not have simple dominances but are still sparse in certain areas. A thorough analysis of these matrices could reveal indexing routines that might significantly shorten runtimes. For finely meshed problems, each time step in CANYON takes appreciable computational effort. A parallel versi on of CANYON could reduce runtimes dramatically for these large problems. There would be a tradeoff between parallel overhead and speed gains, but the potential for improvement bears investigation. Additional efficiency and general flexibility would be gained by incorporating variable meshing into CANYON. Under the uniform mesh parameter, the aspect ratio of a differential surface was not changed by finer meshes. With different parameters in PAGE 92 92 each direction more meshes could be added along the length of the canyon without adding unnecessary horizontal and vertical meshes. It would be even more appropriate to program individual mesh parameters for each of the whole surfaces. It is not necessary to m esh the fireball source as finely as the other surfaces because it has a spatially uniform temperature throughout the calculation. Flexible meshing would also facilitate detailed material maps and cross street gap geometries. Aside from efficiency improve ments, CANYON has room for major functional expansion as well The inclusion of partially specular surfaces and the corresponding specular view factors would tie the model closer to reality. Specular view factors can be cumbersome in 2D, so it is worth investigating alternate methods like ray tracing or Monte Carlo to calculate them in 3D. Lastly, it would be valuable to add the logic for modeling surfaces as partially transparent. Surfaces are currently either completely transparent or completely op aque. With partial transparency, buildings could be modeled with glass surfaces and the transmissive nature of the sky could be more accurately captured. Further Analyses The aforementioned code improvements would lend themselves to corresponding validat ion analyses. In addition to those, the improved CANYON should be used to develop detailed fluencedistance yield relationships. The forensic method this investigation has helped develop ultimately depends on back calculating the yield from a n experiment ally measured fluence and distance. This investigation has focused on the forward calculation of fluence, so further research should focus on the backward calculation of yield. PAGE 93 93 Lastly, a formal propagation of uncertainty should be carried out for the CANY ON result The 25% margin on the normalized power curve was the only directly available uncertainty estimate, so other s would have to be assumed and propagated into a final uncertainty for the fluence. Such an estimate is necessary for comparing this approach to other forensic methods. PAGE 94 94 APPENDIX A TRANSMITTANCE CURVE FIT COEFFICIENTS Atmospheric attenuation was accounted for using previously determined transmittance curve fits as described in Chapter 3 (Leonard, 1985). The fits were 5th order polynomials of the form in Equation A 1. The corresponding coefficients for urban aerosols at zero altitude and different sea level visibilities are shown in Table A 1. The curves are only valid for 1 km < x 25 km. 5 6 4 5 3 4 2 3 2 1) (ln ) (ln ) (ln ) (ln ) (ln x C x C x C x C x C C (A1) Where: = atmospheric transmittance x = distance from ground zero [km] Table A 1. Transmittance curve fit coefficients Visibility [km] C 1 C 2 C 3 C 4 C 5 C 6 5 4.654E 01 2.656E 01 1.320E 01 1.405E 01 3.884E 0 2 3.594E 03 10 6.253E 01 1.728E 01 1.556E 01 6.188E 02 8.687E 04 1.228E 03 15 7.138E 01 2.088E 01 9.459E 03 5.577E 02 3.123E 02 4.118E 03 23 7.778E 01 2.040E 01 5.493E 02 9.103E 02 3.299E 02 3.431E 03 PAGE 95 95 APPENDIX B CANYON SOURCE CODE Canyon21 (serial) DEVELOPMENT NOTES: totsfc=6 hardwired into source code & removed as input option Wall temperature outputs added Black body temperature equations embedded (vice reading BB temperature file Cleaned up for inclusion in Master's Thesis. Steady state & constant T flags added (for easy benchmarking). "Return Factor" a dded to surface 6 emissive power. All TABS removed and line length limited to 72 cols. Copy of Canyon 17 to be parallelized Invisible Surface Option added (jVis). When jVis=0 and eps=1.0, a surface absorbs all incident energy and emits nothing back, making it "invisible". Atmospheric attenuation added (with input file selections). Choice of direct or iterative flux method from input file. Calculate & display surface 2 average flux. Display all output units. Change flux output units to [cal/(cm^2*s)] Unify flux file into general output file. Entire flux profile for Surface 2 output at each write step. Minor changes made f rom Canyon 12 hoping to speed up the code. Input reduced to a single mesh parameter (uniform # of meshes in each direction. Iterative solver for flux in each times step instead of direct matrix inversion Fireball Temp interpolator fixed so that T is constant at first datapoint value for times before first data point. Walls approximated as semi infinite solids using EQ 446 from Cengel pg.243 calculates a temperature distribution along the street canyon in terms of total heat instead of flux ALL CALCULATIONS done in SI UNITS. Gray, diffuse surfaces using MODEST Eqn. 5.38. DIFFERE NTIAL METHOD: This version uses EQN C 15 from Radiation Configuration Factors website to calculate VFs for plates in perpendicular planes MODULE PARAMS PAGE 96 96 Set up Problem Parameters INTEGER :: t otsfc, msfc, mxsfc, numRows, ndiffsfc,numout, & nsfc,icount,iError,iFluxChoice,ivis,iTrans,iTsource, & numtemppoints CHARACTER*5 :: VDC, TMPC CHARACTER*1 :: dum DOUBLE PRECISION :: Yield, tmax, deltaTime, deltaTimeNorm,dij, &tol,Sfc2Ave, Flue2Ave, FluxRow1INavg, tfact, skyfact, Pmax, sig, &Amax, tnormmax CHARACTER*19 :: dtinit, dtfinl END MODULE PARAMS MODULE CANDATA CHARACTER*64 :: filein,tfilein CHARACTER*64, DIMENSION (:), ALLOCATABLE :: Description Set arrays for Surface Dimensions DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE :: &dx, dy, dz, Area, VFt, alpha TH, T, Q, qdotw, ddx, ddy, ddz, sX,sY, &sW, sZ, Cp, eps, diffArea, rk, rho,qdot, error, vertavg, Xloc, &Yloc, Zloc, TPxloc, TPyloc, TPzloc, TPtemp DOUBLE PRECISION :: rdbl, PI DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE::VF,TempTim e,Cij,Aij, &CijINV, Cm, Dm, empwr, qdotNETd, qtemp, qdotNETdO VFdiff is the 6dimensional array containing the view factors from one differential area to every other differential area on all 6 surfaces. It will be completely filled in future v ersion of the code. DOUBLE PRECISION, DIMENSION (:,:,:,:,:,:), ALLOCATABLE :: VFdiff, &dist DOUBLE PRECISION, DIMENSION(:,:,:),ALLOCATABLE :: emitpwr,Tdiff, &qdotINdiff, qdotOUTdiff, qdotNETdiff, qvoldiff, TdiffNEW, &FluenceIN, F luenceNET INTEGER, DIMENSION (:), ALLOCATABLE ::isfcs,jVis,n1TP,n2TP,TPsfc INTEGER, DIMENSION (:,:), ALLOCATABLE :: ixy END MODULE CANDATA  PAGE 97 97 PROGRAM CANYON USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) MAIN PROGRAM PI = 2.D0*DASIN( 1. D0) Print program header CALL HEADR Enter input file !! WRITE(*,*) Enter input file name:' !! READ(*,*) filein filein='canyon.inp' Start timing routine CALL CLKDTO (dtinit, hstart,0.D0) WRITE (1,*) 'Ru n Started:',dtinit Open other files OPEN(8,FILE = filein) OPEN(9,FILE = 'VFdiffDump.put') Read input file CALL READ_INP IF (iError.NE.0)THEN GOTO 999 ENDIF March through view factors CALL VIEW_FACTORS Calculate location of each differential surface & distances between each pair CALL LOCATE CALL DISTANCE Calculate differential surface locations of requested temperature PAGE 98 98 points CA LL TEMP_REQUEST Run Heat Transfer calculation & time steps CALL HEAT_TRANSFER WRITE(*,*) WRITE(*,*) WRITE(*,*) 'Output saved to "CANout.put"' End timing routine CALL CLKDTO (dtfinl, hstop,0.D0) fsec=(hstop hstart)*0.0100 WRITE(*,*) '' WRITE(1,*) WRITE(*,*) 'Run Completed:', dtfinl WRITE(1,*) 'Run Completed:', dtfinl WRITE(*,*) 'CPU time est is ',fsec,' sec.' WRITE(1,*) 'CPU time est is ',fsec,' sec.' WRITE(*,*) WRITE(*,*) 'Enter VFDiff Checker? (y or n)' READ(*,*) VDC IF(VDC .EQ. 'n' .OR. 'N' .OR. 'No' .OR. 'no') THEN GOTO 110 ELSE CALL VF_DIFF_CHECKER ENDIF !110 WRITE(*,*) 'Enter Temperature Integration Checker? (y or n)' READ(*,*) VDC IF(VDC .EQ. 'n' .OR. 'N' .OR. 'No' .OR. 'no') THEN GOTO 111 ELSE !112 WRITE(*,*) 'Enter desired normalized time' READ(*,*) timenorm CALL TIME_TEMP(timenorm,tempint) WRITE(*,*) 'tempint = ', tempint WRITE(*,*) 'Again? (y or n)' READ(*,*) VDC IF(VDC .EQ. 'n' .OR. 'N' .OR. 'No' .OR. 'no') THEN GOTO 111 ELSE PAGE 99 99 GOTO 112 ENDIF ENDIF 999 END PROGRAM CANYON ! SUBROUTINE ALLOCARR USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) Allocate all arrays ALLOCATE(dx(totsfc), dy(totsfc), dz(totsfc), Ar ea(totsfc), &VFt(totsfc),alphaTH(totsfc),isfcs(totsfc), T(totsfc), Q(totsfc), &qdotw(totsfc), ddx(totsfc), ddy(totsfc),ddz(totsfc), Cp(totsfc), &eps(totsfc), diffArea(totsfc),rk(totsfc), &rho(totsfc), qdot(totsfc), vertavg(mxsfc)) ALLOCATE(VF(totsfc,totsfc), Cij(ndiffsfc,ndiffsfc), &Aij(ndiffsfc,ndiffsfc), CijINV(ndiffsfc,ndiffsfc), &Cm(ndiffsfc,2*ndiffsfc), Dm(ndiffsfc,2*ndiffsfc), &empwr(ndiffsfc,1), qdotNETd(ndiffsfc,1), qtemp(ndiffsfc,1), &qdotNETdO(ndiffsfc,1), error(ndiffsfc)) ALLOCATE(sX(2),sY(2),sW(2),sZ(2)) ALLOCATE(VFdiff(totsfc,mxsfc,mxsfc,totsfc,mxsfc,mxsfc), &dist(totsfc,mxsfc,mxsfc,totsfc,mxsfc,mxsfc)) ALLOCATE(emitpwr(totsfc,mxsfc,mxsfc),Tdiff(totsfc,mxsfc,mxsfc), &qdotINdiff(totsfc,mxsfc,mxsfc), qdotOUTdiff(totsfc,mxsfc,mxsfc), &qdotNETdiff(totsfc,mxsfc,mxsfc), qvoldiff(totsfc,mxsfc,mxsfc), &TdiffNEW(totsfc,mxsfc,mxsfc), FluenceIN(totsfc,mxsfc,mxsfc), &FluenceNET(totsfc,mxsfc,mxsfc)) ALLOCA TE(Xloc(ndiffsfc),Yloc(ndiffsfc),Zloc(ndiffsfc)) ALLOCATE(ixy(mxsfc,mxsfc), jVis(totsfc)) ALLOCATE(Description(totsfc)) END SUBROUTINE ALLOCARR ! PAGE 100 100 PURPOSE OF VIEW FACTOR ROUTINE: To calculate all the required view factors for the entire problem NOTES: AUTHOR (date): T.F. Stachitas, March 2009 SUBROUTINE VIEW_FACTORS USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy nsfc=totsfc Set Primary Problem Parameters/Initializations ! Write inputs to the output file & command window WRITE(*,*) WRITE(*,*) 'totsfc=', totsfc WRITE(*,*) 'mxsfc=', mxsfc WRITE(1,52) 'Input File Name:',filein WRITE(1,*) WRITE(1,49) 'totsfc = ',totsfc WRITE(1,47) 'mxsfc =', mxsfc,'(Surface mesh parameter)' WRITE(1,*) 49 FORMAT(2(A9,4X,I2,4X)) 47 FORMAT(A9,4X,I2,3X,A24) msfc=mxsfc WRITE(*,46) 'ndiffsfc = ', ndiffsfc 46 FORMAT(A11,X,I10) Read surface dimensions 48 FORMAT(I5, 4F13.0,3X,A64) WRITE(1,*) 'Surface Dimensions [m]' 50 FORMAT(4(' A5,6X),A9,3X,A11) WRITE(1,50) 'sfc', 'dx', 'dy', 'dz', 'Area[m^2]', 'Description' 51 FORMAT(I5, 100F13.5) 52 FORMAT(A16,3X,A20) PAGE 101 101 Calculate surface areas CALL AREAS DO i=1, nsfc WRITE(1,48) isfcs(i), dx(i), dy(i), dz(i),Area(i), & Description(i) END DO WRITE(1,*) Initialize Surface Temperatures (K) [Steady State] from input file WRITE(1,*) 'Initial Surface Temperatures [K]' DO i=1, nsfc WRITE(1,48) isfcs(i), T(i) END DO WRITE(1,*) Differential View Factors Subroutine Calls Initialize counter & number of differential surfaces iNegVFdiff = 0 iVFdiff = 2*((mxsfc*mxsfc)**2+(mxsfc*mxsfc)**2+(mxsfc*mxsfc)**2) VFdTEST = 0 VFdiff = 0 DO iFROM=1,6 DO iTO=1,6 IF (iFROM .NE. iTO) THEN iPLUS = iFROM + 1 iMINUS = iFROM 1 IF( (MOD(iFROM,2).GT.0).AND.(iTO.EQ.iPLUS)) THEN CALL OPP_SIDE_DIFF(iFROM, iTO, iNegVFdiff, iVFdiff) ELSEIF((MOD(iFROM,2).EQ.0).AND.(iTO.EQ.iMINUS))THEN CALL OPP_SIDE_DIFF(iFROM, iTO, iNegVFdiff, iVFdiff) ELSE CALL COMMON_EDGE_DIFF(iFROM,iTO,iNegVFdiff,iVFdiff) ENDIF ENDIF END DO END DO IF (iNegVFdiff .GT. 0) THEN PRINT*, 'ERROR: negative differential view factors.' END IF PRINT*, iNegVFdiff, 'negative differential VFs calculated' PRINT*, iVFdiff, 'differential VFs calculated' PAGE 102 102 Dump VFdiffs to file called 'VFdiffDump.put' Write index key to 'VFdiffDump.put' WRITE(9,*) 'Differential View Factors' WRITE (9,*) 'Run Started:',dtinit WRITE(9,701) 'Input File Name:',filein WRITE(9,*) WRITE(9,*) WRITE(9,*) 'Index Key:' WRITE(9,*) WRITE(9,*) 'iFROM: Originating whole surface [Range: 1 6]' WRITE(9,*) n1: Leading alphabetical index of differential &surface on iFROM [1 mxsfc]' WRITE(9,*) (i.e. n1=y index if iFROM is in y z plane)' WRITE(9,*) n2: Trailing alphabetical index of di fferential &surface on iFROM [1 mxsfc]' WRITE(9,*) (i.e. n2=z index if iFROM is in y z plane)' WRITE(9,*) iTO: Receiving whole surface [1 6]' WRITE(9,*) m1: Leading alphabetical index of differential &surf ace on iTO [1 mxsfc]' WRITE(9,*) m2: Trailing alphabetical index of differential &surface on iTO [1 mxsfc]' 70 FORMAT(A7,3X,3(I2,A1),X,3(I2,A1),A1,D20.10) 701 FORMAT(X,A16,3X,A20) 702 FORMAT(A29) DO iFROM=1,nsfc DO iTO=1,nsfc WRITE(9,*) WRITE(9,702) 'VFdiff(iFROM,n1,n2;iTO,m1,m2)' DO nc=1, mxsfc DO nf=1,mxsfc DO mc=1, mxsfc DO mt=1, mxsfc WRITE(9,70) 'VFDiff(',iFROM,',',nc,' ,',nf,';',iTO,',',mc,',' & ,mt,')','=',VFDiff(iFROM,nc,nf,iTO,mc,mt) ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO Calculate Whole Surface View Factors PAGE 103 103 Initialize whole surface VF array to zero VF=0 Common Edge View Factor CALL COMMON_EDGE Opposide Ends View Factor CALL OPP_SIDES Self seeing VFs remain zero DO j=1,3,2 VF(j+1,j+2)=VF(j,j+2) END DO VF(1,4)=VF(1,3) rest=1. VF(1,2) VF(1,3) VF(1,4) VF(1,5)=rest/2. VF(1,6)=VF(1,5) Sets opposite end VFs equal (symmetry) DO j=3,totsfc VF(2,j)=VF(1,j) END DO Recipricity on entire VF matrix 40 DO k =1,totsfc DO i=1,totsfc VF(i,k)=VF(k,i)*Area(k)/Area(i) END DO END DO Summation & Symmetry to get VF(#,5) & Vf(#,6) (This section could be condensed with loop development) rest3=1 VF(3,1) VF(3,2) VF(3,4) VF(3,5)=rest3/2. VF(3,6)=VF(3,5) IF(VF(5,3) .EQ. 0) THEN GOTO 40 END IF rest4=1 VF(4,1) VF(4,2) VF(4,3) VF(4,5)=rest4/2.D0 PAGE 104 104 VF(4,6)=VF(4,5) IF(VF(5,4) .EQ. 0) THEN GOTO 40 END IF rest5=1 VF(5,1) VF(5,2) VF(5,3) VF(5,4) VF(5,6)=rest5 IF(VF(6,5) .EQ. 0) THEN GOTO 40 END IF WRITE(*,*) WRITE(*,*) 'Surface Areas [m^2]=' DO i=1, totsfc WRITE(*,61) 'Sfc(',i,')=',Area(i) END DO 61 FORMAT(A4,I1,A2,100F12.2) Summation check leaving each surface (Verify that sum[VFt(i)]=1) WRITE(*,*) WRITE(*,*) 'Verification that View Factors leaving a surface sum & to 1:' DO i=1,totsfc VFt(i)=0 DO j=1,totsfc VFt(i)=VFt(i)+VF(i,j) END DO WRITE(*,58) 'VFt(',i,')= ',VFt(i) END DO 58 FORMAT(A4,I1,A3,F16.14) Print to CMD Window & Output file WRITE(1,*) 'Whole Surface Diffuse View Factor Matrix:' WRITE(1,*) '[ VF(row,column) ]' WRITE(1,53) isfcs WRITE(*,*) WRITE(*,*) 'Whole Surface Diffuse View Factor Matrix:' WRITE(*,*) '[ VF(row,column) ]' WRITE(*,60) isfcs 53 FORMAT(6(I13)) 60 FORMAT(6(I10)) 59 FORMAT (I5, 100F10.5) PAGE 105 105 DO i=1,totsfc WRITE(*,59) isfcs(i), (VF(i,j), j=1,totsfc) WRITE(1,51) isfcs(i), (VF(i,j), j=1,totsfc) END DO WRITE(*,*) WRITE(1,*) WRITE(1,*) 'Differential View Factor Summary Data:' WRITE(1,*) ndiffsfc, 'differential surfaces defined' WRITE(1,*) iVFdiff, 'differential VFs calculated' WRITE(1,*) iNegVFdiff, 'negative differential VFs calculated' 56 FORMAT(I4,3X,16G11.3) END SUBROUTINE VIEW_FACTORS SUBROUTINE HEAT_TRANSFER USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy nsfc=totsfc sig = 5.6704D 8 !Stefan Boltzmann Constant [W/(m^2K^4)] Calculate Thermal Diffusivity UNITS:: rk: [W/mK], rho: [kg/m^3], Cp: [J/kgK] ==> alphaTH: [m^2/s] DO i=1,nsfc alphaTH(i) = rk(i)/(rho(i)*Cp(i)) WRI TE(*,*) 'alphaTH', i, '=', alphaTH(i) ENDDO Change eps to 1.0 for "invisible surface" condition DO iTO=1,nsfc IF(jVis(iTO).EQ.0)THEN eps(iTO)=1.D0 ENDIF ENDDO Write Heat Tranfser properties from input file to output file 62 FORMAT(I5, 6X,I1,1X, 4(F10.4, 5X),E15.5) 63 FORMAT(A3,7X,A4,4X,A3,8X,A9,6X,A11,7X,A8,5X,A14) 620 FORMAT(I5, 6X,I1,1X, F10.4,9X,A3,12X,A3,14X,A3,12X,A3) PAGE 106 106 WRITE(1,*) WRITE(1,*) 'Heat Transfer Propertie s:' WRITE(1,63) 'Sfc','jVis','eps','Cp[J/kgK]','rho[kg/m^3]', & 'rk[W/mK]', & 'alphaTH[m^2/s]' DO i=1, nsfc IF(jVis(i).EQ.0)THEN WRITE(1,620) isfcs(i),jVis(i),eps(i),'N/A', & 'N/A','N/A','N/A' ELSE WRITE(1,62) isfcs(i),jVis(i),eps(i),Cp(i),rho(i),rk(i), & alphaTH(i) ENDIF ENDDO WRITE(1,*) Initialize Surface Temps from input file. DO i=1,ndiffsfc CALL MAT_INDEX(i,iFROM,n,k) Tdiff(iFROM,n,k)=T(iFROM) ENDDO Define coefficient matrices for flux calculation (using Cij & Aij matrices from Eqn 5.38, MODEST Pg. 175) DO i= 1,ndiffsfc CALL MAT_INDEX(i,iFROM,n,k) DO j = 1 ndiffsfc CALL MAT_INDEX (j,iTO,l,m) xt = dist(iFROM,n,k,iTO,l,m) IF((xt.LE.1000).OR.(iVis.EQ.999))THEN Tau = 1.D0 ELSE CALL ATTENUATE(xt,Tau) ENDIF dij = CronDel(i,j) Cij(i,j)=(dij/eps(iTO)) ((1.D0/eps(iTO)) 1) & *Tau*VFdiff(iFROM,n,k,iTO,l,m) Aij(i,j)=dij Tau*VFdiff(iFROM,n,k,iTO,l,m) ENDDO ENDDO Write Cij matrix to 'Cij.out' file OPEN(12, file='Cij.out', A CCESS='SEQUENTIAL') WRITE(12,*) 'Cij Coefficient Matrix' WRITE (12,*) 'Run Started:',dtinit WRITE(12,623) 'Input File Name:',filein PAGE 107 107 WRITE(12,*) WRITE(12,*) 'i & j : Differential surface number (single index)' WRITE(12,*) [Range: 1 ndiffsfc]' WRITE(12,*) ndiffsfc = 6*mxsfc*mxsfc' WRITE(12,*) WRITE(12,622) 'i \ j',(j, j=1,ndiffsfc) WRITE(12,*) DO i = 1,ndiffsfc WRITE(12,621) i,(Cij(i,j), j=1,nd iffsfc) ENDDO 621 FORMAT(I5,3X,9999(ES15.5,3X)) 622 FORMAT(4X,A3,6X,9999(I5,13X)) 623 FORMAT(X,A16,3X,A20) Write Aij matrix to 'Aij.out' file OPEN(13, file='Aij.out', ACCESS='SEQUENTIAL') WRITE(13,*) 'Aij Coefficient Matrix' WRITE (13,*) 'Run Started:',dtinit WRITE(13,623) 'Input File Name:',filein WRITE(13,*) WRITE(13,*) 'i & j : Differential surface number (single index)' WRITE(13,*) [Range: 1 ndiffsfc]' WRITE(13,*) nd iffsfc = 6*mxsfc*mxsfc' WRITE(13,*) WRITE(13,622) 'i \ j',(j, j=1,ndiffsfc) WRITE(13,*) DO i = 1,ndiffsfc WRITE(13,621) i,(Aij(i,j), j=1,ndiffsfc) ENDDO Guess Initial Flux & Fluence profiles (for iterative method) qdotNETdO = 1 qdotNETd = qdotNETdO FluenceIN = 0 FluenceNET = 0 Invert Cij Matrix (if direct inversion is selected) IF(iFluxChoice.EQ.1)THEN CALL MAT_INVERSION ENDIF Calculate Pmax [W] (maxi mum thermal power) Pmax = 1.33D13*(Yield)**0.56 PAGE 108 108 Calculate radius at thermal maximum [m] Rmax = 50*(Yield)**(1.D0/3.D0) Calculate fireball surface area at thermal maximum [m^2] Amax = 4*pi*Rmax**2 Calculate tmax [s] (t ime of maximum thermal power) tmax = 0.0417D0*(Yield)**0.44 Calculate normalized time step deltaTimeNorm = deltaTime/tmax Calculate the number of normalized time steps numsteps = Temptime(numRows,1)/deltaTimeNorm numsteps = tnormmax/deltaTimeNorm WRITE(*,*) 'Temptime(numRows,1) =', Temptime(numRows,1) WRITE(*,*) 'tmax = ', tmax WRITE(*,*) 'deltaTime = ', deltaTime WRITE(*,*) 'numsteps = ', numsteps Write Sky Factor to output f ile WRITE(1,*) WRITE(1,*) 'Sky Factor:' WRITE(1,668) skyfact WRITE(1,*) 668 FORMAT(F6.3) Write Simulation Parameters to output file 65 FORMAT(A9,4X,A8,4X,A8,5X,A3,9X,A13,5X,A12) 66 FORMAT(F5.2,10X,F5.2,3X,ES10.3,2X,ES10.3,8X,I6,6X,I15) WRITE(1,*) 'Simulation Parameters:' WRITE(1,65) 'Yield[kT]','tnormmax','dtime[s]','tol', &'# Print Steps','# Time Steps' WRITE(1,66) Yield, tnormmax, deltaTime, tol, numout, numsteps WRITE(1,*) Write Atmospheric Conditions to output file WRITE(1,*) 'Atmospheric Conditions:' WRITE(1,670) 'iVis[km]','TempFactor', 'Temp Source File Name:' IF(iTsource.EQ.1)THEN WRITE(1,669) iVis,tfact,'none' ELSE WRITE(1,67) iVis, tfact,tfilein ENDIF PAGE 109 109 WRITE(1,*) 67 FORMAT(3X,I3,11X,F7.3,18X,A20) 670 FORMAT(A8,9X,A10,3X,A25) 669 FORMAT(3X,I3,11X,F7.3,18X,A13) WRITE(*,38) 'ntemp','iters','timenorm','tempint','qdotIN2avg', &'Flu enceIN2avg' timenorm = 0 Open file for flux profile storage called 'FluxPf.out' OPEN(6, file='FluxPf.out', ACCESS='SEQUENTIAL') WRITE(1,*) &' WRITE(6,*) 'FLUX RESULTS: Surface 2 Flux & Fluence Profiles' WRITE(1,*) 'SUMMARY RESULTS: Surface 2 Average Flux & Fluence' WRITE(6,*) Write surface 2 flux map locations to 'FluxPF.out'. WRITE(6,*) 'Flux Map Locations:' WRITE(6,37) 'Sfc#','iTO','l','m','Xloc[m]','Yloc[m]', & 'Zloc[m]' iTO = 2 DO m = 1,mxsfc DO l = 1,mxsfc CALL MAT_INDEX_REV(iTO,l,m,j) WRITE(6,45) j, iTO,l, m, Xloc(j), Yloc(j), Zloc(j) ENDDO ENDDO WRITE(6,*) WRITE(6,671) 'PrintStep', 'Time','Timenorm','iterations', &'DiffSfc#','FluxIN','FluenceIN' WRITE(6,672) '[s]','[cal/(c m^2*s)]','[cal/(cm^2)]' WRITE(6,*) 671 FORMAT(A9,4X,A4,5X,A8,3X,A10,2X,A8,4X,A6,13X,A9) 672 FORMAT(14X,A3,40X,A14,5X,A12) Write surface 2 average flux HEADER to 'CANout.put' file WRITE(1,*) WRITE(1,32) 'PrintStep', 'Time', 'T imenorm','FBTemp', &'FluxINavg','FluenceAvg','FluxRow1INavg' WRITE(1,321) '[s]','[K]','[cal/(cm^2*s)]','[W/cm^2]', PAGE 110 110 &'[cal/cm^2]','[cal/(cm^2*s)]','[W/cm^2]' 32 FORMAT(A9,3X,A4,6X,A8,4X,A6,17X,A9,12X,A10,12X,A13) 321 FORMAT(12X,A3,19X ,A3,10X,A14,2X,A8,7X,A10,5X,A14,2X,A8) *******TIME STEPS******* Step through calculation in Time & write results to output files ntemp = 1 nprint = numsteps/numout DO i=1,numsteps timenorm = timenorm + deltaTimeNorm Time = timenorm*tmax icount = 0 Call time step subroutine CALL TEMP_MARCH(timenorm, tempint) Utilize Print Step Parameter IF ((i.GT.nprint).AND.(ntemp.LE.numout)) THEN iTO = 2 DO m = 1,mxsfc DO l = 1,mxsfc CALL MAT_INDEX_REV(iTO,l,m,j) Write flux & profile to 'FluxPF.out' after converting to [cal/(cm^2*s)] & [cal/cm^2*s]. IF((m.EQ.1).AND.(l.EQ.1))THEN WRITE(6,673) ntemp,Time,timen orm,icount,j, & qdotINdiff(iTO,l,m)*2.39D 5,FluenceIN(iTO,l,m)*2.39D 5 ELSE WRITE(6,674) j,qdotINdiff(iTO,l,m)*2.39D 5, & FluenceIN(iTO,l,m)*2.39D 5 ENDIF ENDDO ENDDO Write su rface 2 average fluxes to 'CANout.put' file WRITE(1,31) ntemp,Time,timenorm,tempint,Sfc2Ave,Sfc2Ave*4.184 & ,Flue2Ave, FluxRow1INavg, FluxRow1INavg*4.184 Write surface 2 avg flux and fluence to comand screen WRITE(*,40) ntemp,icount,timenorm,tempint,Sfc2Ave,Flue2Ave Write selected temperatures to 'temps.put' file [For development: Surface 3 (left wall), bottom row] WRITE(11,676,advance='no') ntemp,Time,timenorm,tempint DO iTP=1,numtemppoints PAGE 111 111 iTO=TPsfc(iTP) l=n1TP(iTP) m=n2TP(iTP) WRITE(11,675,advance='no') Tdiff(iTO,l,m) ENDDO WRITE(11,*) Update print counters ntemp = ntemp + 1 nprint = nprint + nu msteps/numout ENDIF IF(iTrans.EQ.0)THEN EXIT ENDIF ENDDO 676 FORMAT(I8,3X,F7.5,3X,F9.5,3X,F9.2,5X) 675 FORMAT(F9.2,3X) 673 FORMAT(I8,3X,F7.5,3X,F9.5,7X,I2,9X,I4,6X,ES14.7,5X,ES14.7) 674 FORMAT(48X,I4,6X,ES14.7,5 X,ES14.7) 33 FORMAT(A39) 31 FORMAT(I8,3X,F8.4,3X,F8.4,3X,F11.4,3X,ES13.7,3X,ES13.7,2X, & ES13.7,2X,ES13.7,3X,ES13.7) WRITE(*,38)'ntemp','iters','timenorm','tempint','qdotIN2avg', &'FluenceIN2avg' 38 FORMAT(2(A8,2X),1X,A8,4X,A8,1X,2(A18,2X)) 39 FORMAT(2(A8,2X),X,A8,5X,3(A18,2X)) 40 FORMAT(I9,I8,E14.5,3X,F8.2,3X,2(F15.5,3X)) 4 FORMAT(I9,I8,E15.5,5X,3(F15.5,5X)) 45 FORMAT(I4,3X,I1,2X,I2,2X,I2,2X,3(F10.3,1X)) 37 FORMAT(A4,2X,A3,2X,A1,3X,A 1,4X,3(A8,3X),2X,A8,5X,A20) 36 FORMAT(A10,X,I4,3X,A9,X,F8.4,4X,A10,X,E10.4,4X,A11,2X,I3) 20 FORMAT(A4,6X,A7,5X,A20,3X,A19,15X,A28) 21 FORMAT(I4,3X,F12.5,2X,F14.7,11X,F14.7) Benchmark output (Cengel pg. 726) IF((iTrans.EQ.0).OR.(iTRA NS.EQ.2))THEN Sfc2Sum = 0 Sfc2Ave = 0 iTO = 2 DO l=1,mxsfc PAGE 112 112 DO m=1,mxsfc Sfc2Sum = Sfc2Sum + qdotNETdiff(iTO,l,m) ENDDO ENDDO Sfc2Ave = Sfc2Sum/(mxsfc*mxsfc) Qdotnet = Sfc2Ave*area(iTO) WRITE(*,*) 'Net Avg Stdy Flux on Sfc 2 = ',Sfc2Ave, '[W/m^2]' WRITE(*,*) 'Net Stdy Energy Trans Rate on Sfc2 = ',Qdotnet,'[W]' ENDIF IF(iTrans.EQ.2)THEN FluNT2Sum = 0 FluNT2Ave = 0 iTO = 2 DO l=1,mxsfc DO m=1,mxs fc FluNT2Sum = FluNT2Sum + FluenceNET(iTO,l,m) ENDDO ENDDO FluNT2Ave = FluNT2Sum/(mxsfc*mxsfc) Qnet = FluNT2Ave*area(iTO) WRITE(*,*) 'Net Average Fluence on Sfc 2 = ',FluNT2Ave,'[J/m^2]' WRITE(*,*) 'Net Energy on Sfc 2 = ', Qnet, '[J]' ENDIF END SUBROUTINE HEAT_TRANSFER SUBROUTINE TEMP_MARCH(timenorm, tempint) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (A H,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy nsfc=totsfc sig = 5.6704D 8 !Stefan Boltzmann Constant [W/(m^2K^4)] Update fireball temperature from ENW data or embedded BB temps IF(iTrans.EQ.1)THEN IF(iTsource.EQ.1)THEN CALL BB_TEMPS(timenorm,tempint) ELSE CALL TIME_TEMP(timenorm,tempint) ENDIF PAGE 113 113 DO j=1,mxsfc DO k=1,mxsfc Tdiff(1,j,k)=tempint ENDDO ENDDO ENDIF Calculate black body emissive power for every differential area at 'old' temp UNITS:: sig: [W/(m^2K^4)], Tdiff: [K] ==> emissive power: [W/m^2] DO iFROM=1,totsfc DO j= 1,mxsf c DO k= 1,mxsfc IF(iFROM.EQ.6)THEN emitpwr(iFROM,j,k)=skyfact*sig*Tdiff(iFROM,j,k)**4 ELSE emitpwr(iFROM,j,k)=sig*Tdiff(iFROM,j,k)**4 ENDIF ENDDO ENDDO ENDDO Calculate steady state flux at each subsurface based on selected method. IF(iFluxChoice.EQ.1)THEN CALL MAT_HEAT_FLUX icount = 0 ELSE CALL FLUXITERATE ENDIF Semi infinite Solid tem perature update method [Cengel pg. 241 eq. 446] IF(iTrans.EQ.1)THEN DO i=1,ndiffsfc CALL MAT_INDEX(i,iTO,l,m) TdiffNEW(iTO,l,m) = Tdiff(iTO,l,m) + ( qdotNETdO(i,1)/rk(iTO)) & *DSQRT(4*alphaTH(iTO)*deltaTime/PI) ENDDO Store new temperatures as old temperatures Tdiff = TdiffNEW ENDIF "Invisible Surface" boundary condition (T maintained at initial PAGE 114 114 temp from input file) DO i=1,ndiffsfc CALL MAT_INDEX(i,iFROM,n,k) IF(jVis(iFROM).EQ.0)THEN Tdiff(iFROM,n,k)=T(iFROM) ENDIF ENDDO Calculate incident flux (net + emitted) & fluence DO iTO=1,nsfc DO l=1,mxsfc DO m=1,mxsfc CALL MAT_INDE X_REV(iTO,l,m,i) qdotNETdiff(iTO,l,m) = qdotNETdO(i,1) qdotINdiff(iTO,l,m) =qdotNETdiff(iTO,l,m)+emitpwr(iTO,l,m) FluenceIN(iTO,l,m)=FluenceIN(iTO,l,m) & +qdotINdiff(iTO,l,m)*deltaTime FluenceNET(iTO,l,m)=FluenceNET(iTO,l,m) & +qdotNETdiff(iTO,l,m)*deltaTime ENDDO ENDDO ENDDO Calculate Surface 2 average flux [cal/(cm^2*s)] & fluence [cal/cm^2] Sfc2Ave = 0 Sfc2Sum = 0 Flu e2Ave = 0 Flue2Sum = 0 iTO=2 DO l=1,mxsfc DO m=1,mxsfc Convert W/m^2 to cal/(cm^2*s) Sfc2Sum = Sfc2Sum + qdotINdiff(iTO,l,m)*2.39D 5 Flue2Sum = Flue2sum + FluenceIN(iTO,l,m)*2.39D 5 ENDDO ENDDO Sfc2Ave = Sfc2Sum/(mxsfc*mxsfc) Flue2Ave = Flue2Sum/(mxsfc*mxsfc) Calculate bottom row average flux FluxRow1INavg = 0 FluxRow1Sum = 0 iTO=2 m=1 DO l=1,mxsfc PAGE 115 115 FluxRow1 Sum = FluxRow1Sum + qdotINdiff(iTO,l,m)*2.39D 5 ENDDO FluxRow1INavg = FluxRow1Sum/(mxsfc) 35 FORMAT(A7,D15.5) 34 FORMAT(A5,I1,A4,D15.5) END SUBROUTINE TEMP_MARCH SUBROUTINE TIME_TEMP(timenorm,tempint) This subroutine reads the timetemperature fireball data from Glasstone and prepares it for scaling to any yield. USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) DO i=1,numRows IF((timenorm .GE. Temptime(1,1)).AND. & (timenorm.LT.Temptime(numRows,1)))THEN IF((timenorm .GE. Temptime(i,1)).AND. & (timenorm.LT.Temptime(i+1,1))) THEN timelow = Temptime(i,1) timehi = Temptime(i+1,1) templow = Temptime(i,2) temphi = Temptime(i+1,2) tempint=temphi ((timehi timenorm)*(temphi templow)/ & (timehi timelow)) ENDIF ELSEIF(timenorm.LT.Temptime(1,1))THEN tempint = Temptime(1,2) ELSEIF(timenorm.GT.Temptime(numRows,1))THEN tempint= Temptime(numRows,2)+(timenorm Temptime(numRows,1)) & *((Temptime(numRows 1,2) Temptime(numRows,2) ) & /(Temptime(numRows 1,1) Temptime(numRows,1))) ENDIF ENDDO END SUBROUTINE TIME_TEMP SUBROUTINE BB_TEMPS(timenorm,tempint) This subroutine calculates the black body fireball temperature based on a constant radius at thermal maximum and the normalized power curve (ENW). PAGE 116 116 USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER ( I N) timenorm=1 Pt = Pmax*((2*timenorm**2)/(1+timenorm**4)) tempint = (Pt/(Amax*sig))**0.25 tempint = tempint*tfact END SUBROUTINE BB_TEMPS SUBROUTINE MAT_ HEAT_FLUX This subroutine uses matrix methods to calculate steady state heat fluxes within a time step using Eqn.5.38 from Modest Pg.175. VARIABLES: qdotNETdO: Singledimension array of heat flux on each differential area CijINV: Inverse of Cij coefficient matrix (VFs, eps,& dij's) Aij: Coefficient matrix with dij's & VFs USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy Recast 3 d 'emitpwr' into 1d 'empwr' DO i=1,ndiffsfc CALL MAT_INDEX(i,iFROM,j,k) empwr(i,1) = emitpwr(iFROM,j,k) ENDDO qtemp = MATMUL(Aij,empwr) qdotNETdO = MATMUL(CijINV, qtemp) END SUBROUTINE MAT_HEAT_FLUX DOUBLE PRECISION FUNCTION CronDel(i,j) The Kronecker Delta Function IF (i .EQ. j) THEN CronDel = 1.D0 ELSE PAGE 117 117 Cr onDel = 0.D0 ENDIF END SUBROUTINE MAT_INDEX_REV(n1,n2,n3,i) PURPOSE: To shift 3 dimensional differential surfaces references to a single dimension surface number. (iFROM,n,k) > i (1,4,1) > 4 VARIABLES USED: i : single dimension surface number n1: Surface number (16) (i.e. iFROM or iTO) n2: 1st subsurface index (alphabetically 1st dim of x,y, or z) n3: 2nd subsurface incex (alphabetically 2nd dim of x,y, or z) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy i = ((n11)*mxsfc*mxs fc)+n2+(n31)*mxsfc END SUBROUTINE MAT_INDEX_REV SUBROUTINE MAT_INDEX(i,n1,n2,n3) PURPOSE: To shift single dimension surface number inputs back to a 3d subsurface index reference i > (iFROM,n,k) 4 > (1,4,1) VARIABLES USED: i : single dimension surface number n1: Surface number (16) (i.e. iFROM or iTO) n2: 1st subsurface index (alphabetically 1st dim of x,y, or z) n3: 2nd subsurface incex (alphabetically 2nd dim of x,y, or z) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy PAGE 118 118 mxsfc2=mxsfc*mxsfc IF(i.LE. mxsfc2)THEN n1 = 1 ELSEIF((i.GT.mxsfc2).AND.(i.LE.2*mxsfc2))THEN n1 = 2 ELSEIF((i.GT.2*mxsfc2).AND.(i.LE.3*mxsfc2))THEN n1 = 3 ELSEIF((i.GT.3*mxsfc2).AND.(i.LE.4*mxsfc2))THEN n1 = 4 ELSEIF((i.GT.4*mxsfc2).AND.(i.LE.5*mxsfc2))THEN n1 = 5 ELSEIF((i.GT.5*mxsfc2).AND.(i.LE.6*mxsfc2))THEN n1 = 6 ENDIF n2 = MOD(i,mxsfc) IF(n2.EQ.0)THEN n2 = mxsfc n3 = (i/mxsfc) (n1 1)*mxsfc ELSE n3 = (i/mxsfc) + 1 (n11)*mxsfc ENDIF END SUBROUTINE MAT_INDEX SUBROUTINE MAT_INVERSION PURPOSE: Di rect inversion method on Cij coefficient matrix via Gauss Jordan elimination USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy Size of matrix to be inverted n=ndiffsfc m= 2*ndiffsfc DO I=1,n DO J=1,m IF (J .GT. n) THEN IF (J .EQ. n+I) THEN Cm(I,J)=1.D0 PAGE 119 119 ELSE Cm(I,J)=0.D0 ENDIF ELSE Cm(I,J) = Cij(I,J) ENDIF ENDDO ENDDO D is a dummy matrix for Row Operations purposes Puts 1's on diagonal and 0'1 below it by sorting a column, then dividing each row by its leading entry ( if non zero), subtracting the "top" row, then stepping down the diagonal DO iP=1,n 2 DO I=iP,n 1 DO K=I+1,n IF (Cm(I,iP) .LT. Cm(K,iP)) THEN Dm(I,1:m) = Cm(I,1:m) Cm(I,1:m) = Cm(K,1:m) Cm(K,1:m) = Dm(I,1:m) ENDIF ENDDO ENDDO DO I=iP,n IF (Cm(I,iP) .NE. 0.D0) THEN Cm(I,1:m) = Cm(I,1:m)/Cm(I,iP) ENDIF ENDDO DO I=iP+1,n IF (Cm(I,iP) .NE. 0.D0) THEN Cm(I,1:m) = 1.D0*Cm(iP,1:m)+Cm(I,1:m) ENDIF ENDDO ENDDO IF (Cm(n 1,n 1) .LT. Cm(n,n1)) THEN Dm(1,1:m) = Cm(n 1,1:m) Cm(n 1 ,1:m) = Cm(n,1:m) Cm(n,1:m) = Dm(1,1:m) ENDIF DO I=n 1,n PAGE 120 120 IF (Cm(I,n 1) .NE. 0.D0) THEN Cm(I,1:m) = Cm(I,1:m)/Cm(I,n 1) ENDIF ENDDO Cm(n,1:m) = 1.D0*Cm(n1,1:m)+Cm(n,1:m) IF (Cm(n,n) .NE. 0) THEN Cm(n,1:m) = Cm(n,1:m)/Cm(n,n) ENDIF Gets 0's above diagonal DO J=1,n 1 DO I=1,n J IF (Cm(I,I+J) .NE. 0.D0) THEN Cm(I,1:m)=Cm(I,1:m)/Cm(I,I+J) Cm(I+J,1:m) ENDIF IF (Cm(I,I) .NE. 0.D0) THEN Cm(I,1:m)=Cm(I,1:m)/Cm(I,I) ENDIF ENDDO ENDDO 40 FORMAT (100F8.3) DO I=1,n DO J=1,n CijINV(I,J)=Cm(I,J+n) ENDDO ENDDO Write C ijINV matrix to 'CijINV.out' file OPEN(5, file='CijINV.out', ACCESS= 'SEQUENTIAL') WRITE(5,*) 'CijINV Coefficient Matrix' WRITE (5,*) 'Run Started:',dtinit WRITE(5,890) 'Input File Name:',filein WRITE(5,*) WRITE(5,*) '[D irect inversion selected for this run]' WRITE(5,*) WRITE(5,*) 'i & j : Differential surface number (single index)' WRITE(5,*) [Range: 1 ndiffsfc]' WRITE(5,*) ndiffsfc = 6*mxsfc*mxsfc' WRITE(5,*) WRITE(5,889) 'i \ j',(j, j=1,ndiffsfc) WRITE(5,*) DO i=1,ndiffsfc PAGE 121 121 WRITE(5,888) i,(CijINV(i,j), j=1,ndiffsfc) ENDDO 888 FORMAT(I5,3X,9999(ES15.5,3X)) 889 FORMAT(4X,A3,6X,9999(I5,13X)) 890 FORMAT(X,A16,3X,A20) END SUBROUTINE MAT_INVERSION SUBROUTINE FLUXITERATE PURPOSE: This subroutine uses Gauss Seidel iteration to resolve the flux in each timestep because there app ear to be numerical issues useing direct inversion methods. USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy Initialize blackbody emissive power DO i=1,ndiffsfc CALL MAT_INDEX(i,iFROM,j,k) empwr(i,1) = emitpwr(iFROM,j,k) ENDDO Set error arrays to zero errorOLD = 0 error = 0 500 DO i=1,ndiffsfc rnewsum = 0 roldsum = 0 EmitTrm = 0 IF (i.EQ.1) THEN rnewsum = 0 ELSE DO j=1,i 1 rnewsum = rnewsum + Cij(i,j)*qdotNETd(j,1) ENDDO ENDIF DO j=i+1,ndiffsfc roldsum = roldsum + Cij(i,j)*qdotN ETdO(j,1) ENDDO PAGE 122 122 DO j=1,ndiffsfc EmitTrm = EmitTrm + Aij(i,j)*empwr(j,1) ENDDO qdotNETd(i,1)=(1/Cij(i,i))*(EmitTrm rnewsum roldsum) ENDDO Calculate error arrays DO i=1,ndiffsfc error(i) = ABS((qdotNETd(i,1)qdotNETdO(i,1))/qdotNETdO(i,1)) ENDDO P=MAXVAL(error) Check maximum error against tolerance from input file IF (P.GE.tol) THEN Store "new" flux as "old" flux qdotNETdO=qdotNETd icount = icount + 1 GOTO 500 ENDIF END SUBROUTINE FLUXITERATE SUBROUTINE DISTANCE PURPOSE: This subroutine calculates the distance between each pair of differential surfaces for use in defining exponential atmospheric attenuation based on the transmittance of air. Distance formula: d = SQRT(dz**2 + dy**2 + dz**2) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy DO iFROM=1,nsfc DO n=1,mxsfc DO k=1,mxsfc CALL COORDDEF(iFROM,n,k,a1,b1,c1) DO iTO=1,nsfc DO l=1,mxsfc DO m=1,mxsfc CALL COORDDEF(iTO,l,m,a2,b2,c2) PAGE 123 123 dist(iFROM,n,k,iTO,l,m) & = DSQRT((a1 a2)**2+(b1b2)**2+(c1 c2)**2) ENDDO ENDDO ENDDO ENDDO END DO ENDDO END SUBROUTINE DISTANCE SUBROUTINE LOCATE Purpose: To define an aray that locates the center of each differential surface in terms of x y z distances in [m]. USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy DO i = 1,ndiffsfc CALL MAT_INDEX(i,iTO,l,m) CALL COORDDEF(iTO,l,m,a1,b1,c1) IF(iTO.EQ.1)THEN Xloc(i) = a1 Zloc(i) = b1 Yloc(i) = 0 ELSEIF(iTO.EQ.2)THEN Xloc(i) = a1 Zloc(i) = b1 Yloc(i) = dy(3) ELSEIF(iTO.E Q.3)THEN Yloc(i) = a1 Zloc(i) = b1 Xloc(i) = 0 ELSEIF(iTO.EQ.4)THEN Yloc(i) = a1 Zloc(i) = b1 Xloc(i) = dx(5) ELSEIF(iTO.EQ.5)THEN Xloc(i) = a1 Yloc(i) = b1 Zloc(i) = 0 ELSEIF(iTO.EQ.6)THEN PAGE 124 124 Xloc(i) = a1 Yloc(i) = b1 Zloc(i) = dz(1) ENDIF ENDDO END SUBROUTINE LOCATE SUBROUTINE COORDDEF(nsurf,n1,n2,a,b,c) PURPOSE: To properly dimension the components of the distance formula based on the iFROM iTO surface pairs. USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy IF((nsurf.EQ.1).OR.(nsurf.EQ.2))THEN a = (n10.5)*ddx(nsurf) b = (n20.5)*ddz(nsurf) c = 0 ELSEIF((nsurf.EQ.3).OR.(nsurf.EQ.4))THEN a = (n10.5)*ddy(nsurf) b = (n20.5)*ddz(nsurf) c = 0 ELSEIF((nsurf.EQ.5).OR.(nsurf.EQ.6))THEN a = (n10.5)*ddx(nsurf) b = (n20.5)*ddy(nsurf) c = 0 ENDIF END SUBROUTINE COORDDEF SUBROUTINE ATTENUATE(xt,Tau) PURPOSE: To calculate the transmittance to each differential surface based on the visibility designated in the input file. Based on J.M. Leonard AFIT Master's Thesis. USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACT ER*180 dummy PAGE 125 125 Coefficients express an equation of the form: Tau = C1 + C2(ln(x)) + C3(ln(x))^2 + C4(ln(x))^3 + C5(ln(x))^4 + C6(ln(x))^5 Where: Tau = Planckianweighted transmittance x = Ground distance (km) b etween receive & target (1 km < x PAGE 126 126 SUBROUTINE READ_INP USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER*180 dummy DO i=1,3 READ(8,'(A)') dummy END DO READ(8,*) mxsfc totsfc hardwired to 6 totsfc=6 DO i=1,3 READ(8,'(A)') dummy END DO ndiffsfc = 2*(mxsfc*mxsfc + mxsfc*mxsfc + mxsfc*mxsfc) CALL ALLOCARR nsfc=to tsfc DO isfc=1, nsfc READ(8,*) isfcs(isfc),dx(isfc),dy(isfc),dz(isfc), & Description(isfc) END DO DO i=1,3 READ(8,'(A)') dummy END DO DO i=1, nsfc READ(8,*), isfcs(i), T(i), dummy END DO DO i=1,4 READ(8,'(A)') dummy END DO DO i=1, nsfc READ(8,*), isfcs(i), jVis(i), eps(i), Cp(i), rho(i), rk(i) END DO DO i=1,2 READ(8,'(A)') dummy PAGE 127 127 END DO Read sky factor from input file READ(8,*) skyfact DO i=1,3 READ(8,'(A)') dummy END DO READ(8,*) Yield, tnormmax, deltaTime, numout, tol, iFluxChoice IF((iFluxChoice.LT.1).OR.(iFluxChoice.GT.2))THEN WRITE(*,*) 'ERROR: Invalid Flux Method Choice. Must be 1 or 2' WRITE(*,*) '1=Direct Matrix Inversion' WRITE(*,*) '2=Gause Seidel Iteration' iERROR = 1 ENDIF DO i=1,2 READ(8,'(A)') dummys END DO READ(8,*) iTrans IF((iTrans.LT.0).OR.(iTrans.GT.2))THEN WRITE(*,*) 'ERROR: Invalid iTrans Choice. Must be 1 or 0' WRITE(*,*) '1=Transient Problem' WRITE(*,*) '0=Steady State Problem' iERROR = 1 ENDIF DO i=1,3 READ(8,'(A)') dummys END DO Read SeaLevel Visibility from input file READ(8,*) iVis, tfact, iTsource IF((iVis.NE.999).AND.(iVis.NE.5).AND.(iVis.NE.10).AND. & (iVis.NE.15).AND.(iVis.NE.23))THEN WRITE(*,*) 'ERROR: Invalid visibility choice in input file. &Must be 5, 10, 15, 23, or 999 km' iERROR = 1 ENDIF Read fireball temperature data file PAGE 128 128 IF(iTsource.NE.1)THEN tfilein = 'tempdata.dat' Open t emperature data file OPEN(10,file=tfilein, status='old') READ(10,*) dummy READ(10,*) numRows ALLOCATE(TempTime(numRows,2)) DO i=1,numRows READ(10,*) dtime, dtemp TempTime(i,1)= dtime Te mpTime(i,2)= dtemp*tfact ENDDO ENDIF DO i=1,3 READ(8,'(A)') dummy ENDDO Read number of temperature requests & allocate required arrays READ(8,*) numtemppoints ALLOCATE(TPsfc(numtemppoints),TPxloc(numtemppoints), &TPyloc(numtemppoints),TPzloc(numtemppoints), &TPtemp(numtemppoints),n1TP(numtemppoints),n2TP(numtemppoints)) DO i=1,2 READ(8,'(A)') dummy ENDDO Read each temperature request DO i=1,numtemppoints READ(8,*) dummy,TPsfc(i),TPxloc(i),TPyloc(i),TPzloc(i) ENDDO END SUBROUTINE READ_INP SUBROUTINE TEMP_REQUEST This subroutine figures out which differential surface contains each of the requested temperature points. USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) PAGE 129 129 DO i=1,numtemppoints IF((T Psfc(i).EQ.3).OR.(TPsfc(i).EQ.4))THEN DO j=1,mxsfc ytest=0 ytest=ddy(3)*j IF(TPyloc(i).LE.ytest)THEN n1TP(i)=j EXIT ENDIF ENDDO DO j=1,mxsfc ztest=0 ztest=ddz(3)*j IF(TPzloc(i).LE.ztest)THEN n2TP(i)=j EXIT ENDIF ENDDO ENDIF IF(TPsfc(i).EQ.5)THEN DO j=1,mxsfc xtest=0 xtest=ddx(5)* j IF(TPxloc(i).LE.xtest)THEN n1TP(i)=j EXIT ENDIF ENDDO DO j=1,mxsfc ytest=0 ytest=ddy(5)*j IF(TPyloc(i).LE.ytest)THEN n2TP(i)=j EXIT ENDIF ENDDO ENDIF ENDDO Open Temperature Profile ouput file OPEN(11,FILE='temps.put') WRITE(11,*) 'Temperature profiles' WRITE(11,*) 'Run Started:',dtinit WRITE(11,*) PAGE 130 130 Write temperature request point map to 'temps.put' WRITE(11,325) 'Requested Locations:','Nearest Node Locations:' WRITE(11,*) 'Point# Sfc# xloc[m] yloc[m] zloc[m] &xloc[m] yloc[m] zloc[m]' DO iTP=1,numtemppoints iTO=TPsf c(iTP) l=n1TP(iTP) m=n2TP(iTP) CALL MAT_INDEX_REV(iTO,l,m,i) WRITE(11,326) iTP,TPsfc(iTP),TPxloc(iTP),TPyloc(iTP), & TPzloc(iTP),Xloc(i),Yloc(i),Zloc(i) ENDDO WRITE(11,*) 326 FORMAT(2X,I2,8X,I1,X,3(F9.2,X),3 X,3(F9.2,X)) Write column headers to 'temps.put' temp request output file WRITE(11,322,advance='no') 'PrintStep', 'Time', 'Timenorm', &'FBTemp' DO i=1,numtemppoints WRITE(11,323,advance='no') 'T(',i,')' ENDDO WRITE(11,*) WRITE(11,324) '[s]' 323 FORMAT(100(A2,I2,A,7X)) 322 FORMAT(A9,3X,A4,6X,A8,6X,A6,9X) 324 FORMAT(12X,A3) 325 FORMAT(16X,A20,13X,A23) END SUBROUTINE TEMP_REQUEST ! SUB AREAS PURPOSE: To Calculate areas of each surface SUBROUTINE AREAS USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) Mesh surfaces into differential areas PAGE 131 131 ddx=dx/DBLE(mxsfc) ddy=dy/DBLE(mxsfc) ddz=dz/DBLE(mxsfc) DO i=1,totsfc IF (dx(i) .EQ. 0.D0) THEN Area(i) = dy(i)*dz(i) diffArea(i) = ddy(i)*ddz(i) ELSE IF (dy(i) .EQ. 0.D0) THEN Area(i) = dx(i)*dz(i) diffArea(i) = ddx(i)*ddz(i) ELSE IF (dz(i) .EQ. 0.D0) THEN Area(i) = dx(i)*dy(i) diffArea(i) = ddx(i)*ddy(i) END IF END DO END SUBROUTINE AREAS SUBROUTINE COMMON_EDGE USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) hr=dy(3) rl=dz(1) wr=dx(1) H=hr/rl W=wr/rl H2=H*H W2=W*W a=W*DATAN(1/W)+H*DATAN(1/H) dsqrt(H2+W2)* & DATAN(dsqrt(1/(H2+W2))) b=0.25*dlog(((1+W2)*(1+H2)/(1+W2+H2))*(W2*(1+W2+H2)/ & ((1+W2)*(W2+H2)))**(W2)*(H2*(1+H2+W2)/((1+H2)* & (H2+W2)))**(H2)) VFside = (1/(W*PI))*(a+b) VF(1,3)=VFside END SUBROUTINE COMMON_EDGE PAGE 132 132 SUBROUTINE OPP_SIDES USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH, O Z) IMPLICIT INTEGER (I N) DO i=1,3,2 IF(i .EQ. 1) THEN c=dx(1) d=dz(1) e=dy(3) ELSE c=dz(3) d=dy(3) e=dx(5) END IF X=c/e Y=d/e X2=X*X Y2=Y*Y f=dlog(((1+X2)*(1+Y2)/(1+X2+Y2))**0.5)+X*dsqrt(1+Y2)* & DATAN(X/dsqrt(1+Y2)) g=Y*dsqrt(1+X2)*DATAN(Y/dsqrt(1+X2)) X*DATAN(X) Y*DATAN(Y) VFends=(2/(PI*X*Y))*(f+g) VF(i,i+1)=VFends END DO 30 FORMAT(100F8.3) END SUBROUTINE OPP_SIDES SUBROUTINE PERP_PLANE(a,b,c,d,G) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) IF (a .EQ. 0) THEN a = 0.000000000000000000001D0 END IF PAGE 133 133 IF (b .EQ. 0) THEN b = 0.000000000000000000001D0 END IF IF (c .EQ. 0) THEN c = 0.000000000000000000001D0 END I F IF (d .EQ. 0) THEN d = 0.000000000000000000001D0 END IF IF (b .EQ. c) THEN b=1.00000000000000001D0*b END IF a2=a*a d2=d*d G = (1D0/(2D0*PI))*((bc)*DSQRT(a2+d2)*DATAN((b c)/DSQRT(a2+d2)) & 0.25D0*(a2+d2(b c)**2)*DLOG(a2+d2+(b c)**2)) END SUBROUTINE PERP_PLANE SUBROUTINE COMMON_EDGE_DIFF(iFROM, iTO, iNegVFdiff, iVFdiff) USE PARAMS US E CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) Differential VF for perpendicular finite rectangles In earlier versions of CANYON these IF statements were necessary for different meshing in each dimension, which was subsequently \ removed. IF (((iFROM.EQ.1).OR.(iFROM.EQ.2)) &.AND.((iTO.EQ.3).OR.(iTO.EQ.4)))THEN mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf=ddx(iFROM) ddc1=ddz(iFROM) ddc2=dd z(iTO) ddt=ddy(iTO) ELSEIF(((iFROM.EQ.1).OR.(iFROM.EQ.2)) &.AND.((iTO.EQ.5).OR.(iTO.EQ.6))) THEN mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc PAGE 134 134 ddf=ddz(iFROM) ddc1=ddx(iFROM) ddc2=ddx(iTO) ddt =ddy(iTO) ELSEIF(((iFROM.EQ.3).OR.(iFROM.EQ.4)) &.AND.((iTO.EQ.1).OR.(iTO.EQ.2))) THEN mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf=ddy(iFROM) ddc1=ddz(iFROM) ddc2=ddz(iTO) ddt=ddx(iTO) ELSEIF(((iFROM.EQ.3).OR.(iFROM.EQ.4)) &.AND.((iTO.EQ.5).OR.(iTO.EQ.6))) THEN mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf=ddz(iFROM) ddc1=ddy(iFROM) ddc2=ddy(iTO) ddt=ddx(iTO) ELSEIF(((iFR OM.EQ.5).OR.(iFROM.EQ.6)) &.AND.((iTO.EQ.1).OR.(iTO.EQ.2))) THEN mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf=ddy(iFROM) ddc1=ddx(iFROM) ddc2=ddx(iTO) ddt=ddz(iTO) ELSEIF(((iFROM.EQ.5).OR.(iFROM .EQ.6)) &.AND.((iTO.EQ.3).OR.(iTO.EQ.4))) THEN mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf=ddx(iFROM) ddc1=ddy(iFROM) ddc2=ddy(iTO) ddt=ddz(iTO) ENDIF Redefine in terms of dummy variables "from", "common", & "to" DO nf=1,mfsfc DO nc=1,mcsfc DO mc=1,mcsfc PAGE 135 135 DO mt=1,mtsfc DO i=1,2 sX(i)=ddf*DBLE(nf+i 2) sY(i)=ddc1*DBLE(nc+i 2) sW(i)=ddc2*DBLE(mc+i 2) sZ(i)=ddt*DBLE(mt+i 2) ENDDO sum=0 DO L=1,2 DO k=1,2 DO j=1,2 DO i=1,2 a = DBLE(sX(i)) b = DBLE(sY(j)) c = DBLE(sW(k)) d = DBLE(sZ(L)) CALL PERP_PLANE(a,b,c,d,G) sign = DBLE(( 1)**(i+j+k+L)) sum = sum +(sign)*DBLE(G) ENDDO ENDDO ENDDO ENDDO rMult = (1/((sX(2)sX(1))*(sY(2)sY(1)))) VFdTEST = rMult*sum iVFdiff = iVFdiff + 1 IF (VFdTEST .LT. 0) THEN iNegVFdiff = iNegVFdiff + 1 END IF Differential VFs are stored in the proper index of the VFdiff array IF ((iFROM.EQ.1).AND.(iTO.EQ.3))THEN VFdiff(iFROM,nf,nc,iTO,mt,mc) = VFdTEST ELSEIF((iF ROM.EQ.1).AND.(iTO.EQ.4))THEN VFdiff(iFROM,(mfsfc nf+1),nc,iTO,mt,mc)=VFdTEST ELSEIF((iFROM.EQ.1).AND.(iTO.EQ.5)) THEN VFdiff(iFROM,nc,nf,iTO,mc,mt) = VFdTEST ELSEIF((iFROM.EQ.1).AND.(iTO.EQ.6))THEN VFdiff(iFROM,nc,(mf sfcnf+1),iTO,mc,mt)=VFdTEST  PAGE 136 136 ELSEIF((iFROM.EQ.2).AND.(iTO.EQ.3))THEN VFdiff(iFROM,nf,nc,iTO,(mtsfc mt+1),mc)=VFdTEST ELSEIF((iFROM.EQ.2).AND.(iTO.EQ.4))THEN VFdiff(iFROM,(mfsfc nf+1),nc,iTO,(mtsfc mt+1),mc)=VFdTEST ELSEIF((iFROM.EQ.2).AND.(iTO.EQ.5))THEN VFdiff(iFROM,nc,nf,iTO,mc,(mtsfc mt+1))=VFdTEST ELSEIF((iFROM.EQ.2).AND.(iTO.EQ.6))THEN VFdiff(iFROM,nc,(mfsfc nf +1),iTO,mc,(mtsfc mt+1))=VFdTEST ELSEIF((iFROM.EQ.3).AND.(iTO.EQ.1)) THEN VFdiff(iFROM,nf,nc,iTO,mt,mc) = VFdTEST ELSEIF((iFROM.EQ.3).AND.(iTO.EQ.2))THEN VFdiff(iFROM ,(mfsfc nf+1),nc,iTO,mt,mc)=VFdTEST ELSEIF((iFROM.EQ.3).AND.(iTO.EQ.5)) THEN VFdiff(iFROM,nc,nf,iTO,mt,mc) = VFdTEST ELSEIF((iFROM.EQ.3).AND.(iTO.EQ.6))THEN VFdiff(iFROM,nc,(mfsfc nf+1),iTO,mt,mc) = VFdTEST ELSEIF((iFROM.EQ.4).AND.(iTO.EQ.1))THEN VFdiff(iFROM,nf,nc,iTO,(mtsfc mt+1),mc)=VFdTEST ELSEIF((iFROM.EQ.4).AND.(iTO.EQ.2))THEN VFdiff(iFROM,(mfsfc nf+1),nc,iTO,(mtsfc mt+1),mc)=VFdTEST ELSEIF((iFROM.EQ.4).AND.(iTO.EQ.5))THEN VFdiff(iFROM,nc,nf,iTO,(mtsfc mt+1),mc) = VFdTEST ELSEIF((iFROM.EQ.4).AND.(iTO.EQ.6))THEN VFdiff(iFROM,nc,(mfsfc nf+1),iTO,(mtsfc mt+1),mc)=VFdTEST ELSEIF((iFROM.EQ.5).AND.(iTO.EQ.1)) THEN VFdiff(iFROM,nc,nf,iTO,mc,mt) = VFdTEST ELSEIF((iFROM.EQ.5).AND.(iTO.EQ.2)) THEN VFdiff(iFROM,nc,(mfsfc nf+1),iTO,mc,m t)=VFdTEST ELSEIF((iFROM.EQ.5).AND.(iTO.EQ.3)) THEN VFdiff(iFROM,nf,nc,iTO,mc,mt) = VFdTEST ELSEIF((iFROM.EQ.5).AND.(iTO.EQ.4))THEN PAGE 137 137 VFdiff(iFROM,(mcsfc nf+1),nc,iTO,mc,mt)=VFdTEST ELSEIF((iFROM.EQ.6).AND.(iTO.EQ.1))THEN VFdiff(iFROM,nc,nf,iTO,mc,(mtsfc mt+1))=VFdTEST ELSEIF((iFROM.EQ.6).AND.(iTO.EQ.2))THEN VFdiff(iFROM,nc,(mfsfc nf+1),iTO,mc,(mtsfc mt+1))=VFdTEST ELSEIF((iFROM.EQ.6) .AND.(iTO.EQ.3))THEN VFdiff(iFROM,nf,nc,iTO,mc,(mtsfc mt+1))=VFdTEST ELSEIF((iFROM.EQ.6).AND.(iTO.EQ.4))THEN VFdiff(iFROM,(mfsfc nf+1),nc,iTO,mc,(mtsfc mt+1))=VFdTEST ENDIF END DO END DO END DO END DO END SUBROUTINE COMMON_EDGE_DIFF SUBROUTINE OPP_PLANE(a,b,c,d,z,G) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) IF (a .EQ. 0) THEN a=1.D 20 END IF IF (b .EQ. 0) THEN b= 1.D 20 END IF IF (c .EQ. 0) THEN c=1.D 20 END IF IF (d .EQ. 0) THEN d=1.D 20 END IF I F (b .EQ. c) THEN b=1.00000000000001D0*b END IF PAGE 138 138 G = (1D0/(2D0*PI))*( (bc)*DSQRT((a d)**2+z**2)*DATAN((b c)/ & DSQRT((a d)**2+z**2))+(ad)*DSQRT((b c)**2+z**2)*DATAN((a d)/ & DSQRT((b c)**2+z**2)) (z**2/2D0)*DLOG((ad)**2+(bc)**2+z**2)) END SUBROUTINE OPP_PLANE SUBROUTINE OPP_SIDE_DIFF(iFROM, iTO, iNegVFdiff, iVFdiff) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) Differential VF for perpendicular finite rectangles In earlier versions of CANYON these IF statements were necessary for different meshing in each di mension, which was subsequently \ removed. IF ((iFROM.EQ.1.AND.iTO.EQ.2).OR.(iFROM.EQ.2.AND.iTO.EQ.1)) THEN z = dy(3) mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf1=ddz(iFROM) ddf2=ddx(iFROM) ddt1=ddx(iTO) ddt2=ddz(iTO) ELSEIF((iFROM.EQ.3.AND.iTO.EQ.4).OR. & (iFROM.EQ.4.AND.iTO.EQ.3))THEN z = dx(1) mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf1=ddy(iFROM) ddf2=ddz(iFROM) ddt1=ddz(iTO) ddt2=ddy(iTO) ELSEIF((iFROM.EQ.5.AND.iTO.EQ.6).OR. & (iFROM.EQ.6.AND.iTO.EQ.5))THEN z = dz(1) mfsfc=mxsfc mcsfc=mxsfc mtsfc=mxsfc ddf1=ddx(iFROM) ddf2=ddy(iFROM) ddt1=ddy(iTO) PAGE 139 139 ddt2=ddx(iTO) ENDIF Redefine in terms of dummy variables "from", "common", & "to" DO nf=1,mfsfc DO nc=1,mcsfc DO mc=1,mcsfc DO mt=1,mtsfc DO i=1,2 sX(i)=ddf1*DBLE(nf+i 2) sY(i)=ddf2*DBLE(nc+i 2) sW(i)=ddt1*DBLE(mc+i 2) sZ(i)=ddt2*DBLE(mt+i 2) END DO sum=0 DO L=1,2 DO k=1,2 DO j=1,2 DO i=1,2 a = DBLE(sX(i)) b = DBLE(sY(j)) c = DBLE(sW(k)) d = DBLE(sZ(L)) CALL OPP_PLANE(a,b,c,d,z,G) sign = DBLE(( 1)**(i+j+k+L)) sum = sum +(sign)*DBLE(G) END DO END DO END DO END DO rMult = (1.D0/((sX(2)sX(1))*(sY(2)sY(1)))) VFdTEST = rMult*sum iVFdiff = iVFdiff + 1 IF (VFdTEST .LT. 0) THEN iNegVFdiff = iNegVFdiff + 1 END IF IF(((iFROM.EQ.1).AND.(iTO.EQ.2)).OR. & ((iFROM.EQ.2).AND.(iTO.EQ.1)))THEN VFdiff(iFROM,nc,nf,iTO,(mcsfc mc+1),(mtsfc mt+1))=VFdTEST ELSEIF(((iFROM.EQ.3).AND.(iTO.EQ.4)).OR. & ((iFROM.EQ.4).AND.(iTO.EQ.3)))THEN VFdiff(iFROM,nf,nc,iTO,mt,mc) = VFdTEST PAGE 140 140 ELSEIF(((iFROM.EQ.5).AND. (iTO.EQ.6)).OR. & ((iFROM.EQ.6).AND.(iTO.EQ.5)))THEN VFdiff(iFROM,nf,nc,iTO,mt,mc) = VFdTEST ENDIF END DO END DO END DO END DO END SUBROUTINE OPP_SIDE_DIFF SUBROUTINE VF_DIFF_CHECKER USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) WRITE(*,*) 'Differential VFs are stored in the format VFDiff &(FROMsfc,nA,nB;TOsfc,mA,mB)' WRITE(*,*) WRITE(*,*) 'Enter desired differential VF as 6digit string &separated by spaces' WRITE(*,*) '(i.e. "1 2 3 4 5 6" to request VFDiff(1,2,3;4,5,6)]:' 109 READ(*,*) iFROM, iFA, iFB, iT O, iTA, iTB WRITE(*,70) 'VFDiff(',iFROM,',',iFA,',',iFB,';',iTO,',',iTA,',' &,iTB,')','=',VFDiff(iFROM,iFA,iFB,iTO,iTA,iTB) 70 FORMAT(A7,6(I2,A1),A1,D20.10) WRITE(*,*) 'Again?' READ(*,*) VDC IF((VDC.EQ.'n').OR.(VDC.EQ .'N').OR.(VDC.EQ.'No').OR. & (VDC.EQ.'no')) THEN GOTO 110 ELSE GOTO 109 ENDIF 110 END SUBROUTINE VF_DIFF_CHECKER ! SUB CLKDTO PAGE 141 141 PURPOSE OF CLKDTO ROUTINE: To get hundreths of seconds from the F90 date and time NOTES: This should be compatible with ANSI F90 on all platforms Note hdsec is a hundredths of sec arbitrary time hack AUTHOR (date): G.E. SJODEN (from CLKDTM kernel) SUBROUTINE CLKDTO (dtm, hdsec, strstp) USE PARAMS USE CANDATA IMPLICIT DOUBLE PRECISION (AH,O Z) IMPLICIT INTEGER (I N) CHARACTER(19) :: dtm INTEGER, DIMENSION(8) :: ndt DOUBLE PRECISION :: strstp Using Intrinsic ANSI F90 Commands CALL DATE_AND_TIME(dtm,dtm,dtm,ndt) hsec=(((ndt(5)*3600.)+(ndt(6) *60.)+ndt(7))*1.D2)+(ndt(8)*0.10) hdsec=(((ndt(2)*30.5)+ndt(3))*8.64D+06) + hsec Write the ndt integer array to character variable dtm dtm=' WRITE(dtm,10) ndt(2),ndt(3),ndt(1),ndt(5),ndt(6),ndt(7) 10 FORMA T(I2.2,'/',I2.2,'/',I4.4,1X,I2.2,':',I2.2,':',I2.2) RETURN END SUBROUTINE CLKDTO ! SUB HEADR PURPOSE OF HEADR ROUTINE: To print header when program is initiated AUTHOR (date): T.F. Stachitas (Mar 2009, etc) PAGE 142 142 SUBROUTINE HEADR Open primary output file OPEN(1, FILE='CANout.put', ACCESS='SEQUENTIAL') WRITE(*,'(A2) ')' WRITE(*,'(2X,A48)')' CANYON View Factor Code WRITE(*,'(2X,A48)')' Version 2.0 DOUBLE PRECISION WRITE(*,'(2X,A48)')' WRITE(*,'(2X,A48)')' T. Stachitas, T.A. Mock, G. Sjoden WRITE(*,'(2X,A48)')' Jun 2009 WRITE(*,'(A2)')' WRITE(1,'(A2) ')' WRITE(1,'(2X,A48)')' CANYON View Factor Code WRITE(1,'(2X,A48)')' Version 2.0 DOUBLE PR ECISION WRITE(1,'(2X,A48)')' WRITE(1,'(2X,A48)')' T. Stachitas, T.A. Mock, G. Sjoden WRITE(1,'(2X,A48)')' Jun 2009 WRITE(1,'(A2)')' END SUBROUTINE HEADR  PAGE 143 143 APPENDIX C CANYON SAMPLE INPUT FILE This is the CANYON input file for the baseline case discussed in Chapters 4 and 5. The box outline emphasizes that it should be copied exactly as formatted in order to be read properly by CANYON. In order to use the file, the contents of the box should be saved to a text file called canyon.inp. When altering input values the spacing within a line ca n be changed (i.e. format free), but the spacing between lines must remain consistent with this example. Input file for CANYON: 3 D Radiant Heat Transfer code Surface mesh [number of meshes each direction] mxsfc 8 Surface Dimensions [m]: s fc dx dy dz Description 1 40 0 100 'Canyon Entrance (Fireball)' 2 40 0 100 'Canyon Exit' 3 0 1610 100 'Left wall (Concrete)' 4 0 1610 100 'Right wall (Concrete)' 5 40 1610 0 'Floor (Road surface)' 6 40 1610 0 'Sky' Initial Surface Temperatures [K]: sfc T 1 2860 Fireball 2 300 Canyon Exit Plane 3 300 Concrete (Left wall) 4 300 Concrete (Right wall) 5 300 Asphalt (Road surface) 6 300 Sky Heat Transfer Surface Properties: Sfc jVis Emittance(eps) Specific Heat (Cp) Density(rho) Conductivity(rk) [J/kgK] [kg/m^3] [W/mK] 1 1 1.0 840.0 1800.0 1.50 2 0 1.0 840.0 1800.0 1.50 3 1 0.7 840.0 1800.0 1.50 4 1 0.7 840.0 1800.0 1.50 5 1 0.7 840.0 1800.0 1.50 PAGE 144 144 6 0 1.0 840.0 1800.0 1.50 Sky Factor: 1.0 Simulation Parameters: Yield[kT] tnormmax deltaTime[s] nu mout tol Flux Method (1=dir 2=iterat e ) 20.0 20.5355 5D 7 1000 1D 4 2 iTrans (0=steady state, 1=transient, 2=constant T's) 1 Atmospheric Conditions: iVis[km] TempFactor T empSource (1=embedded BB temps, ELSE= tempdata.dat ) 999 1.00 1 Temperature Output Points: # Points 6 Point# Sfc# xloc[m] yloc[m] zloc[m] 1 3 0 300 7 2 3 0 900 7 3 3 0 1500 7 4 5 20 300 0 5 5 20 900 0 6 5 20 1500 0 PAGE 145 145 APPENDIX D CANYON SAMPLE OUTPUT FILE (ABRIDGED) This appendix is an abridged version of the primary output file for the baseline case Intermediate time step outputs have been skipped but the file is otherwise complete. The output file is boxed to emphasize its exact appearance. PAGE 146 146 CANYON View Factor Code Version 2.0 DOUBLE PRECISION T. Stachitas, T.A. Mock, G. Sjoden Jun 2009 Run Started:10/18/2009 16:45:49 Input File Name: canyon.inp totsfc = 6 mxsfc = 8 ( Surface mesh parameter) Surface Dimensions [m] sfc dx dy dz Area[m^2] Description 1 40. 0. 100. 4000. Canyon Entrance (Fireball) 2 40. 0. 100. 4000. Canyon Exit 3 0. 1610. 100. 161000. Left wall (Concrete) 4 0. 1610. 100. 161000. Right wall (Concrete) 5 40. 1610. 0. 64400. Floor (Road surface) 6 40. 1610. 0. 64400. Sky Initial Surface Temperatures [K] 1 2860. 2 300. 3 300. 4 300. 5 300. 6 300. WholeSurface Diffuse View Factor Matrix: [ VF(row,column) ] 1 2 3 4 5 6 1 0.00000 0.00049 0.34523 0.34523 0.15453 0.15453 2 0.0 0049 0.00000 0.34523 0.34523 0.15453 0.15453 PAGE 147 147 3 0.00858 0.00858 0.00000 0.66448 0.15918 0.15918 4 0.00858 0.00858 0.66448 0.00000 0.15918 0.15918 5 0.00960 0.00960 0.39796 0.39796 0.00000 0.18488 6 0.00960 0.00960 0.39796 0.39796 0.18488 0.00000 Differential View Factor Summary Data: 384 differential surfaces defined 147456 differential VFs calculated 0 negative differential VFs calculated Heat Transfer Properties: Sfc jVis eps Cp[J/kgK] rho[kg/m^3] rk[W/mK] alphaTH[m^2/s] 1 1 1.0000 840.0000 1800.0000 1.5000 0.99206E 06 2 0 1.0000 N/A N/A N/A N/A 3 1 0.7000 840.0000 1800.0000 1.5000 0.99206 E 06 4 1 0.7000 840.0000 1800.0000 1.5000 0.99206E 06 5 1 0.7000 840.0000 1800.0000 1.5000 0.99206E 06 6 0 1.0000 N/A N/A N/A N/A Sky Factor: 1.000 Simulation Parameters: Yield[kT] tnormmax dtime[s] tol # Print Steps # Time Steps 20.00 20.50 5.000E 07 1.000E 04 1000 6388120 Atmospheric Conditions: iVis[km ] TempFactor Temp Source File Name: 999 1.000 none SUMMARY RESULTS: Surface 2 Average Flux & Fluence P rintStep Time Timenorm FBTemp FluxINavg FluenceAvg FluxRow1INavg [s] [K] [cal/(cm^2*s)] [W/cm^2] [cal/cm^2] [cal/(cm^2*s)] [W/cm^2] PAGE 148 148 1 0.0032 0.0205 1461.3094 1.4043765E 02 5.8759112E 02 3.8321288E 05 1.4041662E 02 5.8750316E 02 2 0.0064 0.0410 2066.5214 2.3261979E 02 9.7328120E 02 9.6263307E 05 2.3252934E 02 9.7290278E 02 3 0.0096 0.0615 2530.9212 3.8668285E 02 1.6178811E 01 1.9350525E 04 3.8635631E 02 1.6165148E 01 4 0.0128 0.0820 2922.4144 6.0641591E 02 2.5372442E 01 3.5023208E 04 6.0454663E 02 2.5294231E 01 5 0.0160 0.1025 3267.2926 9.0300160E 02 3.7781587E 01 5.8916732E 04 8.9493719E 02 3.7444172E 01 6 0.0192 0.1230 3579.0244 1.2693106E 01 5.3107956E 01 9.3449503E 04 1.2524965E 01 5. 2404453E 01 7 0.0224 0.1435 3865.5924 1.6982478E 01 7.1054688E 01 1.4066172E 03 1.6723938E 01 6.9972959E 01 8 0.0256 0.1640 4132.1773 2.2001952E 01 9.2056166E 01 2.0271846E 03 2.1617302E 01 9.0446791E 01 9 0.0287 0.1845 4382.3546 2.7741940E 01 1.1607228E+00 2.8198400E 03 2.7196735E 01 1.1379114E+00 10 0.0319 0.2050 4618.7030 3.4109885E 01 1.4271576E+00 3.8059796E 03 3.3398191E 01 1.3973803E+00 991 3.1653 20.3151 2264.2535 3.5419351E 02 1.4819456E 01 1.5341558E+00 3.5978755E 02 1.5053511E 01 992 3.1684 20.3356 2263.1120 3.5369233E 02 1.4798487E 01 1.5342689E+00 3.5927402E 02 1.5032025E 01 993 3.1716 20.3561 2261.9722 3.5321118E 02 1.4778356E 01 1.5343818E+00 3.5878102E 02 1.5011398E 01 994 3.1748 20.3766 2260.8341 3.5271228E 02 1.4757482E 01 1.5344945E+00 3.5826982E 02 1.4990009E 01 995 3.1780 20.3971 2259.6977 3.5221488E 02 1.4736671E 01 1.5346071E+00 3.5776017E 02 1.4968685E 01 996 3.1812 20.4176 2258.5631 3.5171898E 02 1.4715922E 01 1.5347195E+00 3.5725204E 02 1.4947425E 01 997 3.1844 20.4381 2257.4301 3.5124349E 02 1.4696028E 01 1.5348318E+00 3.5676482E 02 1.4927040E 01 998 3.1876 20.4586 2256.2989 3.5074952E 02 1.4675360E 01 1.5349439E+00 3.5625868E 02 1.4905863E 01 999 3.1908 20.4791 2255.1693 3.5025621E 02 1.4654720E 01 1.5350558E+00 3.5575320E 02 1.4884714E 01 1000 3.1940 20.4996 2254.0415 3.4978313E 02 1.4634926E 01 1.5351676E+00 3.5526844E 02 1.4864432E 01 Run Completed:10/18/2009 19:10: 51 CPU time est is 8702.65880549567 sec. PAGE 149 149 REFERENCES Bridgman, C.J., 2001. Introduction to the Physics of Nuclear Weapons Effects. Defense Threat Reduction Agency and Air Force Institute of Technology. engel, Y.A., 2007. Heat and Mass Transfer: A Practical Approach. McGraw Hill, B oston. Compaq Fortran Engineering Group, 2001. Compaq Visual Fortran Version 6.6. Clarke, J.A., Yaneske, P.P., Pinney, A.A., 1991. The harmonization of thermal properties of building materials. Building Research Establishment, Watford, UK. Dassault Systmes SolidWorks Corp, 2009. SolidWorks 2009 x64 Edition SP4.1. (http://www.solidworks.com ) Oct 2009. Duderstadt, J.J., Hami lton, L.J., 1976. Nuclear Reactor Analysis. John Wiley & Sons, Ann Arbor, MI. Ehlert, J.R., Smith, T.F., 1992. View factors for perpendicular and parallel rectangular plates. Journal of Thermophysics Vol 7, p. 173. Fluent Inc., 2009. ANSYS FLUENT 12.0. ( http://www.ansys.com/products/fluiddynamics/fluent/default.asp ) Oct 2009. Glasstone, S., Dolan, P.J., 1977. The Effects of Nuclear Weapons. Department of Defense, Washington D. C. Howell, J.R., 2001. A Catalog of Radiation Heat Transfer Configuration Factors. 2nd Ed. University of Texas at Austin, Austin, TX. (http://www.me.utexas.edu/~howell/index.html ) Sep 2009. Kaw, A., Kalu, E.E., 2009. Numerical Methods with Applications. University of South Florida, Tampa, FL. (http://numericalmethods.eng.usf.edu/topics/textbook_index.html ) Sep 2009. Leonard, J.M., 1985. Nuclear thermal transmittance in the atmosphere using LOWTRAN 6 computer code. AFIT/GNE/ENP/85M 13. Air Force Institute of Technology, Wright Patterson AFB, OH. Marrs, R.E., Moss, W.C., Whitlock, B., 2007. Thermal radiation from nuclear deton ations in urban environments. UCRLTR 231593. Lawrence Livermore National Laboratory. Medalia, J., 2008. Comprehensive Nuclear Test Ban Treaty: Background and Current Developments. RL33547. Congressional Research Service. (http://fpc.state.gov/documents/organization/106157.pdf ) Sep 2009. PAGE 150 150 Modest, M.F., 2003. Radiative Heat Transfer. Academic Press, San Diego, CA. Sparrow, E.M., 1965 Radiant emission, absorption, and transmission characteristics of cavities and passages. Symposium on Thermal Radiation of Solids San Francisco 1964. Scientific and Technical Information Division, NASA, Washington D.C. Tummers, B., 2006. Data Thief III. ( http://ww w.datathief.org ) Sep 2009. PAGE 151 151 BIOGRAPHICAL SKETCH Tucker Flagg Stachitas was born in 1986 in Philadelphia, Pennsylvania where he lived for 12 years before moving with his family to Jacksonville, Florida. He received his International Baccalaureate diploma from Allen D. Nease High School in 2004. He attended the United States Naval Academy, graduating with distinction in 2008 with a bachelors degree in mechanical engineering and a minor in Spanish. While at the Naval Academy he was an 8time All American on the International Pistol Team, a 2008 National Rifle Association Intercollegiate Team Champion, and a member of the Junior National Team for the 2006 International Shooting Sports Federation World Shooting Championships. His hobbies include playing piano, guitar, and Ultimate Frisbee. 