AN INVESTIGATION OF BLAST WAVES GENERATED BY CONSTANT VELOCITY FLAMES by Robert Thomas Luckritz Dissertation submitted to the Faculty of the Graduate School of the University of Maryland in partial fulfillment of the requirements for the degree of Doctor of Philosophy 1977 Ct ,.) .i. \ APPROVAL SHEET Title of Dissertation: An Investigation of Blast Waves Generated by Constant Velocity Flames Name of Candidate: Robert Thomas Luckritz Doctor of Philosophy, 1977 • 1 I / 1 , Dissertation and Abstract Approved: ; ~:~p~·~{!1M~;~h~ifcf (; . Date Approved: Professor Chemical Engineering University of Maryland ABSTRACT Title of Dissertation: An Investigation of Blast Waves Generated by Constant Velocity Flames Robert Thomas Luckritz, Doctor of Philosophy, 1977 Dissertation directed by: Joseph M. Marchello Professor Chemical Engineering The relevant flow field parameters associated with the generation and propagation of blast waves from constant vel- ocity flames were systematically studied through numerical integrations of the non-steady equations for mass, momentum, and energy. The flow was assumed to be that of an adiabatic inviscid fluid obeying the ideal gas law and the flame was simulated by a working fluid heat addition model. The flame velocity was varied from infinitely fast (bursting sphere) through velocities characterized by the nearly constant pressure deflagration associated with low Mach number laminar flames. The properties noted included peak pressure, positive impulse, energy distrib.ution, and the blast wave flow field. Results were computed for the case of a methane-air mixture assuming an energy density, q = 8.0, an ambient spe- cific heat ratio, y 0 = 1.4 and a specific heat ratio behind the flame, y 4 = 1.2. In the source volume, as the flame velocity decreased to Mach 4.0 the overpressure increased. For flame velocities below Mach 4.0 the overpressure decreased, and approach the acoustic solution originally developed by Taylor. In the far field the overpressure curves for super- sonic flame velocities coalesced to a common curve at approxi- mately 70% of Baker's pentolite correlation. Far field overpressures for subsonic flame velocities decreased as the flame velocity decreased. For the flame velocities investigated the near field impulse was greater than the impulse from Baker's pentolite correlation. In the far field the flame generated impulse decreased to 60 to 75% of the pentolite impulse. In cases where the flow was expected to reduce to a self-similar solution and/or show Rayleigh line behavior it did. The calculations showed that the flow field behaved normally where expected, and for flow velocities where steady state behavior is not expected, non-steady behavior was observed. ACKNOWLEDGEMENTS The author extends a special thank you to Professor Joseph M. Marchello and Professor Roger A. Strehlow for their help and guidance during all phases of this research. To Mrs. Florence Goldsworthy who assisted in the prep- paration and typing of this manuscript, the author expresses his sincere gratitude. Sincere appreciation is extended to the U.S. Coast Guard for the opportunity to pursue this degree and the financial support provided. The computer time for this project was supported in part through the facilities of the Computer Science Center of the University of Maryland. TABLE OF CONTENTS ACKNOWLEDGEMENTS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii LIST OF TABLES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii CHAPTER I. INTRODUCTION............................ ...... 1 A. Ideal (Point-Source) Blast Waves.... 3 B. Non-Ideal Blast Waves............... 7 C. Homogeneous Energy Addition Blast Waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 D. Constant Velocity Flame Blast Waves. 18 E. Problem Definition .................. 20 II. THEORETICAL CONSIDERATIONS ................... 24 A. Governing Equations ................. 24 B. Steady One-Dimensional Flow Discon- tinuity Relationships .............. 27 C. Energy Scaling ...................... 38 D. Damage Equivalence .................. 40 E. Damage Mechanisms ................... 41 III. COMPUTATIONAL PROCEDURE ...................... 47 A. Boundary Conditions ................. 49 B. Dimensionless Variables ............. 50 C. Source Model. . . . . . . . . . . . . . . . . . . . . . . . 52 1. Energy Addition Wave ............ 53 2. Change of Specific Heat Ratio ... 63 D. Numerical Integration ............... 64 iii E. Testing of Program ................... 68 1. Bursting Plane ................... 69 2. Bursting Sphere .................. 69 3. Wave Width................... . . . . 71 IV. RESULTS AND DISCUSSION . . . . .. .. . .. . . . .. .. . . . . . 77 A. Flow Field Properties from One-Dimen- sional Steady State Theory ........... 77 B. The Effects of Energy Addition Wave Velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 1. Flow Field Properties ............ 80 2. Damage Parameters ................ 97 3. Energy Distribution .............. 100 V. COMPARISONS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3 A. Blast Waves Generated from Non-Ideal Energy Sources ........................ 163 B. Some Aspects of Blast from Fuel-Air Explosives ............................ 169 C. Pressure Waves Generated by Steady Flames ................................ 171 D. The Air Wave Surrounding an Expanding Sphere ................................ 174 VI. CONCLUSIONS/RECOMMENDATIONS .. . ................ 194 A. Conclusions. . . . . . . . . . . . . . . . . . . . . . . . . . . 194 B. Recommendations ....................... 197 iv APPENDIX A. Computer Program for the Model ........ 199 B. Computer Program for Analyzing Data . . . 225 NOMENCLATURE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229 BIBLIOGRAPHY. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 234 V LIST OF TABLES 1. Hugoniot Curve-Fit Data .. . ..... . ....... . ............. 54 2. Summary of Parameters Used for Cases Investigated ... . .......................... 55 vi LIST OF ILLUSTRATIONS 1. Pressure-Time Relationship for an Ideal Blast Wave.. 5 2. Motion in a Shock Tube.............................. 15 3. Three-Dimensional Diagram of Three Parameters affecting Non-ideal Blast Wave Behavior ............ 22 4. Pressure-Volume Plot of End States for a One-Dimen- tional Steady Process With Heat Addition ........... 30 5. Scaled P-I Curve for Fixed Level of Damage .......... 43 6. Schematic Diagram of Spring-Mass System to Model the Dynamic Response of a Structural Member.·............................................ 44 7. Energy Deposition Function .......................... 59 8. Wave Width of Energy Deposition Term ................ 60 9. Deposition Time versus Inverse Mach Number as Function of Wave Width ............. ·................ 62 10. Computational Grid for Finite Differencing Technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 11. Pressure Distribution from a Mach 4225. Energy Wave............................................... 70 12. Overpressure versus Energy Scaled Distance .......... 73 13. Overpressure versus Energy Scaled Distance .......... 75 14. Overpressure Distribution through Energy Addition Wave ............................................... 109 15. Pressure Distribution from an Infinite Velocity Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 16. Pressure Distribution from a Mach 8.0 Energy Wave ... 111 17. Pressure Distribution from a Mach 5.2 (Planar Geometry) Energy Wave .............................. 112 18. Pressure Distribution from a Mach 5.2 (Spherical Geometry) Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 19. Pressure Distribution from a Mach 4.0 Energy Wave ... 114 20. Pressure Distribution from a Mach 4.0 Energy Wave ... 115 21. Pressure Distribution from a Mach 3.0 Energy Wave ... 116 vii 22. 117 Pressure Distribution from a Mach 3.0 Energy Wave ... 23. 118 Pressure Distribution from a Mach 2.0 Energy Wave ... 24. 119 Pressure Distribution from a Mach 1.0 Energy Wave ... 25. 120 Pressure Distribution from a Mach 0.5 Energy Wave ... 26. 121 Pressure Distribution from a Mach 0.25 Energy Wave .. 27. 122 Pressure Distribution from a Mach 0.125 Energy Wave. 28. Pressure-Volume Behavior from a Mach 8.0 Energy Wave ............................................... 123 29. Pressure-Volume Behavior from a Mach 5.2 (Planar Geometry) Energy Wave .............................. 124 30. Pressure-Volume Behavior from a Mach 5.2 (Spherical Geometry) Energy Wave .............................. 125 31. Pressure-Volume Behavior from a Mach 4.0 Energy Wave ............................................... 126 32. Pressure-Volume Behavior from a Mach 3.0 Energy Wave ............................................... 127 33. Pressure-Volume Behavior from a Mach 2.0 Energy Wave ............................................... 128 34. Pressure-Volume Behavior from a Mach 1.0 Energy Wave .......... . ................. . ............. . .... 129 35. Pressure-Volume Behavior from a Mach 0.5 Energy Wave ..... . ........... . ............................. 130 36. Pressure-Volume Behavior from a Mach 0.25 Energy Wave ............................................... 131 37. Pressure-Volume Behavior from a Mach 0.125 Energy Wave ............................................... 132 38. Pressure-Time Behavior at Eulerian Radius from an Infinite Velocity Energy Wave ...................... 133 39. Pressure-Time Behavior at Eulerian Radius from a Mach. 8. 0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 40. Pressure-Time Behavior at Eulerian Radius from a Mach 5. 2 Energy Wave ............................... 135 41. Pressure-Time Behavior at Eulerian Radius from a Mach 4. 0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 viii 42. Pressure-Time Behavior at Eulerian Radius from a Mach 2. 0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 7 43. Pressure-Time Behavior at Eulerian Radius from a Mach 1.0 Energy Wave ............................... 138 44. Pressure-Time Behavior at Eulerian Radius from a Mach O. 5 Energy Wave ............................. ~ . 139 45. Pressure-Time Behavior at Eulerian Radius from a Mach 0.25 Energy Wave .............................. 140 46. Pressure-Time ·Behavior at Eulerian Radius from a Mach O. 125. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 47. Particle Position-Time Behavior from an Infinite Velocity Energy Wave ............................... 142 48. Particle Position-Time Behavior from a Mach 8.0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 49. Particle Position-Time Behavior from a Mach 5.2 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 50. Particle Position-Time Behavior from a Mach 4.0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 51. Particle Position-Time Behavior from a Mach 2.0 Energy Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 52. Particle Position-Time Behavior from a Mach 1.0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 7 53. Particle Position-Time Behavior from a Mach 0.5 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 54. Particle Position-Time Behavior from a Ma.ch 0.25 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 55. Particle Position-Time Behavior from a Mach 0.125 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 56. Overpressure versus Energy Scaled Distance .......... 151 57. Impulse versus Energy Scale Distance .. . ............. 152 58. Energy Distribution-Time Behavior from an Infinite Velocity Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 59. Energy Distribution-Time Behavior from a Mach 8.0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 ix 60. Energy Distribution-Time Behavior from a Mach 5.2 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 61. Energy Distribution-Time Behavior from a Mach 4.0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 62. Energy Distribution-Time Behavior from a Mach 2 .. 0 · Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 7 63. Energy Distribution-Time Behavior from a Mach 1.0 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 64. Energy Distribution-Time Behavior from a Mach 0.5 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 65. Energy Distribution-Time Behavior from a Mach 0.25 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 66. Energy Distribution-Time Behavior from a Mach 0.125 Energy Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 67. Energy Remaining in Source Volume versus Reciprocal Mach Number of Energy Wave .............. 162 68. Pressure Distribution from TD=0.2 Homogeneous Energy Addition .................................... 176 69. Pressure Distribution from TD=0.2 Homogeneous Energy Addition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 70. Pressure Distribution form TD=2.0 Homogeneous Energy Addition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 71. Pressure-Volume Behavior from TD=0.2 Homogeneous Energy Addition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179 72. Pressure-Volume Behavior from TD=2.0 Homogeneous Energy Addition .................................... 180 73. Pressure-Time Behavior at Eulerian Radius from TD=0.2 Homogeneous Energy Addition ................. 181 74. Pressure-Time Behavior at Eulerian Radius from TD=2.0 Homogeneous Energy Addition ................. 182 75. Particle Position-Time Behavior from TD=0.2 Homogeneous Energy Addition ........................ 183 76. Particle Position-Time Behavior from TD=2.0 Homogeneous Energy Addition ........................ 184 77. Overpressure versus Energy Scaled Distance .......... 185 X 78. Maximum Overpressure versus Source Volume Deposition Time .... . .............................. 186 79. Energy Distribution-Time Behavior from TD=0.2 Homogeneous Energy Addition ........................ 187 80. Energy Distribution-Time Behavior from TD=2.0 Homogeneous Energy Addition ........................ 188 81. Pressure Distribution from Fishburn's Energy Addition ........................................... 189 82. Pressure-Volume Behavior from Fishburn's Energy Addition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 83. Pressure, Particle Velocity, and Temperature Similarity Plots for Energy Addition of Kuhl, et al ..................................... . ........ 191 84. Pressure-Volume Behavior for Energy Addition of Kuhl, et al ..................................... . .. 192 85. Comparison of Pressure and Particle Velocity Distribution of Mach 0.125 Energy Addition Wave with Results Predicted by Taylor ................... 193 xi I. INTRODUCTION The increasing eriergy nee·ds of the United States and other advanced technology countries have resulted in the hand- ling, transportation, and storage of ever increasing quantities of highly volatile and highly combustible fuels. Present projections of energy needs for the future indicate a con- tinued expansion of energy demands in these countries. As with any technological advance the luxuries provided by the use of large quantities of these energy sources are accom- panied by an increased risk in the event of their accidental release. In addition to the ever increasing need for additional fuel, the government and the public have become cognizant of the necessity for protection of the environment from pol - lution by the contaminants present in many of our more abun- dant fuel supplies. Natural Gas is one energy source which is presently available, easily distributed, and rel a tive l y low in pollution potential. However, the supply of easily accessible Natural Gas in the United States is limited and many existing distribution facilities in large metropolitan areas are unable to meet peak winter demands. To alleviate this situation many utilities are storing the natural gas in a liquefied state and/or providing for the importation of shipload quantities from such areas as Alaska, Algeria, Libya, and Indonesia. The release of natural gas from accident, natural dis- aster, or sabotage could subject personnel and facilities 1 2 near the release to great risk, Among these risks are the danger of fire and/or explosion if the release were ignited. The question which concerns both governmental decision makers and the public at large is precisely what would be the effects of large scale releases of a flammable gas such as Natural Gas. Compounding this difficult question is the conclusions which can be extrapolated from accidents as a result of the release of similar exothermic compounds. A survey of accidental explosions that have occurred over the past 40 years was compiled by Strehlow(l). He noted a sharp increase in annual damage from accidental explosions since 1964 and attributed this increase to larger spills of a variety of chemical substances with many spills occuring in the neighborhood of expensive process equipment. In his paper he recommended an investigation into the effects of the overall flame-propogation rate and the nature of the blast wave produced by the deflagrative combustion of a large unconfined vapor cloud. There are also basic fundamental questions concerning the fluid dynamic flow field developed by an accidental ex- plosion. The flow fields generated by high explosives have been investigated in detail for weapons applications and industrial blast technology. To date there has been only minimal effort directed to investigating the effects of ac- cidental (non-ideal) explosions. This dissertation addresses one aspect of accidental 3 (non-ideal) explosions, namely the consequences of the pro- pagation of constant velocity flames after delayed ignition . That is, what happens when there is a large scale release of flammable gas with widespread dispersion of the vapors, followed by ignition? Other related problems such as the effects of a burning pool of f l arrnnable fuel or the effects of rapid release which of the mixture are not does not involve delay· ed 1.·g · · n1.tion addressed. The problem 1.·s presented in terms of a systematic study of the effects of constant Lagrangian velocity flame through a flammable, compressible mixture . The behavior of the flow is studied in the com- pressible medium surrounding the flammable mixture during and after heat addition. A heat addition-working fluid model is used to replace the combustion process. This model and the equations of mass, momentum, and energy coupled with the equation of state are used to study the effects of heat addition waves . Both the near field and far field effects including peak pressure, impulse, and energy distribution were studied to show systematic trends and effects for an energy density approximating that of a stoichiometric mixture of natural gas in air, a common fuel . A. Ideal (Point Source) Blast Waves A blast wave is a pressure wave of finite amplitude generated by the rapid release of energy, such as an explo- sion o The structure will vary as a function of the energy source which produces it . 4 Nuclear and high explosive explosions generate what are known as ideal or point source blast waves . These explosions are described as a finite amount of energy deposited in an infinitely small increment of time at an infinitesimal point in a uniform atmosphere. They generate a shock wave which monotonically decreases in strength as it propagates from the energy source . The properties of the shock wave and the flow associated with it can be determined by solving the non- steady, non-linear equations of fluid mechanics . The Eulerian pressure-time history at a reference point would show ambient conditions until the shock wave arrived at time t , with an almost discontinuous rise to the peak a + over-pressure of the shock wave, Ps + p 0 , as illustrated in figure 1 from Baker(Z) . This peak overpressure, p+ + p, S 0 would be followed by nearly exponential pressure decay through the ambient pressure, p0 , at time ta+ t+, to a min- imum pressure of less than ambient, p0 -p;, then increasing until the pressure again reaches ambient, p0 , at time ta+ t+ + t- . The time during which the pressure is greater than ambient, t through ta + t+ is know as the positive phase. a , The time during which the pressure is negative, t a+ t+ through ta +t++t-, is known as the negative phase . As an ideal (point source) blast wave propagates away from its source there are three regions of interest: (1) The near field wave where pressures are so large that external pressure can be neglected . In this region 5 Positive Phase Pressure I -----t--·~.._-~ I 0 '--------,-----------+' --------+-------- 0 . Time Figure 1. Pressure-time relationship for ideal blast wave. 6 self-similar solutions and analytical formulations are ade- quate. This is followed by (2) An intermediate region of extremely practical im- portance because both the overpressure and impulse are ·suf- ficiently high to do significant damage. The flow field in this region cannot be solved analytically and must be solved numerically. This in turn is followed by (3) A far field region which yields to an analytical approximation involving extrapolation of overpressure-time curves from one location to another. As the shock wave decays, its Mach number approaches unity and the lead wave nears the acoustic limit. There is theoretical evidence that an "N" wave which propagates as an acoustic level phenomena must form. However, atmospheric non-uniformities prevent the observation of this phenomena. Assuming that the atmospheric counterpressure is small when compared to the shock overpressure, a constant value of specific heat, y, and an instantaneous (over an infinitely small time) energy deposition at a point, Taylor( 3), and Sedov(4) reduced the equations of fluid mechanics to non-linear differential equations with one independent variable. These differential equations were then solved to determine the blast wave behavior in the time-space domain. Their analysis determined the pertinent flow variables between the origin and the lead pressure wave and showed that: (1) the particle velocity and density decrease from a maximum value at the shock front to zero at the origin, (2) the preisure decreises, 7 in a nearly exponent ial manner near the shock front, from a maximum v alue a t the h k f t + 1 f · soc ron , Ps , to a va ue o approxi- mately 36% of Ps+ at the origin (for y=l . 4 gas), and (3) the temperature increases without bound as the origin is approached. While investigating these point source solutions , Bethe(S) observed from the shock relations : P2 = (y+l)Ml2 = Po (y-l)M1 2+2 (y -1)+ 2 ~ 1 I-1 where P2 is the density of the fluid immediately behind the shock, P0 is the ambient density of the fluid, and M1 is the approach fluid Mach number that most of the mass in the ' system is concentrated near the shock. As gannna approaches one and the Mach number of the shock becomes large, the effect becomes more pronounced. Using these same conditions it can be seen that in the limit as y~l, p2/p 0 approaches infinity, i.e . all the mass in the system bounded by the lead shock wave is located in or immediately adjacent to the wave . B. Non-Ideal Blast Waves Actual explosions do not generate ideal blast waves. Because of the explosive configuration, the finite reaction time , and the finite volume of the explosive, the pressure wave generated by a real explosion will not follow exactly the time-pressure distribution of an ideal blast wave . Near the energy source which is driving the pressure wave there may be a more 8 gradual build-up 1.·n pressure than the nearly discontinuous Pressure r~se associated with ideal explosions. Other irreg 1 .. u ar1.t1.es such as fragments, . ground effects, re- flections, etc. may cause pressure-time fluctuations incon- sistent · h wit ideal blast wave theory. In general, non-ideal explosions are those where the source energy density is l ow and/or the energy deposi t ion time is long. There a r e an infi ni te number of non- ideal source behaviors that yield blast waves wi th an i nfinite number of different structures, all non-ideal. In point source analysis for ideal blast waves, the assumption is made that initially the energy is added to an infinitely small mass. Therefore, the total energy from the source is available to the surrounding gas to drive the lead shock wave. However, in a real explosion the energy is di- vided between the source volume and the surrounding atmos- phere. Only the energy in the surrounding atmosphere drives the lead shock. This partitioning of energy causes the curves of overpressure vs. radius to lie below the curves from ideal or point source theory. However, as the energy density is in- creased and/or the time of deposition is decreased, as occurs in nuclear or high explosive explosions, the Ps - RE curves approach the ideal (point source) curves. This is attri- buted to the more efficient transmission of energy to the surrounding gas; thereby making more energy available to the shock and nearby flow field. To model the rate of reduction of shock strength caused 9 by the energy which remains in the source volume a non-simi- lar solution in the form of series expansions of key non- dimensional flow parameters was developed by Sakurai(6). He transformed the dependent and independent variables to an- other set where some of the variable were not as sensitive and then expanded each variable as a function of the Mach number squared. The variables were then incorporated into the conservation equations. Solutions, to various orders of accuracy, were obtained by collecting terms of like orders of magnitude and solving each set of differential equations produced, subject to applicable boundary conditions, and calculating the coefficients to the expansions . In the solu- tion he used an energy source with an instantaneous energy deposition time, but indicated that sources with finite times of energy deposition could be modeled . For the second order approximation Sakurai calculated the shock pressure for a y=l . 4 gas to be: 0 . 69 R -1 E +2 . 33 j=O I-2 Ps 1 . 33 R-2 +2 . 16 j=l I-3 = Po E 1 . 96 R- 3 +2 . 07 j=2 I-4 E where I-5 10 Explosion energy per unit area j=O E. = (Explosion energy per J -1 j=l unit l ine)(2~) (Explosion e~ergy)(4n) -1 j=2 and j is the geometry factor (0,1, and 2 for planar, cylin- drical and spherical flow fields respectively) . Data on shock arrival times were obtained by Oshima( 7) from exploding wire experiments and were extensively compared with the predictions calculated by Sakurai . An increase in the range of validity was shown for the higher order approx- imations o These analyses were performed with the assumption that the energy is deposited instantaneously o The heat release which occurs as a result of chemical reaction associated with a reactive fluid-dynamic process has both spatial and temporal dependence . In many cases this invalidates the simplifying self-similar assumptions and the theoritician must resort to numerical integration techniques to obtain a solution. The conservation equations that describe blast waves are three non-linear partial differential equations . Two numerical techniques which have proven useful in the solu- tion of numerous types of non-linear partial differential equations are the method of characteristics, a procedure from the theory of partial differential equations, and, with the development of high speed computers, finite differences o When the finite differencing technique is used for the 11 st udy of blast waves it is preferred to express the conserva- tion equat1.·ons of flu1.·d d · · h · L · ynam1.cs 1.n t e1.r agrangian form . In this method a fluid particle is followed from its initial position to a later position while its intensive properties vary as a function of time . The principle advantages are the computational grid does not distort with time and new grid points can be added as the lead wave uncovers new material . One of the primary areas of interest on the study of blast waves is the generation and propagation of shock waves contained in the flow field and the deviation of these shock waves from those which would be generated in an ideal (point source) explosion . A shock wave can be described as a non- isentropic region in which the fluid properties rapidly change from their initial equilibrium states to a final state in which the temperature, density, and pressure are greater than ahead of the wave . The change in fluid properties occurs within a few mean free path lengths, the average distance a molecule must travel before it is influenced by the presence of another molecule . Because of the steep gradients in the non-isentropic region, the shock can be replaced by either a discontinuity satisfying the Rankine-Hugoniot "jump" rela- tions(S) or, when using finite differencing procedures, by "spreading" this region to one of large but finite gradients over the length of a few computational cells . When perform- ing numerical integrations using the finite differencing technique, gradients within the boundaries are assumed to the finite . Normally the shock is spread over the computational 12 cells by inco · , h rporating into t e momentum and energy equations a ficticious dissipative term developed by Von Neumann and Richtmyer< 9) for their study of the propagation of plane shock They incorporated a dissipation term which was pro- portional to the absolute value of the velocity gradient and only became significant in the shock region. In a later analysis Lax and Wendroff(lO) restricted the magnitude of gradients in strongly compressive regions by using the inherent dissipative mechanism in a modified central differencing scheme which attenuated the high frequency com- ponents of the solution. The application of either dissipative mechanism to es- tablish finite gradients does not violate the conservation of mass, momentum, or energy, as noted by Richtmyer and Morton(ll). The dissipated energy, which is only a minute amount of the total energy, appears as internal energy of the fluid . Von Neumann<12 ) and Brode<13 ) were two of the first to apply the dissipative technique of Von Neumann and Richtmyer to the numerical solution of propagating spherical blast waves " By numerically integrating the differential equations of gas motion in Lagrangian coordinates, Brode determined the strong shock-point source solutions " He determined that the strong shock-point source solu- tions of overpressure versus radius follows the inverse cube law down to an overpressure of approximately 10 atmospheres at which point actual overpressures are 3% higher than ·predicted . 13 U · (2) sing a form of Sachs' scaling, he proposed that the inverse cube relation be replaced by the following t· equa ion for pres- sures greater than 5 atmospheres: P = 0.1567 R- 3 + 1 . s € I-6 For lower pressures he developed the following empirical fit: p = 0 . 137 + 0 . 11~ + 0 . 269 _ 0 . 019 s R 3 R 2 RE € € I-7 0 . 1

O there was a gradual energy addition of finite power . Adamczyk(lB) performed a systematic study of the fluid dynamic and thermodynamic fields associated with the genera- tion and propagation of blast waves from the homogeneous deposition of energy o Using a Von-Neumann/Richtrnyer-type finite difference integration procedure,numerical solutions of the relevant flow parameters were generated by integrating the one-dimensional non-steady fluid dynamic equations of motion in Lagrangian form. Solutions were calculated for planar, cylindrical and spherical flow fields o Varying both the energy density of the source region and the time of energy deposition over two orders of magnitude he noted tha t they both affect the primary causes of structural damage, shock overpressure and positive phase impulse o A two-order 18 of ma · d gn1.tu e change in the time of energy deposition caused th e near field, peak sho"ck overpressure to vary by a factor of 80 and the near field positive-phase impulse ·to vary by a factor of 6. However, he found that the shock front "for- gets" the influence of source non-idealities as it propagates from the origin. D. Constant Velocity Flame Blast Wav·es In the case of delayed ignition of a large volume of flammable gas the flow field will not be that of a bursting sphere as modeled by Brode(l 3) and Ricker<14) or a homogen- eous reaction as studied by Zajac and Oppenheim(lS) and Adamczyk<18). The flow field will develop from energy re- leased as a flame front propagates from the ignition source through the combustible mixture to the edge of the source volume. Because of the finite source volume and the finite time required for the flame front to propagate from the igni- tion source to the edge of the kernel, the explosion will be non-ideal. Combustion processes and non-steady one-dimensional flow in ducts were investigated by Rudinger(l 9). Assuming the chemical reaction takes place instantaneously as the unburned gas passes through an advancing flame front and the burning velocity is directly proportional to the abso- lute temperature of the unburned gas, he used the method of characteristics to calculate ·the properties of flame fronts with moderate, high,· and detonative flame velocities. The conservation equations were reduced to a manageable 19 form by omitting terms of small magnitude. Flow variables were then assumed to be uniformly distributed over any sec- tion of the duct leaving only time and one space coordinate as independent variables. The propagation of gas particles and pressure waves were then followed graphically in a coordinate system of these two variables on a plot called a wave diagram . Although this solution was strictly for one- dimensional flow, it led to the study of more complex flow fields . A self-similar solution for evaluating the structure of blast waves was developed by Oppenheim( 20), et al. The blast wave was assumed geometrically symmetrical and non- steady . The solution is in terms of two dimensionless independent variables, radius, R, and time, ~. The blast waves were examined in respect to two parameters, one des- cribing the front velocity and the other the variation of the density irrnnediately ahead of the front. The evolution of pressure waves generated by steady flame propagating in an unbounded atmosphere with planar, cylindrical, and spherical geometry was studied by Kuhl Kamel, and Oppenheim( 2l). They considered a self-similar flow field with both the deflagration and shock front pro- pagating at constant velocity and constant gas dynamic para- meters along lines of similarity Y = r/rs . They introduced reduced blast wave parameters as phase-plane coordinates and determined the appropriate integral curves on this plane . A numerical solution for the case of a hydrocarbon-air 20 mixture was developed which showed that the transition be- tween the blast wave solution and the acoustic solution is continuous . Pressure curves were generated as a function of deflagrative burning velocity for an expansion ratio, vf, equal to 7 . A simplified method for calculating blast parameters generated by a propagating deflagration was developed by Strehlow(ZZ) . Assuming that the pressure and density be- tween the shock and the flame is spatially constant, regard- less of geometry, the equations reduce to algebraic form allowing simple iterative solutions . Comparing his results with the exact self-similar solutions of Kuhl, et al . , Strehlow showed his results were identical for the case of planar flow when the pressure between the shock and flame are known constant . However, when the geometry changes to cylindrical or spherical the divergence of the flow field causes the pressure to decrease from the flame to the shock and the results varied from the exact solution but were with- in acceptable limits . E u Problem Definition The classical problem of ideal or point source explos- ions has been extensively studied by many investigators . Ideal blast wave theory is well understood and conveniently summarized by Baker(ZB) . Non-ideal explosions are not well under stood and many of the studies which have been done have not pr ovided com- plete answers to the questions of interest . The solutions 21 of Kuhl, et al o are limited by the self-similar assumption which applies only during the energy addition. There is no solution for the structure of the blast wave after the energy addition o Fishburn only investigated selected cases o There- fore his work did not show any trends. A systematic study of all the parameter affecting the generation and propaga- tion of non-ideal blast waves is needed . In the investigation of non-ideal explosions there are many parameters which affect the structure of the blast wave flow field o These parameters include the energy density of the source volume, the energy deposition time, the heat capacity ratio of the source volume and the surroundings, the flame velocity, and the flame thickness . By considering these parameters as planes or dimensions in an n-dimensional space a convenient tool for visualizing this investigation in relation to other studies is available. Figure 3 illustrates three of the dimensions investigated: (1) Energy density (2) Energy deposition time (3) Flame velocity (Plotted as the reciprocal) An investigation of bursting spheres (infinitely fast energy wave with instantaneous deposition time) was performed by Ricker(l4 ) . His studies are located at various energy densities on the bursting sphere line in figure 3 . Adamczyk(lS) expanded on the studies of Ricker . Add- ing energy uniformly throughout the source volume (infinite velocity, infinitely thick wave) he varied the energy density 22 7' it~~ 'Bursting Sphere Line Infinitely Thick Wave 100% ~kness 50% ----t~ 1/u Reciprocal Flame Velocity Infinitely Thin Wave 30% 20% 10% 5% Figure 3. Three dimensional Diagram of three parameters affecting Non-ideal Blast Wave Behavior. (Not To Scale) 23 and energy deposition time, over two orders of magnitude. This dissertation is part of a systematic study of the parameters affecting non-ideal explosions. In it the inves- tigations of Ricker and Adamczyk are expanded into a third dimension, a study of the effects of a constant velocity flame propagating from the origin to the edge of the source volume . The investigation was done using the energy density of natural gas, a common fuel. Cases were systematically run at selected velocities and the results were then com- pared to the homogeneous energy addition and the common lim- it case of bursting sphere. II. THEORETICAL CONSIDERATIONS A. Governing Equations Blast waves in air are non-steady flow fields propagat- ing through a compressible fluid medium bounded by a gas dynamic discontinuity. To predict the effects of propagat- ing blast waves it is essential to know the time history of the flow field properties at all locations within the med- ium. These properties are determined by the fundamental laws of nature appl~ed to fluid flow. Air, at or near standard temperature and pressure, is considered to be an inviscid fluid. Shock waves that appear in the flow can be treated as discontinuities or by using an artificial viscosity tech- nique. With these conditions the fundamental conservation equations can be expressed as: * + 'i/·(pV) = 0 ~ + V-.nv- =-(V .p)+ f at v p E Ci i i v2 a [ p (e+2 ) l v2 _ at + 'il· [p(e+z-)Vl = (Mass) II-1 (Momentum) II-2 -'v· (pV)+pQ+pEcifiVi (Energy) II-3 · h d · V 1.·s the flow veloc1.'ty vector, Pis where p 1.s t e ens1.ty, 0 energy' Q is the heat the fluid pressure, e is the internal addition rate per unit mass, ci is the mass concentration of P ·es 1.· and f 1.'s the body force act1.'ng on species i. s ec1. , i 24 25 Asswh ing an adiabatic inviscid fluid with no body forces, 0 the terms involving f and q become zero. Applying the ther- modynamic equation of state, p v = mR 6 with: o e e = Ee. (e. + f 0 cv_d e ) il. l. l. II-4 0 where e. is the energy of formation and c is the constant l. V. l. volume heat capacity of species i, internal energy can be linked to temperature, e, and density. There are then four equations to solve for the four prime variables of interest; u (local flow velocity), p (density), p (pressure), and e (internal energy per unit mass) . For simplificatio~ it is desirable to model the actual reactive fluid using a working fluid heat addition model. For a flow process the basic thermodynamic quantity is the enthalpy, h, explicitly defined by: h = e + p v II-5 Enthalpy is used rather than internal energy, e, and equation II-4 is replaced by h = l: c.h . II-6 i l. l. e 0 where hi = J cp. d0 + (thf) i II-7 0 l. 26 An actual flame process is an adiabatic process with no heat transfer to or from the system. However, large tempera- ture changes occur within the system as a result of chemical reactions. If the temperature is held constant during an exothermic chemical reaction heat must be removed from the system. The product enthalpy is then much less than the reactant enthalpy and the difference is ~h, the heat of reaction which was removed from the system. Since the system being modeled is an adiabatic system, the heat of reaction will not be re- moved but will become part of the system. Energy is con- served because the differing bond energies of the different molecules that appear or disappear lead to changes in the thermal energy of the system. With these observations the chemical reactions of the system can be replaced by a simple heat addition to a working fluid. Assuming: 0 h3 = h' = I C d0 3 0 P3 II-8 0 h4 = h4 + ~hf = I CP4 d0 + A 0 II-9 where h 3 and h4 represent the enthalpy before and after heat addition respectively, and a positive value of A represents heat addition to the flow. The derivation of the full equa- tions can be found in many texts on combustion, e.g. Williams< 23 ) and Strehlow< 24 ). 27 B. STEADY ONE-DIMENSIONAL FLOW DISCONTINUITY RELATIONSHIPS Although details on steady one-dimensional flow discon- tinuities are available in most text books on combustion it is desirable to proceed with a brief review of basic princi- ples and concepts for comparisons with non-steady behavior which are to be made in succeeding sections. Using the heat addition working fluid model there are four equations which, because of their complexity cannot be solved without certain assumptions and restriction. For the case under consideration, blast waves, a convenient simplification is that the shock in the blast wave can be approximated as being a one-dimensional phenomenon. Shock waves are extremely thin and fluid proper- ties across the shock adjust within a few mean free path lengths. Thus in the scale under consideration the curvature of the shock approaches that of a planar wave and the one- dimensional relationships apply for the shock in plane-, line-, and point symmetrical blast waves. The basic non-steady, one-dimensional conservation· equa- tions of fluid dynamics can then be expressed as: (Mass) II-10 ' 2. . ('-1) a (purJ) + L(pu rJ+prJ)-jpr J = 0 at ar (Momentum) II-11 . 2 . a . 2 3[purJ(e+! )+purJ] at[prJ(e+z )] + ar = O (Energy) II-12 where e - C 0 = ~ " -y-T (State) II-13 28 and j = 0, 1, and 2 for planar, cylindrical, and spherical symmetry respectively. Because shock waves are so thin the shock wave in blast wave structure can also be approximated as being quasi-steady. The equations of fluid dynamics can then be solved for the case of one-dimensional, constant area, inviscid flow to yield what are generally called the normal shock equations, i.e. the conditions for transition across a shock wave with heat addition: plul = P4U4 (Mass) II-14 2 2 P1 + plul = P4 + P4U4 (Momentum) II-15 I 2 I 2 h1 + u1 /2 = h4 + U4 /2 + >.. (Energy) II-16 Hugoniot Substituting the mass and momentum equation into the energy equation yields the Rankine Hugoniot equation. I hl + >.. = \(P4 - pl)(vl + v4) II-17 With the enthalpy relationship, h = Cp0 =~~~and the equation of state this becomes: II-18 which represents the locus of final states, p4v4 for any 29 initial conditions of p1 and v1 with heat addition A. For the case of no chemical reaction A is zero, y is assumed constant, and equation II-18 becomes the shock Hugo- niot, i.e. the locus of all possible solutions for normal shocks without chemical reactions for one set of upstream con- ditions, p 1 and v 1 . II-19 By algebraically manipulating the shock Hugoniot it can be shown that it will asymptotically approach Pz and v2 as v2 and Pz respectively approach infinity: (Pz) (r-1\ p·l + - \r+1) (~~)+ (f+t) as Vz + 00 as Pz + oo II-20 II-21 Figure 4 is a plot of the shock Hugoniot. However it is physically known that the situation of pressure decrease a- cross a shock wave does not exist. Therefore, in actuality the only physically real solution is the shock Hugoniot for increasing pressure and decreasing specific volume. This can be proven by an examination of the entropy change or by attempting to plot a discontinuous expansion for the shock Hugoniot by the Method of Characteristics. For the case of heat addition to a constant garrrrna, ideal gas working fluid, Strehlow(Z4 ) determined that the reacted end state Hugoniot can be represented by a rectangular hyperbola l p 30 Strong detonation Excluded Weak detonation (supersonic combustion) Shock Hugoniot v--.... Figure 4. Pressure-volume plot of end states for a one-dimensional steady process with heat addition. 31 in the p-v plane. Zajac and Oppenheim(ZS) have shown that this type of hyperbola accurately represents the shape of the real gas hugoniot. The assymptotes of the Reactive Hugoniot are: as v4 -+ 00 II-22 II-23 Two points for plotting the Reactive Hugoniot can be calculated by asstnning a constant pressure expansion and a constant voltnne pressure rise: II-24 II-25 Thus for the reactive Hugoniot the values of p4 and v4 can be determined for plotting the curve. Since isotherms are hyper- bolas that asymptotically approach the p=o, v=o axis, the Hugoniot curves always cross the isotherms such that increas- ing p along a Hugoniot increases temperature. Since the temperature hyperbola asymptotes the axis, p = 0 represents a value of 0 = 0. This represents the hypothetical but im- possible case where all the random kinetic energy of the molecules has been converted to ordered flow velocity. 32 Rayleigh Line Manipulating further the normal shock equations by sub- stituting the mass equation into the momentum equation , an- other relationship between the pressure and specific volume can be developed, the Rayleigh Line. II-26 The equation for the Rayleigh line specifies that the approach mass flow rate squared, (p1u1) 2 , is equal to the negative slope of a line in the p-v plane connecting the initial and final states of the process under consideration. Thus if the initial conditions of p1 ,v1 , and u1 , are known, the final conditions p4 and v4 can be determined by drawing a 2 line through the initial point with a slope equal to -(plul) This straight line intersects the Hugoniot at the final state. The Rayleigh line defines an important characteristic of steady state flow; since the density, p1 , and the flow velocity, ul, are squared, their side of the equation will always be positive. Therefore, the slope of the steady state Rayleigh line must always be negative. Thus for steady state, one- dimensional flow, certain areas of the p-v plane are excluded for final end states as illustrated by figure 4. Equation II-26 can be rewritten in terms of the flow Mach number of the Rayleigh line process both ahead of the shock wave and behind the energy addition: 33 M 2 . . 1 (p4/Pi-) - l II-27 = 1 Y1 (v4) 1 - - vl (~;) (:~) - 1 and M 2 II-28 = 4 Y4 (:~) (~~) 1 - When exothermic chemical reactions occur in a steady-state flow situation the Rayleigh line may intersect the Hugoniot in one, two, or no locations for both supersonic and subsonic incident flow velocity. Above a limiting subsonic velocity and below a limiting supersonic velocity the Rayleigh line does not intersect the Reactive Hugoniot and there are no possi- ble steady-state solutions. Below the limiting subsonic velocity and above the limiting supersonic velocity the Rayleigh line intersects the Reactive Hugoniot twice, indi- cating two possible end states for both subsonic and super- sonic velocities. For very low subsonic veloci t i e s the Ray- leigh line can intersect the Reactive Hugoniot only once because the Hugoniot enters the imaginary region of negative pressure. The very low velocity subsonic solution corresponds to normal flame propagation. Physically, laminar flames are represented by this solution . For ordinary flames the fl ame velocity is very low and therefore there is only a ver y slight pressure drop across the wave . 34 At the tangency point of the Reactive Hugoniot and the Rayleigh line there exists only one propagation velocity . This velocity corresponds to exactly sonic velocity at station 4, and is called Chapman-Jouguet flow or CJ flow . The upper CJ point represents the proper end state for detona- tions. The existence of exactly sonic flow at the tangency point can be shown by differentiating the Hugoniot and equat- ing this to the slope of the Rayleigh line d(:~) dC~) Hugoniot dC~) Rayleigh (~)-1 C~)-1 II-29 II-30 35 The condition for tangency is then: (:~) = II-31 which can be rearranged to: II-32 This is identical to Equation II-28, the equation for the Mach number of the Rayleigh line process at point 4, proving that the Mach number behind the shock at the upper and lower CJ points is sonic. For a strong detonation or a weak deflagration: II-33 and for a weak detonation or strong deflagration: II-34 36 Thus M4 < 1 for a strong detonation or weak deflagration and M4 > 1 for a weak detonation or strong deflagration. Investigating further the characteristics of the upper CJ point, the Hugoniot equation and the Raleigh line can be combined to determine an explicit relationship for the pres- sure and volumetric ratio ahead of the shock and behind the energy addition: = II-35 II-36 The CJ point is the tangency point of the Rayleigh line and the Hugoniot curve. For this point there exists only one solution. Therefore, the expression under the radical sign must equal zero at the tangency point and can be expressed as follows: II-37 37 Knowing A, the approach flow Mach number for the Chapman- Jouguet points can be evaluated: II-38 The pressure and specific volume can also be calculated from equations II-35 and II-36. At the CJ points the quanti- ties within the radical signs of the equations become zero and the equations reduce to : II-39 y 4 2 (½)(ylMCJ+l) (y 4 + l):1ci II-40 38 Even though one-dimensional steady-state heat addition is impossible over velocities which lie between the lower and upper CJ points , this investigation included the addition of energy at these forbidden velocities. This is possible since the calculation is fully non-steady~ Therefore, the flow will follow a solution in accordance with the non-steady equations of mass, momentum, and energy, and will not be re- stricted by one-dimensional steady flow considerations. C. ENERGY SCALING Classical blast studies have been primarily directed to an investigation of the blast waves generated by either high explosives or nuclear weapons. When conducted on a large scale, experimental studies of blast waves are dangerous, expensive, and difficult to control. Large isolated areas are required where access and egress may be closely monitored and controlled to ensure the tests are conducted safely. In addition, the res~lts are subject to the ~ffects of atmospheric and topographical conditions which make the interpretation of data difficult and subject to error. The cost and other problems associated with large scale tests make their use in a systematic study of flow field behavior prohibitive. Energy scaling is a tool which has been used extensively in the comparison and extrapolation of the results of tests involving different quantities and composition materials depositing energy in a source volume. The two most widely used methods of energy scaling involve Hopkinson's scaling law and Sachs' scaling law. 39 Hopkinson or "cube root" scaling is commonly used. This scaling, first formulated by Hopkinson<26 ) states that self- similar blast waves are produced when two similar explosive charges with characteristic dimensions varying by a length scaling factor, a, are detonated in the same atmosphere, an observer whose location from the scaled explosive is a times the dis t ance from the standard, will feel a blast wave of sim- ilar form with amplitude P, duration crt, and impulse , cr l. All characteristic times will be scaled by the same factor as the length scale factor, a. Pressure, temperature, densi- ties, and velocitie s are unchanged at homologous times. Hop- kinson's scaling law requires that the model and prototype energy sources be of similar geometry and the same type of explosive or energy source. A more complete discussion of this scaling is available in Baker(Z). A more general blast scaling law than Hopkinson's was developed by Sachs to account for changes in ambient condi- tions and the effects of altitude. Sachs developed dimen- sionless groups that involve pressure, impulse, time, and ambient parameters as unique functions of the dimensionless distance parameter: tap 1/3) 0 0 E l/3 t II-41 The Sachs' law identifies the blast source only by its total energy, Et, and therefore is not restricted to sj milar 40 geometry and explosive type as Hopkinson's law. However, it would not be expected to be consistent · for scaling of close- in (near field) effects of non-ideal explosions. Although the short comings in the use of these scaling parameters are obvious, they provide a convenient tool for comparing and analyzing theoretical and experimental data. D. DAMAGE EQUIVALENCE The concept of equivalence between non-ideal explosions is not fully understood. With equivalent far field over- pressures, the near field behavior of non-ideal explosions may vary greatly. A means is needed to evaluate the effec- tiveness for blast damage of any particular accidental ex- plosion and how this effectiveness varies with parameters affecting the development of the blast wave. The conn:non procedure in an actual accident is to ob- serve the blast damage pattern to determine the weight of TNT (tri-nitro-toluene) required to develop blast wave over- pressures to do similar damage at the same distance from the explosion center<27 ). Next, the maximum equivalent TNT weight of the fuel or chemical is determined by calculating either the heat of reaction of the mixture or the heat of combustion of the substanced released. The mass equivalence of TNT is expressed as: (WTNT) Equivalent = 6H *m · C C 4.198* 106 where 6Hc is the heat of combustion of the hydrocarbon II-42 41 (cal/kgm), me is the total mass of the reactive mixture (kg), and 4.198 X 106 is the heat of explosion of TNT (joules). The common expression "per cent TNT equivalence" has been developed for comparison with data available from the test- ing of TNT and is determined by: %TNT = ( (WTNT) / (WTNT) . 1 *100. II-43 damage equivalent In an actual hydrocarbon explosion the damage as a function of scaled distance does not agree with that pre- dicted from TNT equivalence. High explosives, such as TNT, contain internally much of the oxygen need for chemical reactions. Once initiated, the explosion proceeds almost instantaneously to completion. Hydrocarbons, on the other hand, must react with the oxygen in the air,making mixing an important parameter. A finite time is required for the flame to propagate through the combustible mixture influencing the development of the blast wave. Also, the calculated heat of combuS t ion is based on reactions to an equilibrium concentration of carbon dioxide and water. In actuality the reaction is not carried to equilibrium and at elevated temperatures the molecules may begin to dissociate, thereby further altering the effec- tive heat release. E. DAMAGE MECHANISMS In the flow field associated with a blast wave there will be transient overpressures and wind induced drag forces. 42 The damage and injuries sustained by people, buildings, animals, and vegetation will vary, depending on the pressure- time history of the blast wave. Large overpressure of short duration may cause ear damage with little physical displace- ment of the body, whereas lower overpressure of longer-dura- tion may cause lung damage and other severe body injuries. Similarly buildings may be constructed to resist overpressure of short duration, but may fail from the impulsive drag associated with lower overpressures of longer duration. Damage and injuries are not restricted to the peak overpressure or impulsive drag alone, but to the combination and interaction of these effects. The exact relationships are quite complex, but a convenient simplification to corre- late blast wave properties to damage effects on a wide variety of targets has been discussed by Baker, et al. (28 ). He states that for any object, levels of constant damage of one type can be plotted on a pressure-impulse (P-I) Diagram, or empirical or analytical equations can be developed to describe the pressure-impulse (P-I) relationship. An example is shown in figure 5. To illustrate this concept, he considered the spring- mass system illustrated in figure 6 and subjected it to a specific time varying force to represent the dynamic res- ponse of a structure. The equations for a curve represent- ing the combinations of scaled force and scaled impulse which cause the same scaled response X of the system were max determined to be: 100. 10. P= 2P xmax K , . ., ., Force Asymptote \ Impulse Asymptote/ 1. Impulsive Loading Realm 43 10. I= [X (KM)½] max Quasi-Static Loading Realm 100. Figure 5. Scaled P-1 Curve for Fixed Level of Damage. p (t)--- M Force p 44 p (t) = p -t/T e I = I: P ( t) dt = PT 0 .__ ___________ _ Time 0 Figure 6. Schematic Diagram of Spring-Mass System to Model the Dynamic Response of a Structural Member. p = I = 2P X K max I 45 2 = ------..--:::---- ---- [ 2- exp (-w2T2 /100) J tanh wT = wT 1 X (KM) ~ max [2-exp( -w 2T2/100)Jtanh wT II-44 II-45 where X is the maximum displacement of the system, K is max the spring rate, mis the mass, w is the natural frequency of the system, and Tis the characteristic loading time. By varying wT in these equations, a scaled response curve or Pressure Impulse (P-I) curve can be determined, similar to the curve in Figure 5. This curve represents the combinations of scaled force and scaled impulse which cause the same scaled response X of the system. This iso-max response curve can be compared to an iso-damage curve of a building or similar structure. For a given structure vary- ing levels of damage can be determined as functions of the pressure and impulse the structure is subject to. Predic- tions can then be made of the level of damage which the building would suffer based on the predicted pressure and impulse of the flow field associated with the blast wave. The causes of damage can be separated into regions on the iso-damage curve, the impulsive losding realm in which overpressure is controlling, the quasi-static loading realm in which impulse is controlling, and the dynamic loading realm in which the combination of overpressure and impulse determine the damage. This technique has generality because once the pressure and impulse are known for any explosion, whether it is ideal 46 or non-ideal, the P-I technique can be used to evaluate dam- age at any location. Sachs scaling and other methods of scaling do not have the flexibility of the P-I technique since they only relate pressure and impulse for high ex- plosive and point source explosions. The P-I technique is a very general technique and more useful for accidental explosions than energy scaling or the TNT equivalency argument. III. COMPUTATIONAL PROCEDURE The computational techniques used are based on the Von Neuman-Richtmyer concept of artificial viscosity as devel- oped by Brode(l3) and Wilkins< 29 ). Using this technique Professor A.K. Oppenheim( 30) of the University of California, Berkley, developed a computer program for studying the flow field of blast waves. The program is written for a one- dimensional, non-steady flow field in planar, cylindrical, and spherical geometry. The system is idealized with several simplifying assumptions: (1) The system is symmetrically one-dimensional. (2) The high energy source volume is separated from the surroundings by a massless barrier and there is no transfer of mass between the high energy gas and the surround- ings. (3) The flow is inviscid with shock wave formation the only dissipative process in the surrounding atmosphere. The computer program was modified by Adamczyk(lB) of the University of Illinois to allow heat addition along particle paths by incorporating a homogeneous energy addi- tion term with temporal and spatial dependence. The program was further modified by the author to incorporate a wave energy addition term and variable gannna, both with temporal and spatial dependence. In the computer program the conser- 47 48 vation equations are expressed in Lagrangian coordinates, since through their inherent conservation of mass they lend themselves more easily to a computational scheme. Partial derivatives are taken along a particle path such that u = ¾f and the equations of mass, momentum, energy, and state are; Mass III-1 d\) \) a (rju) = IT ~ ar Momentum au ~ IT = -\) ar III-2 Energy ae av " IT = -p IT - III-3 e - ~ y-1 State III-4 where vis the specific volume, r is the radial position, j is the geometry coefficient (0, 1, 2 for planar, cylindrical, and spherical flow fields, respectively), pis the pressure, e is the internal energy, "is the heat addition term assumed in the heat addition model, with spatial and temporal depen- dence, and y is the ratio of specific heats, also with spatial and temporal dependence. The f h fl fl.·eld and their variation properties o t e ow Wl..th the 1.·ntegration of the governing time are determined by . h equati'on of state, and the kine- conservation equations, t e matic equation coupled with the energy source term, A, and conditions at r=O, t=O, subject to the appropriate boundary and ahead of the lead wave. E' A. . -- - .-:--':-': -.. . - - - - ~ ,..;,n:"t ~ _...... - ,, 49 Boundary Condition~ For the cases studied the boundary conditions are, 1. At t = 0 and O~r~ 00 u = u(r,o) = 0 p = p(r,o) = p 0 e = e(r,o) = eo \) == v(r,o) = \) 0 2. At r=o and u = u(o,t) = 0 ~ = (~) = 0 yr yr (o,t) ae (ye) = 0 ar = ar (o,t) ~= (~) = 0 ar ar (o' t) 3. Ahead of the lead wave. p = Po \) = \)0 III-Sa III-Sb III-Sc III-Sd 1Ir-6a III-6b III-6c 1Ir-6d III-7a III-7b III-7c III-7d so B. Dimensionless Variables To aid in the computations, all variables are non-dimen- sionalized with respect to the thermodynamic state of the at- mosphere into which the front propagates, mR00 =p0 v 0 , and a reference point at the edge of the energy source volume. The non-dimensional independent variables are defined as: n = r/r 0 t = t/t0 111-8 111-9 where t 0 is a characteristic time proportional to the time it takes an acoustic signal to propagate from the origin to the kernel edge when traveling at the ambient undisturbed sound speed, a 0 , and r 0 is the outermost edge of the source volume at a time t = t = O· r /y t = o o 111-10 o a 0 Using P 0 to represent the ambient atmospheric pressure, v 0 the b . 1 nd a the ambient am ient value of the specific vo ume, a o sp d f d d d t variables can ee o soun, the non-dimensional epen en be expressed as: 111-11 u = = lil-12 111-13 P = p/p (for equation of state) 0 P* = p/p -IT(for conservation equations)Ill-14 0 111-15 A= ~/p v· 0 0 111-16 51 In non-dimensional form the conservation equations are: Mass o\J! - 1jJ a (njU) fi - --r an nJ III-17 Momentum au_ 1jJ aP fi - - ar, III-18 Energy aE _ -P ~+A a-r - aT III-19 E = P\j! (y-1) State III-20 u = an fi where III-21 and the boundary conditions become: 1. At T = 0 and o~n~00 III-22a U(n,o) = 0.0 III-22b P (n, O) = 1.0 III-22c E(n,O) = 2.5 III-22d \J! (n, O) = 1.0 2. At n=O and 05T500 III-23a U(O,T) = 0.0 ap 1 . 0 III-30c Po vo where: s(MD,T) = M~2 T-D w 58 and " ~ /cos(3Tc)-9.0 cos (n o+s.0(/16.0 III-31 The three-dimensional shape of the energy addition function is shown in figure 7. At time r=O the system exhibits am- bient conditions throughout. + At time T , energy addition begins at the center of the kernel in accordance with the energy source term until s=l., when all the energy has been added. At positions of increasing radius the start of the energy addition begins at later times in accordance with dD MD= cfr· The energy addition is done in the energy wave in accor- dance with a selected wave width which can be varied to model the width of the flame. In this model the wave width, W, is the fraction of the source volume to which energy is being added at any time step as shown in figure 8: III-32 This can also be visualized as the fraction of the transit time for the wave to propagate through the source volume, TT, that the energy is being added to a particular cell, Tc, and can also be expressed as: III-33 -C: :, --~ .. c,.: C: w 59 D1 D0 Lagrangian Radius, D Figure 7. Energy Deposition Function . l 60 I I I I I I I I I I I I I I I ~ I I I I I I I I I I I I I I I I I I 0 1,C...__--L. ___ __::.::..._ _______ ....L,___L_ ________ _ 0 D--- Figure 8. Wave Width of Energy Deposition Term. 61 The energy wave propagates at a constant Lagrangian velocity or Mach number,~. where: dD = dT III-34 The transit time of the energy wave through the source volume is inversely proportional to the velocity of the energy wave. For equal wave widths, as the velocity increases both the source volume transit time, TT' and the cell deposition time, Tc, decrease. Figure 9 shows the effects of wave width on cell deposi- tion time. As the wave width increases, the cell deposition time increases for the same energy wave velocity. The source volume deposition time, TD, is the sum of the transit time of the energy wave plus the cell deposition time at the edge of the source volume: III-35 This can also be expressed in terms of the energy wave Mach number: III-36 For an infinitely thin wave W=O and the source volume deposi- tion time equal the energy wave transit time. As the width of the energy wave becomes finite, energy is being added to a, E I= C: 0 ·.:: ·;;; 0 0. ., 0 .. 4.0 3.0 2.0 1.0 .8 .6 .4 Infinitely Thick Wave I I I I 2.4/ I I I I I I I I I 0.5 1.0 Bursting Sphere (r=O, 1/Mw=O) I I 62 I I I I I I / 1.~ }' I / I / .y / I 0/ I I 0.1 2.0 3.0 4.0 1/Mw (Inverse MACH Number) Figure 9. Deposition time vs. Inverse MACH Number as function of Wave Width. / / 63 the last cell after the leading edge of the energy wave reaches the edge of the source volume.· Figure 9 shows that the greater the width of the energy wave the longer the source volume deposition time . 2. Change of Specific Heat Ratio The ratio of specific heats, gamma, for a combustible mixture is known to vary from approximately 1.1 to 1.67 depending on the composition of the mixture and the complex- ity of the molecules in the individual components of the mixture. In addition, the value of gamma can also change as the chemical composition changes to maintain chemical equili- brium or as a function of temperature. As temperature in- creases, the species in air go through various changes including the dissociation and ionization of oxygen and nitrogen. At a temperature of 2500°K the dissociation of oxygen molecules begins. For the combustible mixture being investigated the temperature ratio is 9:1 which corresponds to a temperature of 2700°K behind the energy addition. Thus for the case under consideration the dissociation of oxygen begins, raising the heat capacity and lowering the heat capacity ratio, gamma. An evaluation of an effective value of garrnna and heat release associated with real combustion processes as a func- tion of stoichiometry was performed by Adamczyk(lB). For the case which is being investigated a combustible mixture with an energy density approximating that of me t hane is used. For methane, Adamczyk calculated an effective gamma of 1.202, 64 rounded off to 1 . 200 here. Before a vapor cloud is ignited, uniform ambient conditions exist throughout both the source volume and the surroundings. After ignition the flame front heats the medium through which it propagates and changes the chemical composition, lowering the heat capacity ratio. To model this change a variable gamma was developed in which the heat capacity ratio changes from an ambient condition of 1.4 to 1.2 when energy addition to the cell is completed. A. y = yo - (yo-Y4)(~) III-37 where Ai/A is the fraction of the energy which has been added. D. Numerical Integration The numerical integration was done using a Von Neumann- Richtmyer/type, explicit, finite differencing technique. The equations of motion were integrated for an expanding flow field with constant Lagrangian distance spacing at finite times. The time steps were determined using the Courant Stability criteria as presented by Wilkins(29 ). 1:::, n+½ = . (An+½ 1 4 An-½) T . min LlTR •• uT III-38 where III-39 and: ~(~: -~:: 2 ~ min over all i's 1 2 J III-40 where: z2 = 1 2 64C2 65 - 1/J~~t)2 (n~+½ - n~)2 n-1 n-½ + 1/J·+1 !JT 1. Yi 2 2 = Yn Pn ,,,n 2 i+½ i+½'t'i+½ III-41 III-42 III-43 The computational grid for the finite differencing scheme is shown in figure 10. Velocity is evaluated at full steps in radius, cell boundaries, and half steps in time to maintain the proper relationship between the derivatives as demanded by the conservation equations . Thermodynamic properties, P, 'l/ , and E are evaluated at full steps in time and half steps in radius. Since TI is a relationship between the velocity and effective pressure, it is evaluated at both half steps in space and time. The sequence by which the equations are treated is first the momentum equation, followed by the kinematic equation, continuity equation, and energy equation. Using the nomenclature in figure 10, the conservation equa- tions were written in finite difference form as follows: Momentum Equation n+!.: u. 2 = l. u:1-½ l. III-44 III-45 t T (U) (P) (U) (P) (U) n+1 pn+l n+l i-½ pi+½ (P) n+ ½ Un+½ u~+½ Un+½ i-1 + I + i+l (U) n n n P. ½ pi+½ ,- 2 (P) n-½ n-½ n-½ n-½ u. 1 + u . + ui+l ,- I (U) n-1 n-1 n-1 P. ½ pi+½ ,- (P) i-1 i-½ i+½ i+l o- Figure 10. Computational Grid for finite differencing technique. 67 n n nr:1 n ni+l - n. - n. 1 I k 1 + 1 1-= III-46 2 n n tµi +1-1 tµi-½ and IT is normally zero except in regions of excessive pressure gradients (shock waves), in which case: n-½ = IT . +1 1 -~ [ 1 . 1) 2 ) ] un-~ - un-~ + c2 ~ i+l i (- 1- + 1 o n -:-n=T lJJ 1· +k lJJ1· +k - 2 - 2 III-47 Since the artificial vicosity is required only to smooth out the effects of excessive pressure gradients the condition is introduced that if n tµi ±½ ~ n-1 lj)i ±½ III-48 or LiU ~ 0 III-49 n-1 0 III-50 IT . +1 = 1 -~ KINEMATIC EQUATION III-51 CONTINUITY EQUATION . . n+l n n+½ n+½ J n-½ n+½ J n+ 1 n { T - T } { U i + 1 ( n i+l) -U i ( n i . ) + ~) lJJi+½ = lJJi+½ + n M. +1 J. ~ III-52 where: and: n M-+1 = J_ ~ 68 ·+1 ·+1 n J n J (ni+l) · - .(ni) · (j+l)ij;r:1+1 J_ ~ III-53 3 [U1;+1-1 J ) I II-54 J_ where j is equal to 0, 1, or 2, for planar, cylindrical and spherical flow fields respectively. ENERGY EQUATION and En+l = i+½ En+l i+½ = pn+l i+½ n+l Yi+½ III-55 n+l ij;i+½ III-56 - 1 n+J:a where A.+/ is the energy addition term and y is the local ratio J_ ~ of specific heats. E. Testing of Program To establish credibility of results and ensure that the computer program effectively models the system under 69 evaluation test cases were run in which expected results were available to compare to computer output. 1. BURSTING PLANE To test the computational technique of the program the case of one-dimensional, constant area flow similar to a membrane bursting in a shock tube was run. The initial con- ditions of a high temperature, high pressure constant gamma gas with a step change to ambient conditions at the membrane were used. The calculated results were then compared to result predicted by equation I-10. The results varied by less than 0.01%, establishing the validity of the calculation technique used. 2. BURSTING SPHERE An infinite velocity energy wave propagating through the compressible fluid medium is a constant volume energy addition or bursting sphere. To test this case on the model, a wave velocity was selected at the maximum velocity which could be incorporated into the program, limited by the initial step size (this corresponds to a dimensionless Mach number, MQ = 4225). Figure 11 is a pressure vs radius plot of the energy wave at dimensionless time increments of 0.0001. After the wave has propagated through the source volume, the pres- sure-radius distribution is a bursting sphere. The wave addition of energy yields a pressure difference of less than 0.001% from the energy distribution for a bursting sphere, but imparts a velocity to the particles of approximately l.6xl0- 3 . These differences are considered well within the allowances of 10.0 7.0 p 4.0 1.0 0.0031 0.0024 0.0016 "Z r, r (r ( I 0.0008 0 .0001 2.00 1 .67 1.33 1.00 0.67 0. '33 PRESSURE / PO DISTRIBUTION VS. 015TANCE / DO AND TlME / TO Figure 11. Pres sure dis tribution ver sus Eulerian dis t ance and time f r om a Mach 4225.0 energy wave. 0 . 00 71 the problem under consideration. 2 . Wave Width In a flame the heat of reaction does· not appear instan- taneously but is controlled by the reaction rate of the chem- ical species. A flame propagating through a flammable mixture will have a finite time of deposition of energy to the individual particles as it passes. Therefore it is nec- essary to model the energy addition in the energy wave by adding the energy simultaneously over several cells. The wave width determines the number of cells to which energy is added. In addition, the stability criteria used in determining the time increment relates the time step size used in the calculations to the energy being added to the cells. If the wave width limits energy addition to only one cell at a time, each cell would require a complete time cycle of energy addi- tions and the energy addition would be effectively a series of explosions. If energy is added simultaneously to several cells the time step size is limited only by the most restric- tive energy addition step. Thus, with energy addition simul- taneously in several cells computer time is reduced in proportion to the number of cells within the energy addition wave. The wave width also affects the deposition time of energy addition to each cell. Figure 9 shows that as the wave width increases the time for energy deposition within the individual cell also increase. A series of cases were run at a supersonic energy wave 72 velocity of Mach 4 and a subsonic energy wave velocity of Mach 0.5 to investigate the ·effects of wave width on the model. For the supersonic case (Mach 4) Figure 12 illustrates the effects of wave width on peak overpressure. During the energy addition there are significant fluctuations and dif- ferences in overpressure as the wave propagates through the source volume. For a wave width of 0.2 the energy is added to ten cells simultaneously and as the final energy is added to the last cell in the wave there has been some pressure transfer to adjoining cells during the relatively long deposi- tion time. As the wave width decreases the number of cells in which energy is being added decreases with an accompanying decrease in the cell deposition time. Since the energy is added rapidly the increase in energy of the cell is reflected in a pressure rise with very little pressure transferre.d :to adjoining cells. Also, in the finite differencing scheme all the cell properties are assumed to be concentrated at the cell center. For a narrow wave propagating through the kernel, i.e. containing 1 or 2 cells, the finite differencing scheme may result in large pressure and energy variations in adjoining cells because after energy addition is completed in one cell the energy addition in the adjoining cell may be only starting. During the time of energy addition to the new cell the energy (pressure) in the old cell will be trans- ferred to adjoining cells. Thus there may be successive peaking of the pressure in the cells caused by the wave 50.0 10.0 Lw a::: ::J en en w a::: [L a:: w > 0 1.0 0.4 0.0 0.1 73 DVERPRESSURE VS RRDIUS Q=8 r2J BAKER C) BURSTING SPHERE MAC~ WAVE WIDTH 0.2 + WAVE WIDTH 0. 1 X WAVE WIDTH 0.05 ~ WAVE WIDTH 0.025 f.o RADIUS (ENERGY SCALED) Figure 12. Overpressure versus energy scaled distance. 74 encompassing too few grid points as it propagates through the kernel. This peaking is reflected by the overpressure waves for wave widths of 0.025 and 0.05 (2.5 and 1.25 cells respectively). However, it should be noted that as the pressure wave propa- gates from the source, the peak overpressures coalesce into the same overpressure curve. This implies that one of the effects of wave width is the rate at which the non-steady flow assyrnptotically approaches a . maximum value of peak pressure during the energy addition . For the subsonic wave velocity, Mach 0.5, figure 13 shows similar results, except at much lower overpressures. For a wave width of 0.2 the fluctuations in overpressure are much smaller than the 0.1 and 0.05 case; but all these cases approach similar overpressures at the edge of the kernel. The narrow wave width (0.05) initially has fluctuations in the overpressure, but as the wave propagates to the edge of the kernel the subsonic velocity of the wave allows equaliza- tion of the pressure. Also, the time of energy deposition per cell for the 0.05 wave width at Mach 0.5 is 8 times longer than the Mach 4-0.05 wave width, and twice as long as the Mach 4 - 0.2 wave width. However, in the far field the 0.2 wave width shows a noticeably lower overpressure than the case of a 0.1 and 0.05 wave width. A wave width at 0.1 was chosen because : (1) The solutions assyrnptotically approach the peak value before the energy wave has propagated an excessive 50.0 10.0 L.1...J c.c: =:) en en L.1...J c.c: Q_ cc L.1...J > 0 1.0 0.4 .-- 1 I 11 0.04 0:1 75 ~VERFRESSURE VS RRDIUS O=B ~ B~KER ~ BURSTJNG SPHERE MRCH 0,5 ENERGY RDDJTTDN A WRVt WIDTH 0,2 + W~VE WJDTHO, 1 X W~VE WIDTH 0.05 \ I , 1.0 . AROIUS (ENEAG1 SCRLEDJ Figure 13. Overpressure versus energy scaled distance. 76 distance through the source volume, (2) In the subsonic case , the expansion of the source volume was approximately the same for the 0.1 and smaller wave widths, (3) The calculated results did not require ex- cessive computer time, and (4) This approximation reasonably modeled the physically realistic solution. IV. RESULTS AND DISCUSSIONS A. Flow Field Properties from One-Dimensional Steady State Theory. Fuels at stoichiometric concentrations have an energy density ranging from 5.8 for hydrogen to 9.6 for ethylene oxide. Using an energy density of q=8.0 (approximately that of methane, q=7.93), a Y4 of 1.2, and a y0 of 1.4, the shock Hugoniot and reactive Hugoniot can be plotted and the system constants calculated. From equation III-26: P4IP0 = q + 1 P+lpo = 9.0 IV-1 For a constant volume energy addition equation II-24 can be rearranged to the following: A = IV-2 A= (45. - 2.5) = 42.5 77 78 For a constant pressure energy addition equation II-25 states: and V 1 7.67 The approach flow Mach number for the Chapman-Jouget tangency point can be evaluated from equation II - 38: MCJ [ A (y i-1) (Y{Y42)} [(1+ A(y42-l) = l+ -Y1P1 vl Y1 CY1-l) . Y1P1V1 (Y/-Y42))2 - Y1(Y1-l) -( ~;)T r MCJ = 5.179 & 0.165 IV-3 IV-5 The steady-state, one-dimensional flow properties at the Chapmen-Jouguet points can be evaluated from equations II-39 and II-40. IV-6 79 (~:4) 1 CJ = ( .~4) = 0.56 & 14.78 1 CJ IV-7 These steady state predictions will be compared with the results generated by the non-steady heat addition model. B. The Effects of Energy Wave Velocity. In this analysis, Lagrangian constant velocity energy waves were varied over several orders of magnitude to ascer- tain the flow field properties of the propagating. wave sys- tem. These properties were then compared to those of burst- ing sphere. All cases were run with the same total energy and the variables are summarized in Table 2.:. 1. Flow Field Properties The flow fields of the cases investigated were plotted to illustrate the results. Figure l(is the Lagrangian pressure distribution as the energy wave reaches the edge of the source volume. Figures 15 through 27 show the Eulerian pressure distributions at various times. Figures 28 through 37 show the pressure - specific volume behavior of the indi- vidual particles. Figures 38 through 46 show the pressure versus time history at fixed Eulerian radius. Figures 47 through 55 show the displacement of the particles with time. Note: All figures in this chapter are collected at the end to simplify comparisons. 80 BURSTING SPHERE In case 1 (bursting sphere) there ·is initially a con- stant pressure of 9.0 within the energy source volume, decreasing at the edge to an ambient pressure of 1.0 in the surroundings. Figure 15 shows that following the instant of burst an expansion wave begins to propagate into the high pressure source volume and a shock wave develops, propagating away from the source volume. The expansion · wave propagates into the source volume at the local velocity of sound and reaches the center at a time of 0.257. The center of the sphere is a singularity point and the expansion wave reflects as another expansion wave. The pressure at the center drops to a minimum value of 0.0656 at T=0.625. The system attempts to equilibrate the pressure by returning the mass removed by the expansion wave. The system over compensates and at T=0.680 the pressure peaks at the center and is reflected as a shock wave. This wave behavior can be seen in the particle path plot of figure 47. The initial expansion wave exhibits itself by the outward movement of the particles. Since the conditions within the source volume are initially uniform, the local velocity of sound is uniform and a straight line can be drawn from the source volume edge to the center along the front of the expansion wave. As time progresses the source volume has over expanded and the particles reverse there out- ward movement. At T=0.680 the particle momentum reflects from the center as a shock wave. The second shock wave 81 progression can be seen by the inflections in the particle paths. The decreasing strength of the shock wave ·is shown by the decrease in the inflection of the particle paths as the wave propagates outward. This second shock tranfers mass away from the center generating another expansion wave. This expansion wave generates a third shock at T=l.85. If the calculations had been run to longer times, the reflection of expansion waves and shock waves from the center would have continued, but figure 47 shows that successive shocks become much weaker. Both Boyer, et al. ( 3l), and Huang and Chou<32), have reported similar multiple shock waves propagating away from bursting spheres. The pressure-time behavior at fixed Eulerian radius is shown by figure 38. Inside the source volume (n=0.825) the pressure rises instaneously to 9.0 and remains until the ex- pansion wave propagates from the edge. The pressure decreases to less than ambient at T=0.38. The second shock reflects from the center and passes at T=0.9. At positions outside the source volume there is a rapid pressure rise as the lead shock arrives followed by a nearly exponential pressure decrease to less than ambient. MACH 8.0 In case 2, the energy addition wave propagated at a dimensionless Mach number of 8.0, which for steady-state one-dimensional flow corresponds to supersonic combustion or a weak detonation. The energy wave movement is so rapid relative to the ambient velocity of sound that there is 82 minimal reinforcement of pressure, even during energy addition, Figure 14 shows there is no pressure transferred ahead of the energy addition and the pressure peaks at the end of the energy addition. The peak pressure in the source voltnne is greater than bursting sphere because of the reinforcement of the energy (pressure) propagating with the energy wave. Figure 16 shows the pressure distribution of the flow field. After the energy addition ends, the shock wave propagating from the source volume develops. Comparing this flow field to the flow field in figure 15 (bursting sphere), at equal radii in the far field the shock overpressures are equal and the flow fields behind the shock are similar. The particle behavior during energy addition is shown by figure 28; in cell #1 the pressure initially rises and, due to the non-steady behavior, the cell expands to the Reactive Hugoniot in the excluded region for steady-state solutions. When the energy addition wave has progressed through five cells the energy addition begins to approach the steady-state solution and the exclude region is no longer entered during later energy addition. The p-v behavior of cells 20 through 50 is a straight line which is indicative of a steady-state Rayleigh line. Figure 48 shows that as the energy addition wave over- rides the particles there is a small volumetric expansion during and shortly after the energy addition wave passes the particles. There is no further expansion until the energy addition wave reaches the edge of the source volume and the 83 expansion wave propagates into the center. The pressure-time history at a fixed Eulerian radius is shown by figure 39. Inside the source volume (n=0.825) the pressure changes from ambient to the peak within a time of 0.0106 because the wave velocity is so high that there is no pressure wave propagating ahead of the energy addition wave. After the energy wave passes there is a gradual pressure decrease until the expansion wave propagates through the position. The pressure continues to decrease until the ex- pansion wave is reflected from the center and reaches the position. The pressure drops below ambient at T=0.62, followed by a reflected shock which arrives at T=l.05. As the pressure wave propagates outside the source volume the peak pressure decreases at larger radii. However, at larger radii the pressure decrease behind the shock is not nearly as great as an exponential decrease and approaches a linear decay. MACH 5.2 PLANAR GEOMETRY Steady-state theory is based on the assumption of con- stant area flow. For comparison, the development of the flow field for Chapman-Jouguet conditions was first studied for the case of planar geometry (constant area). Figure 17 shows the development of the blast wave during energy addition. Of particular note is the p-v behavior shown by figure 29. When the energy addition wave passes through the last cell the pressure has reached the predicted steady-state value. The change in cell properties is a straight line from the initial 84 to final conditions, implying Rayleigh line behavior. The cells at the edge of the kernel appear to tangent the isen- tropic behavior behind the energy addition. At the CJ point the Reactive-Hugoniot and isentrope are tangent verifying that the Ma.ch 5. 2 wave exhibits CJ behavior, as it should· MACH 5.2 SPHERICAL GEOMETRY In case 3 an energy addition wave of Mach 5.2, the Chapman-Jouguet value -for steady-state conditions, was run in spherical coordinates. At this velocity the Rayleigh line for steady-state conditions tangents the Reactive Hugoniot. Figure 14 shows very little pressure increase ahead of the energy wave with the pressure peaking at the end of energy addition. The development of the flow field is shown in figure 18. As the energy wave propagates the peak pressure rises and assyrnptotically approaches but does not reach the predicted CJ pressure of 15.08. This can be attributed to the divergence associated with the spherical flow field. The p- v behavior of the individual cells, shown in figure 30, is quite similar to the behavior for the Mach 8.0 addi- tion. The center cells experience a pressure increase and ex- pansion into the excluded region. As the flow field develops the cell behavior approaches Rayleigh line behavior. The cell at the edge of the source volume (cell 50) almost tangents the isentrope. The particle displacement, shown in figure 49, is simi- lar to the other supersonic cases . Before 'the energy wave arrives there is no displacement of the particles. During 85 the energy addition there is some particle displacement caused primarily by the expansion behind the eriergy wave. After the particles expand to nearly equal pressure (P~S.25) behind the addition there is little particle movement until the wave has propagated through the source volume and the expansion wave propagates into the source volume. This is followed by a series of reflected shocks and expansion waves. The pressure-time behavior of the flow field at Eulerian positions is shown in figure 40. Within the source volume (n=0.825) there is an almost discontinuous rise to the peak pressure decreasing to nearly uniform pressure behind the energy wave. The expansion wave propagates from the edge of the source volume, causing a rapid pressure decrease to less than ambient at T=0.67. A reflected shock arrives at T=l.05. At greater radii the sharp peak becomes more and more diffuse. MACH 4.0 In case 4 the energy addition wave propagated at Mach 4.0. This is an impossible velocity according to steady-state theory. At this velocity the Rayleigh line for the steady- state solution does not intersect or tangent the Reactive Hu- goniot. The structure of the blast wave during and after energy addition is shown in figures 19 and 20. The energy addition wave moves supersonic relative to both ambient conditions and conditions behind the energy addition (a4/a0 =2.78). Since the acoustic velocity behind the energy addition approaches the energy addition wave velocity the pressure is reinforced 86 and peaks within the energy addition wave as shown by figure 14 (note: Figure 14 is based on Lagrangian positions. Fluid compression and expansion gives the Eulerian distribu- tion of figure 19). As the energy addition wave propagates through the source volume the peak pressure rises, reaching a maximum pressure of 19.7 at the edge of the source volume when the energy addition ends. The particles are displaced outward by the shock, reaching a particle velocity as great as 3.6 at the peak. When the energy addition reaches the edge of the source volume the pressure decreases and a shock wave is formed. As the shock wave propagates away from the source volume an expansion wave propagates into the source volume. As the pressure peak goes through the transition from an energy addition wave to a shock wave, a "valley" in the pressure distribution can be seen at T=0.25. Since the peak pressure occurs at the middle of the energy addition wave, as the wave propagates through the edge of the source volume the pressure at the leading edges of the energy addition wave continues to propagate. However, in the center of the addi- tion wave (tapered region of the source volume) the energy is less than at the edges of the source volume, resulting in a valley in the pressure distribution curve. The pressure-time distribution at Eulerian radius is shown by figure 41. Within the source volume (n=0.825) there is a high (P=l5.0) but very short pressure peak as the energy addition wave passes . The wave passage is followed by a 87 pressure decrease approaching the uniform pressure (P ~S.45) behind the energy wave. The 'propagation of the expansion wave into the source volume causes a rapid pressure drop with the pressure decreasing to below ambient at T=0.68. Outside the source volume (n=l.15) the shock passage has a peak pressure of P=l0.6 which rapidly decreases to P=5.0 followed by nearly exponential decay through ambient. At greater radii the high peak of short duration disappears and the blast wave structure becomes similar to that of a bursting sphere (Figure 15). From the particle paths in figure 50 it can be seen that the effects of energy addition do not affect the flow field ahead of the energy addition wave. i.e., when the energy addition reaches the edge of the source volume (T=0.21) there has been no movement of the particle. As the energy addition wave propagates through the source volume the shock wave which is formed entraps particles and moves them outward. Behind the wave the particle velocity decreases and a nearly uniform pressure exists. When the energy addition ends an expansion wave propagates into the source volume. However, since the pressure behind the energy wave is lower than for the bursting sphere, the effects at the center singularity point are re- duced. From the pressure-specific volume plot of figure 31 it can be seen that since the approach flow Mach number is less than the Chapman Jouguet velocity, in the late stages of heat addition the pressure does not peak at the end of energy 88 addition but decreases until the Reactive Hugoniot is reached. Examining the energy addition as it begins at the center, the first cell experiences a pressure increase and volumetric expansion until energy addition begins in the second cell. This prevents further expansion of the first cell and further energy addition results in a pressure increase, and specific volume decrease. The behavior of the second and third cells is quite similar. However, in the fourth and fifth cells · there is some compression of the particle during the energy addition. In cells 10, 20 and 30 there is initially com- pression as the pressure rises until the properties reach the Reactive Hugoniot. The particles then experience an expansion and pressure decrease along the Reactive Hugoniot until energy addition ends. Cells 40 and 50 are subjected to a pressure rise before energy addition begins and do not reach the Reactive Hugoniot. At the end of the energy addition there is a specific volume increase to bring the cell properties to the Reactive Hugoniot. These characteristics of the flow field indicate that the flow field remains non-steadv. i.e .. there is no steadv-state solution. The flow approaches a quasi-steady-state, but because the p-v behavior during the energy addition is a curved line the addition is definitely not Rayleigh line behavior. MACH 3.0 The Lagrangian pressure distribution for the Mach 3.0 energy wave has a pressure rise ahead of the energy addition 89 as shown in figure 14. The pressure peaks at the leading edge of the energy addition wave and decreases during energy addition. Figures 21 and 22 show the flow field behavior of the Mach 3.0 is similar to the flow field generated by a Mach 4:.o energy wave, but at lower overpressures. The Mach 3.0 addition is an impossible steady-state solution for the ambient conditions. However, the pressure wave ahead of the energy wave raises the temperature to 0/0 =2.4, changing the 0 properties. The p-v behavior in figure 32 shows the cells at the edge of the source voltnne exhibiting similar behavior with the pressure rise ahead of the energy wave greater as the edge is approached. The p-v behavior during energy addition is not a straight line, indicating non-Rayleigh line behavior. But there is a pressure decrease during the energy addition in- dicating the energy addition is approaching deflagrative be- havior. MACH 2.0 In case 6 the energy addition wave propagated at Mach 2.0. This velocity is supersonic relative to the ambient conditions, but subsonic relative to the properties behind the energy addition wave (a4 /a0 =2.78). This permits energy to be trans- fered ahead of the energy addition wave and the pressure dis- tribution asstm1es the form shown in figures 14 and 23. As the energy addition wave propagates through the source volume a pressure "hump" (P=8.0) develops ahead of the wave . With the arrival of the energy addition the pressure decreases to a 90 nearly uniform pressure (P=4.0) behind the energy addition. Since a pressure decrease ·across the ene~gy addition is a characteristic of a deflagration, an examination of figure 33 will explain the behavior. Initially the acoustic velocity throughout the flow field is the same, ambient. When the energy addition begins in the first cell the energy wave is propagating supersonic relative to the entire flow field. For the first five cells there is no propagation of pressure ahead of the energy wave and during energy addition the cell pro- perties change from nearly ambient to a pressure-specific volume relationship on the Reactive Hugoniot. When the energy wave reaches the tenth cell a pressure "hump" has begun to propagate ahead of the addition wave and the cell properties have been displaced along the shock Hugoniot (P ~l.4) before the energy addition begins. As the energy wave reaches cell 20 the pressure wave ahead of the energy wave has changed the cell properties along the shock Hugoniot (P ~S.7). For cells 30, 40 and 50, the pressure ahead of the energy wave approaches a uniform value of P=8.0, with a pressure drop and specific volume expansion across the energy wave. Since the p-v-line for the energy addition in the final cells approaches a straight line which tangents the isentrope, this case approaches the special case of the lower Chapman-Jouget state for the pressure-specific volume properties ahead of the energy addition. The displacement of succesive plots of the Reactive Hugoniot is caused by transfer of energy away from the cell during the energy addition. Although an energy of 91 42.5 is added to each cell, the cell ene!gy of the cells near the edge of the source volume at the end of energy addition is only 38. The other energy has been transferred into the flow field. Figure 42 illustrates the pressure distribution of the flow field at fixed Eulerian radius. At a location inside the source volume, n = 0.825, there is a rapid pressure rise to P=8.0 at T=0.26 as hhe energy wave approaches. The pres- sure falls through the energy addition to a nearly uniform pressure (P=4.0) behind the energy addition. This pressure is nearly constant until the energy wave propagates past the edge of the source volume and an expansion wave propagates towards the center. The expansion wave causes a pressure decrease through ambient pressure at T=0.85. At the position just outside the source volume, n=l.15, the expansion of the source volume during energy addition results in the energy wave traversing this Eulerian radius. The position is first subjected to the pressure field ahead of the energy wave followed by a pressure decrease during the energy addition. The expansion wave then causes the pressure to decrease to below ambient at T=l.10. At greater radii the peak pressure decreases and the blast wave begins to approach the form of a shock wave. However, the effects of the rapid pressure rise ahead of the energy addition can still be seen at the n=l.6 and n=2.3. The particle displacements can be seen in figure 51. As the energy wave propagates through the source volume the 92 particle movement occurs primarily ahead of the wave. The particle velocity is a maximum at the leading. edge of the wave and decreases to a minimum at the end of the energy addition. After the energy addition is completed the flow field experiences a series of expansion and shock waves reflecting from the center. MACH 1.0 In case 7 the energy addition wave propagated at the ambient velocity of sound. The addition of energy increases the local velocity of sound and energy (pressure) is trans- ferred ahead of the energy addition as shown in figure 14. Figure 24 shows the flow field approaching a self-similar solution. As the energy addition wave propagates from the origin the flow field develops and the peak pressure asymptotes to P=3.5. The leading edge of the flow field ex- periences a rapid pressure rise at the limits of energy transfer. This is followed by a slow pressure rise to the peak pressure at the leading edge of the energy addition wave. Across the energy addition wave the pressure drops to a nearly uniform pressure of P=2.6 behind the wave. This self-similar wave structure continues until T=0.85 when the energy wave has propagated through the source volume. The wave structure changes with the peak moving to the leading edge of the pressure rise as the expansion wave is generated. This can also be seen in figure 43. When energy addi- tion is completed the edge of the source volume has expanded to a radius of 1.5. The positions n=0.825 and n=l.15 are both 93 traversed by the energy wave. At position n=0.825 there is initially a rapid pressure rise when the pres·sure ahead of energy addition arrives·. This is followed by a slow pressure rise to the peak pressure at the beginning of the energy addition wave. The pressure drops through the energy addition to a nearly uniform pressure behind the wave. This uniform pressure continues until the expansion wave forms at the end of the energy addition and propagates back into the source volume. Similar behavior is noted at n=l.15. The position n=l.6 is located just beyond where energy addition ends. The pressure decrease through the energy addition has been replaced by an expansion wave. The leading edge of the blast wave is similar to the pressure profile ahead of the energy addition, however the expansion wave results in the pressure decreasing to below ambient behind the wave. At greater radii the blast wave has a rapid rise to the peak pressure followed by a rapid decrease tapering to a nearly linear decrease through ambient pressure. Most of the particle displacement shown on figure 52 takes place ahead of the energy addition wave. As an example, for the particle initially at D=0.8 the energy addition begins at T=0.76. From figure 34 it can be seen that initially the parti- cle p- v behavior is definitely non-steady. When the energy addition wave has propagated through 20% of the source volume the flow field begins to approach a self-similar solution. 94 Initially the particle goes through a pressure rise along the shock Hugoniot. During energy addition the particle goes through a weak deflagration along a Rayleigh line. MACH 0.5 For case 8 the energy wave is propagating subsonic rela- tive to both ambient conditions and conditions behind the energy wave. Comparing figures 14 and 25, the compression and pressure rise ahead of the energy wave can be seen. As the pressure propagates ahead of the wave there is first a pressure rise along the shock Hugoniot followed by an isentropic compression to the beginning of the energy addition. There is an expansion and pressure decrease through the energy addition with nearly equal pressure behind the energy wave. As the flow field develops the pressure increases and ~symptoticallY approaches a final pressure of P=l.88. In the final stages of energy additions the flow field approaches self-similar behavior. Figure 35 shows the energy addition is a pressure decrease along a straight line in the p-v plane, implying Rayleigh line energy addition as a weak deflagration. The energy wave is propagating much slower than the lower CJ deflagration condition. The low peak pressure associated with this energy addition results in a large expansion through the energy addition wave. This can be seen in the parti cle displacement curves of figure 53. The particles are initially displaced by the pressure rise ahead of the energy wave. As they go through the 95 expansion associated with the ·energy addition their velocity decreases to nearly zero as shOwn by the nearly constant position after the initial displacement. The particle positions remain nearly constant until the expansion wave propagates through the source volume. Since the source volume has experienced considerable expansion during energy addition the secondary shocks are much weaker than for the cases of supersonic addition. This is also shown by figure 44. Inside the source volume (n=0.825) there is initially a rapid pressure rise beginning at T=0.55 followed by a slower rise until energy addition begins (P=l.85). The pressure decreases during energy addition to nearly constant (P=l.69) behind the energy addition, until the expansion wave at the end of energy addition (T=l.69) propagates to the position (T=l.96) causing a rapid pressure decrease to below ambient. A second shock is formed, but the pressure does not exceed ambient. The expansion through the energy addition results in a large expansion of the source volume. When energy addition is completed (T=l.69) the edge is at an Eulerian radius of 1.66. The expansion of the source volume causes the positions n=l.15 and n=l.6 to experience behavior similar to n=0.825, only the initial pressure rise occurs later and the propa- gation of the expansion wave into the source volume occurs earlier. At n=2.3 the pressure rise is similar to the rise ahead of the energy addition, however, at greater radii (n=3.2) the peak appears to be moving to the front of the shock. MACH 0.25 96 The Mach 0.25 case ·is quite similar to the Mach 0.5 case except the lower energy wave velocity allows the solution to approach acoustic behavior. Figure 36 shows a slight com- pression and pressure rise to P=l.32 ahead of the energy wave and a Rayleigh line energy addition with pressure decrease to 1.25. This is indicative of a nearly constant pressure de- flagration. The edge of the source volume has expanded to n=l . 84 when energy addition is completed. Figure 45 shows that at an Eulerian position inside the source volume (n=0.825) the pressure begins to rise at T=0.71, the time required for an acoustic signal to propagate from the center. The pressure rises to P=l.31 ahead of the energy wave and decreases to P=l.25 behind the addition. The expansion of the source volume causes similar behavior at n=l.15 and n=l.6. Outside the source volume the pressure rise is similar to the pressure rise ahead of the energy wave. The overpressure decreases, but since the initial overpressures were low the shock wave decay is slowed. There is a gradual expansion of the flow field as shown in figure 54. MACH 0.125 For the Mach 0.125 the energy wave is propagating so slowly the energy addition approaches a nearly constant pressure deflagration. Figure 27 shows a nearly isentropic 97 pressure rise to P=l.08 ahead of the energy wave. Through the energy addition the pressure de~reases to P=l.075, a nearly constant pressure expansion. Similar behavior is seen in figure 46. At the time required for an acoustic wave to propagate to the Eulerian positions the pressure begins to rise. Figure 55 shows particle movement ahead of the energy wave, with a large expansion through the wave. This case was run only until trends were established because excessive computer time was required. 2. Damage Parameters Experimentally the parameters which are normally observed in blast wave studies are peak pressure, P, and positive im- s pulse, r+·· calculated from the pressure-time history of the blast wave. Using these parameters and the P-I technique described earlier, accurate estimates of structural damage can be made. The peak overpressure as a function of energy scaled dis- tance for cases one through eight and Baker's pentolite data correlation are shown in figure 56. The behavior of the high explosive pentolite does not compare directly with the gas mixture under consideration but is plotted for illustrative comparison. These variables are plotted as they were defined in equations I-8 and I-9. In all cases the overpressures were considerably below the overpressure from an explosion of pentolite with the same total energy. This is caused by the non-ideal structure of the blast wave and the low energy density. 98 Bursting sphere (infinite wave velocity) is the limit case for the wave addition of energy and results in a constant overpressure from the center to the edge of the kernel. After energy addition, a shock front develops, propagating away from the source volume. Beyond the energy source volume the shock overpressure has a maximtun value of P=3.40. In the far field the overpressure of the bursting sphere approaches 70% of the high explosive curve for the same energy scaled radius. As the energy wave velocity decreases through Mach 4.0, the near field overpressure associated with the energy addition increases. Because of the large overpressure associated with the energy wave the shock propagating away from the source vollll'lle initially has a peak pressure greater than the bursting sphere case but decreases to 90% of bursting sphere in the intermediate field. In the far field the overpressure curves coalesce to approximately 70% of the pentolite correlation. As the velocity decreases from Mach 4.0, the near field overpressure decreases. For each 50% decrease in the energy wave velocity the near field overpressure decreases by the following relationship: overpressure (50% velocity)~0.35*[overpressure (100% velocity)] IV-8 In the near and intermediate field all the supersonic cases initially have an overpressure greater than bursting sphere. At an Eulerian radius of n=l.98 the overpressure curves of the supersonic cases intersect and at a radius of n=2.0l their pressures begin to drop below the bursting sphere overpressures. The overpressure in the Mach 2.0 addition 99 decreases to approximately 75% of bursting sphere at a radius of n=2.73. The overpressure then approaches bursting sphere and reached 90% when the calculation was ended. In the case of the energy wave propagating at the ambient velocity of sound, Mach 1.0, the expansion behind the energy addition results in shock wave ahead of the energy addition. When the energy addition ends, this shock wave continues to propagate with only a very gradual decrease in overpressure. Between a radius of n=l.96 and n=2.24 this case has the greatest overpressure. The overpressure then begins to drop rapidly as the expansion waves behind the shock decrease the shock overpressure. If the flow behavior behind the Mach 1.0 addition is similar to the Mach 2.0 the overpressure will begin to approach the bursting sphere in the far field, as it did in the Mach 2.0 case. The subsonic energy additions exhibit expansions of the source volume behind the wave. However, as the velocity de- creases the expansion does not produce the near field and intermediate field overpressures necessary to approach the overpressures from bursting sphere. The Mach 0.5 and Mach 0 . 25 overpressures approach, 84% and 23% of bursting sphere, respectively. Figure 57 is a plot of non-dimensional impulse, I, versus energy-scaled distance, RE.I is defined by Sachs' relationship and is expressed as: I = IV-9 where I+ is the positive (p )2/3(E )1/3 ho . Tl p ase impu se, a0 and p0 are the ambient atmospheric values of sound speed and pressure, 100 respectively, and ET is the total energy deposited within the source volume. For comparison the impulse of a high explo- sive, pentolite, is also plotted. Because impulse is the integral of overpressure with time, the overpressure and impulse plots exhibit similar behavior when plotted against similar parameters. For the supersonic energy addition, the impulse is higher in the near, intermediate and far field than the subsonic cases. As the energy wave velocity decreases the impulse decreases for the entire flow field. I In the near field the impulse from the theoretical energy addition is greater than the experimental correlation for pentolite because of the positive pressure behind the energy addition wave which exist until the end of the energy addition. In the far field the impulse varies from 60 to 75% of that for the high explosive (pentolite). 3. Energy Distribution In an ideal or point source explosion all the energy is transferred to the surroundings and is available to drive the blast wave. In a non-ideal or diffuse explosion the source releases energy relatively slowly over a sizeable volume. In addition, the mass in the source volume retains a portion of the energy, reducing the amount of energy available to drive the blast wave through the surroundings. The energy which remains in the source volume can be used as a measure of the "effectiveness" of the explosive process relative to an ideal (point-source) explosion. 101 The concept of "was.te ~nergy" was introduced by Taylor( 3) who surmised that some energy would rerriain or be "wasted" in the central core region of the blast zone. This energy which remains in the source volume ·after the shock passage and an adiabatic expansion to ambient pressure is unavailable to the pressure wave and has also been called "residual energy" by Strehlow and Baker(Z?). They noted that the energy distri- bution in the system and how it shifts with time are two important properties in determining the behavior of an explo- sive process. Adamczyk(lB) analyzed his non-ideal explosions (produced by homogeneous addition of energy) and noted that the time over which energy is added to the source region determines the structure of the blast wave and the partitioning of energy between the source volume and the surroundings. He considered two idealized limit cases of constant volume energy addition and constant pressure expansion. The first case of constant volume energy addition, bursting sphere, can be visualized as an infinitely fast energy addition wave with an instantaneous deposition time. Initially the source volume is at the ambient temperature and pressure of the surroundings. Energy is instantaneously added, raising the temperature and pressure of the source volume to the initial conditions of the bursting sphere. The energy added is: Ai = n fc (84-0 ) + (C -C )0 1 . L V 4 0 V 4 VO OJ IV-10 102 IV-11 and the energy density is given by: P4 1 q = - - Po IV-12 0 q = 4 1 0 - IV-13 0 where y4 is the constant garrnna of the gas in the source vol- ume after energy addition and y is the initial garrnna through- o out the field. If the initial and final gannna's in the source volume are equal, the second term cancels and equation IV-11, is Brode 1 s< 33 ) formula for the energy stored in a bursting sphere. If the bursting sphere undergoes an idealized isentropic expansion where the sphere expands slowly against a counter pressure equal to its instantaneous pressure, the fraction of the total energy remaining in the source volume is: 1 [ ( 1 +q) 1 / y - 1 ] q IV-14 and the fraction of energy transferred to the surroundings is: 1 [(l+q) - (l+q)l/y] q IV-15 where Es is the energy transferred to the surroundings, EB 103 is the energy remaining in the source volume, and ET is the total energy deposited. Equation IV-15 is Brinkley 1 s< 34) or Baker ' s( 2 ) formula for the effective quantity of energy stored in the sphere, expressed as a fraction of Brode's energy. In the limit as q + 00 (point source), Es/ET + 1 and as q + 0, Es/ET + (y-1)/y. For the conditions being investi - gated: and In the second limit case the energy is added infinitely slowly such that the energy of both the source volume and surroundings remain at p 0 • The fraction of energy which remains in the source volume is : 1 y IV-16 and the fraction of energy transferred to the surroundings is: R = y-1 cP . Y IV-17 this is also the limit case for an infinitely rapid (constant volume), but infinitely small (q+O) energy addition. For the 104 conditions investigated: It should be noted that in both limit cases, q+o for bursting sphere and infinitely slow energy addition, there is no blast wave . In the cases studied all internal properties are initially at their ambient values throughout the system . At the in- stant chemical reaction begins, the heat addition model adds energy to the volume encompassed by the heat addition wave. As time progresses this energy is redistributed as internal and kinetic energy throughout the system, where the system contains all materials out to the lead characteristic or lead shock wave. The energy added to the system can be separated into four classifications: (1) Internal Energy increase in the source volume: r . J £ p00 rJ dr --~ dr - y -1 0 IV-18 0 (2) Kinetic Energy of source volume: r (K~)rn= J £ 2 rj PU dr 2 IV-19 0 105 (3) Internal Energy increise in the surroundings: roo p(G-0 )rj rco p0 rj (IE).s= f -f 0 dr 0 dr y -1 y -1 4 0 r r £ £ (4) Kinetic Energy of surroundings: 2 . pu rJdr 2 IV-20 IV-21 where O is the center of the sphere, r£ is the position of the contact surface of the ball containing the high energy gas, and r 00 is the limits of the flow field. Figures 58 through 66 illustrate the energy distribution for the cases investigated and how it varies with time. Fig- ure 58 shows an instantaneous addition of the total energy to the source volume. Since the instantaneous energy addition is a constant volume energy addition, initially 100% of the energy is internal energy in the source volume. As the flow field develops this internal energy shifts to kinetic energy in both the source volume and the surroundings, and internal energy in the surroundings. As the source volume expands its kinetic energy rises and peaks when the expansion fan reaches the center, followed by an oscillatory decay. The kinetic energy in the surroundings increases until there is a max- imum in the rate of displacement of the source volume at t~0.66. The kinetic energy of the surroundings gradually 106 decreases as the shock wave propagates into the flow field. The internal energy of the air continually rises and asymp- totically approaches a final value of 36%. The internal energy of the source volume appears to asymptotically ap- proach a final value of 66%. In case 2(Mw = 8.0) the movement of the energy wave through the source volume generates kinetic energy of the entraped · particles. Since the energy wave moves supersonic there is no energy transfer to the surroundings until the energy addition wave reaches the edge of the source volume (T=0.116). There is a rapid rise in the internal and kinetic energy of the surroundings as the energy wave propagates into the surroundings and continues as a shock wave. The expansion wave which propagates into the source volume develops a large value of kinetic energy in the source volume. The internal energies approach final values of 63% in the source volume and 37% in the surroundings. In case 4 (MW= 4.0), the large overpressure of the energy wave imparts considerable kinetic energy to the parti- cles in the source volume. This kinetic energy maximizes and decreases abruptly when the shock enters the surroundings (T=0.23). The expansion wave then increases the kinetic energy of the source volume until the wave reflects from the center. Subsequent expansion waves reflecting between the center and the shock have less kinetic energy. The internal energy of the source volume decreases from a value of 98% when the addition wave reaches the edge of the source volume 107 to a final value of 60%. The internal energy of the surround- ings approaches 40% of the energy added. In cases 6(Mw = 2.0), 7(Mw = 1.0), B(Mw = 0.5), and 9(¾ = 0.25), figures 62, 63, 64, and 65 respectively, there is energy transfer ahead of the energy addition wave. This causes a movement (displacement) of the particles resulting in an increase in the kinetic energy. As the energy wave approaches the edge of the source volume the particle move- ment ahead of the wave moves into the surroundings with the kinetic energy abruptly decreasing in the source volume and increasing in the surroundings. As the expansion wave pro- pagates into the source volume the kinetic energy increases, but not to the level reached during the passage of the heat addition wave. At later times the kinetic energy of the source volume decreases as successive expansion waves become weaker. The final distribution of energy is; for case 6, 61% source volume, 39% surroundings; case 7, 66% source volume, 34% surroundings; case 8, 74% source volume, 26% surroundings; and in case 9, 77% source volume, 23% surroundings. As the Mach number of the energy addition wave de- creases, the overpressure also decreases resulting in a weaker shock wave propagating into the surroundings and consequently there is less energy transfer. In a non-steady heat addition the limit case of a con- stant pressure expansion can not be reached since any heat addition, even at very low subsonic velocities will result 108 in a pressure rise ahead of the energy addition wave and pressure decrease through the energy addition. For the cases run the energy distribution approached 77% in the source volume and 23% in the surroundings for very slow flame pro- pagation velocities. The energy distribution for case 10 (MW= 0.125) was not calculated since the complete energy addition was not run. Examining the energy distribution for the cases which were run it can be seen that for a constant energy density the energy distribution is significantly affected by the Mach number of the energy wave. The principle mechanism for transfer of energy to the surroundings is the propaga- tion of the shock wave through the flow field. For the cases of a highly supersonic energy addition wave, there is very little kinetic energy in the flow field as the wave propagates. When the energy addition stops there has been only minimal development of the flow field. The distribution of energy between the source volume and the surroundings and how this distribution shifts with time as a function of the flame velocity is sunnnarized by figure 67. For the limit case of infinite energy wave vel- ocity, bursting sphere, 37% of the energy is transfered to the surroundings by the final time line calculation. The energy transfer to the surroundings increases to 41% as the velocity decreases to Mw = 4.0. As the velocity is decreased further the energy transfer to the surroundings decreases to 23% in case 9(¾ = 0.25), the lowest velocity for which the energy distribution is calculated. It 0. ... Cl) > 0 109 100.0 .--------r----., ----------------, I I : I I I M4.0 10.0 10.0 M 2.0 M0.5 1.0 I 1.0 I I M0.25 I I I .9 1.0 1.1 1.2 Tail Head 0.1 '-----'----___. ____ __. __ .__ _ ___,__---'.__ _ __,_ __ __._._ ____ ~ 0.1 .7 .8 1.3 1.4 Lagrangian Position (Note: Eulerian positions will vary because of compression) Figure 14 , Overpressure Distribution Through Energy Addition Wave r 19.0 13.0 p 7.0 0.05 3.00 2.50 2.00 I.So 1.00 a.so PRESSURE/ PO DlSTAlBUTl~N VS. DISTANCE/ DO AND TIME/ TO Figure 15. Pressure distribution versus Eulerian distance and time for blast system generated by an infinite velocity energy - ,, ~ - = ....__ .! - - - - 'L- - - - ' t--' t--' 0 0.00 19.0 13.0 p 7.0 1.0 2.0 ~J.O 3.00 2.50 2.00 1.50 1.00 ·0.50 PRESSURE/ PO DISTAIBUTI~N VS. DISTANCE I DO AND TIME/ TO Figure 16. Pressure distribution versus Eulerian distance and time for a blast system generate by a Mach 8.0 energy addition wave. 0.00 0 0 . lf) ..... 0 0 . N ..... 0 0 . en 0 0 0 ,- . 0 0 N . 0 0 ("f) . 0 0 Ll) . 0 0 1.0 . 0 0 ,...,. . 0 0 co 0 . O"t 0 0 . 0 0 ,- a-t-----.---- -----r- ----r-----.-----.-------.----- --.----~ 0.00 .so 1.00 1.50 2.00 RADIUS 2.50 3.00 3.50 1,1, 00 Figure 17. Pressure distribution versus Eulerian distance and time from a blast system generated by a Mach 5.2 (CJ) energy wave in planar geometry. 0 0 . l.f) ..... 0 0 . ('\J ..... a 0 . en 0 0 . (I") 0 0 0 . 00 .20 Figure 18. \0 "'T ,- ,- . N . 0 co ,- 0 ,- 0 . . 0 0 N . 0 .40 . 60 . 80 LOO 1.20 1.60 ARDIUS Pressure Distribution Versus Eulerian Distance And Time From a Blast System Generated By A Mach 5.2 (CJ) Energy Wave in Spherical Geometry. p 22.0 1 5. 0 8.0 1.0 0.3 2.00 1.67 1.33 1.00 0.67 0.33 PRESSURE/ PO OISTRIBUTION VS. DISTANCE/ DO AND TIME/ TO Figure 19 . Pressure distribution versus Eulerian distance and time for a blast system generated by a Mach 4.0 energy wave. 0.00 17 .0 13.0 p 7.0 0 3.00 ·2.50 2.00 I.SO 1.00 a.so PRESSURE/ PO DISTRIBUTieN VS. DISTANCE/ DO AND TIME I TO Figure 20. Pressure distribution versus Eulerian distance and time for a blast system generated by a Mach 4.0 energy addition wave. 0.00 p 16.0 11.0 6.0 l.Q 0.4 2.00 1.67 t.00 0.67 0.33 PAESSUAE / PO OISTAIBUT!~N VS. DISTANCE/ 00 ANO TIME/ TO Pressure distribution versus Eulerian distance and time for a blast system generated by a Mach 3.0 energy addition wave. Figure 21. 0.00 19.0 13 .0 p 7.0 1.0 2.0 ·3.00 2.50 2.00 1.50 1.00 -a.so PRESSURE/ pa ·o1STAIBUTlflN vs. DISTANCE/ 00 ANO TIME/ TO Figure 22. Pressure distribution versus Eulerian distance and time from a blast sys~em generated by a Mach 3.0 energy addition wave. O.ClO 0 0 . 0 ...... N N 'd' M . ,- . 0 N 0 0 . 0 0 . co ,- ,- M . LO 0 . 0 CJ 0 . CD LLI a:o :::,o en . cn-::J' LLI a: a... 0 CJ . N 0 0 OJ----r----r-----r----r------r-----r---------r--- -~ . 00 .40 .80 1. 20 1. 60 RADIUS 2.00 2.40 2.80 Figure 23. Pressure distribution versus Eulerian distance and time for a blast system generated by a Mach 2.0 energy addition wave. 3.20 t-' t-' 00 4.0 3.0 2.0 1.0 l. a.a 3.00 2.50 ·2.00 1.50 1.00 0.50 PRESSURE/ PO DISTAIBUTieN VS. DISTANCE/ DO ANO TIME/ TO Figure 24. Pressure distribution versus Eulerian distance and time for a blast system generated by a Mach 1.0 energy addition wave. 0.00 0 lf) . N 0 0 . N 0 lf) . u.J CI:o :::Jo en • cn- u.J a: (L 0 lf) . 0 0 ....- ....- ....- ...... co 0 0 ·~------,,---------.-------r-----r--------.-------,.------~----0 .oo .60 1. 20 1. 80 2. 1.rn 3.00 3.60 LL20 L!.80 RADIUS Figure 25. Pressure distribution versus Eulerian distance and time from a blast system generated by a Mach 0.5 energy wave. 1--' N 0 0 ('l .- LJ1 w rr::o =:Jo if) . if).- w rr:: CL U1 CD 0 M co ,... ._,. . I . \ ' \ \ \ I \ \ ._,. ~ N ~ ,..... ,..... ~ ! ,' \ \ \ co ._,. . N 0 ,..... M ._,. ...... O'\ M ...... ,..... 0 . . . . . N M M ._,. N ' \ \ \., \ \ \ C"--1-4-----~-- D.00 . 75 Figure 26 . 1. so 2.25 3.00 3. 75 LJ .50 5.25 6.00 ARDJUS Pressure distribution versus Eulerian distance and time from a blast system generated by a Mach 0.25 energy addition wave. 0 ...... ...... Ln a ...... 0 0 . w ...... a: ::) (f) (f) w 0:lf) Q_O) . M 1.0 co 1.0 0 ,- co ,,:;t . N .- CV) . rr, 0 ~-.------r------,r--------r------.-- ----,------r----- -....----- 0. 00 .75 1.50 2.25 3.00 3.75 4.50 5.25 6.00 RADIUS Figure 27. Pressure distribution versus Eulerian distance and time for a blast system generated by a Mach 0.125 energy addition wave. t-' N N Cl = lf') r a::J = lf') ID u..J cc :::J= c..no c..n • u..Jln cc (!__ lf') r- = lf') N = = .55 1.00 123 CELL F-V BEH~VIOR M~CH 5.00 ENERGT ~DDITION SFHERIC~L GE~METRT ll BEGIN ENERGT ~DD IT I ON ~ENO ENERGT ~OOITION CJcELL 1 c::icE1=1-.2 ~ L 3 +CE LL U XcELL 5 ~CELL l 0 "t"[ELL 20 XcELL 30 ZcELL UO YcELL SO 1 . 15 1 . '.30 1.l.15 1. 60 SPECIFIC VDLUME 1. 75 Figure 28. Pressure versus specific volume behavior from a Mach 8.0 energy wave (D = 1.0 at cell 50) C) C) . CP __. 0 0 . (,D __. 0 0 . ::r __. 0 0 N __. 0 0 . 0 __. 0 0 r.o 0 0 . ::r 0 0 . N 124 CELL P-V BEHAVJ~R MACH 5.2 ENERGY AOOJTJ~N WAVE ~JOTH 0.1 PLANAR GE(jHETRY [!J BEG 1 N ENERGY AOOJ T 1 ON ~ ENO ENERGY ADDJ T 1 ON CELL l ){ CELL 10 + CELL 2 z CELL 20 X CELL 3 y CELL 30 ~ CELL 4 ):?{ CELL 40 ~ CELL 5 * CELL 50 z CELL 60 0 0 . +------r-- ---.---- - -,-- - - ~-- --.-- - --, o . 2.50 3.00 0.00 .50 1. 00 1. 50 SPECJFJC V~LUME 2.00 Figure 29. Pressure versus specific volume behavior from a Mach 5.2 (CJ) energy wave in planar geometry • "11 ....... ' D l/1 . (") __. D 0 N ...... 0 lf") . D .- 0 0 m 0 I.Ji r-- 0 I.Ji . ::;I' CJ 0 . ('t") 0 l/") - 0 125 CELL P-V BEHRVJOR MACH 5.2 ENERGY ROOJ TJ CIN WAVE WJOTH 0.1 SPHERJCRL GECIMETRY [!) BEGJ N ENERGY RDOJ T HJN ~ END ENERGY ROD 1 T J CIN CELL 1 ~ CELL 10 + CELL 2 z CELL 20 X CELL 3 y CELL 30 ~ CELL 4 ~ CELL 40 ~ CELL 5 * CELL 50 ~ CELL 60 ~_j_----~---- ----r-----.----,-------, 0 0.00 .50 1.00 1.50 SPECIFIC V[jLLJME 2.00 2.50 3.00 Figure 30. Pressure versus specific volume behavior from a Mach 5.2 (CJ) energy wave in spherical geometry I 1 _ _._ ~-'1 I'\' 0 Ll"I 0 0 0 ('\J 0 l/1 r - 0 0 Cl u; 126 CELL F-V BEHRVJOA MRCH W,00 ENEAGr RDDIT SFHEAI[RL GE~METAr ll5EGIN ENEAG1 ~OOITION *END ENEAGt ~oorrrnN C:J[ELL 1 ~CEL_U ~L3 +cELL Ll XcELL 5 ~CELL 10 "l'-[ELL 20 :X:cELL 30 ZcELL ·10 YCELL SO <::::i c::i ~~r· ____ ,-___ __,,---------,----------- ---- 2.00 2.50 3.00 t. OD t. 50 0. OD SPECIFIC VOLUME Figure 31. Pressure versus specif.ic volwne behavior from Mach 4.0 energy wave~ .so 0 i.n ("") - C) Cl . N - 0 ll) 0 - 0 0 0) w a: ::::) 0 cno en . wt0 a: Q_ 0 i.n . :::r 0 0 . CT") 0 i.n ...... 127 CELL F-V BEHRVJ[jR MACH 3.0 ENERGY RODJT1t1N WAVE WIDTH 0. 1 SFHEAJCAL GE~NETAY QJ BEGIN ENERGY RDOJTJ(jN ~ ENO ENERGY AODJTJE!N CELL 1 :><: CELL 10 + CELL 2 Z CELL 20 X CELL 3 Y CELL 30 CELL 4 l:?l'. CELL 40 4> CELL 5 * CELL 50 Z: CELL 60 0 0----!------r------.------r------.-----.-----, 0 0.00 .so 1. 00 1.50 SPECIFIC VOLUME 2.00 2.50 3.00 Figure 32. Pressure versus specific volume behavior from a Mach 3.0 energy wave (D = 1.0 at cell 50). CJ CJ 01 0 0 (D 0 0 r-- 0 0 (D 0 0 lJl w a: =:lo cno en . LLJ :::r a: Q_ 0 0 (1"') 0 0 N 0 0 ...... 0 0 128 CELL F-V BEHRVI~R MACH 2.00 ENERGY ROOIT SPHERICAL GE~METRY ~BEGIN ENERGY RDOITI~N *END ENERGY RDDITICN C!JCELL 1 ~ LL 3 +CELL 4 XcELL 5 ~CELL 10 ~CELL 20 ~CELL 30 ZcELL 4Cl YCELL SO 0 -1--- --...--- ----.---- ----,- ---.-----.-----, 0.00 .so 1. 00 1.50 2.00 2.50 3.00 SPECIFIC VDLUME Figure 33. Pressure versus specific volume behavior from Mach 2.0 energy wave. • 0 lJ") ~l I cl C I ~1 w a= 0 L."': ~, C 0 r,-: Cl Lr. N _! ::Jo j enc cn N w a:: (L C) If) - 0 0 - 0 I.I) . i"-\ I ' \ I I \ fl~' I I . I I I I i I I , , I 'I I I 129 CELL P-V BEHRVJOR MACH 1.0 ENERGT ADOJTJ~N WRVE ~JOTH 0. 1 SPHEFJCAL GE~~ETRT QJ BEGJN ENERGY P.OOJTHJN ~ ENO ENERGT AODJ T 1 (jN D=D. 02 ~ D=O. 20 + X ~ 4' 0=0.04 O=O.OE 0=0.06 D=O. 10 Z 0=0.40 Y D=O. 60 ~ D=O. 80 *'. D= 1 . 00 Z D= 1. 20 0 0 o-+---------,-----r--------.------,,-----r------, o.oo 1.00 2.00 3.oo ij,Oo s.oo s~oo SPECIFIC V~LUME Figure 34. Pressure versus specific volume behavior from a Mach 1.0 energy wave. = C"" "' I.I") CJ <"\I Cl en I.I") r- CJ ID L.J...J cc ::) If) (fl :::1' en....; L.J...J cc (L lf1 - Cl Cl lf1 130 CELL F-V BEHi;iVIc:JA M~CH 0 , 5 ENERGY i;;ioo1r Ic:J f\J WJ;;iVE WIDTH 0,1 SFHERIC~L GE~METRY [!) BEGIN ENERGY RDDITION ~ ENO ENERGY ~DDITic:JN CELL ~ CELL 1 10 + CELL 2 z CELL 20 X CELL 3 y CELL 30 e CELL Ll n CELL LlO 't' CELL 5 * CELL so z CELL 50 '° --+-------.-- - ---..-- ----r--- ---,---- - ,------- o. oo 1. OD 2 . DD J . DD U.OD 5 . 0D Cl. OD SFECJFJC VOLUME Figure 35. Pressure versus specific volume behavior from Mach 0.5 energy wave- C er, Lf) N CJ N - I w i §slf) J <.n- <.n . w- I a::: ' (L C) . - If) 0 0 0 If) ' . / ,j I 131 CELL P-V BEHRVJ~R MRCH 0.25 ENERGY ADDITI~N W.~VE WIDTH 0. 1 SPH~9JCR~ GEr~~TR r QJ BEGJt~ El~ERGY RJJJTICN CD EN:J El·J~R3-;' RJJ J T J Ct~ .l D=0.02 :x: J=0.20 + D=C.D~ Z D=0.40 X D=C. 05 y D=O. 60 D=O. 08 ~ D=O. 80 ~ D=Ci.10 * D=l.00 Z D= 1. 20 ~-t----- ..--------r------,,-------.--------..------. 0.00 1.00 2.00 3.00 ij,00 5.00 6.00 SPECIFIC VCILUME Figure 36. Pressure versus specific volume behavior from a Mach 0.25 energy wave. CD ..... ..... ::r ..- _. N ..... . ..... 0 ..- CCI 0 ..... w a: ::) (.D cno en . w- a: Q_ ::r 0 . ...... N 0 ...... 0 0 . ..... CXJ 132 CELL P-V BEHRVJ~R MACH 0.125 ENERGY RDDJTJ[jN WRYE WJDTH 0.1 SPHERJCRL GE[jMETRY QJ BEG J N ENERGY ADD J T J ~N ~ END ENERGY ADOJTJ~N D=0.01 ;( D=D.11 + D=0.03 z D=D.19 X D=D.05 y D=D.29 D=D.07 ):?( D=D.39 ~ D=D.09 * D=D.49 z D=D.59 ~--+----~-----r------r-------.----~---~ 0.00 1.50 3.00 ~.50 6.00 7.50 9.00 SPECIFIC V(jLLJME Figure 37. Pressure versus specific volume behavior from a Mach 0.125 energy wave. u.J 0 L/') c::i 0 0 ai 0 L/') r.: 0 0 u:i II:o ::JLn en • en-:, uJ a: CL 0 0 '" 0 LI) 0 0 BURSTING Sf'HERE T R= .82S R= 1.15 + R= 1.6 X R= 2.3 <'> R= 3.2 o-+------r-----r--- --.-- ----.--------.------,-- ----r------.-----~ .OD .30 Figure 38. .60 Pressure versus generated by an .90 1.20 I.SO 1.80 2.10 2.ijQ TIME time behavior at fixed Eulerian radius from a blast system infinite velocity energy wave (bursting sphere). 2.70 .... c.., I.,,) w 0 l!) c:i 0 0 ai 0 l!) r: 0 0 u:i OCo :::) l!) (f) • (f)::!' w a: CL Cl 0 (T) Cl l!) 0 0 MACH 8.0 ADDITION T R= .825 R= 1.15 + R= 1.6 X R= 2.3 ~ R::; 3.2 0-i-----~----~ ----~--- - ~ - ---~----~--- - ~ - - - - ~--- - ~ .00 . 15 Figure 39. .30 .llS . 60 TIME .75 .90 1.05 1.20 1.35 Pressure versus time behavior at fixed Eulerian radius from a blast system generated by a Mach 8.0 energy wave. 1--' L,..) .p. w D 0 ::r 0 0 c-i 0 0 c:i 0 0 a:i a:a :=io (!) • (!)(D w a: [L 0 0 3 0 0 N 0 0 MACH 5.2 ADDITION T R= .825 R= 1. 15 + R= 1. 6 X R= 2.3 R= 3.2 a+------.-----.,..------,-------.------,,-------.------.-------.-----, .00 .30 .60 .90 1.20 1.50 1.80 2.10 2.llO 2.'70 TIME Figure 40. Pressure versus time behavior at fixed Eulerian radius from a blast system generated by a Mach 5.2 (CJ) energy wave. I-' (.;.) \Jl UJ 0 l/} . ..... 0 0 l/l 0 I/} f'i 0 0 Cl CCo :::>Ln Cl) • Cl) ..... UJ cc 0... 0 0 l/l 0 I/} ,,.j 0 0 MACH ij,0 ADDITION T R= .825 R= 1.15 + R= 1.6 X R= 2.3 ~ R= 3.2 o-+------r--------.-------,-------,-------.-------.-------.-------.-----~ .00 . 30 Figure 41. .60 .90 1.20 TIME 1.50 1.80 2. 10 2.70 Pressure versus time behavior at fixed Eulerian radius from a blast system generated by a Hach 4.0 energy wave . I-' w O'I __j uJ 0 LI) 0 0 0 ai 0 LI) r-= 0 0 (D CCo ::)LI) en • cn=r uJ cc a... 0 0 ...; 0 In 0 0 HACH 2.0 AODITl~N T R= .825 R= 1. 15 + R= 1.6 X R= 2.3 ~ R= 3,2 o-+-----~----r--------,,--------.-------r--- ---r------.------ .------ - -, .oo .30 .60 .90 l.20 TIME 1.50 1.80 2. 10 2.IIO Figure 42. Pressure versus time behavior at fixed Eulerian radius from a blast system generated by a Mach 2.0 energy wave. 2.70 I-' w -...J 0 lfl ~ 0 0 ~ 0 lfl N 0 0 "" w cr::o ::::) lfl cn • en- w a: a... 0 0 0 lfl 0 0 MACH 1. 0 ROD IT ICIN T R= .625 R= I. 15 + R= 1.6 X R= 2.3 e R= 3.2 o--1------------,.------,.-----,-------,------,------.----~----~ .00 .35 • 70 1.05 l.llO I. 75 2.10 2.115 2.80 3. 15 Figure 43. TIME Pressure versus time behavior at fixed Eulerian radius from a blast system generated by a Mach 1.0 energy wave. ..... w (X) CJ ,n en CJ a cri a lfl "' 0 0 "' w CI:o ~lfl (fl • (f)- w er: Q_ 0 0 0 lfl 0 a 0 .00 .55 1. 10 l.65 2.20 TIME HACH 0.5 ADDITION T A= .825 A= l. 15 + A= l.6 X A= 2.3 A= 3.2 2.75 3.30 3.BS 1L95 Figure 44. Pressure versus time behavior at fixed Eulerian radius from a blast system generated bv a Mach 0 . 5 energy wave. l.f1 r- D CD l.f1 ~ D er, w a:lf) =i- <.n . (!)- w cc (L D c, lf) w D MACH 0.25 ADDITICIN T R= .825 R= I.IS + R= 1.6 X R= 2.3 ~ R= 3.2 t- + -----.-----.--------r------,------,r------r------,--------y-------, .oo .60 Figure 45. 1.20 1.ao 2.1.rn 3.oo TIME Pressure versus time behavior at fixed generated by a Mach 0.25 energy wave. 3.60 11. 20 l{.80 5. l{O Eulerian radius from a blast system I-' +' 0 CJ ,.·,1 l!! 0 UJ u...: lf! :_.1 ,: 1 U) . ~ (J) u .l LL (L. lfJ O"> . C) <.D . C) ("(') . a 144 C)---'1!.w_-L__..L_-.'-t__--1_-;-_ ___.__ _ _,__-,-__.__....__'r-----'--__._-,__.__ _ _.__-+-__ _ o.oo .so 1.00 1.50 RADIUS 2.00 2.50 3.QQ Figure 49. Particle position versus time in a blast system generated by a Mach 5.2 (CJ) energy wave. .... Cl ::::r -N Cl - -N Cl co . - C) lf) . - ~~ ~ . h.-- 0 0) D tO 0 (T') . o.oo I i I I I ' I I ' ' '1 \ \ ', \ I ,, \ ' / / / / I / / ( / / .so 1.00 145 I I I I I \ I l I i / i I I .' I 1.50 RROJUS I I I I I I I I : I I I I ' I I 2.00 2.50 3,00 Figure 50. Particle position versus time in a blast system generated by a Mach 4.0 energy wave, 3. . -. ' ..... 0 w . 0 (Y') . 0 146 O -flUL-1.__..l_ __ .L__ j___ 'r---.JL_-...1.-,__._ _ __.__-\-_ --i...._~-,-~-..___t-__ _ . o.oo .so 1.00 1.50 RADJUS 2. 00 2.50 3,00 Figure 51. Particle position versus time in a blast system generated by a Mach 2.0 energy wave. . N tO lD . 0 C) .... lO lO . CT") CT") . a 147 0-1'-Llll_--1__..J.,._--r-..L_-..____'r----'---.----.-------.-----r----. 0,00 ,50 1.00 1.50 RADIUS 2.00 2.so 3.00 Figure 52. Particle position versus time in a blast system generated by a Mach 1.0 energy wave. LI') <.O . - 0 -. - o.oo 148 .so 1.00 2.00 2.50 3.00 1.50 RADIUS Figure 53. Particle position versus time in a blast system generated by a Mach 0.5 energy wave. . :::r c:, C\J :::r . - . - g . 0 .oo .so 1. 00 149 1. 50 RRDJUS 2.00 2.50 3.00 Figure 54. Particle position versus time in a blast system generated by a Mach 0.25 energy wave. 3.50 lJ') t::) . -::,0 l.n -. 0 t- . N IJ"') N - - - N- c:, co ·- - IJ"') ('f") ·- - 0 en - . I.I') '::J'_ . 0 0 . 0,00 I I I .so 1. 00 150 I I 1.50 RADIUS I I 2.00 2.50 3. 00 Figure 55. Particle position· versus time in a blast system generated by a Mach 0 .125 energy wave. I 3.50 30.0 10.0 l .o 0.04 151 ~VERPRESSURE VS RRDJUS Q=8 ~ BRKER IPENf~LJTEJ ~ BURSTJNG 5PMEA£ & MACH 8.0 RODJTJ~N + MACH 5.2 RDDJTJ~N X MACH 4.0 ADOJfJ~N ~ MACH 2.0 RDOJTJ~N + HACH 1.0 RODJTJ~N ~ MACH 0.5 RDDJTJ~N Z MACH 0.25 RDOJTJ~N I R~OIUS (ENERGY SCALED) l.O Figure 56. Overpressure versus energy scaled distance. 1.0 D w _j IT u en >- ~ O. l LL z w w en _j ::J Q_ L ,_... 152 JNPULSE VS RADJUS SPHEAJCAL GECNErAr BAI'\ ER Ip EN r C LJ r E) B.LJR Sr J NG SPHERE ~CH B.O ADDJfJ~N MA CH 5 . '2 ADO J r J ~ N MACH ~.O ADDJfJ~N M'iCH 2.0 ADDJfJ~N t,1gCH l.O ADOJfJ~N M'iCH 0.5 ADOJfJ~N t,1gCH 0.'25 ADDJfJ~N ~ERNEL ADDJfJ~N fAU=0.2 ~ERNEL ADDJfJ~N fAU=2.0 O.OJ_,__-.---.---.---.-----r-....----.---~---.-------~-~--,--..-,-~----r-~--- 1 I 0:04 0. l RRDIUS(ENERGY .SCRLEDl Figure 57. Impulse versus energy scaled distance. 1.0 bz a 0 51 , I 1/)_J r- ' I I I g l t;~1 5 j :z:o I.J..Jo ~~ )( l( )( l( )( 153 )( l( (IE+KE)Ball + (KE)Surroundings (IE+KE)Ball C, 0 1: L~ I ["- ! I lg Lei Ln I l o I ~ ~ In IN I i I ~~~~--9==8~(sK-E~~~B~1a~~-l~o~eo-~,o~~@-~o~,~o-~o~,,~o~-e~~e~1-eo-~o-~10--o-~o~,~o---~I~ 0 55 l 00 1 33 1 66 2.00 2.33 2.66 3 .CD Cl C, • OD • 33 . • TI ME . 0 0 ir~ ~ ~--i--------------------------+----.-------1----1--+---1>--- I 1-- Cl a :z: w u a:o w~ Q. - 1 0 (IE)B~a~l~l-+---+-r---+-----+--------+---+----+----+----+--+ ( IE) Surroundings (KE)Surroundings (KE)Ball 0 0 0 0 0 0 0 L_: I '::) 0 0 ~---,-----,------,----.----,---- -~-_,_--,-- --~ ~ ~ 4 ---- r-- .66 1.00 1.33 1.66 2.00 2. 33 2.56 3 . 00 . oo .33 TIME 8 Energy distribution versus time behavior in a Figure 5 · f generated by an in inite velocity energy wave blast system (bursting 0 0 0 0 0 0 L/) r- 0 0 0 LJ") >- Cl a: w :z:o wC::: :--.,:u-, 'N 154 (E)Total (IE+KE)Ball + (KE)surroundings (KE)Ball 0 0 0 0 - 0 0 L/) r- 0 c::: 0 L/) 0 0 L/) N 0 0 0 o ..... L--==~==~~:;::::::=~~::::;::===:j;t=::=:;:==Eil---.....;.--....--e----re----,~-- -+ 0 0 >- Cl a: 0 0 N 0 0 Wo zo Wei 1- :z: LU u a:o LJ...J~ 0..7 Cl D __J .00 .15 .30 .LIS . 60 .75 .90 TIME (IE)Ball ~ ~ IE)Surroundings ~ (KE)Surroundings (KE)Ball 1.05 l. 20 l. 35 0 0 0 0 0 0 0 ' 0 0 0 0 ~-.-------,.-- --,-- - - -.----~----.----,-----.------.-- --r~ . 00 .15 .30 .LIS .60 . 75 .90 1.05 1.20 1.35 TIME Figure 59. Energy distribution versus time behavior in a blast system generated by a Mach 8.0 energy wave. >-C) a: u...: Cl Cl Cl C, Cl Cl c:i Vl zo u.JC; ~Lu N 0 C 0 Cl Cl r-.i Cl C >- L'J a: u.Jo zo u.J c:i t- z u.J L ~ a:o u.JC; o....- 1 c..!) 0 ....J 155 .OD .30 . 60 .90 1. 20 TlME (IE+KE)Ball + (KE)surroundings ( IE+KE) Ba 11 (KE)Ball 1.50 1.80 2.10 2.lW ( IE)Ba l (IE)Surroundings (KE)Surroundings (KE)Ball 0 0 0 0 - 0 0 lf") r-- 0 ~ 0 l/"J 0 0 L,; N 0 C 0 2. 70 0 0 N 0 0 0 0 ci 0 0 -I Cl 0 Cl 0 ~ ~;____J__~--- ~ ----,--- --~- ---,----4,--,J-----,--_:,~~-.-----+~ .oo .30 .60 .90 1. 20 TIME 1.50 1.80 2. 10 2.llO 2.70 Figure 60 . Energy distribution versus time behavior in a blast system generated by a Mach 5.2 (CJ) energy wave. >-Q a: 0 0 . 0 0 0 0 \.n r-- 0 0 0 \.n LJ.J zo LJ.J ~ :--.! Ln • N 156 { ( E) total { (IE+KE)Ball + (KE)Surroundings { (IE+KE)Ball c:, c:, 0 0 .... 0 c:, \.n ,-.. 0 0 c:, \.n 0 0 \.n N 0 0 0 c:,~/4.e::B_~f:::;:~~=~~:::8::~-,-e----e-,---e---,.e----a-.--e--~'T"--e--~:l--- -t"o { (KE)Ball 0 .oo 30 .~o • 9o 1.20 1.so 1.80 2. 10 2.1rn 2. 7D . TIME >-Q a: 0 0 <">i 0 0 LJ.Jo zo wa 1- z u.J u a::o u.J t:: CL, Q D _j (IE) Surroundings (KE)Surroundings (KE)Ball 0 0 N 0 0 .... 0 0 ' 0 0 0 0 ~ .... 1..-iM--L--.------,-----,----,-----.----->-r--=---.-----.----t-~ .oo .30 Figure 61. , 60 • 90 l. 20 1. so 1. 80 TIME Energy distribution versus time behavior generated by a Mach 4.0 energy wave. 2.10 2.1rn 2.70 in a blast system 8 8 0 0 . VI r-- 0 0 0 "' 0 157 (E)Total (IE+KE)Ball + (KE)Surroundings (IE+KE)Ball ----- (KE)Bal l g 0 0 0 VI I' 0 0 0 VI 0 0 VI N o.._-.;,/IJ~...d:t;::::.=::::s......e:::~~::::!~~::&,,..[..--1----.:i----,-&----Q--r----€~-Er----e--r3---ro C .OQ ,30 ,50 .90 rhf~ 1.50 1.80 2.10 2.irn 2. 70 0 0 0 0 0 0 ;J (IE)Ball ( IE)Surroundings N 0 0 )- '...'.) X w.Jo zo w.Jc::i 1- z uJ u xo U-1~ G-7 0 0 (KE)Surroundings (KE)Ball 0 0 0 0 0 7 1-.Ju......,..~Hi-.J_---,-----,r----.----,----.------'-~---~-----r~ .oo .30 .50 .90 1.20 TIME 1.50 1.80 2. 10 2,70 Figure 62. Energy distribution versus time behavior in a blast system generated by a Mach 2.0 energy wave. 0 0 c::i 0 0 0 Ln C"" 0 0 c::i u, >- Cl a: lLJ z:o wC; :--,:u; "N >- a a: 0 0 c::i 0 0 c-,i 0 0 We z:o Wo f- z: w w a:::o wC; c..... j' a Cl ~ 158 .oo .33 .66 1.00 1.33 TIME (E)Total (KE)Surroundings (IE+KE)Ball (KE)Ball 1. 66 2.00 2.33 2.66 (IE)Surroundings (KE)Surroundings (KE)Ball 0 0 0 0 0 I.I') C"" 0 0 0 &.rJ 0 0 I.I') N 0 0 0 3.00 0 0 ..; 0 0 0 0 0 0 0 -I 0 0 N ...... -J;.. .... _,'"*°,-;---.+J.L__-r-----,----,------,-----,----,--'----€~,---- ,~ I 00 33 66 1.00 1.33 1.66 2.00 2.33 2.66 3.QQ • • . TIME Figure 63. Energy distribution versus time behavior in a blast system generated by a Mach 1.0 energy wave. >- c...!) a: 0 0 0 0 0 0 I.I) r-- 0 0 0 I.I) 0 0 N 0 0 - Wo zo wc:;; I- z LJ_J u a::o LJ_J~ a...-I c...!) El _J 0 0 N I .00 .55 .oo .55 Figure 64. 159 ~ (E)Total (IE+KE)Ball + (KE)Surroundings (IE+KE)Ball l. 10 1.65 • -q { I Et) Ba~ l & & & • & • t ( IE\urroundings (KE)Surroundings (KE)Ball l. 10 I. 65 2.20 TIME 2.75 3. 30 3.65 LI.Im Energy distribution versus time behavior in a blast generated by a Mach 0.5 energy wave. e 0 0 0 0 U'> r-- 0 0 0 U'> 0 0 U'> N 0 0 N 0 ~ 0 0 0 0 0 ' 0 0 N I ll.95 system >-c..::, C) C) C) C) C) C, U1 r-- C) C) C) U1 a: w zo w~ -...•u-, •'-N >- c..::, a: 0 0 N 0 C) LJ..Jo zo l.J.Jc:i ;..-- z w u cr: o w ~ Q_,. c..::, D _j 160 0 0 0 0 _. (KE)Surroundin s 0 0 (IE+KE)Ball (IE)Ball (IE) Surroundings (KE)Surroundings (KE)Ball l/'J r-- 0 0 0 l/'J 0 0 l/'J N 0 0 N 0 0 _. 0 0 0 0 0 -I C) 0 C) 0 N..__-4-,.4-._L~(,---1.-...-L----,----.----.----.-- ---.--~~ ..... ~--T ~ I QQ 60 1.20 1.80 2.110 3.00 3.60 11 . 20 1'.80 5, 110 . . TIME Figure 65. Energy distribution versus time generated by a Mach 0.25 energy behavior in a blast system wave. g c ::: ..., I I cl c , U') J r--' I >- ~ a: c:, 0 c-.i 0 0 LJ.J 0 zo LJ.Jc::i 1- z LJ.J u a::o LJ.J~ CLj" ~ D _J i 161 (IE+KE)Ball + (KE)Surroundings (IE+KE)Ball (KE)surroundings 0 0 0 0 1/'l ..... 0 0 N 0 0 0 0 0 0 0 ... I 0 O 0 O ~ ·~-E~~--,-,&*-e-+--.---+'~--,--*--*r.L.:...---,------,-----,------,----+~ . DO .45 Figure 66. • 90 1. 35 l. BO 2. 25 2. 70 TIME Energy distribution versus time behavior generated by a Mach 0.125 energy wave. 3. 15 3.60 lj,05 in a blast system 100 90 c, 70 t: 'E 'iii E (l) C: ?;; .... ~ 60 w '?fl. 50 Limit For Constant Pressure Expansion _________________ L::.:=1.;'~~~-: .... _____ -_______ ________ _ .5 1.0 2.0 1/Mw.---• 'Figure 6 7 , Energy Remaining in Source Volume Versus Reciprocal Mach Number of Energy Wave. 3.0 4 V. COMPARISONS A. An Investigation Of Blast Waves Generated From Non- Ideal Energy Sources Adamczyk(lB) systematically studied the flow field of blast waves generated by the homogeneous deposition of energy (infinite velocity wave of infinite thickness with finite deposition time). Using a Von-Neumann/Richtmyer-type finite difference integration procedure he generated numerical solutions of the flow field parameters for planar, cylin- drical, and spherical flow fields. In the analysis, Adamczyk determined the time of energy deposition and the energy density within the source to be the two most critical parameters affecting the flow prop- . perties of the blast system. Since the calculations in this dissertation were all done at one energy density, these comparisons will only address his conclusions concerning deposition time. For comparison with the homogeneous energy deposition inves tigated by Adamczyk, two cases with deposition times l n=0.2 and ln=2.0 were run. A shorter deposition time was not deemed necessary since the flow field approaches burst- ing sphere. 1. Flow Field Properties Kernel Addition ln=0.2 163 164 For the case of an energy deposition time of T=0.2, figures*68 and 69 show the pressure time history of the energy addition. As the energy is added, the flow field develops and an expansion wave propagates in from the kernel edge. When energy addition stops the expansion wave has progressed only about 50% of the distance into the kernel. Therefore, the sturcture of the system closely resembles that generated by a bursting sphere. Comparing figures 69 and 15, it can be seen that the flow fields are similar if the time for energy deposition is considered. The expansion wave has propagated 50% of the way into the bursting sphere at time 0.13 and 50% into the homogeneous energy addition at time 0.2. By adjusting the times for flow field behavior to reflect this time difference, the flow fields are similar. This can also be seen by examining figure 75. The energy addition results in an initial expansion wave followed by a second shock reflected from the origin, similar to the bursting sphere. However, the expansion wave does not pro- pagate into the source volume at constant velocity. The local velocity of sound is a function of the local temperature and gamma, both of which are functions of the energy addi- tion. Figure 71 shows that the center of the source volume experiences a constant volume energy addition, similar to the bursting sphere case. The cells on the edge of the source volume experience both a pressure rise and specific volume increase. -;..:c----- Note: f~gures in thia ch~~tcr arc collected at the end to sinolify comparison. 165 The blast wave structure at fixed Eulerian radii is shown in figure 73. Inside the source volume n = 0.825 the pressure rises during energy addition, peaking when energy addition ends. The pressure decreases below ambient at T=0.58. The blast wave behavior outside the source volume is similar to bursting sphere, figure 38. Kernel Addition T=2.0 In run 20 a homogeneous energy addition was done with a deposition time of T=2.0 which is quite long in relation to the characteristic times of the system. For the case of no energy addition an ambient temperature acoustic wave would take a time of T=0.85 to propagate from the edge of the source volume to the center. Since the energy addition takes place over a much larger time, the system distributes the energy as it is added. Figure 70 shows that during the energy addition an expansion wave forms which reaches the center at T=0.68 with a maximum pressure of P=2.7. (note: the travel time of the expansion wave is decreased by the effects of temperature and gamma on the speed of sound.) The expansion wave reflects from the center and an outward propagating pressure "hump" develops. Since the energy is being added slowly there is primarily a low pres- sure expansion of the source volume with the pressure wave propagating into the surroundings. When energy addition has been completed the specific volume of the cells in the source volume is approximately 7. 0 which approaches the spec.ific volume expected from a constant pressure expansion. 166 Figure 76 shows that since the energy addition continues during and after the arrival of the expansion wave at the center there is no second shock generated. The expansion of the flow field is a smooth continuous process. In the very slow homogeneous addition of energy, T=2.0, figure 72 shows very unique behavior is the p-v plane. Initially the energy addition results in a pressure rise in all cells. As the expansion wave propagates into the source volume the energy addition changes from a pressure increase to a specific volume increase. Expansion waves propagating through the source volume tend to equalize the pressure and the intersections of the p-v curves indicate equal temperatures in the source volume. At the end of energy addition the indi- vidual cells have expanded to a specific volume of approxi- mately 6. 75 at P = 1.1. The blast wave develops as a relatively slow pressure rise both inside and outside the source volume as shown in figure 74. The pressure remains greater than ambient throughout the energy addition. An interesting behavior in this case is that the pressure drops below ambient first in the source volume and then propagates outward. 167 2. Damage Parameters Figure 77 shows the overpressure from homogeneous addi- tion of energy. For the rapid energy addition (Tn=0.2) the pressure peak progresses from the edge of the source volume towards the center until a shock waves forms. As expected, the overpressure of shock approximates overpressure from the bursting sphere shock. For the very slow energy addition (T=2.0) the peak pressure propagates from the edge of the source volume to the center and out as a shock is formed. However for this case the overpressures are significantly lower. These overpressureslie between those of the Mach 0.5 and Mach 0.25 cases plotted on figure 56. This would be expected since the times for energy addition in these cases are 1.859 and 3.719, respectively. In Adamczyk's investigation the instantaneous deposition time produced the highest overpressures, whereas in the wave addition of energy the overpressures increase to a maximtnn at a finite time of deposition, Tn=0.28. Figure 78 presents comparisons of the overpressures developed in the wave addi- tion of energy and the homogeneous addition of energy. For the cases investigated the overpressure outside the source volume was greater than the overpressure from the homogeneous energy addition. However, in the source volume, as the deposition time becomes greater than T=0.6 the overpressure 168 in the source volume was greater for the homogeneous addition of energy than the wave addition of energy. This is not con- sidered to be controlling since the overpressure is low (P~S.O) and the area would be subjected to extensive fire damage. The impulse in the cases involving the kernel (homogen- eous) addition of energy is shown in figure 57. For the rapid deposition of energy T=0.2 the impulse is slightly less than bursting sphere in the near field and slightly greater in the intermediate to far field. In the near field the impulse is lower because of the time required for the energy deposition. In the intermediate and far field the impulse is greater because the finite time of deposition extends the positive phase of the blast wave. For the long kernel deposition time (T=2.0) the impulse is much lower because of the low peak pressures developed. The impulse curve lies between the Mach 0.5 and Mach 0.25 energy addition wave curves with Tn=l.86 and TD=3.72, res- pectively. 3. Energy Distribution With consideration given for the time of energy addition, the kernel addition with TD=0.2, shown in figure 79, is similar to the case of bursting sphere. However, during energy addition the energy appears as kinetic and internal energy in both the source volume and the surroundings. The internal energies approach final values of 65% in the source 169 volume and 36% in the surroundings. As expected kernel energy deposition of long duration, Tn=2.0, results in very inefficient energy transfer to the surroundings as shown in figure 80. The final distribution is 74.7% in the source volume and 25.3% in the surroundings. This can be compared to the energy distribution from a Mach 0.5 wave with 74.1% of the energy remaining in the source volume. This behavior appears reasonable, since for the Mach 0.5 energy wave the total time of energy deposition in the source volume is T=l.859. B. Some Aspects of Blast from Fuel-Air Explosives Beginning with the finite differencing technique of the "Cloud" program written by Oppenheim (30), Fishburn <35) added a burn routine similar to that of Wilkins(Z9) to simulate the detonation process. Using the program he studied blast waves generated by (1) centrally initiated, self-similar Chapman-Jouguet detonation, (2) edge initiated spherical implosion, and (3) constant volume energy release followed by sudden venting to the environment. Selecting MAPP gas,methyl-acetylene propadiene mixture, as a representative hydrocarbon, Fishburn used the "TIGER" program to calculate thermodynamic equilibrium for MAPP gas in the CJ plane. Using the calculated detonation pressure, the energy to be added and the detonation ·Mach number were calculated from the steady-state conservation equations (Equations II-36, II-37). 170 The energy was added linearly and garrnna changed propor- tional to the energy release through the front. Several runs were done varying the front thickness and a final wave thick- ness of 10% of the energy addition zone was selected. In figure 2 of Fishburn's paper the plot of a centrally initiated detonation has a constant pressure from the center to the edge of the source volume. This plot was based on known detonative behavior and not program calculations. Cal- culated pressures started near zero at the center and approach- ed the CJ pressure as the energy addi tion app1:i.oached the edge of the source volume. ( 36 ) This behavior is consistent with the results noted in this dissertation. Fishburn noted that the constant volume energy release produced lower peak pres- sures near the charge but slightly higher peak pressures than the centrally initiated detonation to radii greater than R/Rc=2. This behavior was also noted in . this dissertation. Fishburn also did an analysis of the energy distribution by determining the net work done by the detonation products on their environment by the following relationship; Work y~-1] Where Re is the initial radius of the change and Rf is the final radius of the source volume. His calculations showed the fraction of energy deposited transferred to the 171 surroundings to be: explosion = 0.378 high pressure= 0.336 In t his dissertation the fraction of energy released which is transferred to the surroundings as kinetic and internal energy was: Chapman-Jouguet Bursting Sphere = 0.385 = 0.361 The differences in the results may be attributed in part to the different technique used in the calculation. However, the results are comparable. The conditions calculated by Fishburn were used as input parameters for a run using the program modified by the author. Figure 81 shows the development of the blast wave with time. Figure 82 is a pressure-specific volume plot. The particles near the edge of the source volume exhibit Rayleigh line behavior during the energy addition and appear to tangent the insentrope. This indicates that for the specified con- ditions the results approach the expected results from a CJ detonation . C. Pressure Waves Generated by Steady Flames Kuhl, Kamel, and Oppenheim( 2l) studied the self-similar behavior of the flow field associated with flames traveling at constant velocity. Their study was directed to the steady-state condition the system attains when the flame propagation velocity attains a constant velocity. They did 172 not consider ignition, initial flame acceleration, or the pressure wave decay after the source volume is consumed by the flame. Introducing reduced blast wave parameters as phase-plane coordinates, they determined appropriate integral curves on this plane. For one of their calculations they assumed a combustible mixture with a specific heat ratio of y =1.3 0 ahead of the flame and y 4=1.2 behind, a volumetric expansion ratio of 7 for a constant pressure deflagration, and an am- bient sound speed of 345 m/sec .. For comparison these parameters were used as input variables in the program used for this dissertation. The results are plotted as figure 83. In their analysis the flame was treated as a steady deflagration and a piston expanding at constant velocity was used as a representative case. Using subscript p to denote parameters corresponding to the locus of states at the piston face, solutions were obtained in terms of s = r1 zp as the parameter. By integration of the governing differential equa- tions, the solution for a spherical flow field is: zP = z F2/3 2 2/3 = [(t/rµ)a] [(t/rµ)u] zP = 5.67 s = ylZP 1: = (1.3) (5.67) = 7.37 173 An examination of figure 83ashows the blast wave ap- proaching a self-similar solution with a nearly linear de- crease in pressure from a pressure of 1.26 at leading edge of the flame front (X=0.42) to 1.02near the shock front (X=0.95). Comparing this to figure 7 of Kuhl, et al., the ~=4 curve has a nearly linear pressure decrease from a pressure of 1.28 at X=0.45 to 1.02 at the shock front. Thus the finite differencing technique assymptotically approaches p0 but the similarity solution appears to begin to develop a shock front at the leading edge. In figure 83b the energy transfer ahead of the flame can be seen. The calculations appear to be approaching a self-similar solution ahead of the flame front with a near linear decrease to 0/0 =1.005 at X=0.95. Figure 83c shows 0 the particle velocity in the blast wave. 'From a maximum velocity at the flame front it asymptotically, decreases to zero at the shock front. Through the flame front it decreases rapidly and remains at nearly zero. Comparing these results to the results of Kuhl, et al., the maximum values calculated with the finite difference technique at the flame front approach the values calculated 174 by Kuhl, et al., for the ~=4.0 case. However, the blast wave structure is more closely approximated by the ~=7.0 case. D. The Air Wave Surrounding an Expanding Sphere The properties of the flow field generated by a sphere expanding at a velocity, slow relative to the ambient velo- city of sound, were determined by Taylor<3). He integrated the velocity potential equation and developed the following relationships for the pressure and particle velocity dis- tribution outside the expanding sphere: p-p 0 = u = 2 2 M 3 P a s (l-Ms2) (at_ l) r 2t2 (~ - 1) r where Ms is the Mach ntnnber of the surface of the expanding sphere: and R(t) is the Eulerian position of the sphere. The results calculated in case 10 (~ = 0.125) were analyzed and compared to predicted results from Taylor's formulas. The leading edge of the energy wave was used to represent the surface of the expanding sphere. After the self similar flow field developed the Eulerian velocity and Mach number of the energy wave were calculated to be Ms=0.24. 175 Using this velocity the pressure and particle velocity dis- tribution were plotted in figure 85 for comparison with the R results calculated by Taylor for the case at= 0.2. The distributions are nearly identical indicating close agreement of the results calculated with theoretical predictions. 10 .0 7.0 p 4.0 1.0 3.33 2.67 2.00 1.33 0.67 PRESSURE/ PO OISTA!BUTI~N VS. DISTANCE/ 00 AND T!ME / TO Figure 68. Pressure distribution versus Eulerian distance and time from a blast system generated by a TD= 0.2 homogeneous energy addition. a.ao 3.0 2.33 p l.67 1. 0 2.0 -~-__ / 6.0 -5.0 4.0 3.0 2.0 1.0 PAESSUAE / PO DI5TAIBUTI~N VS. DISTANCE/ DO AND TIME/ TO Figure 70. Pressure distribution versus Eulerian distance and time from a blast system generated by a TD= 2.0 homogeneous energy addition. a.a LJ = = . en C) C) CIJ C) C) c-- C) D (.D D D If) C: :JC) no n · LJ :::r C: L 0 D (1") D D N D D 0 D D 0 . 00 1.00 Figure 71. 2.00 179 CELL F-V BEHRVI~R KERNEL TRU=0.2 SPHERICRL GE~METRY 12JBEGIN ENERGY ROOITI~N (!)ENO ENERGY ROOITIDN A CELL 1 + CELL 10 X CELL 20 ~ CELL 30 ~ CELL 40 )<: CELL 50 Z CELL 60 3.00 4.00 5 . 00 SPECIFIC VOLUME Pressure versus specific volume behavior h ( = J:I 6.00 from T 0.2 ,::,, 1.1") r- D Ln Ln N N 0 0 N Ln r- w 0::: :=lo (j") If) (j") • w ...... [C (L If) N 0 0 If) r- . 0 180 CELL P-V BEHAVIOR KERNEL TRU=2.0 SPHERICAL GE~METRY ~BEGIN ENERGY RDDITJ(jN C9END ENERGY ROOITI~N &- CELL 1 + CELL 10 X CE LL 20 ~ CELL 30 ~ CELL 40 :><'. CELL 50 Z CELL 60 Lr; --t------.-----.--------,------,------~ ---~ .50 1. 75 3.00 LL25 5.50 6.75 8.00 SPECIFIC V(jLLJME Figure 72. Pressure versus specific voltnne behavior = 2.0 0 In c:i C) 0 oi 0 In r-= c:, C) to w CCc, ::::Jin en . (n-::r w a: a... c:, c=: er, c:, In c:, c:, C> .00 - -·-·-- ------- -,----·r--- ·--.,- .15 .30 .~s .so TIME KERNEL, TAU= 0.2 T R= .825 R= 1.15 + R= 1,6 X R= 2.3 ~ "= 3.2 .90 1.05 1.20 ' 1.35 Figure 73 . Pressure versus time behavior at fixed Eulerian radius from a blast system generated by a TD= 0.2 homogeneous energy addition. t-' CX) t-' 0 1/) ~ 0 0 ~ 0 U) N 0 0 N w CI:o ~U) (f) • (f)- w a: n.. 0 0 . 0 U) 0 0 KERNEL. TAU= 2.0 ? R= .825 R= 1.15 + ff= 1. 6 X R= 2.3 ~ ff= 3.2 o--t---- - , --- - ---r---- --,- -- ---,---- - -,- -----,---- -,-----,- - --~ .oo .35 Figure 74. .70 1.05 t.tW TIME I. 75 2. 10 2.45 2.BO Pressure versus time behavior at fixed Eulerian radius from a blast system generated by to= 2.0 homogeneous energy addition. 3. 15 ..... Q) N lJ") ('I) . - - L.n -. o.oo .so 1.00 183 1.50 ARDJUS 2.00 2.so 3,00 Figure 75. Particle position versus time in a blast system generated by a TD= 0.2 homogeneous energy addition -- ' - • o.oo .so 1.00 183 1.50 RAD JUS 2.00 2.50 3,00 Figure 75. Particle position versus time in a blast system generated by a t 0 = 0.2 homogeneous energy addition. 3.50 184 . ....... ~-llLUl_...i__-L-,-Ji___--1._-+---L._...1._-r-L.._--t.._ -+--l_-.L.-,--,..____,__-'r __ _ . 0,00 . so 1.00 1.50 RRD1US 2.00 2.50 3,00 Figure 76. Particle position versus time in a blast system generated by a t 0 = 2.0 homogeneous energy addition. so.a 10.0 w cc ::J c.n c.n w a:: Q_ cc w > 0 1. 0 .0 . 4 0.04 185 BVERPRESSURE VS RADIUS O=B ~ BRKER C) BURSTING SPHERE KERNEL ENERGY RDD1IIDN & OEP~SlTJON TIME o 2 + DEP05IT10N TIME 2:o R~J1us (ENERGY SCALED) 1.0 Figure 77. Overpressure versus energy scaled distance. 10.0 1 1.0 .1 .1 186 Constant velocity Wave 1.0 To Homogeneous Energy addition [ Maximum overpressure in Source Volume [ Maximum overpressure in Surroundings 10.0 Figure 78. Maximum Overpressure vs Source Volume Deposition Time. 8 § 0 0 .,; ..... 0 0 0 Lil >- e:, CI: w zo w"! :,...:~ 0 0 0 g "' 0 "! - >- e:, a: We zo We I- z w u a:::o w"! a..-I e:, C ...J 8 "' I 187 .00 .15 .30 .115 .60 • 75 TIME (E )Total (IE+KE)Ball + (KE) Surroundings ( KE >ea 1·1 .90 I.OS 1.20 (IElsurroundfngs (KE)Surroundings B ~ 0 0 .,; ..... B 0 In 0 0 .,; N 0 0 0 1.35 0 0 - g i g · ... 4-.--- .,.-----.----r----,-----,-----,-------.---~--- +-1 .oo .IS Figure 79. .30 .60 TIME • 75 .90 1.05 1.20 1.35 Energy distribution versus time behavior in a blast system generated by a • 0=0.2 homogeneous energy addition. >- L'.) Cl Cl Cl Cl Cl Cl in r- Cl Cl Cl in a: w zo w ": ..... • in ''-N 188 0 0 0 0 -(E)Total (IE+KE)Ball + (KElsurroundi 98 (IE+KE)Ball 0 0 0 0 U'l 0 0 U'l N 0 (KE) Ba 11 ° Cl 0 0 --.--~~~-~-....... --.t.---....----li,-......,..--....,..-~4._-e--,-e--€>--Gr-S----€.-~-e----ro >- L'.) a: Cl Cl ,_; 0 Cl Luc, zc, Wei 1- z w u a:o w": CL- L'.) 0 _j ' .00 .35 . 70 1.05 1.1.lO 1.75 2.10 2.llS 2.80 3.15 .00 .35 Figure 80 . TIME (IE)Ball (IE)Surroundings (KE)Surroundings .7o 1.os 1.1rn 1.75 2.10 2.llS 2.80 TIME 0 0 N 0 0 - 0 0 0 0 0 -I Energy distribution versus time behavior generated by a Tn=2.0 homogeneous energy in a blast system addition. 0 0 . 0 N 0 0 . (.D - 0 CJ . N ...... w CC:CJ =:)o Cf) • Cf) co w cc Q_ Cl 0 . =1' CJ Cl 0 \0 ,- N -=:t . ,- ,- . . co ,...... . I \ 0 N \ . ' 0-+-- ----,-------.-------.-------,--------.--------,,--------,--------, .00 .20 . • l!Q .60 .80 1.00 1.20 1.1.lO 1. 60 RADIUS Figure 81. Pressure distribution versus Eulerian distance and timt Lt~ a 01ast system generated by a Mach 5.55 energy addi t ion wave (Fishburn's case study). ~ 00 '° D Cl . co ...... G 0 lD ...... 0 D ::I' ...... 0 D N .-I 0 0 0 ..... w a: =:Jo (f)o (f) • wro a: o._ 0 0 . (D 0 0 . :::1' 0 0 N / 190 CELL P-V BEHRVJ~A FJSHBUAN ENERGY RDDJTJ~N WAVE vU D f H O. l SPHEAJCAL GECINEfAT [!] BEGJN ENER GY FlODJfJLlN OJ EN D ENER GY AODJ f J LIN Ji CE LL l ~ CELL l 0 + CELL 2 Z CELL 20 X CELL 3 y CELL 30 ~ CELL 4 J:!l'. CELL LlO ~ CELL 5 * CELL 50 & CELL 60 C) c=: -L--- -.-----,----~-----r-----r----.:i 0 50 3. 00 0. 00 , 50 1. OD 1. 50 2. 00 2 • SPECIFIC VOLUME Figure 82. Pressure versus specific volume behavior Mach 5. 55 ener . wave (D = l ,. 0. at eel from a 0 ,.., 0 ~ 0.0:J ~ ~l :2 - u, wo er:~ ::::l ..... a: er: i..LJ CLc:, 40 W· ..... - u, I , 12 191 .25 .37 .50 . 62 .75 . 67 I.DO RROIUS (a) \ I ai,-----,-----.----....-----,---~---...------r----, >- ..... 0 I/") u, "' --o uo :30 i..LJ > LU d~ ;: , .. er: a: CL 0 u, 0.00 . 12 .25 .37 .50 .62 . 75 ,S7 I. OD RROIUS (b) ,· +--- - -,--- -...-----....- ---.--- --.-----,---- --,-----, 0.00 .1 2 Figure 83. .25 .37( ) .so c RADIUS .62 .75 .67 I.OD Pressure, temperature, and particle velocity versus similarity radius from energy addition 0 :::r • _j - 1./) ('l") i • _J - 1./) N . ...., ...... LJ....J r ~ o CFJN en · -J u...1--' .r ,)_ If) ....... ' . _j ....... 0 ....... i • _J _ ... tJ) 0 ·~ ......... I \ j 192 CELL P-V BEHRVJ~A I KUHL. ET AL. ENERGY ROOJTJ~N WRVE WJDTH 0. 1 SPHEAJCAL GE~METRY QJ BEGJN ENERGY ADDJTJ~N cp END ENERGY ROD J T 1 ON ~ CELL 1 :X: CELL l 0 + CELL 2 z CELL 20 X CELL 3 y CELL 30 <'> CELL 4 ~ CELL 40 Lf' CEL L S * CELL SD Z CELL 60 g . . . J ______ _ 0.00 1.00 2.00 3.00 tLOO 5.00 6.00 SPECIFIC VOLUME Figure 84. Pressure versus specific volume behavior from Kuhl. et al., case energy wave (D = 1.0 at cell 50). u - a 193 r---------------------- 0.1 0.08 Taylor R - =020 at · 0.06 0.20 Mw=0.125 0.04 0.02 0.10 0.00 0.00 L-___.JL--_J,.--'--..J--.L-_.JL-~-...:::C==-....... __.J 0.4 0.6 0.8 1.0 0.0 0.2 r at P·Po Po Figure 85. Comparison of pressure and particle velocity distribution of Mach 0.125 energy addition wave with results predicted by Taylor. VI CONCLUSIONS AND RECOMMENDATIONS In conclusion, this dissertation presents a systematic theoretical study of both the near and far field effects of constant velocity flames. Earlier studies included only the development of a self-similar solution during energy addition. None of the previous studies included blast wave behavior after the end of energy addition. In this dissertation, the non-steady, one-dimensional fluid dynamic equations of motion with divergence and energy source terms, subject to the appropriate boundary conditions, were integrated using a Von Neumann/Richtmyer - type finite difference numerical integration procedure. The calculations yielded the thermodynamic changes and fluid-dynamic be- havior associated with the propagation of the blast wave. Particular attention was directed to changes in peak pres- sure, positive impulse and energy distribution. In parti- cular the relationship of non-steady behavior to steady- state behavior was noted. A. Conclusions On the basis of this investigation the following conclusions have been reached: 194 195 1. Near Field Behavior For Methane a. In assessing potential damage from non-ideal ex- plosions, preliminary estimates can be made from the Values predicted by steady-state theory. For the cases of super- sonic combustion from CJ velocity through infinite velocity (bursting sphere) the overpressures symptotically approach the values predicted by steady-state theory. b. As the flame velocity decreases from infinite vel- ocity, even through velocities impossible by steady-state theory, the pressure increases to a maximwn of P~20.0 at a Mach number of 4.0. 1.) As the energy wave velocity increases above 4.0 the flame is moving so fast relative to the expansion behind the flame that the reinforcement of the pressure decreases. 2.) For flame velocities below Mach 4.0 a signi- ficant amount of the energy is taken up in the expansion through the flame front. This results in a decrease in the peak pressure as the Mach number decreases. For a 50% decrease in the wave velocity the following relationship holds: overpressure(50% velocity)= 0.35 {overpressure(l00% velocity)} 3.) For flame velocities much less than the am- bient velocity of sound the pressure ·and particle velocity distribution closely match the results originally predicted by Taylor (l 3) 196 2. Far Field Behavior For Methane a. In the far field, the overpressure for all super- sonic flame velocities approach 65% of high explosive at equivalent energy scaled radius. b. At subsonic flame velocities the overpressure is significantly less than either the high explosive or the supersonic energy addition. When calculations were terminated, the Mach 0.5 case had reached 84% of the supersonic overpres- sure and the Mach 0.25 case had reached only 23% of the supersonic overpressure at n=lO.O. 3. General Observations a. For equal source volume deposition times the wave addition of energy produced greater overpressures than the homogeneous energy addition. This is attributed to the propagation of the energy addition wave interacting with the fluid dynamics of the flow field to develop greater over- pressures. In the homogeneous energy addition there is no reinforcement of pressure. b. In cases where the flow should reduce to a self- similar solution and/or show Rayleigh line behavior it did. The calculations showed that the flow field behaved normally where expected, and in the forbidden region, where steady- state behavior is not expected, non-steady behavior was observed. c. Maximum energy transfer to the surroundings from 197 the blas t process occurs at a flame velocity of Mach 4.0, corresponding to the maximum overpressure in the flame. 1.) At flame velocities greater than Mach 4.0 the energy transfer to the surroundings decreased to the energy transfer associated with constant volt.nne energy addi- tion (bursting sphere) in the limit. 2.) At flame velocities less than Mach 4.0 the energy transfer to the surroundings decreased, approaching the energy transfer predicted for constant pressure deflagra- tion. d. For the energy density investigated, q = 8.0, the use of ideal (point source) theory results in an overestima- tion of the damage potential of these explosions. e. In as much as the energy density, q, of most hy- drocarbons are all approximately equal, the conclusions reached can be applied with reasonable confidence to other gases and flammable liquids having energy densities in the range of 6 000027 000 30 FORMATC lHO,•NON-PnSITIVE TIME STEP•J 000028 000 31 FORMAT(• '•'AT TIME 't1PE12,&,• THE ENERGY INTF'GRAL F'QIJALS •, gg~g~z ggg 321 FORMAT(• ~!~1~t:,:1~~s~iA~,; CELLS•) 000031 goo 3j FoRMAT('l'•//) 000032 00 34 FoRMAT,lHO,• LW wCM WAA WSLSOR WSREXP WCMt•,1I5,5D20.tOJ 000033 000 35 FORMAT 1HO,• THE TA TAP TERMIN 1 ,/ 4E?0.10I 000034 000 3b FORMAT 7A4J 000035 000 37 FORMAT(• '•'224'•T5•11E10o41 000036 000 38 FORMAT(• •,• 205' 0 I5,11E10•41 000037 000 39 FORMAT(' '•'M nT OTMIN 000038 000 l 'OEM R2MP R2M OOOOJ9 000 40 FORMAT(I5•11E10o4) 000040 000 41 FORMAT(• •,•234••~15,7E15,9/8Et5o9J 030041 000 42 FoRMAl(' '•'10A'•~I5,2E10t4I O 0042 000 4j FoRMA (' '•'186'•T~•5E20. 4/6E20ol4) 000043 000 44 FORMAT<• '•'204'•T5•6E18.12) 000044 000 4~ FORMAT(• '•'247••~1~,4E20ol4) OTTfS AS2 GRADTA AS1 AooTN 000045 000 4b FORMAT(' '•'300'•~I5,6E151ql 000046 000 47 FORMAT(' •,•305•,,5,sr20. 41 000047 000 48 FORMAT(• t, 'LAGRAt..1GIAN POSITION OF HEAD IS ',Ft2o6• •• AS2' I 000048 000 l • ANO TArL IS •,F12.61 000049 000 49 FORMAT(• •,•THE 5PECIFIC VOLU,iiE OF CFLL•,15,, IS NEGATIVE•! OOO U50 000 50 FORMAT(• •,3I5,7E14o8l O 51 FORMAT(' •,•STORAr.E FILE OVERFLOW') gggg5} 808 52 FORMAT(• •,••••••••••GAMMA CHANGES AT WAVE FRONT•••••••'I gggg~~ ggg ~~ ~g~~:+i: ::~~~Yt1~6~r 11~sl~~rL1TY 1N TIME sTEP•, 000055 000 5~ FORMAT(' ',215,7E15o8 000056 ooo 56 FoRMATC' •,•142•,,s,6E15.717Ets.71 000057 000 WRITE 6,33J 000058 000 C•••• DETERMINE INITIAL OR RESTART CONDITIONS 000059 000 CALL MTIMECXCP11,XMEMI 000060 000 TA:XMEM 000061 000 INnEX: 0 000062 000 WUN: 1.00 000063 000 ZERO: OoDO 000064 ooo Two= 2.00 ogoo65 000 THREE: 3.Do o oo&6 ooo Elr.HT = a.uo 27 N 0 V, •••••• AMAIN 000067 000 000068 000 000069 000 00001~ 00007 goo 00 000072 000 000073 000 00007'+ goo 000075 00 000076 000 000077 000 000078 000 000019 000 ogoo8r 888 0 001:1 000082 000 00008J 000 00008'+ 000 000085 000 000086 000 000087 000 00001:!8 000 00001:19 000 000090 000 000091 000 000092 000 000093 000 000094 000 000095 000 00Q09b 000 000097 000 000098 000 080099 300 0 0100 00 00010! 000 00010 000 000103 000 000104 000 000105 000 000106 000 000107 000 000108 000 0001~9 0 00 0001 0 goo 000111 00 000112 000 000113 000 000114 0 00 00011i 000 0001 · 000 gsg1u 000 000 000119 000 000120 000 000121 000 000122 000 000123 000 OOn124 000 00·0125 000 000126 000 000127 000 000128 000 000129 000 000130 000 OOOlj! 000 0001 000 OO.o13J 000 000134 goo 000135 00 0100136 000 I 000137 000 000138 000 ogo1J9 000 0 OllfO 000 00014! 000 00011f 000 •••••• C C C SIXTEE: 16,DO ENIN(: 9,00 1c;r=o Pl: DACOS(-WUNI R~AD15,24)LSTART•nTRACE,OTAPE OSTART: LSTART •FG, 0 IFIOSTART) GO TO 57 CALL RESTAH Go TO 513 !,7 RFAD( 5,?b ) NrYCLE, NPUNCH, NSTORE",NS RFAUC5,J6) LABEL CALYNI~lUk11 IF(,NOT, OSTART)GO TO 59 T:Zl.:RO NPEAK:t CALL PUOAl T:OT IsT=IST+N 59 WRITE 16,52) LSTART: LSTART + 1 NFRPR:LSTART IF(OTAPE)NFRPR-20 TAP: TA+ 30, - THF = TEHMIN- Two IFIOTRACE)WRil~(16,35JTHE,TA,lAP,lFRMIN IF(THE ,LT, TAI TME: lAP OSPHER: J ,EO, 2 OPLANE: J ,f.Q, 0 OPRINT: ,FALSE, ONOVIS: CL ,LF., Z~RO ,AND, CO ,LE, ZF.RO OPIJN: ,FALSE, OPIIT I : TI PUN , GT, ZERO OSKIP: NS ,LE, 0 OAODEN: ,FALSE, OFE=,FALSE, OF.ENO: ,FALSE, OEXIT: ,FALSE, OVFR: ,TRUE, NPR=O NNFIN:201-11 IFIN ,GE, NNFIN)GO TO 3)5 lF(NFINAL ,GT, NNFIN)NF NAL: NNFIN IF I ,NOT. OESnR ) GO TO 84 ENSTEP: ENMAX / ln0,00 DO 81 L: 1,10 AA2: AAl AA!: AA CM2 : CJ'Jll CMl: CM IF ( L - 2) 7~, 74, 76 72 AA: OLOG I SLOSOR +WUN+ WUN I SOREXP Go TO 77 71f AA: .qsoo • AA GO TO 77 76 AA =AA1+(SLOSOR-C~ll•IAA2-AA1)/(CM2 77 CM=DEXPIAA>-WlJN+lnEXP(-AA/SOREXP)-WUN) / SOREXP • IF ( OARS (( CM - SLOSOR ) / SLOSOR) ,LT, ,.oo-7 • 60 TO f\2 81 CONTINUE 82 IF(,NOT, OSTART)GO TO 83 E2CL:ZERO 83 SORSPA: MINCOS - MAXCOS 84 IF( .NOT. OEWAvl GO TO ql Do 89 LW: 1,,0 WAA2:WAA1 WAAl: WAA wcM2: WCMl wcM1: WCM IF(LW - 2) 85•A6,87 85 WAA: OLOG(WSLSOR+WUN+WUN/WSREXP) GO TO 88 8& WAA: ,9~DO•WAA nHE 011477 -CMl) 1 N 0 (j\ ~ -··-·· A~AIN ogot43 soo 0 0144 00 O·Oo145 000 00014~ 888 00014 000148 000 000149 000 000150 000 000151 000 000152 ooo 00015J 000 0001s" 888 ogo1 5 o 0150 000 000157 000 000158 000 000159 000 000160 000 000161 000 000162 000 00016J 000 000164 ooo 000165 000 000166 000 000167 000 000168 000 ogo1~9 888 0 01 0 000171 000 oooI72 000 00017J 000 &88Hi goo 00 000176 000 OOOlH 000 ogo1 000 o 011z 000 00018 000 000181 000 000182 000 000183 000 000184 000 000185 000 000186 000 000187 000 000188 000 000189 000 ggg1~r 000 000 000192 000 000193 000 000194 000 000195 000 000190 000 0001~7 000 0001 8 000 000199 000 000200 000 000201 ooo 000202 000 00020.3 000 000204 000 0002oi 000 00020 000 000201 000 000208 ooo ggg~~z 888 000211 000 000212 000 00021.3 000 000214 000 000215 000 0002H, 000 0002u 000 0002 ooo •••••• 87 68 89 90 C••• C 91 CC•• ~ ... C GO TO 88 WAA: WAA1+CWSLS0Q-WCM1l•(WAA2-WAA1I/IWCM2-WCM1) WCM : Of XP(WAA )-W11N+ (Of:XP (-WAA/lliSREXPI-WUNI IWSR[XP lF(OAAS((WCM-W~LSOR)/WSLSOR) .LT. 1.00-11 Go ·ro 90 coNTJNUE WSRSPA = MNWCOS-Mxwcos lF(OlRACE)WR1TFl16,J41LW,WCM,WAA,WSLS0R,WSRFXP,WCM1 TwlO = WIOWAV/wVEt . ENMAX: F.NWAV - SLOSOR: WSLSOR TMAXE : TWID SoRExP = WSREXP ENST~P: ENWAV/10n,DO s1wuwv = w10wAv AA : WAA PHI: lERO lF(NCYCLE ,GT, OIGO TO 91 WT1 : T • WVEL IF(WTl ,LT, R(1,911N: 10 •••REMAINDER nF THESE WAVE CALCULATIONS AT 204••••• IF ( .NOT, oPUTI I GO TO 9~ M: T/TIPUN TNFX: M+ 1 TLINE: TlPUN•TNEx SET INDEX NUMAER ••• CALCULaTIONS FOR NEW TIME STFP •*•••••• C 92 INDEX: INDEX+ 1 C ••••• ~Al~~EslA~l[~~v t~1 TIME STEP ••••• IF(,NOT, OEXITl GO TO 95 C 91J 9~ wRllE(6,54) GO TO 312. rHECK CENTRAL PROCESSOR ELAPSFO TIME CALL MTIME(XCP11,XMEM) TR= )(MYM l~N;T!~TC,5, lo: THE - TB A: TB IF ( TCN,LT. Tn) GO TO 99 WRITE (6,271 GO TO 312 C ••••••••• CHECK FOR STnRAGE OVERFLOW••••••••••• 99 1F(0TAPE)GO TO 10~ lFllST ,LT• 10n00)60 TO 103 C••• 103 WRITE l 6, 51) Go TO 312 CHECK NUMRER OF STEPS IF ( INDEX ,LE. NSTEPS I GO TO 107 WRITTE(6,29) GO O 312 C 106 CHECw FOR DIMENSION LIMIT STOP OR MESH FXPANSION 107 NO: N NM2=N-2 108 109 110 111 OFLP: IP(t,NM2)-p(1,N))/TWO - DO 111 J: NM2•NFINAL PAB5: DABS(P(t,I1-WUN) IF(OTRACE)WRIT~(16,42ll•N,PABS,V(1,I) IF ( PARS ,LE• 9,0D-14) GO TO 115 N : I + It NL: N-1 PNL=oAAS(P(l,NL)-wUNI PN= OARS(P(l,N)-W11NJ IF(PNL .LE• 9,no-17)GO TO 108 GO TO 109 IF(PN •LE• 9o0n-17JGO TO Ill DO 110 J:1,2 UO 110 M:1•201 WR JTE ( 6, 55J M, J ,P ( ,Jt M), V (J, M), R ( J• MI ,U(JoMJ, 000276 000217 00027t1 000279 000280 000281 000282 00028.3 000284 000285 000286 000287 000288 000269 030290 0 0291 000292 000293 000294 000 000 0 00 0 0 0 0 0 0 00 0 000 00 0 000 000 000 0 00 0 00 888 00 0 0 00 0 0 0 00 0 000 000 000 000 000 000 000 000 000 000 00 0 000 000 888 000 000 000 000 000 888 000 000 000 000 (100 000 000 000 000 000 000 000 000 000 0 0 0 000 000 000 000 000 000 000 000 0 00 000 000 000 000 000 000 000 000 0 0 0 000 0 0 0 C 1114 C 1114 11~ llb 117 C••• C••• C••• ~ C C 127 C••• C••• C••• Go TO .312 DFTEpMINE PROPERTlfS AT N[W TIME NL= N-1 IF(N .EQ. NOi r.0 TO 117 DO llf, 1 : NO• M P(l•II = P(l,I-11-DELP IF(P(l,I) .r,r. WUN) GO To 116 P( 1 •II : WUN GO TO 117 CON r INUE DTN: ( OT+DTLI/TWn MOMENTUM CnNSERVATION Efrru~M~ ~'h~ ~b~8rHA~~ ANO TRAJECT0RIF'l IF(OT1RACE)WRITE(lg.~2)NCYCLE,1NOFX,NL,N,I,N~lEP5,T,nT,nTL,DTN,PARS Ul2• I : UL FIRST SJGMa FORWARD SJGMAF: -P 1,11-nll,1) FIRST Fl Fl: (R(l,21-R(l•tll/V(l,1) R(2•11: H(l,11+~ •UT RlMP: RIJ,21 Do llt2 M: 2tNL VELO,-IT'( AT GENERAL POHIT INTERMEOIATE SIGMA VALUES SJGMAB: SIGMAF SJGMAF: -P(l•Ml-~(1,M) PHI PARAMETER FO: Fl HlM: RlMP RlMP : RU,M+l) Fl: (RlMP-RlM)/V1l,MI PHI: ( Fl+FO I/ TWO VELOCITY U2M: Ull,Ml+(OTN1PHl)•(SIGMAF-SIGMAB) U(2•MI : U2M TRAJECTORY R(2•M): RlM+u2M.oT IF (OTRACE l WR I h· ( 16, 56)M•R ( 2,M) ,R1'4,lllM,U( 1 ,M) ,PHI, 1 .. 2• T SIGMAF,SyGMAB,Fl,FO,RlMP,RlM,V(l,M),R(2•1) CON INUE C••• RIGHT BOUNnARY CONDITIONS U(2•NI : UR R(2•N): R(l,Nl+yR•Ol ~::: ~g~[.~t~{E y ,J~N~~~~f1fgNVOLU'4ES AASEn ON CONTINUITY C 152 CONTyNUITY u13 = u,2,1,.u,2,,,.u,2.11 IF(.NOT. OESOR)GO 10 160 IF(OENDJGO TO 1~9 ElCL: E2CL PHI: T / TMAXE• aA E?CL: ( ENMAX / ~LOSOR) • (( DEXP ( PHI ) - WUN) + • (UEXP( -PHI* SOREXP) - WUN)/ SORFXP) OFND:T oGT. TMAXE OAODEN: .TRUE. IF(OENOIE2CL=E~MAX DFLlNG: E2CL - EtCL GO TO 160 1590 OFE=.TRUEo 6 . IF ( .NOT. O~WAV) GO TO 163 IFIOfENDl GO Tn 163 WIDWAV: STWDWV RRHlAO:T•WV~L + R~F IF( RRHEAO .LT. WIDWAV) WIDWAV: RRHfAD RRTAIL: RRHEAO - WIOWAV PHIW: ZERO HLHEAD: RRHEAn/R~F MwHEAO:RLHEAO RLTAIL: RRTAIL/R~F MwTAIL: RLTAIL MWTMl: MWTAIL - 1 11ATE 011477 Pf\Gf N 0 (X) •••••• 00(1295 000296 000297 000298 000299 000300 000301 000302 00030j 000304 000305 00030t> 000307 000308 000309 000310 000311 000312 000313 0003h 000315 00031b 000317 888~1~ 000320 030321 0 0322 ogo323 0 0324 000325 88032b 0327 000328 000329 000330 000331 000332 000333 000334 ogo335 0 033b 000337 000338 000339 000340 000341 080342 0 0343 000344 000345 ooo3'4b 000347 000348 000349 000350 000351 000352 000353 030354 0 OJ!>!, 000356 000357 000358 000359 000360 0003bl 001)362 001)363 000364 000365 000366 000367 000368 sssin AMAIN 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 800 00 000 888 800 00 000 888 000 00 0 000 000 000 000 000 000 000 000 000 000 000 000 800 1) 0 000 000 000 000 0 0 0 000 00 0 000 000 0 00 888 000 000 000 000 000 00 0 000 000 000 000 000 000 000 888 •••••• IF(MWTAIL •EQ• MWHEADJ MWTAIL: MWTMl lF(MWTMl oLFo MXWCOSJGO TO 163 MWHl::AD:o MWTAIL:O 16, OfENo: oTRUEo ,J Fl = U(2,ll IF I .NOT. OPLaNE J Fl: Fl•(l(R12,t)+Rf1,11)/TWOl••JI GAMW: GFMW U2MP: Ul2,l) GAAR: ZERO NX: (NCYCLE/NNl*NN OPRfNT: NX .Eof MCYCLE F(OTRACEJWR TF(16,25JMWHEAO,MWTAIL,MWTM1,RRHEAO, RRTA IL,R1_HEAD, RL TAIL, IJ2MP, Fl ,1Jt3, RI;,, NI ,GFMW, GAMW C C C••• •••DO LOOP TO CftLCULATE CELL PROPERTIES••••• C 165 ••• C C C••• 178 1A4 185 C••• C••• C DO 235 14: 1,NL U2M: U2MP U?MP: y12,M+l) li~ ~ ~EM)MI CHI PARAMETER IF I .NOT. OSP~tR GO TO 178 U03 : Ul3 Ul3: U2MP••3 CHI: Dl•DT•IU13~~3)/ 12.DO FO : Fl Fl: U2MP lF IOPLANEI GO TO 184 Fl: U2MP•IIIR(2•M+l)+Rll,M+11)/TWOJ••JJ V?M 1: VlM+DT•IFl-~O+CHII/XIMJ FIV2M oGEo ZEQO)GO TO 185 WRflEl6,49)M WR TE 6,!iOIM,NCYC1F.dNDEX,V2M,V1M,F1,FO,CHJ•ll1:'l,U03 GO TO 31?. Vl2•MI: V2M APTIFICIAL VISCOSITY CALCULATE ft OlSSIPATlON TERM IF ( ONOVIS I GO TO 204 EXISTENCE rRITERIA IFF I U2MP •GE• U?.M I GO TO 202 I V2M oGf. vlM I GO TO 20?. COMPLETELY CENTERED PAHAMETfRS IFIOTRACE)WRIT~(16,43JM•U2M,IJ2MP,CHl,UO:'l•lll~,FO, 1 FJ,y2M,VtM•XIMJ,OT AC: OSQRT(P l,Ml1TWO•IV2M+V1MII HfTAC: ((WUN/V2M1+(WUN/VlM)I/TWO C C C C C C C C 202 203 204 l . LINEAP TERM ur)IF : OAA<;(U2MP-ll2MI QL: CL•AC•HfTAC•,~IF QUAORA TIC TERM Q~: CO,..CO•HETAC•rrOIF•UOIF TOTAL ARTl~ICIAL VISCOSITY FOR 1/2 POSITION FORWARD Q2M: QL+QQ GO TO 203 8f~,~) ~E~~M 1FIOTRACEIWRIT~(16,441M•AC,HFTAC,P11,Ml,UOIF,QL,OQ END OF vlSCOSITY CALCULATIONS f.NF.RGY CONSERVATIOM OF PERFECT GAS CALCULATE NEW SPECIFIC ENERGY ANO PPESSURE QRAR: I Q2M+Q(l•Ml)/TWO ENUM: E(l,Ml-(Plt,MJ/TWO+QBARJ•IV2M-V1M) IF(OTRACEJWRIT~(16t3RIM•ENUM,E11,Ml,Pll,MJ,~12,MJ, Q(l,MJ,V12,Ml,Vll•MJ,QQ,OL,UOJF,AC •*•WAVE ENERGY 10DITION••••• DATE 011477 PI\GE N 0 1.0 •••••• AMAIN 000:571 000 000372 000 000373 000 000374 000 000375 000 ooo.376 000 o8o37J 000 0 037 000 ooo.579 000 000380 000 000.381 000 000382 000 080383 000 0 0.384 000 888~~6 888 000387 000 000388 000 000.389 000 000390 000 000391 000 000392 000 000393 000 000394 000 000.395 000 000396 000 000397 000 000398 000 000399 000 000400 000 000401 000 000402 000 000403 000 000404 000 000405 000 000406 000 000407 000 000406 000 gg~~n ODO 000 888~H goo 00 000413 000 000414 000 ogo41g 003 0 041 00 000417 000 000416 goo 000419 00 000420 000 000421 000 000422 000 080423 000 0 0424 000 ogo42i 000 0 042 000 000427 000 000428 000 000429 000 000430 000 000431 000 000432 000 000433 000 000434 000 000435 000 000436 000 000437 000 000431:1 000 000439 000 000440 000 000441 000 000442 oog 000443 00 000444 000 000445 000 00044b 000 •••••• 210 IF C .NOT. OEWftV) GO TO 225 IFC.NOT. OEf.NOJ GO TO 210 GAMW: GFMW MC: M - MNWCO<; IFCM .GE. MXWCoS)GAMW: GMW IFCM oGTo MNWCoS .ANO. M oLT. MXWCOS)GAMW: GCOS(MC) Go TO 225 OWVENO: .FALSF"• GAMW: GMW lFCM .GT. MXWCnS) GO TO 220 lFCM .GT. MWHEaO)GO TO 220 GAMW:<;FMW lF(M .LE. MWTM1) GO TO 220 IFCM tLTo MWlAyJ) OWVEND: .TRUE, ••• pH CALCIJLAr ONS FROM 154 ••••• WF1CLCM):Wf2Cl(M) UR: CRLHEAO-M>•R~F PHIW: DR• WAA / STWOWV IFCPHIW oGT. PHI) PHI: PHIW WF2CLCMI: CENWAV1WSLSOR)•ICOEXPIPHIW)-WUN)+ 1 (OEXP(-PHIW•WSREXp)-WUNI/WSREXPI IF(WE2CL(MI .Gr. ENWAVI WE2CLIM): ENWAV IF (OWVENOI wE,CL(M) : ENWAV 219 WOLENG: Wf2CL(M) - WE1CLCMI DG = CG-GF>•CWF.:2C1 CMJ/ENWAVI GAMW: GMW - Ur, IFCM.GT,MNWCOS, GO TO 219 WADENG : WDLENG GO TO 223 WSPAN : M - MXWCOc, ' wsPAN = WSPANIWSRc,PA•Pl WSF : (DCOSC THRE~•WSPANt-ENINE•DCOSCWSPANl+ 1 FIGHTJ/SIXTEE WAOENG: WOLENG * WSF MC=M-MNWCOS 224 t C••••• ~··· 225 GcOS(MCI = GMW - nG•WSF GAMW: GCOS(MCJ GO TO 223 IIIA8ENG: ZERO EN M: ENUM + WAD~NG lF(M .E~. MwH~aDIRHEAO:IRC2,M))+OR•CR(2,M+lJ-RC2,M)) IF (M oNE.MWTAlt I GO TO 224 ooR=CRLTAlL-Ml•RE~ RTA1L:CR(2,M)) +UnR•CRC2oM+1)-Rt2,M)J IFCOTRACEIWRJT~(l6t37)M•RC2,M),OR,PH1W,ENUM,WE?CLCM)• WE1CLCMl,WAOENGtPCl,M),O~AR,VC2tM)•WSPAN ••••• HOMOGENE011S ENERGY SOURCE••••• IFCeNOT, OESORJGO TO 230 IFC,NOT. OEEIGn TO 226 GAMW:GFMW Mc:M-MJNCOS IF(M .GE, MAtCnS)GAMw:GMW IFCM .GT. MJNCnS .AND• M oLT, MAXCOS)GAMW:Gr.OS(MC) GO TO 230 22b G11MW:GMW IF C ,NOT. OADnEN) GO TO 230 OADUEN = M .LT, Maxcos IF ( ,NOT, OADnEN) GO TO 230 or,:(G-GF)t(E2CL/lNMAXI GAMW:GMW-DG IF CM ,GT• MlNCOS I GO TO 227 AOOlNG: DELENG GO TO 2~6 227 SPAN: M - MAXcos SPIIN: SPAN/ SORSpA • Pl SF=locos ( THRFE. SPAN - ENINE. nco~ (SPAN) • + EIGHTI / SlyTEE AODENG=llELENG•SF Mc:M-MINCOS GCOSCMCl:G~W-IDG•~F) GAMW:GCOSIMC) ENUM: ENUM+ AnOENG 228 C••••• 230 EOEM: WUN+ GAMW•CWUN-VlM/V2M)/TWO DATE 011477 N t-' 0 •••••• 000447 000448 000449 888~~~ 000452 00045~ 03045 0 045 00045b 000457 888~~8 000460 000461 000462 000463 0004b4 000465 0004&6 000'+67 0001+68 000469 000470 000471 000472 ggg~H 000475 000476 000477 000478 000479 000480 ggg~g~ 000483 00048'+ 000485 00048b 000'487 000488 000'489 000490 000491 000'492 000'493 UOn'494 000'495 000'496 000497 000'496 000499 000500 000501 000502 000503 000504 000505 000506 000507 000508 030509 0 0510 000511 000512 000513 000514 000515 000516 000517 888~1~ 000520 000521 000522 AMAIN 000 000 000 888 000 000 888 000 000 888 000 000 000 000 000 000 000 000 000 000 000 000 000 888 000 000 000 000 000 000 800 00 000 000 000 000 000 000 000 000 000 OCtO 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 888 000 000 ooo. •••••• C • 236 237 C••••• C 238 239 C1238 C 240 246 C1242 247 C C C C C 248 ENEW: ENUM/EOEM E12•MI: ENF.W P(2•MI = GA~w•rNEw/V2M A(Ml:OSQRTIIWUN+GaMWl•Pl2,Ml•Vl2•MII IF(OTRACE)WRI1F(l6,41)M,NL,ACM),GAMW,Pl?•M),V2M,VtM,ENEW, EOEM•AD0FNG,OELENG,SPAN•SORSPA,RTAILtRC2,Ml,RH[A0•00R ••• END OF LOOP FOR CALCULATION C£LL PROPFRTIES••••••• ••• RACK TO 1&5 ••••• WRITEC18,236JT,NCvCLE,MWHEAO,MWTAIL,IPl?•Il•Vl2,IJ,I:1,5J, IPl2,111,Vt2•Tll,fI=10•60,101 FORMAT(' '•El0.3•~15, OF9.5/12F9o51 OP[AK: NSAM oNE• 0 OfNT: .TRUE, IF I OSKIP I r.0 TO 237 N~X: INCYCLE/NSl•NS OPUN: NSX ,EQ. NrYCLE IFl,NOT, OPRINTI GO TO 238 CALt SAMPLE CAL INT CALL FIOIF ••• FIDIF - 433 •••••••• f~lif~~v~~~TttTt~he,RRHEAO•RRTAIL !Fl.NOT, OPUNJ60 TO 240 lFIOPRlNTJGO Tn 239 ~:ct 1~~PLE CALL PUOAT 1ST: 1ST + N PUOAT: 7q1 DETERMINE NEXT TSTEP,TIME AND RFINITIAL PROPERTlfS lF I OVER I Gn TO 248 lFIOPRINTJGO Tn 241 IF(OPllN)GO TO ,41 CALL SAMPLE CALL INT CALL PUDAT PUOAT: 7ql IST=IST+N NX=INPR/NFRPR)•NFpPR IFINPR ,EQ• NX1GO TO 246 GO JO 247 CALL F OIF FIOIF: 411J NPR=NPR+l WR1TEC6,31JT,ET,~ OVER: ,THUE, DTL: OT OT: OTHOL TLAS: T OPRINT: ,TRUE, IFIOTRACEJWRITFl16,45)NCYCL£,IN0£X,T,OT,DTHOLtDTL GO TO 289 TLAS: T OTL: OT STABILITY CRITERIA R?MP : R 12, 1) Gf: GF OFLIT: ,FALSE, IFIOTAACE)WAITFl16,391 Do 281 M: 1,NL A2M: R2r,IP R2MP: RC2,M+ll ROIF: R2MP-R2M V?: V(2,MI Vl : VU,M) OATE 011477 PAGf. 6 N t-' t-' •••••• AMAIN oon52.3 000 000524 000 OOn525 000 000526 000 OOn 5 27 0 00 oon528 000 oons29 000 0005.30 0 0 0 0005.31 00 0 0005.32 0 00 0005.33 000 OOn5.54 00 0 oon5.3; 000 88Rs~ 888 0005.38 000 0005.39 000 000540 000 000541 000 000542 000 00054.3 000 000544 000 00054i 000 00054 000 000547 000 000548 000 000549 000 000550 000 000551 000 000552 000 000553 000 000554 000 ooo55i 000 00055 0 0 0 000557 0 00 000558 0 0 0 000559 000 000560 000 000561 1)00 0005&2 000 00056.3 000 000564 000 oon565 000 000566 000 0011567 000 000566 000 0005~9 000 oon57 f goo 0005 00 000572 000 00057.3 000 0005~4 000 ooo57i 000 0005 0 0 0 000577 000 000578 000 000579 000 00056~ 000 00058 000 000582 000 000583 000 000584 000 000565 000 ggo586 888 0587 000568 000 000569 000 000590 000 000591 000 000592 000 00059.5 000 00059'+ 0 0 0 ooo5ii 000 0005 000 ooosii 000 0005 000 •••••• 265 267 AS2 : P{2,Ml•V?•G~ VOOTN: TWO•{V~-V1l/(V2+Vll/DTL lF ( VOOTN ,LT. ZERO l GO TO 265 AS2 : ZERO GO TO 267 8~1: CO•ROIF•VDOTN HS2: 64,DO•HSl•ASl S2: AS2 + B52 1F!S2 ,GT. 1ERn)GO TO 268 S:,: WUN WRITE(6,531DEM,S2,AS2,RS1,AS2,vooTN,Vl•V2,RDtF,R2MP,R?M, • CO,OTL C C 268 270 • 276 1 2AO 2A9 l C 290 C 297 C••• C ~ 309 l OF:XlT: ,TRUJ• DEM: 3,00•D QRT(~~) OT: TWO•RD1 IOEM IF( OFLIT I Gn TO 280 OFLIT: ,TRUE, DTMIN: OT IF (OEWAVI GO TO 270 IF ( ,NOT, OESnR I GO TO 2AO GRADEN: I ENMAX I SLOSOR •AA/ TMAXE) • ( DEXP ( PHI I - DEXP ( - PHI• SOREXP) ) IFIGRAOfN •NE• ZERO) GO TO 276 GO TO 280 GRADTA: WUN/ GRADEN UTTES: ENSTEP • ~RADTA 1F(OTRACEIWRIT~!l6,401M•DT,DTMIN,DTT~S,GRAOTA•OEM,R2MP, R2~,HS2•n5l,VD0TN•GE IF ( oTTES .GT, DTMJN) GO TO ?80 OTMIN: nTTES IF I OT ,LT. OrMIN) DTMIN: OT IF(OTRACE)WqJT~(16•40)M•OT,OTMIN,OTTFS,GRAOTA•DEM,R2MP, I R2M,HS2•RS ,VOOTN•GE F ( M oF.O• NLT GE: G OT::: DTMIN IF ( OT oLE, Z~qO J 60 TO 309 LIMITING CONSTRAINT PIF (OTL oLE• 7ERO) GO TO 289 DTU = 1 • 4 ao.oT, IF ( OT oGT, DTUP) OT: DTUP T: T + OT 1F(OTRACEIWRIT~(l~•40 1M• DT,oTUP,OTL,DTMIN,OTTES• GRAOEN,D~M,PHI,VOOTN,GE•RDIF REINITlAL PROP~RTIES 00 297 M: l,NL U(l•MI : U(2,M) R(l•M>: R(2tM) V(l•M>: V(2,M) Q(l•M>: Q(2,M) P(l•M>: P12,M) E , KLBl?.)1 KPA(2), KR0(2), KLOC2) EQUIVALENCE (OJSTrII, Mil)), COIL l 4061 407 411 1'13 1115 lllb 420 lt22 !RUFF XOR ( IUUFF, LYR l lrEAR OR(IY~AR•rRUFF) lC OM M AND(LCOMM,INSTR) IoAY = OR(lCOMM,lnAY) 17.EHO: ANO(LZFRO,ISHOV) lZERo = IZfRO•ICHaR l5THIP = AND (JlEpO,ISHOV) lZEHo: XOR(JZERO,JSTRJPI IoAY: IOAY•JCHAR 1STRIP = ANO(lOAY,ISHOVI ITEST: XOR(ISTRlp,JOAY) Oc;I=ITEST oEG• IlJ:-RO ILEFFT: ANOILBLAN,ISHOV) ILE T: XORIILEFl,LALANI !STRIP= YR(Il[FT,~~TRIP) IOA~F~O?JAY/Jg~x~f RIP IoAY: ANO(IOAY,I~HOV) IoAY: ORIIOAY,IL~FT) InAY: 1DAY•ICHAH IM: 1 DO 361 I: l, 12 IF I IMON ,E(h UAII) ) IM: I CONTINUE K: 4o ( IM-I I DO 364 I: 1,3 NLLIIJ: MONTH(K+y) NLL(lt) : HlAY NLL(SJ : lYEAR•ICHAR FoRMATC' '•'••••••••• GAMMA CHANGES Q ENERGY AODITION ••••'I FORMAT( ' '• 4A4•a5, ' TOTOEOTOEOEOTOOTOEOEOT DF.TONAT•, 'ON wAy~s oTO~oToEoTgo ,, 3A4, 'EOE ',A6,• •,A6, I FORMAT T 9 ,4A4,A5,t TOOE i,2A4,14X,315,9X,2A1t, 3A4,'EOT '•A6,• '•A61 FORMAT(' '•4A4,A5,• TOTOE'•2A4,4X,2Fl5.A,4X•2A4, 3A4,'FOT '•A6,• '•A6) FORMAT C •1•,/1 I FORMAT Ct •, 4A4•ft5J 'ETOEO• , 2A4, 3BX, ?.A4,:'IA4, •EOE ' ,A6,' '•16 FORMAT(' '•4Alt,A5,• T000E'•2A4,4X,6IS,4X,2A4,3A4• 'EOT •,A~,, ',A6) FORMAT C ' •, 4A4,A5, ' T00EO•, 2A4, 6X, 7A4, 4Xt2A4, 3A4, 'EOT •, A~,, '•A61 FORMAT(/) OFOH: .NOT, OFOR IF ( OFOR) GO TO 413 DO 1111 I: 1,2 KLOl 1 I> : KUH I> KRO II : KRIH I) r.o TO 416 00 4lo:; KLOl II : KLFCl) KROlJ): KRFCl) ILIM : 15 DO 4:,0 WRITE C 6,402) WRITEC6t398) I : 1,2 K: 1•3 DO 420 T: 1,8 WRITE C 6, 3991 NLL, KRUN•Mf1),MC2) 00422 1:1,9 WRITE C 6, 40J) NLL, KLO• KRO, KRUN•MC t) t"1(2) WHITE( 6, 405 I N1 L• KLO, LABEL• KRO, KRUN,MC1),MC2) WRITE(6,40'+)NLL,KLOtLSTART•NCYCLE,NPUNCH,NSlORF,NS,NSTEP5,KRO, 1 KRUN•Mfl1,M(2) WRITE(6,404)NLL,K10,NFJNALtNN,NNN,NBUFF,NFRE~,NWAVE, l KRO,KRUN,M(1),M(2) WRITE(6,'+0l )NLL,K1 O, TERMIN•TIPUN,KRO,KRIIN,Mf I) ,MC21 WRITEC6,4011NLL,KLO•G•GF,KR0,KRUN,~(11,M(2) WRITE (6,400) NLL, Kt O, .J• NOP, NL I, KRO, KRUN, MC 11 , MC 2) DO 426 I: 1,9 WRITE C 6,4031 Ni.1!..•KLO, KRO,KRUN,M(l h!W!f2) no 42A I: 1,ILIM WRITE I ~,399)NLL 0 KRUN,~C\J,MC2> IFCK eEQ, 2JGO 10 4.30 IF C OFOR I 11 M: 2 CONTINUE 01\TE 011471 Pr.Gt •••••• BURST •••••• 000126 000 WR1TE(6,406) 000129 000 R£TUHN 000130 000 p431 IN IT 1 L : ,, 75 000131 000 2431 MAIN PROGQAM = 319 000132 000 ENO lilHDG•P •••••• FIOIF •••••• WELT•L FIOIF ELT 6:-01/14-14:42 FIDIF oonoo ooo 00000 000 C 433 SUBROUTINE FIDrF FI01F 000003 000 000004 000 IMPLICIT DOUBLE PQECISJON (A-H,P-l),LOGJCAL (01 DOURLE PRECISION ~ACH ooooos ooo 11 COMMON/ ARRAYS I Ul2,20l>, R(2•201),V(2,?01 1,Q(2t20I>• P12,2011• Xf20ll, E(2,20J), NCELLf20tl, WE2CL(201J,WElCLl?Ol)tA(20l),GCOS(lOl COMMON/ TIMF. / TERMIN • TJPUN, T• OT, OTL, 000006 000 U 000007 000 2 000008 000 13 000009 000 oonolo ooo 000011 000 000012 ooo 000013 000 OOOOh 000 88881~ 888 000017 000 000016 000 88~8~3 888 000021 000 000022 000 ogoo2J ~og 0 0024 0 000025 00 00002t, 000 000027 000 000028 000 000029 gog OOOOJO O 000031 000 000032 000 000033 goo 000034 00 000035 000 000036 000 000037 000 OOoOJ8 000 oonoJ9 ooo 000040 000 000041 000 000042 000 00004J 000 000041f 000 00001f5 000 000046 000 000047 000 000048 000 00001+9 000 000050 000 000051 000 000052 goo 000053 (IQ 000054 000 000055 000 000056 000 000057 000 000058 000 000059 goo 0000&0 00 000061 000 0000&2 000 0000&3 000 000064 000 • KRUN(J), LAAELC7) COMMON / PARAM / c, • co, G, GF. UL• UR , GMW• GFMW, f"NMI\X, rr, ~ PJt,SLOSOR, SOREXP, TMAXF• Two, WUN, ZERO' 4 MINCOS, MAXCOS• J, JPl, N• NIL NLl,NOP, ~ NPART, NSTEPS, OENT, OESOR, OPEAK, OPLANE, 6 OPRJNT• OPUTI, OSKJP, oSPHFR,oTRACF.• 7 MWTAIL,MwHEAD,WVEL1WlnWAV1ENWAV,Rl~AD,RTAIL, A WSLS0R,W~REXP,MNWCo~,MXWCo5,0fWAV,RFF,f2CL 22 COMMON/ ARGJNT / rNDEX, LSTART, NCYCLE• NFINAL, NSTORE• NS, 7 NN, NNN, NPEI\K• NSAM• NSHIF,N~UFF,NFRF.O,NWAVE 2 ,NPUNCH 451 FORMAT t TIME=', 1POl2o6• 3X, 'DT =·· D1?o6t 5X, 'INOEX =•, • lo;, ax, •vTIMD ='• 012.6, 7X, •CYCLE=·· 15 / 453 FORMAT( • CELL rENTER DIST CENTER PRESS PRESS OJFF •, • •ciLL SP VnL cINTER VFL CF.LL Vtsc ,, • •c LL ENER~Y CON INUITY CELL•) 454 FORMAT(' ,•TRE E~ERGY W VE BEGINS AT•,Fl,.9, 1 ' ANU ENnS AT •,F14o9) 45& FoRMATI' •,1s,2F1~.9,1PE14•4,0PF14o9,1P2E14o4,0PF1-.9,F9.4,t~) 457 FORMAT •o•, 20X• 'THE LEADING MACH NUMAER :, , F7.4 ) 458 FORMAT C lHl) 45':I FORMAT C •+•, 6Xt A( •••, l:'IX >, •• Pf:AI( ••' 460 FORMAT('+', 6X, 8rls•,1Jx),T~ ~EAD s~T) 461 FORMAT(•+•, &X, 8r•!•,t3X>••! TAIL r!•> 462 FORMAT(lHO,•NO LE~O SHOCK WAVE YET•) VTIME: T - oT1TWn C461 PRINT HfAU~RS WRllE ( 6, 4581 WRI E l 6, 451 > T,OT,INOEX,VTIME, NCYCLE WRITE C 6,453 ) C ••••• CALCULATE CURRENT PROPERTIES MM: NFRF.Q IF(NFREQ,LE,O) MM: NL/25 IF ( ~M,EQoO) ~M: l lFCoNOT. OEWAVtGO TO 470 MWHP9: MWHEAD + Q MwTM5: MWTAIL - o; -70 R2: R(2,l) U.?. : U(2, 1) RPJ2: R2••JP1 IJ: 0 DO Rl: R2 Ul - U2 R2 RC2,M+1) U2 - UC2,M+1) PO P(2,M)-WUN OMID: (Rl+R2J/TWn UMUD: CU1+U2)/TWn RPJl: RPJ2 M: 1,NL RPJ2: R2••JPl IFCM oLTo 9)GO TO 485 IF(M ,EOo NPEAw)GO TO 4A5 IF ( M oLTo MWTM5) GO TO 484 IF CM oLEo MWHP9) GO TO 485 F C MoEQ, NL l GO TO 485 DA1'E 01141'1 l'AGE •••••• FIDIF' •••••• 000065 000 lt84 MX: I I M-NSHIF 1/MM I •MM+ NSHJF 000066 000 IF I MX .NE.MI GO TO 492 000067 000 1485 CM: CRPJ2-RPJJl/v(2,MI/PJl/XCMl ogoo68 000 C 48b PRINT CURRENT PROPERTfJS 0 00b9 000 WRITE I 6, 456 I NCELLCM), DMID• P ,M1,Pn,vc,,M1, 000070 000 • UMUO, QC2,MI, EC2,MI, CM,M 000011 000 IJ = IJ + 1 000072 000 IF IM oEQ• MWTAIL I WRITE (6,461) oooo7J 000 IF CM oEQ. MWH~AO I WRIT~Cb,4601 000074 000 IF ( M oN~.NP~AK I GO O 491 000075 000 WRITE C 6, 45 I 000076 000 GA: y 000077 000 lF M .LT. MWHEftO)GA: GF 000076 000 MACH: DSQRT I r PC2,MI - WUN l•CGA+WUNI I TWO/GA+WUN) 8888J~ 888 :~l C fFCI( .GT. 601r.O TO 49J ON INU 000081 000 q9J IFCP 2,NPEAK) 0 GT. WUN) WRITEC&,457) MACH 000082 000 IF\PC2,NPEAK\ ·kE· wvNI WRIT~16,462) 00008J 000 F OEWAVIWHI E1 ,454 RHEAO,R IL 000084 000 RFTURN 000085 000 Cl494 INITIL : 730 000086 000 C2494 MAIN PROGi:,AM 238 000087 000 C34lLOGfCAL co, OOOOOJ 000 000004 000 11 COMMON/ ARRAYS/ U 2,2011, R(2•20 J,V 2,201 l,QC2,?0JI, U PC2,20ll• xc2011, EC2,20Jl, NCELLC20ll, 000005 000 000006 000 000007 000 000008 000 000009 000 000010 000 000011 000 000012 000 000013 000 000014 000 8888½~ ggg 000017 000 000018 000 000019 000 000020 000 ogoo21 ooo 0 0022 000 000023 000 00002i. 000 000025 000 000026 000 000027 000 000028 000 000029 000 000030 000 000031 000 000032 000 000033 000 ooooJi. ooo ooooJs ooo 000036 000 000037 000 000038 000 000039 000 000040 000 C C 000041 000 8888:~ 888 000044 000 000045 000 2 WE2CLC20l),WF:1CLC2011,AC201),GC05C101 COMMON / PARAM / C1, C01 G, GF, UL, UR f GMW• GFMW, F:NMAX, E'T, 3 PJt,~LOSOR, SOHEXP, MAXF• Two, WUN, ZERO, 4 MINCOS, MAXCOS• J, JPl, ~, NL, NLI,NOP, 5 MPART, NSTEPS1. OENT, OESOR, OPE•K, OPLANf:, 6 OPRINT• OPUTI, OSKIP, OSPHFR•OJRACE• 7 MWTAIL,MwHEAD,WVEL,WinwAV,ENWAV,RHEAO,RTA L, 7 8 WSLSOR,W~REXP,MNWCOS,MXWCOS,OEWAV,REF,E2CL 50 FORMAT( l - 508 FORMAT('O'•lsx,•K~RNEL PRESSURES P(KERol/PCAMA 0 1:•,FS~2•/ l X,• TCKER.J/T(AMAol: 1 ,F~.2,/ 2 lOX,'THEQF ARE•,IS,• CELLS WITH•,I5,• FAIRJNG CELLS',/) 511 FORMAT C '4', ~X• 1EN~RGY SOURCE!' ENERGYCSoR. MAX. / ENERGY•, • •CINT. AMA.t: •, F8.l• t F NAL F.NF.RGY DEPOSITION TIME'• • •(TAU MAXol: •, F5.2, / , 2UX, • SHAPJNr, CONSTANT 1: '• 3 F6•2• • SHAPfNG rONSTANT 2: •,F5.2, • crLL NO OF RO:,, 4 F5.2, /, 2 X• t(MIN.I CELL NO, OF COS. DtSTo: '• 14, • t (MAX.I CELI . NO. OF COS. OISTo : •, 14 I RF.AD CS,5071 PRE~S, TEMP, N, NOEC WRITE ( 6, 508) PRESS, TEMP, N, NDEC IF(N .Gr.201>Gn To 609 Gr;MW: GFMW IF(OEWAVJGGMW - GMW F(OESORIGGMW=r.MW NL: N-1 LIMIT:N-2•NDEC HUFF: N-NDEC IF C OESOR I Wi:,ITEC 6, 511) ENMAX, TMAXE, SLOSOR, • SOREXP, ~UFF, MINCOS, MAXCOS REF: WUN/RUFF VOL: TEMP/PRESS Rl: ZERO RlJ: ZERO RCl•lJ : Rl ~H:JI ~ ~sps NCELLCU : 1 ACMl:oSORTCCWUN+G~MWl*PCl,Ml•VCl•M)l flATE 01h77 PAf',F l ) •••••• 000046 000047 1!)00048 000049 000050 000051 OOn052 0Uo05J 000054 000055 000056 OOnU57 000058 000059 00110~0 OOnObl 000062 00006J 0Qf)064 000065 000066 00()067 000066 000069 001)070 000011 000072 oooo7J 000074 000075 000076 ogoo77 0 0078 000079 000080 000081 000082 000063 000084 000085 88~8~~ 001)086 OOOCJ69 000090 000091 000092 001)093 000094 , oono9s 001)09b 000097 oon096· 000099 000100 000101 000102 0100103 000104 ogo1os 0 0106 000107 888188 000110 000111 000112 ooollJ 0001111 000115 ssgn~ 000118 000119 000120 000121 GENOAT 000 000 000 000 000 000 000 888 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 800 IJO 000 000 000 000 000 000 000 888 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 800 00 000 000 000 000 000 000 888 000 000 000 000 •••••• 568 587 Ell•ll : TEMPIGGM~ IF ( LIMIT .tQ. 1 I GO TO 545 no 544 M: 2•LIMIT R::tJ: RlJ HUFF: M-1 fH : BUFF • REF RlJ : Rl .. JP1 R!l•M) : Rl PI 1 • ~O : PHF:SS Vll•MI: VOL X(M-11: (RlJ-R2Jl/YOL/PJl NcELL(MI : M lF(M .GTf NLl,r.GMW: GMW A(Ml:QSQRf( WUN+Gr.MWl•Pll,Ml•V(ltMI, E(l•M): TEMP/ Gr.MW LL :: LIMIT XRMXO:: NUEC XRMXo: TWO•XRMXO•REF PS:: PRf.SS - WUN TS:: TEMP-WUN DO 56~ I: l, NDEC M: LL+ 1 BUFF: M-1 R;>J: RlJ Rl: BUFF•REF RtJ : Rh•JPl Rll•M•) : Rl XCM-11: IRlJ-R2Jl/YOL/PJl NCELL(M): M X"'1X0 : I XMXO: XMXO•REF sc~LE: TWO•([XMXn/XRMXOl••2J PV: PRESS-SCALE*PS TV: TEMP-SCALF:•T~ VOL: TV/PV VI I •M) : VOL P( I •Ml : PV lFIM .GT\ NLI)r.GMW: GMW A(Ml:OSQRT( WUN+Gr.MWl•Pll,Ml•V(ltMIJ E(l•MI : TV/GGMW LL: LL+~SEC 5A7 I: l,NDEC M: LL+I HUFF: M-1 R2J: RlJ Rl: AUFF•REF RlJ : Rh•JPl X(M-11 : (RlJ-R2Jl/YOL/PJl R(l•M) : Rl NCELL(MI : M XRMX: NOEC -I XRMX: XRMX•REF ScALE: TWO•l(XRMv/XRMXOl••2) PV: WUN+SCALE•PS TV: WVN+SCALE•TS VOL: TV/PV PC l •M) : PV V(l•MI: VOL IF(M oGT. NLllGGMW: GMW A(M1l:OSQRl ( (WUN+Gr.MWJ •Pll ,MhVU ,M)) E(l•MI: . V/GGMW LL: N+l F no 599 M: LL, 201 BIIF : M-1 Rl: BUFF• REF NCELL(MI : M Pl l •Ml : WUN VCl•M>: WVN lF(M eGT, Nt.llr.GMW: GMW E(l•MI: WUN/GMW R::tJ: RlJ RlJ: R ••JPl R(l•M): Rl A(M):OSQRT((WUN+Gr.MWJ•P(1,MJ•V(1•MJ) XCM-11: (RlJ-R2Jt/PJ1 no 60A M = 1,201 DATE Olh77 P11G1:: t r •••••• GENOAT •••••• 000122 0011123 000124 000125 000126 0001 27 000128 000129 001)130 000131 000132 000133 000134 000 ODO 000 000 000 000 000 000 000 000 000 000 000 608 609 C1609 U(l,p.l) : U(2•P.0 : QC 1 •M) : 0(2,M) - WF2CL(M) R(2•Ml : Pl2•Ml : l/t2•M> : El2•M) - CONTINUE RF.:TURN END ZERO ZF.HO ZERO ZERO : ZERO R ! l ,M) P l,M) VlltMI El1'MI JNITJL : gHOG•P •••••• INITIL ****** lilELT•L INITIL ELT 68-01/14-14:42 INITIL 000001 000 SUBROUTI~E INITIL INITTL 000002 000 C 611 000003 000 001'1004 000 000005 000 001)006 000 000007 000 000008 000 000009 000 000010 000 000011 000 000012 000 000013 000 000014 000 001)015 000 88~81~ 888 000018 000 000019 000 000020 000 000021 000 030022 000 0 0023 000 000024 000 000025 000 000026, 000· 000027 000 000028 000 00!)029 000 000030 000 8888i1 888 0000·33 000 000034 000 00(),0'35 000 OQ.0036 000 0000.57 000 000036 000 000039 000 000040 000 000041 000 000042 000 000043 000 8888:~ 888 000046 000 000047 000 000046 000 000049 000 000050 000 000051 000 000052 000 00()053 000 000054 000 00·0055 000 000056 000 IMPLICIT DOU~LE P~ECISION IA-H,P-Zl,LOGICAL 101 DOUULE PRECISION uACH 11 COMMON/ AARAYS I Ul2•2011, Rl2•201),V12,?0! 1,Gl2,20JI, U Pl2,2011• Xl20ll, EC2,2011, NCELLl20JI, 132 WE2CLl20ll,wF.1CLl20ll,A(~Ot1,r,CoS1tOI COMMON/ TIME/ TE~MIN, TIPUN, T• OT, DTL, e KRUN(31, LABfL(71 COMMON/ PARAM / C1, CO, G, GF, UL, UR , GMW• GFMW, FNMAX, ET, 3 PJ1,SLOS0R, SOHF.XP, TMAXF• Two, WUN, ZERO, 4 MINCOS, MAXCOS• J, JP!, N• NL~ NLI,NDP, 5 NPART, NST[PSf OF.NT, OESOR, OPEAK, OPLANF, 6 OPRINT• OPlJ I, OSKIP, oSPHFR•OJRACE• 7 MWTAIL,MwHEAD,WVEL,WIDWAV,ENWAV,RHEAD;RTA L, 22~C0MM0N / An~~~YRiW~~~~~:M~~~~~t~X~f9~t~~w~~f~~E:E~~T0RF., NS, 7 NN, NNN, NPEAK• MSAM, NSHlF,NAUFF,NFRE(hNWAVE 2 ,NPUNCH 627 DtMEANSION KINIT(3),KREST(2),KNtJMf9) 628 UAT KINIT I 'INlT'•'IAL '•'RUN /, • KHEST / •REST••'ART ' /, • KNUM / 'l TE ,'2 TE•,•3 TE•, 1 4 lE•,•5 TE•, 6 ,t• '6 TE'•'7 TE','8 TF.•,•9 TE• / JJ DATA ALLOW / l,no-5 / 634 FORMAT() 635 FORMAT!) 63b FORMAT I 637 FORMAT() 638 FORMAT(/, 40X, 'FIRST TIME STEP :',F13,6/ 2 .. gx, •GAMMAl :• ,F13,6/ 3 4 X, •GAuMAIJ :• ,F 3.6) 639 FORMAT(30X,'PLANAn GEOMF.TRY•) 640 FoRMATC30X,•CYLINnRICAL GEOMETRY') 641 FoRMAT(30X,•SPHERrCAL GEOMETRY•) 642 FORMAT(3ox,•oES1GNATED MAX TIMF.: •,14,, SF.CoNns•, 1 30Xt 1 NUMBEn OF TIME STEPS :•,15/1 6 .. 3 FORMATJlOX,•RESULTS WILL BE SlOREO FOR RESTART•) 644 FORMAT l5X, 1 LlNEAR ARTIFICIAL VISCOSITY COFFFICIENT :•,013.6/ l 15X, 1 QUA0RftTIC ARTIFICIAL VISCOSITY COEFFICIENT :•,nt3,6/) 645 FORMAT(lOX,•RESllLTS WILL NOT HE SlORF.O FOR RFSTARTING•) 646 FoRMAT(lOX,•RESULTS OF EVERY•,15,• CYCLFS ARF ~TOREO ON TAPF') 647 FORMAT(70x,• Er =',012,61 648 FORMAT I '0'• ~ox, ' THE SHOCK FRONT MACH NUMBFR ='•F7,4) 649 FoRMAT(lOX,•THE MaYIMUM NUMH[R OF CELLS 15 ',15,/ 1 lOX•'RES\JLTS AHF. PRINTF.:n EVERY•,15,, CYCLES'•/ l lox, 'THE C11RRENT NU"1BER OF flATA POH,rTs ts•, 15,/ f ig:;:vA~M~,S~A~~fictfyC~~LT~~Mrf~;·iau~nAPY 1S•,Fl5.tO,/ t fox, 'THE F,-ow VELOCITY AT THE RIGHT BOUNDARY IS' ,Fl~. 10,/1 6~0 FoRMAT(lOX,•THE E~ERGY FUNCTION SLOPF CONSTANT EOUALS•,Fto ... ,, I lox•' THE Ei>.1ER6Y SLOPING CONSTANT EQUAi s•, FlO ,4 • / 2 lOX,•THE MaXfMUM TIME OF ENfR6Y ADDITfON IS',Fl0,4,/ 3 lOX,•THE Max MUM ENERGY ADDFD IS•,Ft0.4,/ ~ lOX,•THE SpAlIAL ROUNDING FUNCTION AEGINS AT CFLL',l!'i,/ !'i lOX,•THF.: 011TERMOST [OGE OF THF: EMF:R«;Y FUNCTION IS AT' DATE 01 h77 PI\GF: 2 N ,.... 00 •••••• INITIL •••••• 000 000 000 000 000 0 00 000 000 000 000 , 000 000 000 000 ono 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 000 6 651 652 2 3 .. !i 6 6537 651t 655 l C C651t C C C 'CELL'•T51 FORMAT() FORMAT(lOX,•THF WAVE VELOCITY 1S',1P[10 0 4,/ lox,•THE WAVE FRONT WIDTH IS•,F.10.4,/ lOX, 'THE EMERGY ADDED IS• •El0,4,/ lOX,•THE WAVE SLOPE CONSTANT EOUALS'•OPFl0,4,/ lOX,•THE ENERGY SLOPE CONSTANT fOUALS•,FI0,4,t . lox, 'THF. SpATJAL ROUNDING FIJNCTION pEGIN'i AT CELL•. Jc;,, lOX,•THF. LftST ENERGY CELL JS •,15/) FORMAT(lOX,•NO EN~RGY ADDITION•,/ FoRMAT(•O•,JoX,•lHERE IS NO SHOCK WAVE ',Eln 4) FoRMAT(lOX,•RF.SULTS WJLL BE STORED AT TIME INfFRVALS•, ' OF •, F 1 0 • 6 OPASS: LSTART .N~. 0 IF(OPASSI GO Tn 665 NPAHT: 2000 OFNA[): .FALSE. T: ZERO OTL: ZERO READ TNPUT RF.AD(5,631tJNSTEPS 0 NFINAL,NN,NNN,TERMIN,TIPUN,N~UFF,NFREQ,NWAVF RFAD ( 5, 635)NDP,J,NL1,CL•C0,G,GF,ut.,UR - IF(OTRACE)WRIT~(6,634JNSTEPS,NFINAL,NN,NNN,TERMIN, TIPUN,N~~F,NFREQ•NWAVE 1F(OTRACF.>WRIT~(6,6J5JNDP,J,NLl,CL,CO,G,GF,tL,UR 0F50R: NBUFF ,NE. 0 OFWAV: NWAVE ,NE. 0 000057 000058 000059 000060 000061 000062 000063 00006«. 000065 000066 000067 000068 000069 000070 oono71 000012 000073 000074 000075 000070 000077 000078 000079 000080 000081 000082 000083 000084 000085 000086 000087 030088 0 0089 000090 000091 000092 000093 000094 000095 000096 030097 0 0098 000099 000100 000101 000102 800103 00010«. 0001005 0001 6 000107 000108 000109 000110 000111 000112 000113 0100114 000115 0001 6 000117 000118 000119 000120 000121 000122 000123 000124 000125 000126 000127 000128 000129 000130 0001331 0001 2 800 00 OQO 000 000 000 000 000 000 IF ( ,NOT OEWaV) GO TO 660 REAO (5,651) wVEL•WIDWAV,ENWAV, t WSL~OR,WSREXP,MNWCOS,MXWCOS TISTO:WUN/(WVEL•5n0,0DOI 800 OU 000 000 000 000 000 000 800 00 000 000 800 00 000 000 000 000 008 00 0 00 000 000 000 00 0 000 000 000 000 000 000 000 000 000 000 000 660 • 661 665 670 C••••• 67!> C1667 C C lF(TIPUN .E~. 7EROIT1PUN:T1STO OF'N n: ,TRUE, - IF( ,NOTl OE~OR) GO TO 661 READ 15,6~7) ~LOSOR, SOREXP,TMAXE,ENMAX• MINrOS,MAXCOS OFNAO: ,TRlJE, GMW: G - WUN GFMW: GF - WUN NL=NOP-1 JPl: J+l PJl: JPl DEFb1E INITIAL MESH POINTS ANO CEtl PARAMETERS OALT: NDP,EO• 0 IFCINOEX .NF., nJ GO TO 675 KRUN(l) : KINIT(l) KRUN(2): KINITl2l KRUN(3): KIN1T(3) GO JO 670 F(INDEX .N~. O) GO TO 728 KRUN(l): KREST(ll KRUN(2J: KRESTl2l KRUN(31: KNUM(LSTART) CALL BURST ••••••BURST: 3~1 •••••••••• IF(OPASSJ Gn TO 715 WRITE 16,64q) NF JNA1 ,NN,NDP,NLI ,uL,UR IF( ,NOTf OENAnJ WRITE(6,6531 IF ( OAL) CA1L GENDAT GENDAT: o;OO IF ( OALT I GO Tn 695 Rl = ZERO GAMW: GFMW M: 1,201 00 694 NCELL(MJ: M lF(M.GTeNJ GO TO 679 R2: Rl A~A0(5,636)K,Rl,Url•M),P(1•M),V(l,M),9(l•MI R(ltMJ : Rt RolF: R1-R2 DATE Olh7T PAGE •••••• JNITIL •••••• 000209 (100 PHIGH: PEST 000210 000 747 GO TO 741 000211 000 PCAL: PEST/IARGU~••POW J 000212 000 fF < BAn~1PcAL1PAASE-WUNJ •Lf'.'•ALLOW1 GO TO 753 0302p goo 0SW T: CA .r.T. PAASf'.' 0 02 4 00 IF (OSWITJ PHyGH :-PEST 000215 000 IF ( .NOT. osw,r I PLOW= PEST 88H~l~ 888 753 10 !n 741 MACH: ~SQRT( PE T-WUNl•GPW/Twg/G+WUN 000218 000 fF I EST .Gr. wlJNJ WRITF:16, 481 MA~H 000219 000 F (PEST .LE, wUN) WRITE16,t,54 PES 000220 000 755 T: T+DT 000221 000 RETURN 000222 000 C1756 MAIN PR0611AM = 58,317 00022J 000 END 000221f 000 C WHOG•P •••••• INT •••••• QELT•L INT ELT 6B-Ol/llf-14:42 000001 000 000002 000 C INT 758 OOoOOJ 000 C 000004 000 000005 000 88888~ 888 000008 000 000009 000 8888H 888 000012 000 oooolJ ooo 000014 000 000015 000 000016 000 000017 000 000018 000 000019 000 000020 000 000021 000 000022 000 OOn02J 000 000024 000 000025 000 000026 000 000027 000 000028 000 000029 000 ooooJo ooo 0000J1 ooo 000032 000 ooooJJ ooo 000034 000 ooooJs ooo oonoJ0 ooo 000037 000 oonoJa ooo C C 11 u • 3 4 s 6 7 8 776 171 • 781 SUBROUTINE INT INT IMPLICIT DOUBLE P1:1ECISION (A-H,P-Z1,t.OGtCAL (01 COMMON I ARRAYS/ Ul2,20ll, Kl2•201),V(2,20l ),GC2,20l), COMMON Pl2,20l)• x12011, El2,2011• NCELLl?Oll, WE2CL(20t)•WE1CL(20ll,A 2011,r.cos 0) / PARM'1 / Ct , CO, G, GF", UL• IIR , GMW• GFMW, ENMAXo f"T, PJl,SLOSOR, SOKEXP, TMAXf• Two, WUN, ZERO. ~~~~i~·N~•~sg~·o~NT~PAfs~~.Nb~E~k!·~BLANf"o OPRJNT• OPUTI, OSKIPt OSPHFRtOTRACE• MWTAlL•MWHFAO,WVfL,WIDWAV,F"NWAV,RHEAO,PTAIL, WSLSOR•WSHEXP,MNWCOS,MXWCOS,OFWAV,REF•E?CL FORMATl2I5,4E12o5) FORMAT(' '•I5•10Et2o5l ET: -IR(l,N)••JP1)/GMW R?.P: Rll•l>••JPl UlMP: UI l, 1J 1J2MP: U12,t) IF(OTRACE)WRIT~(16,776)N,NL,U2MP,U1MP,R2P,R(l•N) DO 781 M: 1,NL RIP: R2P R~P: RlloM+t)••Jpl UlM: UlMP U2M: U2MP UlMP: UlltM+l) U?.MP: Ul2,M+l) U?. : IU1MP•U2MP+111M•U2M)/TWO . IF ( OT RACE I WR ITd 16,777 )M,E T ,F. I 1 ,M) ,ll2,R2P, R!P• VU, M), UlMP,U2MP 0 U1M,U2M ET : ET +IE'1•M)+112/TWOJ•(R2P-R1P)/Yll,M) RETURN ENO INJTIL : 716 MAN PROGpAM 238 QHOG•P •••••• PUOAT •••••• '1ELT•L pUDAT ~~~otr-ol/l 4oAi:-2 PUOAT SUBROUTINE PUOaT 000002 000 C 78~ PUOAT 000003 000 IMPLICIT OOURLE PpECISION(A-H,P-ZJ,LOGICAL(O) 000004 000 RF.AL TTT 000005 000 11 COMMON I ARRAYS/ ¥<2,201), R(2•201),VC2,?01 l,GC2,20lJ, 8888080' 8880 1 ~~ pw~2~flio1f~~~lll,iA¥J~ilio1~~t~b~1~AI· 000 8 00 ~ COMMON/ TIME/ TEpMIN, TIPUN, T• OT, OTL, OATE 011477 PAGf'.' :'I •••••• PUOAT 000009 000 000010 000 000011 000 000012 000 oooolJ 000 000014 000 000015 000 00()016 000 000017 ono 000018 000 000019 000 000020 000 000021 000 8888~~ 888 00()024 000 00002g 000 00002 000 000027 000 000028 000 000029 (100 000030 000 000031 000 000032 000 000033 000 000034 000 000035 000 030030 oog 0 0037 00 000038 000 000039 000 0000140 000 00004! 000 00004 000 000043 000 8888:i goo 00 000046 000 000047 000 000048 000 8888~3 888 · 000051 000 8888~i 888 000054 000 000055 000 000056 000 000057 000 000058 000 030053 888 o noo 000061 000 oono62 000 000063 000 000064 000 000065 000 00006& 000 000067 000 QHnG,P •••••• •••••• KRUN(3J, LABEL(7J I COMMON 3 / PARAM / Ct , CO, G, GF, UL, UR , GMW• GF'MW, ENMAX, ET, PJ1,SLOS0R, SOREXP, TMAXF• Two, WUN, ZERO, MINCOS, MAXCOS• J, JPI, N• NL, NLl,NDP, NPART, NSTEPS, OfNT, OESOR, OPEAK, OPLANf, 4 5 6 7 oPRINT• OPUTI, OSKJP, oSPHFR,oTRACE• . MWTAIL,MwHE~D,WVEL,WIOWAV,ENWAV,RHEAD,RTA L, WSLSOR,W~REXP,MNWCOS,MXWCOS,OEWAV,RFF,F2CL 8 22 COMMON / ARGINT / rNDEX, LSTART, NCYCLEt NFJNAL, · NsTORE, NS, NN, NNN, NPEAK, NSAM, NSHIF,NAUFF,NF'RE~,NWAVF' 7 8042 805 88~ ~ AU C 835 847 B849 849 C C RESTAR • • • • ,NPUNCH FoRMAT(l5•El5,g,2,S,El5,9,15,2r15.g) FORMATll5,El5,Q,I~,5E18ol3) . f8~~:f1: ::flg:1:Jl~:t~~a:~{l2F9 • 5 > OsAMP:NSTORE ,GT, l STORE CURPENT PROPERTIES MM: NF'RE9 IF ( MM,EQo O 1 MM: l LCAR: l WRIT[(l9,804)LCAR,T,NL,NCYCLE,ET•NPEAK,r,,GF LCAR: 2 R2 = IH2,l) u2 = 012,1, A?.:A 11) Do 847 M: l,NL Rl: R2 Ul: U2 At:A2 R?.: R(2,M+l) U?.: Ul2,M+l) A?.:A(M+l) DMID: (Rl+R2)/TWn UMUU: (Ul+u2>1TWn AMI0:(Al+A2)/TWO IF l,NOTo OSAMp) 60 TO 835 IF (M ,EQ. NL 1 GO TO 835 MX: (( M-NSHIF >1MMt•MM+NSHIF fF lMX,NEJM) r.0 8 847 WRI E( 9,805 LCAR,AM ,M,OMIO,UMUO,Pl2tM),V(?.,M)•El2,M) LCAR:LCAR+l TTT~~nl~8~ Ic:TTT IcY=Ic+l WRITEC2U,807>T,1Cy,NCYCLE,NL,NPEAK,P(2tNPEAK),V(2,NPEAK>• R(2•MPEAK),El2,NPEAK>,Ul2tNPEAK),A(NPEAK) WRITE(18,806)T,1Cy,MWHEAD,MWTAJL•P(2,1),Vl2J1),P(~,2), vi12~lo)!~ti!i~l~P~l:~~f:~l2~i~;:~1;!~oi,vYi~3ol:P(2,40), ENO v12,40>,P<2,sn,,v12,so>,P12,&o>,v<2,60> RETURN MAIN PROGpAM: 239 MAIN PROGQAM: 241 •••••• gELT,L RESTAR ~&~08r-01114oA3= 42 RESTAR SUBROUTINE REST~R 000002 000 C 851 REST AR 000003 000 IMPLICI~UOUBLE PPECISION (A-H,P-Z>,LOAJCAL (0) 000004 000 11 COMMON/ ARRAYS I U(2,20l>, Rl2•201),V 2•~01 ),0(2,201), 000005 000 IJ Pl2,201)• X1201), E(2,20t>, NCELL(20l>, 000036 003 2 WE2CL(201),WE1Cb(?.01),Al201),GC0S(l0) 0000 7 00 13 COMMON/ TIME / TE11MlN , TIPON, T• T, OTL, oso33s 808 l KRUN(3), LAAELi7> 0 O 9 O COMMON/ PARAM / C1, CO, G, GF, UL, UR t GMW• GFMW• ENMAX, ET, 000010 000 3 PJt,SLOSOR, SOREXP, TMAXF.• Two, WUN, ZERO• OAT!: 011477 PAGF: N N N •••••• RESTAH 000011 000 000012 000 8888H 888 000015 000 000010 000 000017 000 ououls ()00 000019 000 000020 000 000021 000 000022 000 000023 000 0000~1+ 000 0000 !> 000 000020 000 000027 000 000028 000 000029 000 ooooJo 000 000031 000 000032 ()00 OOn03J 000 000031+ 000 ooooJ~ 000 00003 000 oono~A 000 0000 000 000039 000 000040 000 00004i 000 00004 000 00004J 000 000044 000 oooo"i 000 00001+ 000 ooooitA 000 000016 000 000049 goo 0000~0 00 0000 ! 000 oooos 000 000053 000 gggg~~ goo 00 000056 000 000057 000 030058 goo 0 0059 00 000060 000 000061 000 00006 000 000063 000 000064 000 og1Jo6i 000 0 006 000 000067 000 000068 000 000069 000 000070 000 000011 000 0000 2 000 000073 000 00007'+ 000 000075 000 000076 000 000077 000 000078 000 000079 000 00008y goo ogoos 00 0 0082 000 ogoo83 000 0 008'+ 000 000085 000 · 000086 000 •••••• 4 ~INCos, MAXCOS• J, JP!, N• NL, NLl,NOP, 5 NPART, N5TEPS, DENT, OESOR, OPEAK• OPLANf, 6 OPRINT• OPuTI, OSKJP, OSPHF'"RtOJRACf• . 1 MWTAIL,MwHEAD,WVEL,WIOWAY,ENWAV,RHEAD,R AL, - 8 WSLSORtW~REXP,MNWCOS,MXWcOS,ofWAV,REFtF2CL 22 COMMON I AHGINl / rNDEX, LSTART, NCYCLF:,- NFJNAL, "NSTORE• NS, . 1 MNrNNN,NPEAK,NSAM,NSllJF,NRUFF,NFREO,NWAVE 2 ,NPUNCH 868 FORMAT(4IS,20J~.2At3l5> 8&9 FoRMAT(515,7A4) 870 FORMAT(3I5,JOJ~,2A/3035,28•I5) 871 FORMAT(t4•4028,161?D261 16,14,2D2811A) 872 FoR"IAT JOJ5,28/nJ5.2i,,215,DJ5,2H) 87J FORMAT 2035.2A• I5/2035.28,2L2> 874 FoRMAT(50216.18/4D~4.18) 875 FORMAT t3035,2A/2n35,28,415) IF(INOEX,EO.O) GO TO 892 C••• C •••• STOH~ DATA FOR LATF.R RESTARTJ~G RUN•••••••• C 886 NPUNCH: 1 NCYCLE: NCYCLE-1 T: T - OT WRITE (17,86Q)LST4RT,NCYCLF.,NPUNCll,NSTORE,NS,LAAEL WRlTE(17,86tt)NST~PS,NFINAL,NN,NNN,TERMIN,TIPUN,N8UFF,NFREn,NwAVE WtH E U 7,870) N• ,1, NLI, i, OT, OTL, UL t UR ,RF.:F ,NPF.AK WRITE(17,87J) CL•rO,NPART,r.F,r.,oESOR,OEWAV IF ( .NOT. OEW4V) GO TO A86 WH1TE(17,87~) wVEL•WIOWAV,FNWAV,WSLSOR, T W~REXP,MNWCOS,MXWCOS,MWTAIL•MWHEAO WRI E(17,874)(r.CO~(MC),MC:1,9) wRl~li~~!A7iT5g~~~gR!~o~~~P, • TMAXE, ENMAX,MJNC0S,MAXCOS,F2CL WRITF.(l7,874J(GCO~(MCJ,MC:1,9) . 887 DO 889 M: 1,N WRITE 117,871) M, Rtl,M), P(l,M)tV(l,M)t • UC1,M>,Oll,M),X(M),NCELL(M),WF.2CL(M),fll•M) 889 CONT IN U E ENDFILE 17 ENOFFILE 17 ENDFILE 19 END LE 19 GO To 932 C ••••• REaD RESTART CARDS 892 READ(15r8691LSTARTtNCYCLE,NPUNCHtNSTORE,NS,LARFL · RF.AD(15,668)NSTEP~,NFINAL,NN,NNN•TERMIN,TIPUN,NRUFF,NFREQ,NWAVE RFAD11s1s101 N,J•NLI,T,oT,nTL,uL,UR,REF,NPEAK NL= N- JPl: J+l PJl: JPl RF.AD(lS,~73) CL,Cn,NPART,GF,G,OESOR,OFWAV GMW: G-WUN GFMW: GF-WUN Rl: ZERO GAMW: GMW lF(,NOT, OEWAVJGO TO 895 HEAO(JS,875> WvEL,WIOWAV,ENWAVtWSLSORt 1 wSREXP,MNWCOS,MXWCOS,MWTAJL,MWHEAD RFAD(15,A74)(GC051MC)•MC=1•9> 895 IF( ,NOT. 0F.S0Q)G0 TO 896 RFAO (15, 872) SLnSOR,SOREXP, 896 • . TMaXE, ENMAX,MINCOS,MAXCOS,r.2cL READ(15r87•)CGC051MC),MC=1•9) • - DO 928 M: 1,201 IF (M.GT,N) GO TO 913 R?: Rl READ(lS,871) K,R1,P(1,M),V(1,M)•U(1,MJ, G(l,M), X(M)tNCELL(M),WF.2CL(M),E(ltMI RDfF: R1-R2 RI rM) : Rl GO TO 921 B2: M-N RCl•MI : Rl+R2•R0TF PllrM) : WUN V U(2•M) : U( ,Ml P(2•M) : P ~?~Jl~U~ ZERO R E T U R N M~IN PROG~AM: 53 ENO QHOG•P •••••• SAMPLE ••••*• lilELT,L SAMPLE ELT 6U-Ol/14-14:42 SAMPLE 000001 000 SUBROUTINE SAMPLE SIIMP1E 030002 goo c 93J o 0003 00 C 000004 000 ooooo~ ooo 88R88~ 888 000008 000 000009 000 000010 000 iiiiH iii 000013 000 00fl014 000 000015 000 000010 000 000017 000 000018 000 000019 000 000020 000 000021 000 000022 000 000023 goo 000024 00 000025 000 000026 000 8888jA 888 000029 000 ooooJo ooo 000031 000 000032 000 000033 000 000034 000 ooooJ~ ooo 00003b 000 000037 000 0000J6 000 000'039 000 ogoo40 goo 0 0041 00 000042 000 000043 000 000044 000 000045 000 000046 000 000047 goo 000048 00 0'00049 000 ooooso 000 000051 000 IMPLICIT DOUBLE P~ECISION <11-li,P-Zl ,LOIHCAL fOJ 11 COMMON/ ARRAYS/ U(2,201l, R(2•201J,V{2,?0t 1,0(2,2011, U P(2,20ll• X{20ll, E(2,20tl, NCFLL<201I• 2 WE2CL(201),WE1CL(201J,A(201J,GCOS( 01 COMMON / PARAM / C1 , CO, G, GF, UL, UP , GMW• GFMW, F:NMAX, ET, 3 PJl,SLOSOR, SOREXP, TMAXE• Two, WUN, ZERD, 4 MINCOS• r,tiAXCOS• J, ,JPl, N• NL, Nll,NOP, i MP~~i iN~~Tn~it l~E~!K Ia~ssiPH~~~~~kAit~ANF:' 7 MWTAIL,MWHEAD,WVEL,WIDWAV•FNWAV,RHF:AO,RTAIL, A WSLSOR,WSREXP,MNWCOS,MXWC0S,OEWAV,REF,E2CL 22 COMMON/ ARGINT / TNDEX, LSTART, NCYCLE, NFINAL, -NSTORF • NS, 7 NN, NNN, NPEAK• NSAM, NSHIF,NRUFF,NFREQ,NWAVE 2 ,NPUNCH 947 DATA GAIN /1.0010n / UATA TEST/1.0000lnO/ 949 FORMAT(' '•'964'•bl5,LS,5E15.5) 952 FORMAT{• •,•968••15,715,5E12o6J URE: ZERO 962 96J 964 965 ~BlF==w~~Ro osET = .FALSE. DO 964 ~R~ ~-YGE I: 5,NL Uc;E: Pf2,KJ IF (OSETJ GO TO 962 r81f ~ ~Bl~PRE IF{UGE .LT. TE~TJGO TO 96J IF IFDIF .GE.PnlFJ GO TO 963 osET = • TRUE. IF IUGE •LE.URFI GOTO 965 URE: UGE•GAIN IF(OTRACE)WRITF(16,949)N,I,K,NL,OSET,UGF•IJR~,FOIF,PDIF,PRE K: N-2 K: K+1 NPEAK: K lF(NPEAK eEA. NL) NPEAK: 1 NFREQ: K/NSAM NSHlF: K-NFREA•N~AM IF (NFREQoGT• n J GO TO 96A NFRIQ: 1 N<;H F: 0 968 F(OTRACEJWRITi:-(16,952)0SET,N,NSHIF,NFRF.Q,K,NSAM,NPEAK,I, 1 URE,UGE•r.AIN,FDIF•PDIF RF.TURN FlOIF:465 C1969 PAGF 2 N N ~ Appendix B Computer Pr~gram for Analyzing Data The calculation of the impulse and the eriergy integrals were performed by the following program which read and analyzed data stored on tape by the model. The tape is read from unit 10 and the input variables are read from unit : 5. The following unit 5 input variables must be specified: FIRST CARD ILINE: Number of time lines to be calculated MXWCOS: Cell number corresponding to the outermost cell of the source volume J: Geometry factor (0) Planar (1) Cylindrical (2) Spherical TSCALE: Dummy variable not used in this edition of program. RMAX: Maximum dimensionless radius at which impulse is calculated. TO: Value of last time line.Set to 0.0 for first data set. In addition to the printed output from unit 6 there are four other output units in which the output data is stored. Output unit 11 is for the impulse calculations, unit 12 is for pressure-time behavior at fixed Eulerian radius, unit 13 is for the energy distribution calculations and unit 14 stores the positions of selected particles for plotting of particle dis- placement. 225 226 AMAHJ(lq) I~PLICIT LOGICAL (0) 0IMENSION P(401) •R(401) ,U(401) ,V(401) ,1:(401) ,A(401), * RALl<-f( 102) ,THALE( 102) ,HPAKE(lO?.) ,ETOTAt ( 10?) ,TT(t02) * RALlf:(102l ,AlRKEq02) ,l'\IRIF(l02) ,RPRT(i02,24), ' * RR(~),RL(1n5),Al~P(l05),0IMP(l05),PP(404,5) TIME(404) RFAD(5,9)1LINE,MXWCOS,J,TSCALE,RMAX,TO ' WR I TE ( 6 , 9 ) I L I 1-.J F , M X l'I COS , J , TS CALF.: , R ~AX , TO Y FORMAT() EMAX:O. Pt=Acos<-1. > EsCAL:(PI*l60•/3•>**(1./3.> AySCAL:SQkT<1•4)/ESCAL MX wcp l=tv· x wcos+ 1 JPl=J+l Ii,.,LT:o lrviAX:102-2 lTEST:ILI~E/I~AX+l uo 5!:> I:t,4n1 P:T TMAX:f oT=T-To To=T It->= l IF(IO .EQ. 1)60 TO g6 JTEST:(10/ITEST)*ITEST lF(IO .N£. JTEST)OPLOT: .FALSE. 9b WRITE(6,98)IO,JTEST,NL,NCYCLE,1T,OPLOT,T,or,ro,FT 9b FoRMAT(~I~,L5•4Fl0.5) L)Q 199 I=l•ML RF AO ( 1 0 , 1 !S 9 , E 1•,.in = 991 • ERR: 9 8 9 ) LC AR • A ( I ) , M , R ( l ) , U ( I ) , p ( I ) , * v,E • 159 FoRMAT(I5,El~.q,I5,5ElA.13> C WR I TE < o, 1 !:>9) LC AR• A ( I ) , M, R ( I ) , U ( I ) , P ( I ) , V ( I ) , E ( I ) 19lJ CoNrlNUE !R:1 oo 249 1=1•401 C WRITE(6,218)I•IR•IP,Il,OIMP(IR),R(I>,RL(IR)•RO,RR(lP), C • P(I),PO 227 2Ib FoRMAT(4I5,L5,6Fl0.5) C ***** STORE UATA FOR P-T CURVES***** lF(IP .GT. S)GO TO 219 IF(Rll> .LT. RR(IP))GO TO 219 DR:H(I)-RO URl=RL(lR)-RO up:P(I)-PO PR=Po+(DP•nR1/nR> PP(lT,IP):-PR lF(PR .GT. PMAX)PMAX:PR Ip:lP+l lF(I .GT. NL)GO TO 249 C ***** CALCULATE I ~~ULS~ ***** 21':;l lf-:(R(I) .LT. HL(IR))GO TO 241 22l1 IF( .HOT. 0 tMP(IR))GO TO 239 OR=k - RO URl=RL(lR)-RO lJ.:,=P(I)-PO P~=PO+(DP•~Hl/nR>-1. 1F(PH .LT. n.)GO TO 229 AIM~(lR):AIMP(IR)+PH•DT 22~ lF(~(l) .LT. l.)OIMP(IR) : .FALSE. 2~~ If.~Ml:IR lR=lR+l 1F(IH .GE. IHMAX)GO TO 259 241 CoNfINtJE C ~RIIE(6,244)I,rH,IRM1,0IMP(IHMt),HO,PO,AIMP(IRM1),PR,OP,OR1,nR, C * RL(IR~l),k(l),P(I) · lFlRL(lP) •LT• R(I))GO TO 220 Po=~(I) Ro=~ . _ 244 ~OH ~Af(3I~,L5•10fl0.5) 24~ CoNT1f',1Uf. 2 'i ':;l 1 F ( • I\JO T • OPL. OT) <;O TO 9 79 IµLl:IPLT+1 TT(lf-'LT):T C ***** SlORE OATA FOR PARTICLE PATHS***** .J:0 uo 31 O I= 1 , s:; C w RI IE ( 6, 3 U o > IO, I , J, I PL T, R ( I ) , TT ( I PL T) .3 0 '.1 t,: () R ,·,1 A r ( 4 I 5 , :? F 1 n • ~ ) 3IU KPHf(IPLT,I):K(I) J:5 Do 410 I:lU,50,5 .J:J+I C ~RI1E(6,40q)t,J,RPRT(l) 41U HPHf(IPLT,J)=K(I) 40':I FoR1·1AT(2I!::>,F10.':>) UO 510 I=6U,l!::>0,10 J::J+l C WRITE(n,509)1,,.J,R(I) 50~ FoRMAT(2I5,F10.5) 510 HPRl(IPLT,J):H(I) C ***** CALCULATE ENEPGY INTEGRAL***** RoP=o. bTMASs::o. ATMASs=o. TMASs:o. HALKf.(IPLT):Q. bALlE(IPLT):Q. AikK.F.(IPLT):O. AIRlE.(IPLT):O. AIRMJIH:l.l(G-1.) RnP=o. ~'-'=O • Uo 599 MC:1,MXWCOS RNEW:2.*R(~C)-ROP RNP=RNE~o*JPl RMASS=CRNP-RJ)/VCMC) Uc:;Q:lJ(MC)**2 GAMW:P(MC)•V(MC)/E(MC) lF(IO .LT. 16 .OR. IO .GT. 20>GO TO 555 C WRITE(6,549)10,MC,RMASS,US0,AALKE(IPLT),GA~W,RNFW,VCMC), C 54 u* R(MC) ,tl(MC) ,ECMC> ,PCMC> · 7 FORMATC215,10Etl•3) C C C C C C C 8q';I 97'7 98'~ Ygu 991 9gt:'. 100b 100':I 1014 1019 102Y 1051:, 106'::I 107"=' * * * * 228 bAL AMlj: 1 1 /GAM1v CFLLKE:((USG*KMASS)/2,) BAL~E(IPLT):HAI.KE(IPLT)+CELLKF CFLLIE:(E(~C)-RALA~8)*RMAS5 ti AL IE ( IPLT) =~lAI. IE< IPL T) +CELL IE HTMASS:bTMASS+RMASS H.J:Kl' IP HoP=I-WEW wRirE<6,~HA)IO,MC,IPLT,ROP•RJ,RTMASS,R~ASS,RALIE**2 GtV"1;~=G-1, CFLLKf= ,CELLKF.,GA1"1w,IJSl~ ~ CUl'JTIMJE ETOfALCIPLT):rlALKE(IPLT)+~ALltCIPLT)+AlRKE(IPLT)+AIRJF(IPLT) lF(ETOTAL(lPLT) .GT, EMAX)FMAX=fTOTAL(IPLT) . T~AS5:A1MASS+8TMASS r~ALE ( IPL. T > =H~ I_KE:: ( IPL T) +RALIE ( IPL T) tJPAKE ( IPL T) =Tri ALI:. ( IPL T) +A IPKE ( TPL 1) WH I i F ( 6, f'\Yll) IO, I PL T, ~,c, NL, TT< I PL T) , FTOT AL ( T P1 Tl , RPAKE ( IPL T) TdALE(IPL T),HALIE(IPLT>,RALKE(IPLT),AIPIF(TPLT), · ' AIRKE(ll-'LT),TMASS,ATMASS,RTMASS F oR ·"tA T ( 4 I 5, 11 F q • ~) Cof\lT I r~lJE w R 1 r[ ( 6, g~ I') ) IO, LC AR , M •NL, T • A ( I ) , ~ l I ) , P ( I ) , \/ ( I ) , ll ( I ) FnR ~AT('FILE EPRUH•,4I5,6F10.5) wRI1E(6,qY?)LCAR•T,NL,NCYCLE,fT,NLT,G,GF F o F~ '·.,AT ( • f N n OF F l LE • , I 5 , F 1 0 • 5 , 2 I 5 , F 1 n • i:; • I 5 , 2 F 1 O , 5 ) wHIIE(11,1008)IHMAX,AI5CAL FoH MA T(I~,FlU,~) w R I TE ( 11 , 1 0 O 9 ) ( I , O IMP ( I ) , R L ( I ) , A HlP ( I > , I= 1 , IRMA X ) FoRMAT(I5,L5,2F2U,10) WR I It ( 12, 10 14 > IT, ( RH ( I ) , I: 1, 5 > FoR~AT(I~,5FlU,5) ~JR I TE ( 12, 10 19 > ( I •TI ME ( I ) , ( J • PP C I • .J) • J: 1 , 5) , I: 1, J T) FoH~AT(6(1~,Fl0.5)) wRITE(13,1029)IPLT,EMAX FnR~AT(l~•Fl5e10) WRITE(13,1059>Cl•TT,BALKF(I),AALIE,ATRIECI), AIRKE(I),TRALE(l),APAKE,ETOTAL(I),I:1,IPLT) FoRMAT(I5,8F15,10) wRITE(14,1069>I0•1PLT FoR1VtAT(2I5) ~RITE(14,1079)(1I,TT(lt),(RPRT(Il,I), 1:1,2~>,Il:l•IPLT) FORMAT(I5,12F9.5/13F9,5) WR1TE(6,1Ub9)I0,1PLT WRITE(6,107g)(Il•TT,I:1,24)•JI:1,IPLT) sroP END ...... , Roman D e eC? 1. eo E EB ES ET f h fl.ho f fl.H C h. 1. h' NOMENCLATURE Speed of sound--ambient Speed of sound--ahead of(before) shock wave Speed of sound--behind(after) energy addition Concentration - mass fraction Constant pressure heat capacity Constant volume heat capacity Newtonian speed of sound--ambient Chapman Jouguet condition Lagrangian distance Beginning of rounding term in energy source volume-- Lagrangian distance Extent of energy source volume--Lagrangian distance Width energy addition wave--Lagrangian distance Internal energy Energy of formation--species i Internal energy--ambient Non-dimensional internal energy Energy remaining within the source Energy transmitted to the surrounding gas Total amount of energy deposited at the source Body force vector Enthalpy Effective zero point energy Heat of combustion Enthalpy--species i Enthalpy-working fluid heat addition model 229 230 h 1 Enthalpy-ahead of shock front h4 Enthalpy-behind energy addition i Species I+ Positive phase impulse I Non-dimensional positive phase impulse j Geometry factor K Linear spring constant me Mass M Mach number M1 Mach number-approach flow Ms Mach number-expanding sphere M Mach number-energy addition wave w n moles of gas within the source volume p Pressure p8 Shock pressure P0 Pressure--ambient pl Pressure--ahead of shock P2 Pressure--behind shock P3 Pressure--behind (ahead. of (before) energy addition p4 Pressure--behind(after) energy addition P Non-dimensional pressure P* Non-dimensional dissipative pressure Ps Non-dimensional shock overpressure q Source energy density Q Heat transfer rate Q Heat release during a constant gamma process Qc Heat release/unit mass of fuel r R E R s t a T u u 231 Non-dimensional amount of energy deposited at the origin · · · · Ene!gy/unit mass deposited at the origin Radial distance coordinate Initial source radius Gas constant Energy-scaled shock position Shock position Energy scaling distance Time Time of shock arrival Source deposition time Time at which maximt.nn structural displacement occurs Characteristic acoustic propagating time Time end of positive phase Time--end of negative phase Characteristic loading time Particle velocity Non-dimensional particle velocity Non-dimensional energy wave velocity Volume of the source Initial source volume Flow velocity vector Wave width-energy addition wave Comparable weight of tri-nitro-toulene Weight of explosive material Weight of hydrocarbon explosive Similarity lines 232 Greek y Specific heat ratio Yo Specific heat ratio--ambient Yi Specific heat ratio--ahead of(before) energy addition Y4 Specific heat ratio--behind(after) energy addition n Non-dimensional distance coordinate 0 Temperature 00 Temperature--ambient 01 Temperature--ahead of shock front 0 2 Temperature--behind shock front 84 Temperature--behind energy addition A Energy source term A Non-dimensional energy source term v Specific volume vf Sp'ecific volume expansion ratio v0 Specific volume--ambient ~ Energy addition wave parameter M Energy wave structure parameter IT Artificial viscosity term p Density Po Density--ambient P1 Density--ahead of shock front P2 Density behind the shock front cr Length scaling factor T Non-dimensional time 233 Tc -Non-dimensional cell deposition time TD Non-dimensional source volume deposition time TT Non-dimensional energy wave source volume transit time ~ Non-dimensional specific volume w Natural frequency ~ Non-dimension energy wave Mach number 1. 2. 3. 4. 5. 6. 7. 8 . 9. 10. 11. 12. 13 . 14 , BIBLIOGRAPHY Strehlow, R.A., "Unconfined Vapor-Cloud Explosions-- An Overview," 14th Symposium (International) on Combustion, The Combustion Institute, pp 1189-1200(1973), Baker, W.E., Explosions In Air, University of Texas Press, Austin, Texas (1973). Taylor, G.I., "The Air Wave Surrounding an Expanding Sphere," Proc. Royal Society, Al86, pp. 273-292 (1946). Sedov, L . I., Similarity and Dimensional Methods in Mechanics, Academic Press, New York, N.Y. (1959) . Bethe, H.A., Fuchs, K., Hirschfelder, H.O., Magee, J.L . , Peierls, R.E., and Von Neumann, J . , "Blast Waves," LASL 2000, Los Alamos Scientific Laboratory (1947). Sakurai, A., "Blast Wave Theory," in Basic Developments in Fluid Mechanics, Vol I, (Morris Holt, ed.), Academic Press, New York, N.Y. (1965). Oshima, K., On Exploding Wires (W.B. Chance and H.K. Moore, eds.) Vol . 2, Plenum Press, New York, N.Y, (1967) . Liepmann, H.W. and Roshko, A., Elements of Gas Df:amics, John Wiley and Sons, Inc., New York, N.Y. (1967 . Von Neumann, J. and Richtmyer, R.D., "A Method for the Numerical Calculation of Hydrodynamic Shocks", J. of Appl. Phys., 21, pp 233-237 (1950). Lax, P.D. and Wendroff, B., "Systems of Conservation Laws," Comm of Pure and Applied Math., 13, p 217 (1960). Richtmyer, R.D. and Morton, K.W., Difference Methods for Initial Value Problems, Interscience Publishers, Inc. , New York, N.Y. (1967) . Von Neumann, J. and Goldstine, H., "Blast Wave Calculation," Comm. Pure and Appl. Math., ~. pp 327-353 (1955) . Brode , H.L., "Numerical Solutions of Spherical Blast Waves," J. Appl . Phys., 26, pp 766-775 (1955). Ricker, R. ! "Blast Wav7s from Burs t ing Pres·surized Spher es," M. S. Thesis·, Aeronautical and Astronautical Engineering · Dept . , Univ. of Illinois, Urbana-Champaign, Ill . (1975) . 234 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 235 Zajac, L.J. and Oppenheim, A.K., "The Dynamics of an Explosive Reaction Center," AIAA Journal, 9, pp 545- 553 (1971). - Freeman, R.A., "Variable Energy Blast Waves·," British Journal of Appl. Phys., Ser. 2, l, pp 1697-1710 (1968) . Dabora, E.K., "Variable Energy Blast Waves," AIAA Journal, 10, pp 1384-5 (1972). Adamczyk, A.A., "An Investigation of Blast Waves from Non-ideal Energy Sources," Ph.D. Thesis, Aeronautical and Astronautical Engineering Dept., Univ. of Ill, Urbana Champaign, Ill. (1975) . Rudinger, G., Nonsteady Duct Flow, Dover Publications, Inc., New York, N.Y. (1969). Oppenheim, A.K., Kuhl, A.L., Lundstrom, E.A., and Kamel, M.M., "A Parametric Study of Self-similar Blast Waves," J. Fluid Mech., 52-4, pp 657-682 (1972). Kuhl, A. L. , Kamel, M. M. , and Oppenheim, A. K. , "Pressure Waves Generated by Steady Flames," XIV Symposium (International) on Combustion, The Combustion Institute, pp 1201-1215 (1973). Strehlow, R.A., "Blast Waves Generated by Constant Velocity Flames: A Simplified Approach," Combustion and Flame, 24, pp 257-261 (1975). - Williams, F.A., Combustion Theort, Addison-Wesley Publish- ing Company, Inc., Reading, Ma. 1965). Strehlow, R.A., Fundamentals of Combustion, International Textbook Company, Scranton, Pa . (1968) . Zajac, L.J. and Oppenheim, A.K., "Thermodynamic Computations for the Gasdynamic Analysis of Explosion Phenomena," Combustion and Flame, 13, pp 537-550 (1969). Hopkinson, B. , British Ordnance Board Minutes, 13565 (1915). Strehlow, R.A., and Baker, W.E., "The Characterization and Evaluation of Accidental Explosions," Report NASA CR-134779 (1975). Baker, W.E., Westine , P.S. and Dodge, F.T., Similarity Methods in En ineerin D namics: Theor and Practice of Scale Modelling, Spartan Boos, Rochel e Par , N. J. 3). Wilkins, M.L., "Calculation of Elastic-Plastic Flow" UCRL-7322, I ,Rev. 1, Lawrence Radiation Laboratory 'Livermore CA ( 19 6 9 ) . ' ' ... 236 30. Oppenheim, A.K., "Elementary Blast Wave Theory and Computations," Proceedings of the Conference on Mechanisms of Explosion and Blast Waves, Paper No. 1, Yorktown, Va (1973). 31. Boyer, D.W., Brode, H.L., Glass, I.I., and Hall, J.G., "Blast From a Pressurized Sphere," UTIA Report No. 48 (1958). 32. Huang, S.L., and Chou, P.C., "Calculations of Expanding Shock Waves and Late-Stage Equivalence," Final Report 125-12, Drexel Institute of Technology, Philadelphia, Pa. (April 1968). 33. Brode, H.L., "Blast Wave From a Spherical Charge," Physics of Fluids,~. pp 217-229 (1959). 34. Brinkley, S. R., "Determination of Explosive Yields," AIChE Loss Prevention, 1, pp 79-82 (1969). 35. Fishburn, B.D., "Some Aspects From Fuel-Air Explosions," Acta Astronautica (in press (1977)). 36. Fishburn, B.D., Private Corrnnunication (1977). CURRICULUM VITAE Name: Robert Thomas Luckritz. Permanent address: 90-36 Avenue North Clinton, Iowa 52732. Degree and date to be conferred: Ph.D., 1977. Date of birth: February 24, 1943. Place of birth: Clinton, Iowa. Secondary education: Clinton High School Clinton, Iowa 11 June 1961. Collegiate institutions attended Dates United States Coast Guard 1961-1965 Academy University of Maryland 1969-1971 University of Maryland 1971-1976 Major: Chemical Engineering, Degree B.S. M.S. Ph.D. Professional positions held: June 1973-Dec. 1975 Lieutenant Commander U.S. Coast Guard Date of Degree June 9, 1965 June 5, 1971 May 14, 1977 Cargo and Hazardous Materials Division U.S. Coast Guard Headquarters Washington, D.C. 20590.