University of Southern California Dissertations and Theses

Numerical study of focusing effects generated by the coalescence of multiple shock waves

(USC Thesis Other)

Numerical study of focusing effects generated by the coalescence of multiple shock waves

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Numerical Study of Focusing Effects Generated by the Coalescence of Multiple Shock Waves Dissertation by Shi Qiu In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy (MECHANICAL ENGINEERING) USC Graduate School University of Southern California Los Angeles, California August 2017 c Shi Qiu 2017 All rights reserved Acknowledgments I would like to express my sincere gratitude to my advisor, Dr. Veronica Eliasson, who has brought me into the scientific world. Until now, I still remember the first time when Iwenttoherofficehoursandthatwasthebeginningofmyshockwavejourney. Without her encouragement and guidance, none of this work would be possible. I am also very grateful to my committee member, Dr. Aiichiro Nakano, for his extensive help in high-performance computing technique. His patience, motivation and immense knowledge havebeeninvaluabletome. Heistrulythedefinitionofanexcellent professor. Many thanks go to Dr. William D. Henshaw for his guidance on the software named Overture. Without his help, none of the numerical simulations would be possible. I would like to thank high-performance computing center in University of Southern California for providing free computing resources for my project. I also thank the other committee member, Dr. Ivan Bermejo-Mereno, for his insightful comments, which has incented me to widen my research from various perspectives. I am indebted to all my labmates, colleagues and friends forgiving me advice on this work and helping me practice with my presentations. I will keep from naming anyone specifically for the risk of leaving someone out. You know who you are and I hope you understand how deeply I treasure the experiences we have shared during the past six years. Last but not the least, I would like to thank my parents for their unconditional support for my overseas study. Without their love and constant support, I will never be able to accomplish this. ii Abstract Field of shock focusing has been studied for a long time, and its application has been expanded to a lot of areas such as industrial and civil engineering, biomedical, anti- terrorist and aerospace industry. Due to the nonlinear nature of shock waves and their interactions with nonuniform or moving media, solid or porous surfaces or other shock waves, predictions of shock dynamics and shock focusing events are far from trivial. Our current study aims to deepen the understanding of shock wave focusing processes by studying the interaction of multiple cylindrical shock fronts. The study has been performed using a combination of numerical, theoretical and analytical tools. For the numerical part, simulations of multiple point energy sources were performed, and the resulting shock interaction and coalescence behavior were explored. Three to ten energy sources were placed concentrically around the target, and conditions at the target area were monitored and compared to those obtained using a single source. For each simula- tion, the energy summed over all sources was kept constant, while the radial distances between target and energy sources and the energy source initiation times were varied. Each energy source was modeled as a point source explosion. The resulting shock wave propagation and shock front coalescence were solved using the inviscid Euler equations of gas dynamics on overlapping grids employing a finite difference scheme. Results show that multiple energy sources can be beneficial for creating extreme conditions at the in- tended target area; over 20 times higher peak pressure is obtained for ten simultaneous sources compared to a single source. Moreover, peak pressure at a point away from the target area is reduced by more than a factor of three. The ultimate goal is to develop general strategies that can serve for various types of design objectives on the coalescence of multiple shock waves. Thus, a reduced-order model has been developed, and an optimization framework has been proposed to help with various decision-making processes. To be specific, the model consists of two parts, iii which can predict the early stage (expansion) and the late stage (focusing) of the in- teraction of multiple shock waves, respectively. An analytical method was adopted to describe the expansion of cylindrical or spherical shock waves. The sonic criterion was theoretically studied to predict the transition from the expansion stage to the focusing stage. To predict the focusing stage, Geometrical shock dynamics (GSD) theory was numerically implemented on parallel computers using a spatial decomposition method. In addition, an efficient tridiagonal system solver for massively parallel computers was applied to resolve the most expensive function in this implementation, resulting in an efficiency of 0.93 while using 32 HPCC cores. Furthermore, a modification was made to the GSD theory to improve its accuracy. Comparisons of shock front Mach number history with both the Euler equations and existing experimental data show good agree- ment and the computational cost of the new model is significantly reduced compared to the Euler equations. iv Contents Acknowledgments ii Abstract iii 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Category of Shock Waves . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Previous Studies onShockWaves Followed byanExponential Decay 4 1.1.3 Shock-Shock Interaction . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 Outline and Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2 InteractionandCoalescenceofMultipleSimultaneousandNon-simultaneous Cylindrical Shock Waves 13 2.1 Numerical Scheme and Initial Conditions . . . . . . . . . . . . . . . . . . 13 2.1.1 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.2 Verification of Initial Conditions . . . . . . . . . . . . . . . . . . . 17 2.1.3 Convergence Rate of Numerical Solutions . . . . . . . . . . . . . . 17 2.2 Problem Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Results and Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3.1 Case 1 - Simultaneous Initialization . . . . . . . . . . . . . . . . . 20 2.3.2 Case 2 - Delayed Source Initialization . . . . . . . . . . . . . . . . 22 2.3.3 Case 3 - Semicircle Placement . . . . . . . . . . . . . . . . . . . . 28 2.3.4 Center of Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.5 Energy Dissipation and Shock Reconfiguration . . . . . . . . . . . 31 2.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 v 3 Interaction of Two Cylindrical Shock Waves with Equal Strength 35 3.1 Division on Interaction of Multiple Shock Waves . . . . . . . . . . . . . . 35 3.2 Transition Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.3 Analytical Prediction on the Expansion of Shock Waves . . . . . . . . . . 40 3.4 Optimization on the Number of Shock Waves . . . . . . . . . . . . . . . 44 3.5 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4 A Reduced-order Model for Shock Wave Focusing 48 4.1 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1.1 Serial Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.1.2 Parallel Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Hybrid MPI Combined with OpenMP Approach . . . . . . . . . . 54 MPI Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1.3 Comparison with Analytical Solution . . . . . . . . . . . . . . . . 56 4.1.4 Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.1.5 Symmetric Boundary Condition . . . . . . . . . . . . . . . . . . . 61 4.2 Comparison with the Euler Equations . . . . . . . . . . . . . . . . . . . . 62 4.3 Modified GSD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3.1 Area-Mach Number Relation . . . . . . . . . . . . . . . . . . . . . 67 4.3.2 Modification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.4 Comparison with Previous Experiments . . . . . . . . . . . . . . . . . . . 73 4.5 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5 Conclusions and Future Directions 80 5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Bibliography 84 vi List of Figures 1.1 Chronology of Shock Focusing Research. . . . . . . . . . . . . . . . . . . 3 1.2 Profiles of Different Shock Waves . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Interaction of Two Cylindrical Shock Waves with Equal Strength . . . . 7 1.4 Reflection Patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.5 Schlieren Visualizations of Two Cylindrical Shock Waves Interaction. . . 11 2.1 Normalized Initial Conditions Based on Taylor’s Similarity Law . . . . . 16 2.2 Comparison of Pressure History 5 mm Away from the Shock Center for Our Results and Those of [42]. . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Schlieren Visualizations of Cases with Different Number of Sources. . . . 21 2.4 Schlieren Visualizations of Ten Sources Reconfiguration. . . . . . . . . . 23 2.5 Schlieren Visualizations from the Case with 3 Sources with One Source Delayed by Three Delay Time Units. . . . . . . . . . . . . . . . . . . . . 25 2.6 Schlieren Visualizations from the Case with 5 Sources with One Source Delayed by Three Delay Time Units. . . . . . . . . . . . . . . . . . . . . 26 2.7 Schlieren Visualizations of Shock Front Propagating from 10 Sources Set in Semicircle (10SC). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.8 Pressure Plot along Horizontal Line Centered at Target T1 for Delayed and Non-delayed Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.9 Pressure History at Point S1 for 5C, 10C and 10SC. . . . . . . . . . . . . 31 2.10 Comparison of Non-dimensional Peak Pressure for Cases with Different Radii. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1 Schlieren Visualizations of Expansion Phase with Six Sources . . . . . . . 36 3.2 Schlieren Visualizations of Focusing Phase with Six Sources . . . . . . . . 36 vii 3.3 Schematic Illustration of the Wave Configuration in the Frame of Reflec- tion Point R. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Comparison between the Sonic Condition and Actual Transition from Numerical Results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.5 Schlieren Visualization Illustrating the First Numerical Setup. . . . . . . 41 3.6 Schlieren Visualization Illustrating the Second Numerical Setup. . . . . . 42 3.7 Comparison between the Analytical Method [4] and the Euler Equations for a Cylindrical Shock Wave Expansion Progress. . . . . . . . . . . . . . 43 3.8 Schematic Illustration of Shock-Shock Interaction Configuration between Two Expanding Cylindrical Shock Waves . . . . . . . . . . . . . . . . . . 43 3.9 An example of Shock Front Mach Number versus Shock Interaction Angle. 44 3.10 Transition Conditions on Different Numbers of Sources. . . . . . . . . . . 45 3.11 Schlieren Visualizations of A Case with Three Unsymmetrical Sources . . 47 4.1 Flow Chart for Serial Version. . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Flow Chart for Hybrid MPI and OpenMP Scheme. . . . . . . . . . . . . 55 4.3 Example of Spatial Decomposition of a 5-sided Polygonal Shock Front Split into 21 Points Distributed among Five Processors. . . . . . . . . . . 56 4.4 Flow Chart for MPI Scheme. . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.5 Converging Shock Propagation Showing Cuccessive Positions for an En- neagon and a Decagon. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.6 Wall-clock Time of Strong Scaling for Hybrid Scheme Using 3,242,911- point Shock Front on β Nodes and β Processors per Node. . . . . . . . . 59 4.7 Wall-clock Time of Strong Scaling for MPI Scheme Using 3,242,911-point Shock Front on N node Nodes. . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.8 Wall-clockTimeofScaledWorkloadsUsing25,991N node -pointShockFront on N node Nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.9 Dodecagon Shock Front Propagation at Successive Time Instants. . . . . 62 4.10 Pressure Contour Illustrating Reconfiguration Process from the Cases with 5 (left) and 10 (right) Sources. . . . . . . . . . . . . . . . . . . . . . 63 viii 4.11 Schlieren Visualizations of Initial Shock Geometry from the Cases with 5 (left) and 10 (right) Sources. . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.12 Mach Number History Comparison for 5 Sources Case. . . . . . . . . . . 66 4.13 Mach Number History Comparison for 10 Sources Case. . . . . . . . . . . 67 4.14 Comparison of Mach Number as Function of Shock Radius for Regular GSD, Extended GSD and Analytical Result [4] for an Expanding Cylin- drical Shock Wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.15 A Binary Tree Example to Illustrate Different Expansion Types During the Interaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.16 Overlay of Euler Simulation and GSD Results. . . . . . . . . . . . . . . . 77 4.17 RatioofMaximum Pressure P m, atthe Mach Stem to Ambient Pressure P a , versus Time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.1 Optimization Framework. . . . . . . . . . . . . . . . . . . . . . . . . . . 83 ix List of Tables 2.1 Estimated Convergence Rates Using Second-order Godunov Scheme. . . . 18 2.2 Comparison of Non-dimensional Peak Pressure and Arrival Time at T1 and S1 for Cases with Multiple Simultaneous Cylindrical Shock Waves . 22 2.3 Comparison of Non-dimensional Peak Pressure and Arrival Time at T1 and S1 for Cases with Multiple Non-simultaneous Cylindrical Shock Waves 27 2.4 Comparison of Non-dimensional Peak Pressure and Arrival Time at T1 and S1 for Semicircle Case. . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.1 Computational Cost of Each Function for Three Different Tests. . . . . . 54 4.2 Comparison of Numerical Results with Analytical Solutions from [72]. . . 57 4.3 Computational Cost and Speedup of 1/2, 1/4, 1/24 of the Full Geomet- rical Simulation Using Symmetric Boundary Conditions. . . . . . . . . . 62 4.4 A Comparison of Shock Arrival Times of the Current Euler Simulation versus the Experiment Presented in [35]. . . . . . . . . . . . . . . . . . . 76 x Chapter 1 Introduction The focusing of shock waves is a phenomenon known to generate extreme conditions in the focal region. In nature, the curved shock fronts that appear in a variety of circumstances can lead to focusing effect during propagation. In the ocean, the collapse of cavitation bubbles is one of the most prevalent ways for shock focusing to occur. For large swimmers at shallow depth, the maximum speed of the swimmers is limited by the physical pain caused by cavitation collapse near the caudal fin [40]. In outer space, the interaction of shock waves propagating in the interstellar medium (ISM) that are generated by supernova explosions with interstellar clouds plays a significant role to the evolution of the ISM [48]. On the other hand, in human society, many fields such as industrial and civil engineering, chemical, nuclear, biomedical and aerospace industry can benefit from shock focusing. One of the most common applications is extracorporeal shock wave lithotripsy (ESWL) [54]. It is a minimally invasive therapy to pulverize kidney stones into smaller fragments. During the procedure, shock waves are generated and focused to the localized area in a patient’s body, a process that will lead to stone fragmentation. However, in some applications, focusing of shock waves can be detrimental. One example of this is shock loading of marine structures from underwater explosions. There exist converging geometries on naval vessels such as rudder-hulljunction,propellershaft,andbowthruster,andhenceshockwavesgenerated from underwater explosions can be focused by these structures and cause severe damage [89]. Thus, studying the physics behind the mechanisms of shock focusing can not only help us to understand the nature better but also improve many applications valuable to human life. There exists a large body of research on shock wave dynamics, in two and three dimensions, using analytical, numerical and experimental tools. Shock wave interaction 1 may result in shock focusing, in which multiple wave fronts coalesce and interact to cre- ate extreme conditions. These conditions depend upon the shock media, the speed of the shocks andthe geometrical pattern ofthe coalescent shocks. An overview ofthe chronol- ogyofsome ofthestudies doneinthe fieldofshock focusing isrepresented inFigure1.1, with simulations and experiments listed on the left and right side, respectively. Most studies have focused on stability issues of converging shocks in two dimensions. This research area was initiated by the analytical work of Guderley in the early 1940’s [28] where he predicted a power law increase in the Mach number as the shock propagates towards the focal region. About a decade later, Perry and Kantrowitz [60] studied con- vergent shocks experimentally using a shock tube setup. They noticed light emission due to high temperatures as the shocks converged. These results sparked an interest among researchers to learn more about converging shocks using analytical, numerical and experimental approaches. For the majority of the experimental work, shock tubes were used toinitiateconverging polygonalshocks andstudy theirbehavior astheyprop- agate towards a focal region [24, 45, 84, 88]. Numerical simulations of similar scenarios have been performed using the Navier-Stokes equations, the inviscid Euler equations of gas dynamics or Whitham’s geometrical shock dynamics [3, 5, 21, 72, 91, 92]. 1.1 Background 1.1.1 Category of Shock Waves Based on the flow properties behind the wave front, there exist two types of shock waves. The first type is followed by constant properties lasting for an extended period. A tremendous number of experimental studies on shock focusing have utilized this kind of shock wave, most likely because it can be easily produced from a constant-area shock tube with high repeatability. However, in reality, the second type, which is followed by an exponential decay in properties, occurs in most circumstances, such as lightning, cavitation, and explosion. Typically, any situations that can release a certain amount of energy in an open area will result in the second type of shock. Figure 1.2 shows the pressure-time curve at a fixed location for these two types of shock waves. The former 2 1940 1950 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 Experiments Simulations Perry (1951) [60] Knystautas (1969) [49] Roig (1977) [66] Wu (1981) [93] Saito (1982) [68] Takayama (1984) [84] Keefer (1984) [44] Matsuo (1985) [53] Neemeh (1986) [56] Watanabe (1995) [88] Hosseini (2005) [39] Kandula (2008) [43] Eliasson (2007) [24] Kjellander (2011) [45] Wang (2014) [89] Liverts (2016) [51] Guderley (1942) [28] Brode (1977) [16] Hikida (1981) [36] Schwendeman (1987) [72] Book (1990) [15] Demmig (1990) [20] Aki (1990) [1] Starkenberg (1994) [81] Apazidis (1996) [3] Betelu (2001) [9] Schwendeman (2002) [76] Dimotakis (2006) [21] Balasubramanian (2013) [5] Qiu (2015) [63] Qiu (2016) [64] Figure 1.1 Chronology of shock focusing research in gases. Blue lines: shock fronts followed by constant properties, Red lines: shock fronts followed by exponentially decaying properties, Dashed lines: 3D, Solid lines: 2D. Only first author is shown. one is a step function shown as a red line. At time t 0 , pressure jumps from ambient pressure p 0 to p s + , then remains constant. However, pressure for the latter one which is created by an explosion first jumps to a peak value p s + then decreases to p s − , which is smaller than ambient pressure p 0 . Due to the initial charge and the distance to the explosion center, p s − will return to p 0 either gradually which is shown as dotted line in Figure 1.2 or by a secondary shock. The time duration when the pressure behind the shock is over p 0 is called the positive phase (from t 0 to t + ), while the negative phase is fromt + tot − ,meaningthedurationwhenpressureislowerthanp 0 . Sincetheshockfront 3 expands in three dimensions (spherical) or two dimensions (cylindrical), the intensity of theshock wave andtheflow propertiesbehinditdecay asa functionofpropagationtime and distance. Therefore, coalescence and focusing mechanisms of shock fronts initiated by sudden energy-releasing in open area will result in different conditions at the focal region than those where the shock front is followed by constant flow properties. The present work will focus on the second type of waves that are more generally applicable to everyday life. Figure 1.2 A typical pressure-time curve for a shock wave. Red line: shock wave generated by shock tube; Black line: shock wave generated by an explosion. This figure is reproduced from [58]. 1.1.2 Previous Studies on Shock Waves Followed by an Expo- nential Decay One of the earliest work to model the spherical shock waves was done by Taylor (1950) [85]. Asimilarity law wasderived when theshock frontremains strongwith theassump- tion that energy is released instantaneously in an infinitely small region (point source) and perfect gas condition remains all the time. Later, Lin (1954) [78] modified Taylor’s work to reproduce a cylindrical shock wave. Other contributions on theoretical study on isolated spherical or cylindrical shock wave propagation include those of Sedov (1946) 4 [77], Sakurai (1953, 1954) [69, 70], Goldstine (1955) [27], von Neumann (1958) [10], Os- hima (1960) [59], Bach and Lee (1970) [4]. Among those, Taylor [85], Sedov [77] and von Neumann’s self similar solutions [10] are valid at the early time regime when the shock front remains strong, while Sakurai [69, 70] and Oshima [59] conducted a theoret- ical study of shock waves at intermediate times when the shock strength is moderate. In addition, an exact solution was numerically studied by Goldstine and von Neuman [27]. Later, analytical methods, which can describe the entire propagation regime of an expanding shock wave accurately, were derived independently by Sakurai [71] and Bach [4]. Although shock focusing has been studied for decades, coalescence of multiple shock frontsproducedbypointenergysourcesisnotaswellunderstoodasthoseofshockswith constant flow properties behind. There exist only a few theoretical and experimental studies on the coalescence and interaction of multiple cylindrical or spherical shock waves. Thewavepropagationhistoryatapointinspacewhereshockwavesfrommultiple explosions coincidewasapproximatedwithnonlinearadditionruleslabeledLAMB(Low Altitude Multiple Burst) rules [36, 55]. This method can also be used to predict the triple point location if a Mach stem is generated. Despite the fact that LAMB model has a good agreement when compared to the shock reflections on high-explosive field tests and nuclear tests, Brode (1977) [16] stated that the arbitrary assumptions in the LAMBmodelarenotphysicallycorrect. Inhisstudy,fivesimplifiedapproximationswere developedandcomparedtotheLAMBrules,toconstructupperandlowerboundsonthe over-pressures at a point due to two simultaneous spherical shock waves. Comparisons showedthattheLAMBmodelmayunder-predictthepeakpressuresforinteractingshock waves. Computations of coalescence patterns for shock waves initiated by sequentially detonatingdistributedammunitionstackswereperformedbyStarkenberg andBenjamin [81]. Results showed thatredistributing theammunition stacks intosmaller subdivisions was advantageous to reduce peak overpressures. Implosion waves from eight symmetrically placed and simultaneously detonated ex- plosives were produced by the US Navy and resulted in glowing star-like gas regions where shocks fronts collided [26]. In 1984, equilateral triangular patterns of bare ex- plosives and their effect on overpressure was investigated, and results were compared to 5 those obtained by LAMB rules [44]. The tests showed that the overpressure from multi- ple detonations was increased. Furthermore, this technique could cause greater damage at distances further away from the target than compared to one charge with the same total weight. In [94], the effect on aircraft subjected to multiple simultaneous nuclear bursts arranged in hexagonal close-packed patterns, all at sea-level, was investigated. It was concluded that thermal effects contributed to damage more than gust winds or overpressure. 1.1.3 Shock-Shock Interaction When cylindrical or spherical shock waves with equal strength interact with each other, high pressures occur at the contact region, which can be several times greater than their original pressure. A schematic description of two interacting cylindrical shock waves is shown in Figure 1.3(a). Figure 1.3(b) shows an example of the pressure field from a computation of two interacting cylindrical shock waves, using the techniques outlined in Section 2.1. The pressure along the centerline between the two explosives, also shown by the white line in Figure 1.3(b), is plotted in Figure 1.3(c). This exemplifies that the nonlinear interaction of the coalescent shock waves gives rise to a non-uniform pressure inthe regionbehind thereflected waves, which suggests thatinteraction ofweaker shock waves may result in a higher pressure region than a single strong shock wave. Although the effect is desired and impressive, there are several limitations. Such high-pressure area will not expand far beyond the initial contact point. Besides, if the two shock waves are not of equal strength, the resulting pressure in the contact region will not be as high. Therefore, to better take advantage of this high-pressure region, a deeper understanding of the flow properties and how they change with the expansion of the region is necessary. Based on the concept of image bursts [55], the iteration of two shock waves with equal strength can be treated as a single shock wave reflects off a smooth, flat, perfectly reflecting surface. A vast majority of studies have been performed in an analogous situation, in which a planar shock impinges a wedge surface. Although planar shock waves, followed by constant flow properties, are different from shock waves followed by 6 (a) 34.98 0.86 28.16 21.33 14.51 7.68 (b) 0 0.2 0.4 0.6 0.8 1 2 4 6 8 10 Normalized distance Non-dimensional pressure (c) Figure 1.3 (a) A cartoon of shock interaction from two point sources. (b) Pressure distribution from current simulations during shock interaction from two point sources. The gray scale bar shows non-dimensional pressure. (c) Normalized pressure along a line connecting the centers of the explosions in (b). 7 exponentially decaying properties, the interaction patterns on the surface between the two types of shock waves are the same. Two kinds of reflections exist, regular reflection (RR) and irregular reflection (IR), see Figure 1.4. When a planar shock wave travels overasolidwedge, theflowbehindtheshockisforcedtopropagateparalleltothewedge surface, and a reflected shock is generated. This is usually called regular reflection. As the angle between the wedge surface and the incident shock wave increases, to keep the flow direction parallel to the wedge surface, the angle between the reflected shock and the surface will also increase. However, when the angle between the wedge surface and the incident shock wave is sufficiently large that a single reflected shock is not able to turn the flow direction parallel to the wedge surface, a third shock wave is developed to ensuretheflowdirectionisalongthewedgesurface. Thisiscalledanirregularreflection, and the third shock wave is called Mach stem. The point which connects all three shock waves is called triple point. Due to the constant flow properties behind the incident shock, the flow is locally steady if the triple point (or reflection point in the RR case) is treated as a reference. Thus, this configuration is usually referred as a pseudosteady flow. The first researcher who reported the shock reflection phenomena was Ernst Mach in 1878 [14]. Surprisingly, not until around 60 years later, this topic started to capture otherscientists’ interest. In1943,von Neumann derived thefundamental equations that describe regular and Mach reflection [87]. (Noted that Mach reflection (MR) is only one of four subdomains divided from IR. Details are summarized in Ben-Dor’s shock wave reflection book [13].) Since then, intensive investigations on reflection phenomena were carried out. Numerous studies have been engaged in solving the fundamental problem regarding shock reflection, which is the transition condition from RR to IR. There exist many criteria to determine the transition rule from RR to IR. The most common ones are detachment criterion, mechanical-equilibrium criterion and sonic cri- terion. In1943,vonNeumann reportedthatinordertochangetheflowdirectionbehind theincidentshocktobeequaltothewedgeangle, areflectedshockisrequired. However, there is a maximum turning angle for the reflected shock. Beyond this angle, the reflec- tion point, which connects the reflection shock and incident shock, can no longer attach to the wedge surface. This condition is later known as the detachment condition. Later in 1975, Henderson and Lozzi revisited detachment condition and argued that there is 8 (RR) Incident shock Reected shock (a) (IR) Mach stem Incident shock Reected shock (b) Figure 1.4Shock wave propagates over a solid wedge: (a) Regular reflection (b) Mach reflection. a sudden pressure change at the transition from RR to IR if detachment condition is applied, which means a compression wave or an expansion wave is required to support this sudden pressure change [32]. However, in the experiments, neither of the waves were detected. Therefore, a mechanical equilibrium condition was proposed that the transition from RR to IR should be continuous. This condition has a good agreement with the steady flow experiments by Hornung and Kychaakoff [37]. The sonic criterion wasalsofirstintroducedbyvonNeumannin1943[87]. Hestatedthatinordertotransit from RR to IR, the corner-generated signals need to catch up with the reflection point. In the pseudosteady flow, as long as the flow behind the reflected shock is supersonic if viewed intheframeofthereflection point, thecorner-generatedsignalscannotreach the reflection point and IR is not possible. Noting thatthe sonic and detachment conditions are very close to each other in a practical sense, however, based on Lock and Dewey’s experiments [52],inpseudosteady flows, thetransitionfromRRtoIRoccursatthesonic condition rather than the detachment condition. In 1979, Hornung et al. introduced a new criterion called the length-scale criterion [38]. The authors stated that to develop a shock wave with a finite length (i.e. Mach stem), pressure signals must be able to reach the reflection point. Under their statement, there exist two different transition condi- tions depending on the flow conditions. In steady flows, this statement leads to the mechanical equilibrium criterion, while in unsteady flows, it meets the sonic criterion. Based on the previous experimental results in both steady and unsteady flows, the me- chanical equilibrium criterion is more suitable for steady flows, while in unsteady flows, results have a better agreement with the detachment criterion or the sonic criterion. 9 Thus, this finite length scale criterion seems to be more acceptable than other criteria. However, until recently, the general criterion is still under debate since the agreement between the proposed theories and the experimental data cannot cover up the entire span of incident shock Mach numbers and wedge angles. Compared to shock reflection in pseudosteady flows, the interaction between two cylindrical or spherical shock waves with equal strength can result in the same reflection patterns, see Figure1.5. Aregularreflection occurswhenthese two shockfrontsinitially interact with each other. As they propagate further, at certain conditions, the regular reflection may transit into an irregular reflection depending on the strength of the shock waves and their geometry. Higher pressure and density are formed behind the Mach stem, which indicates that more energy is focused on this direction. This phenomenon hasmotivatedthepresentinvestigationtogaininsightsofdirectingtheenergyofmultiple explosions towards a target, e.g. maximizing pressure, while simultaneously minimizing the collateral damage of the surrounding area. Unlike the pseudosteady flows, the flow behind the cylindrical or spherical shock wave cannot be treated as locally steady flow. Instead, it changes as the wave front expands. Thus, the transition criterion that is suitable to pseudosteady flows may not be applicable for unsteady flows. More details are presented in Chapter 3. 1.2 Outline and Contributions The outline of this thesis is presented as follows: In Chapter two, simulations of multiple energy sources are performed, and the re- sulting shock interaction and coalescence behavior are explored. A number of energy sources are placed around the target, and conditions at the target area are monitored and compared to those obtained using a single source. For each simulation, the en- ergy summed over all sources is kept constant while the radial distances between target and energy sources, and the source initiation times, are varied. Each energy source is modeled as a point source explosion. The resulting shock wave propagation and shock front coalescence are then solved using the inviscid Euler equations of gas dynamics on overlapping grids employing a second-order Godunov scheme. 10 (a) (b) (c) (d) (e) (f) (g) (h) Figure 1.5 Schlieren visualizations of two cylindrical shock waves interaction. (b) and (c) show the regular reflection pattern, while (d)-(h) show the irregular reflection pattern. In Chapter three, based on the results from Chapter two, the interaction of multiple simultaneous cylindrical shock waves is divided into two phases, expanding phase and focusing phase. Phase transition condition is analyzed and compared to several existing RR to IR transition criteria. The resulting new transition curve based on the simulation results, together with an analytical cylindrical shock wave solution, is then applied to search the optimal number of energy sources placed around a target while keeping the total energy to initialize the energy sources as constant. In Chapter four, the theory of geometrical shock dynamics (GSD) is explored as an alternative to a full Euler solution for this particular project. Details on numerical implementation and parallelization for 2D computations using a spatial decomposition method coupled with a front tracking algorithm is presented. A symmetric boundary condition to further speed up the code is given. In order to increase the accuracy of this approximate theory, similitude equations obtained from analytical solutions or experimental data are coupled to the numerical code. Results from the modified code 11 are compared with existing analytical expression and experimental data. The potential of using small energy sources in an optimized pattern to direct and focusdesiredshockwavestowardatargetwhilereducingcollateraldamageistheprimary interest of the current study. There only exists a few studies in the open literature on thecoalescence ofmultiple shockwaves. Mostrecently, Stowe developed areduced order modelforoptimizationonarectangulararrayofsphericalexplosives[80]. However, none of these works focused on the effect of the geometry of the initial placement of energy sources. In addition, when dealing with the nonlinear interactions between shock waves, their analyses are either insufficient or physically incorrect. Thus, to extend previous understanding and further investigate focusing effects generated by the coalescence of multiple shock waves, the following contributions are made: 1. Represent a comprehensive study on the interaction of multiple shock waves, in- cluding the initial placement, the timing and the number of energy sources. 2. Investigate the transition condition from RR to IR on the interaction of two cylin- drical shock waves with equal strength. 3. Explore the optimal number of shock waves placed around a target while keeping the energy summed over all sources as constant. 4. Extend the computational advantage of the theory of geometrical shock dynamics (GSD) by using a spatial decomposition method coupled with a front tracking approach. Also, symmetric boundary conditions are developed to reduce the com- putational cost further. 5. Modify the GSD theory to increase its accuracy to simulate expanding waves. 6. Developacomputationallyfastmodelthatcanpredictandoptimizetheconditions for both the target area and the collateral area. 12 Chapter 2 Interaction and Coalescence of Multiple Simultaneous and Non-simultaneous Cylindrical Shock Waves There exist only a few studies that reveal the potential of generating higher overpres- sure from the coalescence of multiple simultaneous shock waves [26, 44]. To further understand the physical mechanisms behind this phenomenon, a comprehensive nu- merical investigation on interaction and coalescence of multiple simultaneous and non- simultaneous cylindrical shock waves is presented. This chapter describes a numerical approach that is used in the current study to assess the results of shock focusing from multiple cylindrical shock waves. A specific problem consisting of converging shock fronts emerging from cylindrical shock waves, which are arranged in both circular pat- tern and semicircular pattern around a target, is set up to help understand shock front interaction during the focusing process [63]. 2.1 Numerical Scheme and Initial Conditions TheinviscidandadiabaticflowsunderconsiderationaregovernedbytheEulerequations of gas dynamics, i.e. conservation of mass, momentum and energy, ∂ρ ∂t +∇·(ρu) = 0, (2.1) ∂(ρu) ∂t +∇·(ρuu)+∇p =0, (2.2) ∂(ρE) ∂t +∇·(u(ρE +p)) = 0, (2.3) E =e+ 1 2 (u 2 +v 2 +w 2 ), (2.4) 13 where ρ,u,p are the density, velocity vector and pressure. In Equation (2.4), E and e are the total energy and the internal energy per unit mass, and u,v,w are the velocity components in three directions. The equation of state (EOS) for a perfect gas is used to close the system of equations. The EOS is given by e = p ρ(γ−1) . (2.5) The variable e is the internal energy, p is the pressure, ρ is the density; and the ratio of specific heats γ = 1.4 with the assumption that the ambient air is an ideal gas, more details can be found in Anderson [2]. The simulations are performed using Overture, a research tool to solve partial dif- ferential equations using finite difference methods on overlapping grids [17]. The Euler equations of gas dynamics with a shock-capturing second-order Godunov scheme are utilized to compute the shock wave dynamics [19]. Parallel computing based onmessage passing interface (MPI) and adaptive mesh refinement (AMR) are employed to improve the efficiency of the computation [34]. 2.1.1 Initial Conditions Theinitialconditionsforthepointsourceexplosionin3DarebasedonTaylor’ssimilarity law [85]. For 2D computations, the 3D initial conditions were modified to reproduce a cylindrical shock wave [78]. Taylor’s similarity law for pressure, p, density, ρ, and radial velocity, u, can be written as p p 0 =R 0 −3 F(η), (2.6) ρ ρ 0 =ψ(η), (2.7) u =R 0 −3 2 Φ(η). (2.8) Here,p 0 andρ 0 represent the ambient pressure and density ahead of the shock wave. R 0 is the chosen radius of the shock wave front at time zero,η = r R 0 withr representing the 14 radialcoordinatemeasuredfromtheexplosioncenter, andF,ψ andΦarefunctionsofη. Applying the similarity law to the equations of motion, continuity, and equation of state for a perfect gas leads to the following three differential equations in non-dimensional form, ˙ φ(η−ψ) = 1 γ ˙ f ψ − 3 2 φ, (2.9) ˙ ψ ψ = ˙ ψ+( 2φ η ) η−φ , (2.10) 3f +η ˙ f + γ ˙ ψ ψ f(φ−η)−φ ˙ f =0, (2.11) wheref = Fa 0 3 A 2 and Φ =Aφ. The sound speed, a 0 , is the value in ambient air, andA is a coefficient that can be determined from the total energy and the shock front radius. The total energy, E, can be separated into kinetic energy and heat energy as follows: Kinetic energy = 4π Z R 0 0 1 2 ρu 2 r 2 dr, (2.12) Heat energy = 4π Z R 0 0 pr 2 γ−1 dr. (2.13) Expressing equations (2.12) and (2.13) in terms of the variablesf,φ,ψ andη, the total energy E can be written as: E = 4πA 2 ρ 0 2 Z 1 0 ψφ 2 η 2 dη+ p 0 a 0 2 (γ−1) Z 1 0 fη 2 dη . (2.14) In the early stages of the explosion, the shock wave is strong, and one can assume that p s >> p 0 . Consequently, the boundary conditions at η = 1 can be obtained from 15 0 0.2 0.4 0.6 0.8 1 η Pressure, density and radial particle velocity Pressure Density Velocity 0 0.2 0.4 0.6 0.8 1 Figure 2.1 Normalized initial conditions based on Taylor’s similarity law [85]. Figure from [63] . the Rankine-Hugoniot relations ρ s ρ 0 ∼ = γ +1 γ−1 , (2.15) U s 2 a 0 2 ∼ = (γ +1) 2γ p s p 0 , (2.16) u s U s ∼ = 2 γ +1 , (2.17) where the subscript s represents the state behind the shock and U s is the shock wave speed. Given the total initial energy and shock radius, the three differential equations (2.9)-(2.11) can be solved numerically to obtain the initial conditions. These initial conditions result in a sphere of a given radius R 0 with a moderate shock jump suitable for implementation in numerical simulations. The two advantages of computing initial conditions compared to simulating a condensed energy source directly are: (1) First, it eliminates the need to generate an extremely fine mesh for the energy source. (2) Also, sharp discontinuities at the wave front can be avoided. Figure 2.1 shows line plots of non-dimensional initial conditions for pressure, density, and radial velocity. The values are made non-dimensional by scaling with the peak values at the shock wave front. 16 0 5 10 15 20 25 0 0.5 1 1.5 Time [µs] Pressure [MPa] Result from Jiang (1998) Current simulation Figure 2.2 Comparison of pressure history 5 mm away from the shock center for our results and those of [42]. Figure from [63] . 2.1.2 Verification of Initial Conditions In order to verify our initial conditions, a simulation of a single spherical shock wave is performed to compare against the work presented by Jiang et al. [42] who conducted both numerical and experimental studies of a micro-blast wave. The initial conditions (total energy 1.38 J and initial spherical shock radius 1.5 mm) are set to be the same as in Jiang et al., but the grid resolution is slightly different. In Jiang et al., the grid size was set to Δx×Δy×Δz = 0.05mm×0.05mm×0.035mm. In this study, the original gridsize is kept uniformwith 0.125mm×0.125mm×0.125mm. However, AMR is utilized with a grid refinement factor of four; thus the smallest grid size is 0.031 mm. Figure 2.2 showsacomparisonoftheresultingpressureatapoint5mmawayfromtheshockcenter for the results from [42] and the current model. 2.1.3 Convergence Rate of Numerical Solutions To choose a proper grid size, a Richardson extrapolation method is used to estimate the convergence rates for a 2D problem. A single source with initial energy 5000 J and 17 Table 2.1 Estimated convergence rates using second-order Godunov scheme. Components P p P u P v P ρ ¯ t = 19.2 0.80 1.16 1.16 1.01 ¯ t = 23.2 0.88 1.09 1.09 1.00 shock radius 1.5 mm is added to the center of a square domain. The computation is repeated for three different grid sizes (0.1 mm, 0.05 mm, and 0.025 mm). Following the algorithm presented by Roy [67], the convergence rate,P, is computed via the following expression, ku h 1 (x)−u h 2 (x)k ku h 2 (x)−u h 3 (x)k = |h 1 P −h 2 P | |h 2 P −h 3 P | . (2.18) Here, u h (x) represents the numerical solution, h 1 , h 2 , and h 3 denote the different mesh sizes, k·k indicates a discrete approximation to the L 1 norm, which is most often used for hyperbolic equations with discontinuities [6]. The results of the convergence rates at two arbitrary non-dimensional time instants, ¯ t = 19.2 and ¯ t = 23.2, for pressure, p, horizontal and vertical velocity components, u and v, and density, ρ, are presented in Table 2.1. Here, time is non-dimensional according to equation (2.19). The computed convergence rates at both time instants for velocity and density are higher than the expected asymptotic values of 1 for the L 1 norm. The reason may be that the overall error at those resolutions is not entirely dominated by the error associated with dis- continuities in the flow (see [34]). Based on these test results, a uniform grid size of 0.1 mm is used in the simulations presented in the following sections. The time step in the numerical computations in this study is based on the Courant-Friedrichs-Lewy (CFL) number, which is set to 0.75 for all the simulations. 18 2.2 Problem Setup To study shock-shock interaction and shock focusing effects, a problem is set up as follows: the total energy summed over all sources is kept constant at 5000 J and the initial shock wave radius for each source is set to 1.5 mm, with the number of shock wavesvariesindifferentsimulations. Computationsareperformedona120mm×120mm domain, see Figure 2.3(a). Two locations in the computational domain are monitored: (1) the target (T1) - a region where the shock waves are focused and extreme conditions are formed; and (2) a point (S1) to evaluate the collateral damage. The point S1 is located atthe same radialdistance away fromthe center source as targetT1, but placed on the outside of the concentric source arrangement, see Figure 2.3(a). Previous experimental results on shock focusing have shown that the generation of a converging shock with a shape close to that of a circular cylinder is advantageous for generating extreme conditions at the focal region [5, 22, 24, 45, 46, 60, 82, 84]. Based on these studies, the sources in the current study are placed at equal radial distances of r = 30 mm from target T1. The coalescent shock front approaches a circular shape when additional sources are used in the simulations. 2.3 Results and Discussions The purpose of the simulations is to investigate how source arrangement and timing of individual sources affect the conditions in the surrounding areas. Three types of simulations are performed using Euler equations: (1) all energy sources are placed in a concentric pattern around the target with simultaneous initiation, (2) all sources are placed in a concentric pattern around the targetwith simultaneous initiation except one source that is delayed, and (3) all sources are placed in a half-circular pattern around the target with simultaneous initiation. The non-dimensional properties for density ¯ ρ, pressure ¯ p, velocity ¯ u from the Euler 19 simulations, are obtained using the following relations ¯ ρ = ρ ρ ∗ ;¯ p= p p ∗ ; ¯ u = u a 0 ; ¯ t= t t ∗ ; (2.19) where ρ ∗ = 1 kg/m 3 , p ∗ = ρ ∗ a 0 2 , a 0 = √ γRT 0 , R = 287.06 J/kg·K, γ = 1.4, T 0 = 293.15 K,t ∗ =0.728μs. Here,t ∗ is obtained by taking a unit length,L = 0.001 m, divided by the speed of sound, a 0 , times a factor of 0.25, which is equal to the delay time unit Δt used in later sections, t ∗ = L a 0 Δt. 2.3.1 Case 1 - Simultaneous Initialization Four different simulations are performed featuring 1, 3, 5, and 10 sources initialized simultaneously. (Simulations are referred to as 1C, 3C, 5C and 10C.) Figure 2.3 shows numerical schlieren visualizations of all cases at an early time instant and a later time instant just before the converging shock fronts reach the focal region. For each case, probes are placed at the two points of interest, T1 and S1, to collect pressure data as a function of time. Table 2.2 shows the non-dimensional peak pressure, ¯ p, and time, ¯ t, at the time instant when the maximum peak pressure at T1 and S1 were obtained. Forthecasewithasinglesource,itshouldbenotedthatthepeakpressureandarrival time for T1 and S1 differ slightly because of the grid size and the fact that perfect grid symmetry is not obtained. The results from the cases with 3 and 5 sources show that peak pressure at S1 is achieved when the shock front generated by the nearest source (the source at the center of the computational domain) reaches S1. Peak pressures (Table 2.2) show that the case with 5 sources resulted in lower peak pressures at S1 when compared to the cases of a single and 3 sources. This is because each source has a lower energy than in the case of a single or three sources. In the case of 3 sources, the converging shock front remains a triangular shape throughout the focusing process. There is a 1.3% difference in its arrival time at T1 and S1, indicating that the Mach number of the converging shock front does not increase. This is in agreement with previous studies [5, 9]. 20 120 mm S1 T1 120 mm Blast wave (a) (b) (c) (d) (e) (f) (g) (h) Figure 2.3 Schlieren visualizations of the cases with (a) one (1C), (b) three (3C), (c) five (5C) and (d) ten (10C) point explosions in their early stages (top row) and later stages (bottom row). Figure from [63]. 21 Table 2.2 Comparison of non-dimensional peak pressure, ¯ p, and arrival time, ¯ t, at T1 and S1 for cases with 1, 3, 5 and 10 sources. The total energy is kept constant for all cases. T1 S1 Case ¯ p ¯ t ¯ p ¯ t 1C 11.43 17.8 11.43 17.9 3C 36.85 29.3 4.87 28.9 5C 86.95 34.2 3.50 35.3 10C 241.8 35.5 3.21 46.0 As the number of sources is increased to 10, the peak pressure at T1 is 178% higher compared to the case with 5 sources. The reason is that the individual shock fronts coalesce and a reconfiguration process starts as the shock front propagates towards the center, see Figure 2.4. During each reconfiguration process the vertices of the polygonal shock front form Mach shocks that propagate with a higher Mach number than the original shock fronts, and thus overtake them. This confirms the results reported in pre- vious studies [23, 24, 72, 76]. This process repeats itself numerous times until the shock reaches the center. For the case with 10 sources, the Mach number of the convergent shock front is higher than that of 3 or 5 sources, resulting in a higher peak pressure at the focal point. For the cases with 5 and 10 sources, the time when the highest peak pressure isrecorded atT1 isearlierthanthatofS1,indicating thattheconverging shock accelerates during the focusing process. 2.3.2 Case 2 - Delayed Source Initialization A Gaussian function is implemented as a source term and added to the energy equation to control the delayed initiation time ofsources. Equation (2.20)shows the general form of a forcing function, f(x)=g(x)e −α(t−Δt) 2 (2.20) where g(x) represents the energy term in the initial conditions. The width of the pulse is controlled by α, and the Δt determines the time to initialize the shock wave. In this study, α = 5·10 13 . For the cases with 3, 5 and 10 total sources, a single source is retarded by one, two 22 (a) (b) (c) (d) Figure 2.4 Schlieren visualizations of ten sources reconfiguration. a) Initial decagonal shock front; b)-c) Reconfiguration process; d) Reconfigured decagonal shock front. 23 or three delay time units, represented by Δt in equation (2.20). Here, one delay time unit Δt is chosen to be the same for all cases, which is equal to 0.25 non-dimensional time units, as mentioned in section 2.2. For the different cases, one delay time unit means 3.4% of the total time for the shock wave to reach the target for the case with 3 sources, 2.9% for 5 sources case and 2.8% for 10 sources case, since the total time for shock wave to arrive at the target is different. In each case, the source located at the center of the computational domain is delayed. Figure 2.5 shows the numerical schlieren visualizations for the case with three sources, of which one is delayed by three delay time units. Figure 2.6 shows the numerical schlieren visualizations for the case with 5 sources, with the center source delayed by three delay time units. Figures2.5and2.6showthatthepeakpressureatthepointS1iscausedbytheshock frontgeneratedby thesource placedatthecenter ofthecomputationaldomain. Around target T1, a non-symmetrical geometry of the converging shock is created. A summary ofpeak pressures andthetimeinstants atwhich they occurredis shown inTable2.3. As expected,forcaseswith3sourcesand5sources,thepeakpressureamplitudeattargetS1 remains constant but the time instant when it occurs is increased by the corresponding delay time units. However, for the case with 10 sources, the peak pressure at S1 is caused by the interaction of the three nearest sources. Therefore, the increased time at which the peak pressure occurs is smaller than the prescribed delay time units. For the case with 3 sources with one delay time unit, the peak pressure at target T1 shows a modest decrease of 5.6% compared to the simultaneous case. When the delay- timeisincreasedfurther,thepeakpressureexperiencesalargerdecreaseof21.5%fortwo delaytimeunitsand33.2%forthreedelaytimeunits. Asimilartrendisobtainedforthe peakpressure forthecasewith5sourcesattargetT1. Thepeakpressure decreases 4.6% for one delay time unit, 29.4% for a two delay time units, and 44.3% for a three delay time units, see Table 2.3. The time instants when the peak pressures are obtained at T1 occur earlier than the ones at S1, which indicates that the convergent pentagon-shaped shock front still accelerates during the focusing process, despite a lack of symmetry. For the case with ten sources, delaying a single source by 1Δt,2Δt, and 3Δt results in 5.9%, 21.1% and 34% reduction of the peak pressure amplitude at target T1, respectively, as compared to the case with no delay, which is similar to the case with 3 and 5 sources. 24 S1 T1 (a) (b) (c) (d) Figure 2.5 Schlieren visualizations from the case with 3 sources with one source delayed by three delay time units. From (a) to (d), four sequential time instants are plotted and the source in the center of the computational domain is delayed. Figure from [63]. 25 S1 T1 (a) (b) (c) (d) Figure 2.6 Schlieren visualizations from the case with 5 sources with one source delayed by three delay time units. From (a) to (d), four sequential time instants are plotted and the source in the center of the computational domain is delayed. Figure from [63]. 26 Table 2.3 Comparison of non-dimensional peak pressure, ¯ p, and arrival time, ¯ t, at T1 and S1 for cases with 3, 5 and 10 initial sources where a single source has been delayed up to 3Δt. The total energy is kept constant (5000 J) as in the previous cases. 3 sources T1 S1 Delay time ¯ p ¯ t ¯ p ¯ t Simultaneous 36.85 29.3 4.87 28.9 1Δt 34.78 29.7 4.83 30.1 2Δt 28.91 31.2 4.83 31.1 3Δt 24.61 32.6 4.83 32.1 5 sources T1 S1 Delay time ¯ p ¯ t ¯ p ¯ t Simultaneous 86.95 34.2 3.5 35.3 1Δt 82.99 34.4 3.47 36.6 2Δt 61.41 35.2 3.47 37.6 3Δt 48.45 36.3 3.47 38.6 10 sources T1 S1 Delay time ¯ p ¯ t ¯ p ¯ t Simultaneous 241.8 35.5 3.21 46 1Δt 227.6 35.6 3.27 46 2Δt 190.8 35.9 3.32 46.4 3Δt 158.3 36.3 3.35 46.8 27 Table 2.4 Comparison of non-dimensional peak pressure, ¯ p, and arrival time, ¯ t, at target T1 and points S1–S3 for the case with 10 sources placed along a semicircle (10SC) compared with 10 sources in a circular pattern (10C). The total energy is kept constant (5000 J) as in the previous cases. T1 S1 S2 S3 Case ¯ p ¯ t ¯ p ¯ t ¯ p ¯ t ¯ p ¯ t 10C 241.8 35.5 3.21 46 3.21 46 3.21 46 10SC 42.8 26.5 4.67 34.8 7.87 58.5 3.84 38.8 2.3.3 Case 3 - Semicircle Placement Forthe case of10 sources placed in a semicircle (10SC),two additionalpoints ofinterest labeled S2 and S3 are monitored. These new points are located 30 mm to the right of T1 and 60 mm above T1, respectively, as shown in Figure 2.7. Figure 2.7 shows the numerical schlieren visualizations for four different time instants. The conditions at the points of interest are summarized in Table 2.4. The peak pressure at S1 is higher than for the case of 10 sources placed in a full circle (10C). Also, the arrival time is smaller since the total energy at the left side of the half circle is higher. The combined divergent shock front has a higher Mach number due to multiple shock interactions. On the right side, for target T1, the peak pressure is lower than for the 10C case due to the asymmetry of the incident shock fronts. The shock speed on the right side is higher because the reconfiguration process occurs earlier and more often as shock fronts are generated from sources located closer to each other. This results in an initially stronger convergent shock in the shape of a semicircle. Point S2 is on the path of the convergent shock front resulting in a larger peak pressure than at S1 and S3. The lowest peak pressure is at S3 because it is initially influenced only by a shock front combined from the three nearest sources, which are placed asymmetrically with respect toS3. At latertimes, shock waves fromthe othersources reach S3 but they areattenuated andweaker thantheinitialshock. Thepeakpressure atS1ishigherthan that at S3 because of source-symmetric placement and because more sources contribute to the initial shock front. 28 T1 S1 S2 S3 (a) (b) (c) (d) Figure 2.7 Schlieren visualizations of shock front propagating from 10 sources set in semicircle (10SC). Target T1 and points S1–S3 are labelled in (a). The shocks coalesce, and in (b)-(d) the shocks first reach T1 and then S1, S3 and later S2. Figure from [63]. 29 − 5 − 4 − 3 − 2 − 1 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 Horizontal position around T1 [mm] Normalized Pressure 10C 5C 3C 10C 5C 3C Figure 2.8 Pressure plot along horizontal line centered at target T1 for 3, 5, and 10 sources with one source delayed 3Δt (solid lines) and no delay (dashed lines). T1 is located at zero. Figure from [63]. 2.3.4 Center of Energy A center of energy is defined as a unique point where the entire energy is assumed to be concentrated in a similar way as the center of mass is defined. For the 3C, 5C and 10C simulations performed in case 1 with simultaneous initiation, the center of energy is located exactly at target T1. For the 3C, 5C and 10C simulations presented in case 2, when one of the sources is delayed, the center of energy is shifted to the left of target T1. This results in a lower peak pressure at T1. This is shown in Figure 2.8, where the pressure along a 10-mm horizontal direction centered at T1 is plotted for no delay and with a 3Δt delay for cases 3C, 5C and 10C. All peak pressures have been normalized to unity. It can be seen that from 10C to 3C, the center of energy, which is the focal point of the convergent shocks, is shifted further to the left from target T1. This is because the energy of each source in the 3C case is larger than that of each source in the 10C case. For the 10SC simulation with sources placed in a semicircle presented in case 3, the center of energy is shifted to the left of target T1. The resulting peak pressure at T1 is reduced. Given the distance between targets and sources, the total energy, and the number of 30 Figure 2.9 Pressure history at point S1 for 5C, 10C and 10SC. Figure from [63]. sources, the closer the center of energy of all sources is to the target, the better focusing of the convergent shock will be achieved resulting in higher pressures and temperatures. This principle does not apply if a point of interest is initially influenced by a subset of nearest sources (the initial shock front is the strongest even if additional shocks follow at later times). For example, see the pressure history at point S1 for the cases with 5C, 10C and 10SC in Figure 2.9. 2.3.5 Energy Dissipation and Shock Reconfiguration There is an infinite number of ways to arrange multiple sources with a fixed position of the center of energy. Here, the effects of the converged shock front reconfiguration have to be taken into account. Forexample, in the cases with 5 and 10 sources, itcan be seen thata regularpolygonal-shapedshock wave keeps reformingatsuccessive intervals, with the Mach number increasing at each successive interval. This phenomenon is explained byWhitham’sgeometricalshockdynamictheory[92],andhasbeenobserved inprevious numerical and experimental studies [24, 72, 82]. When two shock fronts intersect, a reflection is created. The reflected wave then travels into a region with higher density, pressure and temperature created by the former shocks, see Figure 1.3. Under certain conditions, the reflected waves generate an irregular reflection in which a Mach stem 31 25 26 27 28 29 30 31 32 33 34 35 150 200 250 300 350 Radius [mm] Non-dimensional pressure Figure 2.10 Comparison of non-dimensional peak pressure, ¯ p, at target T1 with 10 initial sources placed along a circular pattern at 5 different radii. The total energy is kept constant (5000 J). Figure from [63]. is created. The Mach stem overtakes and merges with the initial shock front. After each successive reconfiguration, in which a Mach stem reflects off another Mach stem, the Mach number is increased. This process keeps repeating until the waves reach the center. For higher numbers of sources, a larger number of shock intersections occur. Thus, as the distance between the source and target T1 is increased, the converging polygonal shock has additional time to reform and accelerate more often. However, when the energy is kept constant, increasing the shock wave propagation distance from the source to the point of interest means that the shock front dissipation increases. To study which one of these competing effects is dominant, five additional simulations were performed. These simulations featured 10 sources with the same total energy placed at varying radial distances from T1. Results are shown in Figure 2.10. It can be seen that as the radius increases, the peak pressure at target T1 decreases (at least in the range from 25 to 35 mm radii), suggesting that energy dissipation plays a more important role than shock reconfiguration. 32 2.4 Chapter Summary The presented work in this chapter summarizes results from numerical simulations of coalescent shock waves from multiple point sources used to achieve the following objec- tive for this thesis: (1) to increase the extreme conditions at a designated target area by using multiple sources and utilize shock focusing mechanisms from the combined shock fronts, and (2) to substantially reduce collateral damage away from the target area. From the results with 1, 3, and 5 sources, we can see that the peak pressure at target S1 (representing collateral damage) is mainly influenced by the energy of the nearest shock, however the peak pressure at target T1 (at the intended focal region of all coalescent shocks) is influenced by the shock focusing process. In the case of 3 sources (case 3C), the convergent shock front formed a triangular shape that remained triangularduringthefocusingprocess. Forthecasewith5sources(case5C),apentagon is obtained. It only had time to reconfigure once from a pentagon (5-sided polygon) to a decagon (10-sided polygon) during the focusing process. In the case with 10 sources (case10C),thecoalescentshockwavesfirstformadecagon,thenchangeintoanicosagon (20-sided polygon) and back again to a decagon. Subsequently, this process is repeated multipletimesthroughoutthefocusingprocess. Ateachreconfigurationstage, theshock front Mach number is increased. Considering that in reality, it is difficult to initialize all the sources at exactly the sametime, scenarios whereonesourceisdelayed aresimulated. Itisfoundthatforcases with 3, 5, and 10 sources, with one source delayed by Δt, there is not much influence on the peak pressure at target T1. But as the delay time is increased to 3Δt, the pressure recorded at target T1 was reduced dramatically. This is also related to the position of the center of energy. With a longer time delay, the center of energy was shifted further away from target T1. A semicircle source placement case (10SC) is also conducted to compare results with the previous case of 10 sources placed in a circular pattern (10C). It is found that the conditionsofthetargetarealocatedatthefocalpointexperiencedalowerpeakpressure, butthepeakpressures measuredatthreelocationssimulating collateraldamage(S1–S3) show higher and lower values depending on their location. Therefore, it is important to 33 take into account the combined effect of coalescent shock fronts also for regions away from the intended target area. To understand the competing effects of shock reconfiguration during the focusing phaseandattenuationoftheflowvariablesbehindtheshockfront,aseriesofsimulations with different initial radial positions of the sources are undertaken. Results show that for the current geometrical spacing of sources and a total energy of 5000 J, an increase in radius yields lower peak pressures at the focal region. Thus, for these simulations, the energy attenuation is dominant. Despite the fact that comprehensive information can be obtained through entire space and time by using Euler equations, it is computationally expensive. For each case (2D) in the current study, with more than 40 processors used on high-performance computing clusters, it still takes more than 12 hours to complete. For 3D problems, it may take weeks or even months. However, during decision-making processes (e.g. optimize the number of shock waves placed around a target while keeping the total energy summed over all the shocks as constant), due to the highly iterative nature of optimization algorithms, the current numerical method is not applicable. Investigation on the interaction of two identical cylindrical shock waves will be presented to seek an alternative method for optimization purposes. 34 Chapter 3 Interaction of Two Cylindrical Shock Waves with Equal Strength Fromthe previous chapter, itcan be seen thatnon-simultaneous initializationcan result in an unsymmetrical placement, which plays an important role in optimizing the focus- ing effect for the interaction of multiple shock waves. Thus, for optimization purposes, the rest of work will focus on simultaneous initialization. In the case with simultane- ous initialization from Section 2.3.1, four different number of cylindrical waves are used, resultinginamonotonicbehavior. However, asthenumberofshockwaves keepsincreas- ing, will this behavior change and an optimal solution can be obtained? Investigation of this question is presented in this chapter. 3.1 Division on Interaction of Multiple Shock Waves Based on the gradient of shock front Mach number, the interaction of multiple shock waves can be divided into two phases, expansion phase and focusing phase. In the expansion phase, the strength of each shock front is always decreasing due to two facts, the nonlinear interaction with the flow behind it and the expansion in space, see Figure 3.1. As they start to interact with each other, a regular reflection (RR) occurs, and a reflection point is formed to connect neighboring shock fronts. As the expansion phase continues, under certain conditions, based on the geometry and the strength of the interacting shock fronts, an irregular reflection (IR) takes place, resulting in a Mach stem between each shock front. At this moment, focusing phase takes over, and a polygonal-shaped converging shock wave towards the focal point is developed. As the convergingshockwavepropagatesinward,Machstemstakeovertheadjacentportionsof the shock front, which will later result in a reconfigured and enhanced converging shock 35 front. This process keeps repeating and the curvature of the converging shock front decreases for each reconfiguration, eventually resulting in a regular stabilized polygonal converging shock front, see Figure 3.2. (a) (b) Figure 3.1 Schlieren visualizations of expansion phase with six sources. a) Before interaction. b) Regular reflection. (a) (b) Figure 3.2 Schlieren visualizations of focusing phase with six sources. a) Before first reconfig- uration. b) After several reconfigurations. For the focusing phase, the polygonal-shaped converging shock wave is similar to a 36 perturbed polygonal shock wave, which has been numerically studied by Schwendeman [72, 76] and experimentally investigated by Eliasson [22, 23, 24]. Schwendeman utilized an approximate theory called Geometrical shock dynamics to simulate the converging process and same reconfiguration procedure has been reported. This theory is compu- tationally efficient because it only considers the change in shock geometry and neglects the interaction with the flow behind the shock front. For the expansion phase, there exists many analytical studies for a single cylindrical or spherical shock wave, which is included in Section 1.1.2. Thus, if the transition condition can be predicted between two phases, there is hope to describe the whole interaction process without using Euler equations. 3.2 Transition Condition Numerous studies have been conducted on the shock reflection phenomena to seek the transition condition from RR to IR. In 2007, Bendor published a book regarding this phenomena, in which he summarized most of the studies in this area [13]. Based on the book, the comparison between all the existing theories with the experimental data was never satisfactory enough to lead to a general criterion. However, when taking external effects into account (e.g. viscous effects, real gas effects), a general agreement between theories and experimental data can be improved. From previous studies, it is suggested that transition condition in steady flows can generally be predicted by the mechanical equilibrium criterion, while results from experiments inunsteady flows tend to approach to the sonic criterion. In the present study, because the flow behind the shock front is changing with time, the sonic criterion and reflection phenomena in unsteady flows are reviewed briefly. The underlying assumption for the sonic criterion is that to develop a Mach stem, pressure-signals must be able to communicate with the reflection point. In the pseu- dosteady flow case, these pressure signals are generated when the incident shock hits the leading edge of the reflecting wedge (i.e. corner-generated signals). A schematic figure is presented in Figure 3.3 to further explain the termination condition for RR. Viewed in the frame of the reflection point, which is denoted as R in Figure 3.3, the 37 α! β! α# β# IS RS R M! M# M$ RR IS RS R θ M& Figure 3.3 Schematic illustration of the wave configuration in the frame of the reflection point labeled R. flow is locally steady. A supersonic flow with Mach number M 1 propagates along the wedge to the ground. The incident shock (IS) deflects the flow by an angle ofβ 1 and the original angle between the flow and the incident shock is α 1 . The Mach number of the flow between the incident shock and the reflected shock (RS) is represented by M 2 and α 2 denotes the angle between the reflected shock and the flow deflected by the incident shock. Next, the reflected shock turns the flow back to a direction parallel to the wedge, resulting in another deflected angle β 2 and a different Mach number of the flow M 3 . In order to let the corner-generated signals catch up with reflection point R, M 3 must be subsonic. Thus, the transition condition from RR to IR is M 3 = 1. To calculate M 3 , the following oblique shock relations need to be used, cotβ = (γ +1)M i 2 2(sinα 2 −1) tanα, (3.1) 38 and M d = s (γ +1) 2 M i 4 sin 2 α−4(M i 2 sin 2 α−1)(γM i 2 sin 2 α+1) [2γM i 2 sin 2 α−(γ−1)][(γ−1)M i 2 sin 2 α+2] (3.2) whereα represents the angle between the oblique shock and the incident flow,β denotes the flow deflection angle. The Mach numbers for the incident flow and the deflected flow are labeled as M i and M d , respectively. Given the Mach number of incident shock M 0 and the wedge angle θ, M 1 can be calculated as M 1 =M 0 /cosθ and α 1 =π/2−θ. Then, the deflected flow conditions can be computed accordingly based on equation 3.1 and 3.2. In the pseudosteady case, the flow can be treated as locally steady, while in the truly non-stationary flow (e.g. planar shock propagation over a cylinder), this treatment is no longer valid. Maybe this is the reason that a growing number of studies have been focused on this area over the past four decades. Most of the experiments have been conducted on the interaction between an incident shock with constant velocity and circular surfaces [7, 11, 12, 31, 41, 83] and a smaller transition angle compared to the pseudosteady case was reported. In other words, the transition from RR to IR in the truly non-stationary flow was delayed. However, recently in 2014, Kleine et al. [47] argued that the delayed transition angle observed from previous experiments may be mostly caused by insufficient optical resolution due to the fact that unlike numerical predictions,noneoftheexperimentalapproacheswerecapableofdetectinganearlystage of Mach stem with a characteristic length below approximately 0.05 mm. Most recently, an experimental study with high spatial resolution was conducted by Ram et al. (2015) [65]. Compared to previous experimental studies, closer transition angles to the sonic criterion were obtained by using a single-lens reflex camera with the ability to detect a characteristiclengthofaround0.06mm. Meanwhile, recentnumericalresults[29,30,86] indicated that the transition occurs at the sonic condition in both pseudosteady flows and truly non-stationary flows if viscosity is not included. In the current study, one of the objectives is to predict the transition condition from RR to IR. Unlike most of the previous studies on the truly stationary flows, there is no 39 interaction between shock waves and circular surfaces involved. Thus, the viscous effect on the surface, which is known to delay the transition angle [13], does not apply to our case. Based on the conclusions from the above-listed studies, it is most likely that the sonic criterion is the answer to our setup (i.e. interaction between two identical shock waves). In addition, several numerical experiments based on two different setups have been conducted and compared to the sonic transition line, see Figure 3.4. In the first numerical setup, an incident shock passing over a cylinder is studied, which is shown in Figure 3.5. Transition angles for multiple Mach numbers and cylinder radii are plotted as blue dots in Figure 3.4. In the second setup, shock wave interactions are performed using the same initial conditions as those introduced in Section 2.1.1. To reduce the computational cost, a symmetric condition is set on the bottom wall, which is shown in Figure 3.6. The initial energy and height of the shock wave are varied from case to case, and transition conditions for three different cases are plotted as red dots in Figure 3.4. All the simulations are computed using the same code as introduced in Section 2.1 with a uniform grid size of 0.04 mm in both horizontal and vertical directions. All the boundaries in the computational domain are set to slip walls. The transition angle is determined by the location at the first occurrence of a visible Mach stem. As can be seen in Figure 3.4, good agreement is obtained compared to the sonic condition (within 3.5% difference). The small discrepancy at high Mach numbers can be addressed to two reasons; First of all, those interactions occur at the early stage when the shock front is still strong. But the length scales of the shock front are smaller compared to those at the intermediate stage. Thus, finer grid sizes are needed to resolve all the small length scales during the early interaction process. Secondly, due to the time difference between each screenshot when running the simulation is constant, it is harder to capture the right transition moment when the Mach number is higher. 3.3 Analytical Prediction on the Expansion of Shock Waves As can be seen from Figure 3.4, in order to predict the transition condition, both Mach number and interaction angle are needed. However, during the interaction process be- tween two identical shock waves, both the wave Mach number and the interaction angle 40 1 2 3 4 5 0 10 20 30 40 50 60 Sonic condition Data from setup 1 Data from setup 2 M θ [deg.] Figure 3.4 Comparison between the sonic condition and actual transition from numerical re- sults. (a) (b) (c) Figure 3.5 Schlieren visualization illustrating the first numerical setup. An incident shock is passing over a cylinder from left to right. (b) Regular reflection; (c) Irregular reflection. 41 (a) (b) (c) Figure 3.6 Schlieren visualization illustrating the second numerical setup. A cylindrical shock wave is expanding and interacting with the bottom wall. (b) Regular reflection; (c) Irregular reflection. are changing with time. To acquire such information without using numerical methods, analytical prediction on the expansion of shock waves are explored. As mentioned in Section 1.1.2, there exist numerous analytical studies on the propagation of a single cylindrical or spherical shock wave. However, only two of them can describe the entire shock wave propagation regime with good accuracy [4, 71]. Here, the analytical method from Bach (1970) [4] is utilized to describe the expansion of shock waves. In addition, a cylindrical shock wave expansion progress with a 5000 J initial energy is conducted by both the analytical method and the Euler equations. Results are plotted on Figure 3.7 and an excellent agreement is achieved (within 1% difference). From the analytical method and given an initial energy E, a relation, denoted as f, between the shock front Mach number M and the shock radius R can be obtained, which can be represented as M = f(R). In addition, given the distance between two shock centers S, the shock interaction angle θ can be calculated as sinθ = S/2R. A schematic illustration of shock-shock interaction configuration is shown in Figure 3.8. Based on this geometrical relation and M = f(R), a shock front Mach number and shock interaction angle relation can be established. An example of such relation given a specific S and E is illustrated in Figure 3.9. When two cylindrical shock waves initially contact with each other, they are normal to each other. As they keep interacting, due to the expansion of the wave fronts, θ decreases and so does M. Combining this with 42 10 20 30 40 50 60 70 0 2 4 6 8 10 Analytical Numerical M R [mm] Figure 3.7 Comparison between the analytical method [4] and the Euler equations for a cylin- drical shock wave expansion progress with a 5000 J initial energy. S θ Figure 3.8 Schematic illustration of shock-shock interaction configuration between two expand- ing cylindrical shock waves. 43 2 2.2 2.4 2.6 2.8 40 50 60 70 80 90 θ [deg.] M θ = 90 ◦ θ Figure 3.9 An example of shock front Mach number versus shock interaction angle. At the beginning of the interaction, θ is 90 ◦ . As shock fronts keep expanding, both θ and M decreases. the sonic condition, the transition condition can be predicted. 3.4 Optimization on the Number of Shock Waves At the beginning of this Chapter, an optimization problem is raised. Here a detailed de- scription is presented as follows: a total energy summed over all sources is kept constant at 5000 J with the number of shock waves varies. All sources are placed in a concentric pattern around a target with simultaneous initiation. The distance between each source and the target is fixed as 30 mm. The objective is to maximize the peak pressure at the target. Perhaps this is one of the simplest optimization problems since it has only one design variable, which is the number of sources N. However, due to the computational cost of the Euler equations, it can be unfeasible if thousands of iterations are needed. Based on the results from previous sections, an alternative way to solve this problem can be achieved. By using the method introduced in the previous section, transition conditions on different numbers of sources are predicted and shown in Figure 3.10. It 44 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 80 90 Sonic condition N = 5 N = 10 N = 20 N = 40 N = 60 N = 80 N = 200 N = 1000 N increases IR RR θ [deg.] M Figure 3.10 Transition conditions on different numbers of sources. can be seen that as the number of sources increases, even though weaker cylindrical shock waves are initialized, bigger Mach numbers are obtained at the transition condi- tion, which indicates that stronger converging shock waves are developed at transition. Therefore, a monotonic behavior on the peak pressure at the target can be predicted. However, it needs to be mentioned that in reality, this is certainly not true. As the converging shock gets stronger, thermodynamic conditions that are generated behind the shock are higher. At some point, real gas effects can no longer be neglected. 3.5 Chapter Summary In this Chapter, coalescence of multiple shock waves is classified into two phases, ex- pansion phase and focusing phase. Theoretical and analytical methods are explored to describe the expansion phase and the phase transition condition. From the intensive studies in the past decades on shock reflection phenomenon, they imply that the sonic criterion is capable of describing the transition condition from RR to IR in unsteady flows. Inaddition, numerical simulations in two different setups areconducted andgood 45 agreement is achieved when compared with the theoretical work. An analytical method to predict a single cylindrical shock wave expansion process is implemented and compared with the Euler equations. Excellent agreement is obtained. Combining with the theoretical study and some geometrical manipulations, the phase transition condition is predicted successfully. Furthermore, investigation on the optimal numberofcylindricalshockwavestomaximizethepeakpressureatthefocalregionwhile keeping the total energy as constant is performed. By using the technique developed in this Chapter, a monotonic behavior is predicted without using the Euler equations. However, inreality, symmetrical placement cannotalways beachieved. Forproblems withunsymmetricalshockwaves,seeFigure3.11,theshockfrontstrengthinthefocusing phase isnotalways increasing. Onlythe wave in theMach stem area, which is inthe red circle in Figure 3.11(c), is enhanced. Thus, to be capable of optimizing such problems, a method that can predict the focusing phase with high efficiency is required. In the next Chapter, an alternative method to describe shock wave propagation is explored. 46 Target Cylindrical Shocks (a) (b) (c) (d) Figure 3.11 Schlieren visualizations of a case with three unsymmetrical sources, (a) setup demonstration, (b) expansion phase, (c) transition, (d) focusing phase (Mach stems are shown in the red circles). 47 Chapter 4 A Reduced-order Model for Shock Wave Focusing Basically, two major numerical methods are used to study shock wave propagation, one is the Navier-Stokes equations, and the other is the inviscid Euler equations if viscosity in the shocked medium can be neglected. The advantage of these two methods is that a fullflowfieldcanbeaccuratelyobtained. However, thestrengthoftheshockfrontinthe focusingareacanbemuch higherthanthatoftheinitialshockfront,resulting insmaller time scales to maintain the Courant-Friedrichs-Lewy (CFL) condition. In addition, the length of the shock front in the focusing area can be much smaller than that of the initial shock front. In order to resolve all small scales close to, and in, the focal region, high resolution in both time and space are required, which can make the computational task very expensive. Whithamproposed analternative method[91], namedGeometrical Shock Dynamics (GSD), that describes the motion of shock waves in a different way. Unlike the Navier-Stokes equations and the Euler equations, this theory avoids dealing with the flow field around the shock andonly focuses onthe curvature ofthe shock wave itself. As a result, solving a shock focusing event with GSD instead of the Navier-Stokes equations or the Euler equations, the computational complexity is reduced by solving a lower dimensional problem. In addition, the actual computational cost may be reduced by more than an order of magnitude depending on the required grid resolution when dealing with higher dimensional problems. In GSD, the shock front is discretized into small elements. Between each element, orthogonal trajectories are introduced as rays so that each shock front element can be approximated to propagate down a tube whose boundariesaredefinedbytherays, aso-calledraytube. ThemainassumptioninGSDis thatthe motionofthe shock only changes withthe variationofthe ray tubearea. Then, instead of solving the full Euler equations, the motion of the shock can be predicted by deriving the relation between shock strength, which can be represented by the Mach 48 number, M, and the area of the ray tube, A. This is the so-called Area-Mach number (A-M) relation, 1 A(x) dA dM =−g(M), (4.1) where g(M)= M M 2 −1 (1+ 2 γ +1 1−μ 2 μ )(1+2μ+ 1 M 2 ), (4.2) μ 2 = (γ−1)M 2 +2 2γM 2 −(γ−1) . (4.3) Here, γ represents the adiabatic index, x denotes the distance in the ray tube and the cross sectional area A(x) is a function of x, see reference [92] for additional details. In this Chapter, numerical implementation of GSD in clusters is introduced. Com- parison with the Euler equations and some existing experimental studies is investigated. ModificationoftheoriginalGSDtheoryisalsopresentedtoimprovetheresultsforshock expansion progress. 4.1 Numerical Methods There exists a variety ofalgorithms to implement GSD numerically. Here, only a partof those works arelisted. The method ofcharacteristics was used by Bryson and Gross[18] to analyze shock diffractions. Decades later, a front tracking approach was developed by Henshaw [33] in two dimensions and Schwendeman [73] in three dimensions. A con- servative finite difference method was formulated by Schwendeman [75]. Most recently, a fast-marching like algorithm was developed by Noumir et al. [57]. Among all these algorithms, the front tracking method is the most used one due to its computational accuracy and simplicity of implementation. For example, Schwendeman applied this method to study shock motion in non-uniform media [74], Best modified it to investi- gate underwater explosions [8], and Apazidis and Lesser utilized it to compare with the shock converging experiments [3]. However, inorderto utilize the fronttracking method 49 to study shock motion for large-scale applications, the computational cost is still a large bottleneck, which can be addressed through parallel computing techniques. In this study, the front tracking method is implemented on parallel computers and symmetric boundary conditions are developed in order to reduce the computational expense dramatically. 4.1.1 Serial Scheme The numerical implementation of GSD is based on previous front tracking methods [8, 33]. A short review of implementation of the numerical scheme is given as follows. In a two-dimensional condition, the shock front is discretized into N points denoted by x i (t), wherei= 1,2,...,N. The shock front propagates along the direction of its normal vector, and the speed is determined by the A-M relation. The motion of the shock is described by dx i (t) dt =a 0 M i (t)n i (t), i =1,...,N, (4.4) where M i (t) and n i (t) denote the Mach number and shock front normal at x i (t), re- spectively. The speed of sound in ambient air is fixed asa 0 = 1 in all ofour calculations. A fourth-order Runge-Kutta scheme is adopted to numerically integrate equation (4.4). TheMachnumber,M(t), iscalculatedbycombining theA-M relationfromequation (4.1)withtheshockfrontrelationd t A =a 0 Md x A. Thusequation(4.1)canbeconverted into dM dt = M −g(M) A ′ A . (4.5) In equation (4.5), A ′ denotes d x A and A ′ /A can be expressed explicitly as [8], A ′ A = ∂x(s(t),t) ∂s(t) ∂n(s(t),t) ∂s(t) , (4.6) where, as is shown below, the arclength s(t) represents the geometry of the shock front 50 andn(s(t),t) indicates the norm vector of the shock front, s i (t) =    0, if i= 1; s i−1 (t)+|x i (t)−x i−1 (t)|, if i= 2,...,N; (4.7) n(s(t),t) = ∂y(s(t),t) ∂s(t) ,− ∂x(s(t),t) ∂s(t) . (4.8) Figure 4.1Flow chart for serial version. The equation to maintain the CFL condition is marked by ∗. Figure from [64]. Following [33], there are two extra procedures. One is a point insertion and deletion approach that maintains the resolution of the shock front and the CFL condition. Ad- ditionally, a time-step reduction scheme [72] should be applied to replace this approach for a particular converging shock scenario. The other scheme is a two-step smoothing procedure to reduce the high-frequency errors. The main focus of this study is to imple- ment the algorithm for converging shocks. Thus, here the time-step reduction scheme is adopted. However, the parallel implementation can be coupled with the point insertion 51 and deletion approach to tackle other shock wave simulations. The time-step reduction scheme can be expressed by the following equation, Δt<f(R(t),M(t)), (4.9) wheref(R(t),M(t)) is a function of shock radiusR(t) andM(t) at timet. More details ofthis scheme can be foundin [72]. Ifthis condition fails, Δt will bereduced as Δt new = 0.75Δt and equations (4.4) and (4.5) are restarted by using Δt new . The two-step smoothing procedure is given as follows, x i (t)= 1 2 (x i−1 (t)+x i+1 (t)), (4.10) and it is applied every n s iterations, where n s depends on the time increment, Δt, and the average shock arclength between each discrete point, Δs avg . When Δs avg =0.01,n s is usually set between 10 to 50. A schematic flow chart shown in Figure 4.1 illustrates how this algorithm works. In general, the algorithm consists of four functions during the time marching process. At timet,A ′ i /A i iscalculated first by equation (4.6). The two components∂x(s(t),t)/∂s(t) and ∂n(s(t),t)/∂s(t) in equation (4.6) can be obtained by applying a cubic spline in- terpolation method using discrete data s i (t),x i (t) and s i (t),n i (t),i = 1,...,N under adequate boundary conditions, which vary according to the task settings. This step is named updateAda. The advantage of setting updateAda initially is that it can also be used to generate the norm vector of the shock front, see equation (4.8), which is required for updating the shock location in the third step. Next, a fourth-order Runge- Kuttaschemeisemployed, seeequation(4.5),toupdatetheMachnumberatthecurrent time iteration. This step is named updateM. In the next function, called updateX, the shock location for the next time iteration is calculated by using equation (4.4) so that x(t+Δt) is obtained. Later, the time-step reduction scheme is applied to maintain the CFL condition, see equation (4.9). The next function, updateS, computes s(t + Δt) through equation (4.7). In addition, the smoothing procedure, see equation (4.10), is implemented afterupdateS. At the end of each iteration, a terminate condition is given 52 to determine whether the program should be aborted. In order to run this algorithm, a set of coordinates,x i (t 0 ), and Mach number, M i (t 0 ), representing the initial location and strength of the shock frontare given as initial conditions. In addition,updateX and updateS need to be called before the time marching to obtains i (t 0 +Δ(t)), which is the arclength for the first iteration. These initial steps are implemented in the parameter initialization, shown in Figure 4.1. 4.1.2 Parallel Scheme Performance profiling has been done to detect the computational cost for each function describedintheprevioussection. Thetestsetupisa10-sidedpolygonal(decagon)shape shock wave propagatingtowards its center point. Three different cases have been tested, see Table 4.1. The Mach number, smooth step interval and the initial shock arc length between each discrete point are kept constant as M i (t 0 ) = 15,n s = 200,Δs(t 0 ) = 0.002 for all cases while the apothem of the decagon (a line segment from the center to the midpoint of one of its sides) varies from 1 to 100. All three cases are terminated after 10,000 iterations. Figure 4.5 shows an example when the apothem equals 1 and the program ends after 55,000 iterations. The blue solid polygons represent the initial and the first reconfiguration steps, which will be explained in the validation section. The red solid polygons represent the successive shock fronts at each iteration sequence. Table 4.1 shows the computational costs of each function in different cases. It can be seen that most of the computations are conducted by the first four functions, occu- pying over 98% of the total running time. And updateAda is the most expensive one, which occupies over 64% of the total running time. The reason is that to solve equation (4.6),∂x(s(t),t)/∂s(t) needs to be obtained first by applying four cubic spline functions since periodic boundary conditions are used. Then, n(s(t),t) can be achieved through equation (4.8) automatically. Later, another four cubic spline functions are employed to access ∂n/(s(t),t)∂s(t), where n = (n x ,n y ). Thus, there are eight tridiagonal systems of linear equations in total that have to be solved because each cubic spline function can be considered as a tridiagonal system. Commonly, Gaussian elimination can solve linear systems in a serial manner. However, for parallelization purpose, it is not efficient 53 Table 4.1 Computational cost of each function for three different tests. A* and N* denote the apothem and the total number of points of the polygon. A*=1 N*=3,242 A*=10 N*=32,491 A*=100N*=324,911 time (s) percentage time (s) percentage time (s) percentage updateAda 21.5 64.56% 223.19 65.01% 2539.98 67.82 % updateM 9.36 28.11% 94.52 27.53% 937.56 25.03% updateX 1.13 3.39% 11.99 3.49% 115.32 3.08% updateS 0.75 2.25% 8.38 2.44% 92.01 2.46% time-step reduction 0.11 0.33% 1.05 0.3% 10.78 0.29% smooth function 0.188 0.56% 2.12 0.60% 20.51 0.55% compared to other schemes [25, 61, 62, 79, 90]. In this study, a recent approach [25] for a parallel tridiagonal solver has been applied with the following features: the new solver, under the model of the SPIKE algorithm [61, 62], utilizes message passing inter- face (MPI) forinterprocessor communication. In addition, the collective communication is not required, and the machine-zero accuracy is guaranteed. The solver can compute a number of distinct tridiagonal systems together by being executed only once which is more efficient than being called multiple times to solve each single system. Multiple boundary conditions such as periodic boundary conditions and natural boundary condi- tions can be achieved without modifying the solver. Therefore, the solver is used in the function updateAda. For the rest of the functions, both the shared-memory multiprocessing interface OpenMP and MPI have been implemented for parallelization. Thus, combined with the function updateAda under MPI, the whole algorithm is parallelized using two ap- proaches: hybrid MPI combined with OpenMP, and pure MPI. Hybrid MPI Combined with OpenMP Approach In the hybrid scheme, the MPI section and the OpenMP section are designed to be independent of each other to reduce the implementation difficulty. A flow chart is shown in Figure 4.2 to represent the architecture. At first, the input data is partitioned 54 Figure 4.2 Flow chart for hybrid MPI and OpenMP scheme. The equation to maintain the CFL condition is marked by ∗. Figure from [64]. into multiple subsets, each of which is assigned to a particular node, to satisfy the precondition of the MPI tridiagonal solver in theupdateAda function. Figure 4.3 shows an example of the spatial decomposition for the case with a 5-sided polygonal shock in which the shock front has been split into 21 points that are distributed among five processors. It can be seen that the first point and the last point on the shock front share the same location because of the closed geometry. After calculation by the MPI tridiagonal solver, the resulting data is gathered into one node as a whole, and then, broadcasted to all nodes. The purpose of this broadcast is to guarantee that every node has the entire dataset. The gather and broadcast are referred as allGather in MPI section. In the OpenMP section, each node conducts the same computation by using its own processors with the same spatial decomposition method as shown in Figure 4.3. At the endofeach iteration, each nodewill contain thewhole resulting data. Consequently, there is no cost due to communication between nodes in the data partitioning step, so in other words, the communication only happens in the MPI tridiagonal solver and the allGather steps. The main advantage of this approach is that all processors share the same memory, and as a result, parallelization can be easily achieved for all the loop- independent functions with little communication cost. The implementation process is 55 considerably simpler compared to the MPI approach because of the limited amount of data communication. • • • • • • • • • • • • • • • • • • • • Processor 0 Processor 1 Processor 2 Processor 3 Processor 4 4 2 3 0 1 19 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Figure 4.3 Example of spatial decomposition of a 5-sided polygonal shock front split into 21 points distributed among five processors. Figure from [64]. MPI Approach The pure MPI structure is presented in Figure 4.4. In the beginning, all the data is partitioned and assigned to multiple nodes. Compared to the previous hybrid approach, the data-gathering step is eliminated by rescheduling data allocation for loop-carried functions so that each node can conduct its own work without requiring additional information from the other nodes. To be specific, in the functions updateAda and updateS, point-to-point communications have been used to allocate the data to avoid data dependency across adjacent nodes. Although this approach increases the difficulty in both the implementation and debugging process, it is expected to achieve a better parallel performance than the previous approach. 4.1.3 Comparison with Analytical Solution To compare the two parallel schemes, numerical simulations of regular polygonal con- verging shocks with 9 sides (enneagon) and 10 sides (decagon) have been conducted on the above schemes with the same initial conditions as indicated in section 2.2 and 56 Figure 4.4 Flow chart for MPI scheme. The equation to maintain the CFL condition is marked by ∗. Figure from [64]. Table 4.2 Comparison of numerical results with analytical solutions from [72]. Enneagon simulation (MPI) Decagon simulation (hybrid) Analytical Numerical Difference Analytical Numerical Difference M 1 /M 0 1.175 1.161 1.2% 1.155 1.159 0.3% r 1 /r 0 0.442 0.446 0.8% 0.482 0.486 0.8% the apothem is one. Figure 4.5 shows the numerical results, in which the successive shock fronts of each iteration sequence are plotted. As the polygonal converging shock propagates towards the center, the original shape keeps reshaping, and a Mach stem is generated at each corner. In Figure 4.5, r 0 and M 0 represent the distance from the center to the initial side of the polygon and the initial Mach number, respectively, while r 1 andM 1 denote the distance from the center to the first repeated polygon of the same shape as the initial polygon and its Mach number. From Table 4.2, it can be seen that the numerical results from the parallel schemes agree well with the analytical solutions. 57 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 M, r M!, r! −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 M, r M!, r! Figure 4.5 Converging shock propagation showing successive positions for an enneagon (MPI scheme with 4 cores) and a decagon (hybrid scheme with 2 nodes and 2 processors per node). Figures from [64]. 4.1.4 Benchmark Performance analysis has been conducted for the two parallel schemes on the high- performance computing facility (HPCC) at University ofSouthern California. The node set, HPCC sl250s, consists of 256 Dual-Octacore Intel Xeon CPUs operating at 2.4 GHz with 64 GB memory. A strong-scaling test has been performed for both schemes by simulating a decagon converging shock with a total number of 3,242,911 discrete points. The initial conditions arethe same asthose described in section 2.2. As previously men- tioned, in the hybrid scheme the MPI section and the OpenMP section are independent of each other. To perform the strong-scaling test for both sections, the number of nodes is defined as N node and the number of processors per node is denoted as P pn . In addi- tion, the granularity in the MPI section can be defined as N/N node , where N denotes the number of discrete points on the shock front, while in the OpenMP section, the granularity is defined as N/P pn . Then, the strong-scaling test for the hybrid scheme is performed as follows: the granularity in both sections is kept consistent with each other, which means N node is always set to be equal to P pn , while N is kept constant. In this test, we defineβ =N node =P pn , which ranges from 1 to 16. Consequently, the efficiency of the hybrid scheme is defined as the ratio of the speedup of the parallel code divided by β. Figure 4.6 shows the wall-clock time and the efficiency of the hybrid scheme. As 58 a baseline, the ideal case is calculated by the ratio of the wall-clock time for the case with N node =P pn =1 over the number of different β. Figure 4.6 Wall-clock time of strong scaling for hybrid scheme using 3,242,911-point shock front on β nodes and β processors per node. Figure from [64]. For the MPI scheme, the strong-scaling test is performed with the same setup as above,exceptthatN node isvariedfrom1to64andP ppn issetto1forallcases. Therefore, the efficiency of the MPI scheme is defined as the ratio of the speedup of the parallel code divided by N node , and the ideal case is calculated by the ratio of the wall-clock time on a single node over the number of different nodes. As shown in Figure 4.7, the hybrid scheme has a comparable performance with the MPI scheme up until the case of β = 8. The efficiency of the hybrid scheme for the case of β = 2 is 0.92, and for β = 4, it is 0.83. However, as β increases, the efficiency drops from 0.72 with β = 8 to 0.6 with β = 16 cores. In contrast, the efficiency of the MPI scheme stays above 0.9 until 64 nodes are used. The main reason for the difference between the two schemes is that there is a data gather-broadcast step in the hybrid scheme to connect the MPI section with the OpenMP section, which is denoted asallGather in Figure 4.2. This additional, albeit necessary step, makes the communication/computation ratio in the hybrid scheme higher than that of MPI scheme. With the increment of N node and P pn in use, the communication/computation ratio increases, and results in a loss of 59 performance. In the MPI scheme, without this extraallGather step, the wall-clock time decreases nearly linearly with the granularity (N/N node ). The reason that the efficiency in the MPI scheme drops slightly for the case with N node = 2 is likely due to that the communication costisintroduced inthecomputationcomparedtothecase withasingle node. When the number ofnodes increases to 64, the communication cost starts to have a notable effect on the performance. Figure 4.7 Wall-clock time of strong scaling for MPI scheme using 3,242,911-point shock front on N node nodes. Figure from [64]. Additionally, an isogranular-scaling test, where the number of discrete points on the shock front in each node, N/N node , is fixed, has been performed for the MPI scheme. The test case is a decagon shock with 25,991N node discrete points while N node ranges from 2 to 128 and P pn is set to be 1 for all cases. The weak scaling efficiency for the case ofN node nodes is calculated ast 2 /t n , wheret 2 andt n denote the wall-clock time for a 2-node test and a N node -node test, respectively. Figure 4.8 shows the wall-clock time versus the number of nodes scaled with the number of points. The efficiency remains nearlyconstantaround0.95untilthecaseof64nodes,whichindicatesagoodscalability. 60 Figure 4.8 Wall-clock time of scaled workloads using 25,991N node -point shock front on N node nodes. Figure from [64]. 4.1.5 Symmetric Boundary Condition A symmetric boundary condition is developed to further improve the speedup for con- verging shock configurations. From Figure 4.5, it can be observed that as the polygonal converging shock propagates towards the center, the line of symmetry of the shock front is independent of its reconfiguration process. Thus, it is possible to compute only one section of the shock front instead of the full geometry. This symmetry feature indicates thatthenormvectorateachpointwherethelineofsymmetrymeetswiththeshockfront is always directed to the center. This is the boundary condition for each symmetric part of the shock front. For example, since the shock front features a polygonal shape, geo- metrical symmetry is used such that only a triangular portion of the converging shock, from a corner where two shocks meet to the mid-point of a neighboring planar side, is considered. The shock velocity oneither side ofthis triangle isalways directed along the sides of the triangle, i.e. the symmetry lines. With this condition, the problem size can be reduced fromN down toN/2n l , whereN denotes the number of discrete points and n l denotes the number of symmetric lines. A test on 12-sided polygonal (dodecagon) converging shock has been performed to verify the boundary condition. Three different 61 Table 4.3 Computational cost and speedup of 1/2, 1/4, 1/24 of the full geometrical simula- tion using symmetric boundary conditions. All the simulations are terminated after 10,000 iterations. Problem size full size 1/2 size 1/4 size 1/24 size Computational cost 17.15 8.65 4.44 0.89 Speedup 1.0 1.98 3.86 19.26 symmetric parts (1/2, 1/4 and 1/24) have been computed independently and compared with the full size, see Figure 4.9. It can be seen that all of the symmetric parts match the full-size result very well. For performance analysis, the speedup of the symmetric boundary condition is defined as wall-clock time ratio of the full geometry versus the reduced geometry using a single core. Results are summarized in Table 4.3, and it is shown that for the dodecagon shock front, the speedup can be obtained up to 19.26. −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 (a) 1/2 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 (b) 1/4 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 (c) 1/24 Figure 4.9 Dodecagon shock front propagation at successive time instants. Red dashed line represents the full geometrical simulation, blue solid line represents, from left to right, 1/2, 1/4, 1/24 of the full size simulation using symmetric boundary conditions. Figures from [64]. 4.2 Comparison with the Euler Equations From previous sections, it can be seen that tremendous computational cost can be re- duced by using the parallel computing techniques introduced above. However, the ac- curacy of this approximation theory for our specific problem is difficult to determine in advance. In this section, the simulation results from geometrical shock dynamics are compared with the results from Euler equations. 62 A comparison of simulation results from GSD and Euler equations for cases with 5 and 10 energy sources simultaneous initialization is performed. In these two cases, it can be observed that as the shock fronts coalesce and interact with each other, they form a polygonal shape converging shock as shown in Figure 4.10. As discussed in previous chapters, a reconfiguration process starts when Mach stems are generated on the corner of the converging shock where the irregular reflection occurs. Since these Mach stems propagate faster than the nearby portion of the shock front, they overtake them and reconfigure the original shape. This process repeats itself until the shock reaches the center. Figure 4.10 shows that if the converging shock front is presumed to be a regular polygon, by following Schwendeman’s [72] work, the converging shock propagation process can be calculated using the theory of GSD. (a) (b) Figure 4.10 Pressure contour illustrating reconfiguration process from the cases with 5 (left) and 10 (right) sources. Red represents high pressure region and blue shows low pressure region. Outer polygons indicate the initial converging shock front while inner polygons illustrate shock front geometry after one reconfiguration. TosimulatethisprocessbyusingthetheoryofGSD,thephasetransitionconditionis pre-calculatedbasedonthemethodsdevelopedinChapter3. Thus, theexpansionphase 63 does not need to be computed anymore. At the transition condition, Mach stems are formed and a polygonal-like converging shock starts to propagate. The initialization for GSD simulations are set at this time and assumptions are made that the strength of the Mach stems at the very early stage are the same as the other portion of the converging shockfront. Inaddition,duetothefactthatthecurvatureofthisconverging shockfront decreases for each reconfiguration, eventually resulting in a regular stabilized polygonal converging shockfront,theinitialcurvedconverging shockfrontissimplified toaregular polygonshapeforcurrent comparison. Figure4.11shows initialshock frontpositionsfor cases with 5 and 10 sources. They are obtained from the pre-calculated phase transition conditions and the shape of the shock front at this moment is approximated to a regular polygon, which can be seen as the solid lines in Figure 4.11. All the sides of regular polygons are chosen to be tangent to the original curved structure. As can be seen from Figure 4.11, a good agreement between the initial simplified curve for GSD and the curvature at the transition from the Euler equations is achieved. Then, by applying the numerical scheme introduced from previous section, Mach number at successive shock positions are computed. For the case with 5 sources, Mach number of the shock front along the x-axis (shown asdashed line inFigure4.11(a))is shown in Figure4.12. Onthe leftofthex-axis, itcan be seen that the shock wave front propagates inward, while on the right side, the wave frontcoalesces (twoshock waves intersection) andtravelstowards thecenter. Thecenter is located atx = 0, however, for GSD simulations, when the converging wave focuses at thecenter,asingularityoccurs. Thus, forbothcases,thesimulationsarestoppedaround 2 mm away from the center. The red line and the blue line represent the Mach number from the Euler equations and GSD respectively. As the shock wave front propagates from left to right, before the Mach stems overtake it, its Mach number decreases, which can be seen as the red line on the left of Figure 4.12. This is because this portion of the converging shock front is still expanding, see the convex curvature on the left side of converging shock front in Figure 4.11(a). This behavior can be predicted by the Euler equations. However, since the current theory of GSD is based on the assumption that the motion of the shock changes only with the variation of the tube area, if the tube area does not change, neither does the shock strength. This is the reason that the blue 64 (a) (b) Figure 4.11 Schlieren visualizations of initial shock geometry from the cases with 5 (left) and 10 (right) sources. line on the left stays the same as it moves inward. When the coalesced wave front travels from right to left, a Mach stem is generated and the Mach number is increased. Shown in the right of Figure 4.12, the Mach number of both methods is increased at the beginning, indicating the generation of the Mach stem. After the Mach stem is formed and overtakes the neighboring wave fronts, the results show a similar trend as observed in the left of Figure 4.12: For the GSD simulation, the Mach number that is denoted by the blue line stays at the same level while for the Euler equations, the Mach number, which is represented by the red line, drops. For the case with 10 sources, the Mach number of the shock front along the x- axis and the y-axis (shown as dashed line in Figure 4.11(b)) is shown in Figure 4.13. Figure 4.13(a) shows the shock front propagates from right to left along x-axis. In the beginning, fromx = 19 mm tox =10 mm, bothlines show a similar trend as the results observed on the left of Figure 4.12. However, when the shock front reaches y = 10 mm, the Mach stems that are generated by the neighboring coalesced waves start to overtake 65 −8 −6 −4 −2 0 2 4 6 8 10 1.8 2 2.2 2.4 2.6 2.8 3 Euler simulation GSD simulation M Location [mm] Figure 4.12 Mach number history along x-axis for the case with 5 sources, left, shock wave front propagation from left to right, right, coalesced wave front propagation from right to left, x = 0 is the center of the converging wave. the shock front at the x-axis and create a new Mach stem. Then this whole process repeats at y = 2 mm again. Figure 4.13(b) displays the coalesce wave transmits from right to leftalong the y-axis. In thebeginning ataroundx = 20 mm, bothlines show an increase in Mach number that indicates the generation of Mach stem. Later, at around x = 8 mm, this process repeats again, which implies the second reconfiguration of the converging shock. Although Mach numbers from both methods do not agree well with each other, the shock transition locations where Mach stems are generated match well. The reasons for the differences between Mach numbers obtained from these two methods are twofold: first, the originalconverging shock frontisa curved structure, while inGSDsimulations, the arc structure of shock front is simply approximated to a regular polygon. Second, the A-M relation is based on the assumption that the flow behind the shock is uniform, which is not valid to for the shock wave expansion case. Thus, to improve the accuracy of current GSD theory for our particular cases, the above shortcomings are addressed. 66 0 2 4 6 8 10 12 14 16 18 20 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 Euler simulation GSD simulation M Horizontal location [mm] (a) 0 5 10 15 20 25 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 Euler simulation GSD simulation M Vertical location [mm] (b) Figure4.13Mach number history for the case with 10 sources, (a), shock wave front propagation from right to left along x-axis, (b), coalesced wave front propagation from right to left along y-axis. The center of the converging wave is at (x,y) = (0,0). 4.3 Modified GSD When a cylindrical or spherical shock front is expanding outwards, its strength is de- creasing due to two facts: the expanding geometry and the continuing interaction with the nonlinear flow behind it. In the original theory of GSD from Whitham (1957) [91], only problems where the nonlinear geometrical effects are the primary reason for the changes in the shock motion are considered. However, for our cases, the interactions with the flow behind are as important. Thus, the base of GSD theory, which is theA-M relation, is revisited without taking the original assumption into account. 4.3.1 Area-Mach Number Relation ThederivationoftheA-Mrelationisshortlyreviewed, however,moredetailsareoutlined in Whitham (1957) [91] and Best (1991) [8]. Consider that a shock propagates along a tube of which the cross section changes slowly. If x is used to represent the distance in the tube, the cross-sectional area A(x) which is a function of x is uniform at the 67 beginning. A(x) =A 0 =constant, x<0. (4.11) As the cross section starts to change slowly, the derivatives of A(x) are small, A(x)−A 0 A 0 ≪ 1, x≥ 0. (4.12) Then, by averaging the equations of one-dimensional inviscid compressible flow for con- servation of mass, momentum and energy across the tube, a set of partial differential equations can be obtained as follows, ∂ t ρ+u∂ x ρ+ρ∂ x u+ρu A ′ (x) A(x) = 0, (4.13) ∂ t u+u∂ x u+ ∂ x p ρ =0, (4.14) ∂ t p+u∂ x p−a 2 (∂ t ρ+u∂ x ρ) =0, (4.15) where ρ, u, p, a denote the density, particle velocity, pressure and local speed of sound respectively. A ′ (x) represents dA dx , the sound speed is given bya= γp ρ , where for an ideal gas γ is the ratio of specific heats. To solve equations (4.13) to (4.15) forx> 0, they are first linearized about the state u 1 , p 1 , ρ 1 behind the shock as follows, ∂ t (ρ−ρ 1 )+u 1 ∂ x (ρ−ρ 1 )+ρ 1 ∂ x (u−u 1 )+ρ 1 u 1 (A(x)−A 0 ) ′ A 0 = 0, (4.16) ∂ t (u−u 1 )+u 1 ∂ x (u−u 1 )+ ∂ x (p−p 1 ) ρ 1 = 0, (4.17) ∂ t (p−p 1 )+u 1 ∂ x (p−p 1 )−a 2 1 (∂ x (ρ−ρ 1 )+u 1 ∂ x (ρ−ρ 1 )) =0. (4.18) 68 Then, their characteristic form can be derived as, C + :{ ∂ ∂ t +(u 1 +a 1 ) ∂ ∂ x }((p−p 1 )+ρ 1 a 1 (u−u 1 ))+ρ 1 a 2 1 u 1 (A(x)−A 0 ) ′ A 0 = 0, (4.19) C − :{ ∂ ∂ t +(u 1 −a 1 ) ∂ ∂ x }((p−p 1 )−ρ 1 a 1 (u−u 1 ))+ρ 1 a 2 1 u 1 (A(x)−A 0 ) ′ A 0 = 0, (4.20) P :{ ∂ ∂ t +u 1 ∂ ∂ x }((p−p 1 )−a 2 1 (ρ−ρ 1 )) =0. (4.21) The general solution can be obtained from integrating equations (4.19)-(4.21), (p−p 1 )+ρ 1 a 1 (u−u 1 ) =− ρ 1 a 2 1 u 1 (u 1 +a 1 ) A(x)−A 0 A 0 +F{x−(u 1 +a 1 )t}, (4.22) (p−p 1 )−ρ 1 a 1 (u−u 1 ) =− ρ 1 a 2 1 u 1 (u 1 −a 1 ) A(x)−A 0 A 0 +G{x−(u 1 −a 1 )t}, (4.23) (p−p 1 )−a 2 1 (ρ−ρ 1 )=H(x−u 1 t), (4.24) where F, G and H are arbitrary functions which can be determined from the initial conditions of the problem and the boundary conditions at the shock. However, because the C + characteristics behind the shock originate in the uniform region where u = u 1 , p = p 1 , ρ = ρ 1 , A = A 0 , F must be zero. In addition, based on the Rankine-Hugoniot relations, the uniform state behind the shock where u = u 1 , p = p 1 , ρ = ρ 1 can be derived in terms of the undisturbed state ahead of the shock where u = 0, p = p 0 , ρ =ρ 0 and the Mach number M, u 1 = 2a 0 γ +1 (M− 1 M ), (4.25) p 1 = ρ 0 a 2 0 γ(γ +1) (2(γ)M 2 −(γ−1)), (4.26) ρ 1 = ρ 0 (γ +1)M 2 (γ−1)M 2 +2 , (4.27) with M = U a 0 , where U is the speed of the shock. Combine equations (4.22), (4.25), (4.26) and (4.27) together and apply M 0 to represent the Mach number of the shock at 69 A 0 , the A-M relation can be obtained as, A(x)−A 0 A 0 =−g(M 0 )(M−M 0 ), (4.28) where g(M)= M M 2 −1 (1+ 2 γ +1 1−μ 2 μ )(1+2μ+ 1 M 2 ), (4.29) μ 2 = (γ−1)M 2 +2 2γM 2 −(γ−1) . (4.30) Considering that during wave propagation through a tube with a sufficiently large length, massive changes can be accumulated inA(x). Thus, to make the change inA(x) smallenough,thetubeisdividedintosmallelementsandtheA-Mrelationcanbeapplied to each element. However, the uniform state behind the shock when it propagates into each element will no longer be strictly valid. By neglecting this drawback, application of equation (4.28) in the limit of infinite subdivision yields the differential equation 1 A(x) dA dM =−g(M). (4.31) 4.3.2 Modification Without applying the uniform flow assumption, we can manipulate equation (4.19) to get ∂p ∂x +ρ 1 a 1 ∂u ∂x =− 1 a 1 +u 1 ∂p ∂t +ρ 1 a 1 ∂u ∂t +ρ 1 a 1 2 u 1 (A(x) ′ A 0 . (4.32) Notethatphysicalinformationalongtheshockischangingwithbothtimetandlocation x(t). Thus, the total derivative of p and u with respect to t are dp dt = ∂p ∂t +a 0 M ∂p ∂x , (4.33) du dt = ∂u ∂t +a 0 M ∂u ∂x . (4.34) 70 Then, the following equations can be obtained by combining equation (4.33) and (4.34), dp dt +ρ 1 a 1 du dt = ∂p ∂t +ρ 1 a 1 ∂u ∂t +a 0 M ∂p ∂x +ρ 1 a 1 ∂u ∂x . (4.35) Substitutingequation(4.32)into(4.35)withsomemathematicalmanipulationscanyield the following equation, dp dM +ρ 1 a 1 du dM dM dt =− a 0 M ρ 1 a 1 2 u 1 a 1 +u 1 A(x) ′ A 0 + a 0 M a 1 +u 1 −1 ∂p ∂t +ρ 1 a 1 ∂u ∂t . (4.36) In equation (4.36), the last term on the right side, ∂p ∂t +ρ 1 a 1 ∂u ∂t , represents the flow conditionbehindtheshock. Withtheoriginalassumptionapplied,thistermisextremely small, neglecting it will result in the original A-M relation as shown in equation (4.1). In 1991, Best [8] deduced expressions for∂p/∂t and ∂u/∂t based on the assumption thatρ, p,u and a have continuous derivatives of all orders. Equation 4.36 is considered as the zeroth order of a system coupled with a first order term, which is ∂p ∂t +ρ 1 a 1 ∂u ∂t . He noted that Whitham solved this system by truncating the first order term. He investigated and improved the theory by truncating at the second order term. Detailed expressions can be found in [8]. Here, anew approachisproposedto coupletheflowconditions behindtheshock into equation (4.36), where similitude equations, which are obtained from analytic solutions, experimental data or other existing numerical data, are developed and implemented into the last term. Thus, no truncation is needed, and the accuracy is related to the development ofthesimilitudeequations. Tocomparethenewapproachwiththeoriginal theory, propagation of a cylindrical shock wave is computed by both methods. The analyticalmethodintroducedinSection3.3iscoupledintothenewapproach. Compared to the original theory, excellent improvement is achieved, and good agreement between the modified GSD and the analytical solution is shown in Figure 4.14. Thedrawbackofthisnewapproachisthatfortheinteractionofmultipleshockwaves, various similitude curves are needed to account for different types of shock expansions. Despite the expansion of the initial shock wave, the expansion of the Mach stem also 71 10 12 14 16 18 20 3 3.5 4 4.5 5 5.5 6 Analytical Solution GSD Modified GSD M R [mm] Figure 4.14 Comparison of Mach number as function of shock radius for regular GSD, extended GSD and analytical result [4] for an expanding cylindrical shock wave. needs to be included. For interaction between three or more shock waves, Mach stems interaction can occur so that the expansion of the second Mach stem needs to be con- sidered as well. Perhaps it is more clearly to be explained as a perfect binary tree. Take an example of the case with a three asymmetrical sources shown in Figure 3.11. An inverted binary tree is generated in Figure 4.15 to represent different shock expansion types. Each node represents an expansion of shock front and each level indicates the formation of a new Mach stem from the interaction between the nodes in the previous level. Due to the differences between the initialization of the shock wave in each node, to simulate such problems with the modified method, similitude equations designed for each node are required to couple the non-linear flow properties behind that shock. For- tunately, foroptimization purposes, symmetrical placement is the key, which will reduce thecomplexity ofsuch problem tremendously. Forcasesthatsources canbeplacedsym- metrically around the target, a binary tree structure is not suitable anymore. Instead, there will be infinite levels with each level represents a reconfiguration of the converging shock front. At each level, the initialization condition for each new Mach stem is the same due to the symmetrical setup. Thus, only one similitude curve is needed for each level. Furthermore, because the shape of the converging shock front tends to approach 72 a regular polygon as the reconfiguration continues (i.e. as the level increases), the flow behind the shock front will have a smaller impact on reducing the shock front strength, which means that the similitude curve can be neglected after several levels of interac- tions. From the experience gained from the simulations introduced in Section 2.3.1, the curvature of the converging shock front disappears after level two. In addition, for cases that sources cannot be symmetrically placed around the target, local symmetry can still be achieved to ensure the optimal focusing performance. For example, in the case of the three sources asymmetrically placed around a target shown in Figure 3.11, three shock waves are placed in a pattern that two sources are symmetrically placed along the third one for optimization purposes. Thus, at the same level, the similitude curves for each node are the same (same as the symmetrical case). 4.4 Comparison with Previous Experiments Combining the modified GSD with the analytical method for cylindrical or spherical shock waves and the sonic condition, a new model can be built to predict the inter- action of shock waves with equal strength. To further test the new model, a previous experiment, which was done by Higashino et al. (1991)[35], is computed numerically by both the Euler equations and the new model. For simulations conducted by the Euler equations, the computational domain was established based on the experimental setup, however, contained only a single shock wave to take advantage of the line of symmetry along the reflection plane. The computational domain was 105 mm x 150 mm with grid size 0.04 mm for both the horizontal and vertical directions. All boundaries are set as a rigid wall with a slip wall condition. The center of the cylindrical shock wave is located at (−30,0), which is the same as in the experiment [35]. From the concept of image bursts [55], when the shock hits the wall, it is equivalent to the shock interaction with another shock wave of the same strength at the opposite side of the wall. By using this setup, the computational complexity is reduced to half of the original. The energy to initialize the shock was estimated based on the total electrical energy reported by Higashino et al. (1991), in which exploding wires were used to generate the cylindrical shock waves. Theoretically, the initial energy of the shock wave, E, can 73 Level 0: Three shock fronts Level 1: Two Mach stems Level 2: One Mach stem (a) Level 0 (b) Level 1 (c) Level 2 (d) Figure 4.15 A binary tree example to illustrate different expansion types during the interaction. (b), (c), (d) show schlieren visualizations of three different expansion types. 74 be calculated as capacitive energy, E c , calculated as using the capacitance, C, and the voltage, V, expressed as E c = 1 2 CV 2 (4.37) and scaled to two-dimensional space using E = E c L(2π) , (4.38) where L is the length of the wire, as defined by [35] to account for the pseudo-two- dimensionality oftheexperimental work. Usingequation(4.37),thecomputed explosion energiesforthecopperandnichromecasesare25Jand64J,respectively. Thisenergyis thenusedinequation(4.38),wherethecalculatedenergy,E,isusedinthe2Dsimulation. The calculated explosion energy for the copper case is close to that of the experiment of 23.12 J. However, a discrepancy is evident for the nichrome case which is 169.28 J. This is due to the uncertainties derived from the experiment. Because of the finite length of wire used (34 mm), it is likely that the entire wire does not detonate evenly or at all. It is assumed, however, that the entire length is detonated simultaneously. It is also uncertain whether the entire capacitive energy is used to facilitate the shock. While this is assumed, it is possible that the wire does not require the full energy and will detonate before the capacitors have been emptied. Because of these uncertainties, multiple simulations for nichrome case at different energy levels are run. At an energy of 64 J, a good agreement with the experimental results is achieved (pressure data from the simulation is inside the experimental uncertainty). Two pressure gauges were used in the experiments to record the pressure history at two different locations along the centerline of two cylindrical shock waves. According to the experiments, two probes are placed in the same locations in the Euler simulations, and a good quantitative agreement is obtained between experiments and simulations. Table 4.4 compares the shock arrival times that were quoted by Higashino et al. to the Eulersimulations. Despitetheuncertaintyintheshock-initiationtimeoftheexperiment, the arrival times differ by only a maximum of 2μs for copper wires and by a maximum 75 Table 4.4 A comparison of shock arrival times of the current Euler simulation versus the experiment presented in [35]. Experiment (μs) Euler simulation (μs) Copper Gauge 1 83.7 83.8 Copper Gauge 2 181.2 179.2 Nichrome Gauge 1 71.2 66.5 Nichrome Gauge 2 128.7 130.6 of 4.7 μs in the case of nichrome wires. For simulations using the new model, a transition condition is first predicted by combiningthetheoreticaltransitioncurveandtheanalyticalsolutionofcylindricalshock waves basedontheinitialenergyobtainedfromtheEulerequations. Then, themodified GSD starts to kick in with the information transformed from transition condition. It should be mentioned that two similitude relations are generated for the modified GSD method due to the two different propagationstatuses in the focusing phase. Ifexplained by binary tree concept, there are two levels of propagation status, one is the expansion of the original shock wave, while the other one is the growing Mach stem. For the first level, the analytical method of cylindrical shock wave is applied to design the similitude curve. For the second level, numerical data from the Euler simulations is collected to describe thenon-linearflow properties behindthe Mach stem. Because ofthe agreement between Euler simulations and experiments, the GSD results are only compared to the Euler results. Figure 4.16 shows the qualitative comparison of the shock fronts of GSD and Euler simulations via simulated schlieren at the same time instant. The most notable difference in Figure 4.16 is that the Mach stem in the GSD result appears to move faster compared to the Euler simulation result. For the nichrome wire case, the same results were found. This is due to the modified GSDcode, which, while it takes the decaying properties behind the shock into account during the shock expansion process,itisnotabletotakeadvantageofthisduringtheformationoftheMachstem. In this case, a stronger Mach stem is formed compared to the Euler simulation results. To show this in more detail, the pressure history of the Mach stem is shown in Figure 4.17 for GSD, Euler simulations, and experiments. Here, the dotted brown and black curves correspond to the shock-shock interaction 76 (a) (b) Figure 4.16 Overlay of Euler simulation (schlieren background) and GSD (red lines) results. Two instants in time for the expansion of a cylindrical shock wave generated by an exploding copper wire by [35] are shown. (a) Regular reflection. (b) Irregular reflection. 50 100 150 200 1 1.5 2 2.5 3 3.5 4 Cu−Ex Ni−Cr−Ex Cu−Euler Ni−Cr−Euler Cu−GSD Ni−Cr−GSD t [μs] P m /P a Figure 4.17 Ratio of maximum pressure P m , at the Mach stem to ambient pressure P a , versus time t. (Experimental data recreated from [35]). 77 and resulting Mach stem by the modified GSD, the solid green and purple curves are from Euler simulations, and the red and blue data points are from experiments. For both copper and nichrome, the Euler curves can be seen to exist almost entirely inside the experimental uncertainty, while both GSD curves overshoot (faster Mach stem) at first, then converge over time to the Euler and experimental results. The computational cost for each case of the Euler simulations is more than 36 hours with 48 cores, while for the reduced-order model, it takes ten minutes with only one core. 4.5 Chapter Summary In this Chapter, a spatial decomposition method to implement the GSD simulation on parallel computers is adopted. At first, performance profiling is conducted to analyze computational hot-spot in the serial scheme. Two schemes are designed to parallelize the simulation, in which MPI and OpenMP are employed. An efficient tridiagonal solver [25] based on the SPIKE algorithm [61, 62] is also incorporated into the parallel implementation to handle the most computationally expensive function. Performance analysis and comparison between the two schemes are investigated. The hybrid scheme shows its advantage of implementation ease, and the running result demonstrates its speedup, 9.7, on 16 HPCC node-processors. On the other hand, the MPI scheme leads to an improved performance in efficiency and scalability. The strong- scaling experiment proves a high efficiency of 0.93 using 32 HPCC cores, while the isogranular-scaling experiment shows thattheefficiency holdsupto0.83using 64HPCC cores. To further improve the speedup for symmetric converging shock simulations, sym- metric boundary conditions were developed to reduce the problem size considerably. Results have shown that for a dodecagonal converging shock front, a speedup of up to 19.26 can be achieved. The comparison is performed between the Euler equations and GSD for cases with 5 sourcesand10sources. Duetothesimplifiedapproximationthattheinitialcurvedshock frontistreatedasaregularpolygonandtheassumptionthattheflowbehindtheshockis uniform, Machnumbers obtainedfromGSDcodedonotagreewell withtheresults from 78 Euler equations. However, shock transition positions from RR to IR matches well. To make GSD as a proper model for future optimization purposes, modifications of initial shape and A-M relation is made to improve the accuracy of shock expansion prediction. In addition, Euler and modified GSD simulations have been performed to match the experiments of Higashino et al. (1991) [35]. The results of the simulations showed bothgoodqualitativeandquantitative agreementforthedetonationofcopperwiresand nichrome wires. The GSD showed an overshoot in the Mach stem location compared to both Euler and experiments. This is because the GSD code does not account for the flow properties behind the shock front during the generation of the Mach stem. Despite the overshoot, the overall accuracy of GSD methods in shock-shock interaction has been shown. More importantly, the computational cost of the reduced-order model is much lower than the Euler equations (minutes versus days). 79 Chapter 5 Conclusions and Future Directions 5.1 Conclusions A comprehensive study of the interaction of multiple blast waves, which can be used to direct energy towards a target while simultaneously reducing collateral damage away from the target area has been investigated. Two major points can be implied from the numerical results of interaction of multiple simultaneous and non-simultaneous blast waves. First of all, compared to a single strong blast wave, multiple weaker ones can generatehigherpressure region,even thoughtheenergysummed overalltheblastwaves is equal or less than that of the strong one. Secondly, to achieve better-focusing effect, symmetrical placement of the energy sources is crucial. Results have suggested that a closer distance between the center of energy and the target can result in a higher pressure at the target. These findings have great potential to improve the design of many applications that involve shock induced high-pressure regions. However, great potential comes with strict limitations. To turn the potential into reality, many parameters and constraints need to be considered, such as the placement, the initial energy, the timing and the total cost. Extensive optimization work will be included in these decision making processes. Due to the highly iterative nature of optimization algorithms, the transitional computational technique is not desired. Thus, a reduced-order model is developed, involving all theoretical, analytical and numerical analyses. For the theoretical study, transition condition from RR to IR in the interaction of two blast waves with equal strength is investigated. From both previous experimental results and our numerical results, the sonic criterion is chosen to serve for the transition prediction. 80 For the analytical study, an analytical method, which can describe the entire prop- agation region of a single blast wave accurately, has been implemented and tested. Excellent agreement is found when compared to the results from the Euler equations. Together with the sonic criterion, transition condition at the shock-shock interaction can be predicted. In addition, exploration on the number of blast waves to maximize the peakpressure atthefocalregion while keeping the totalenergy asconstant has been performed and a monotonic behavior has been reported. For the numerical study, an approximation theory to predict shock propagation has been investigated, especially on two aspects: the computational power and the accu- racy. Numerical implementation on parallel computers has been conducted to improve thegeneralperformance. Symmetricboundaryconditionshavebeendevelopedtoreduce the problem size for symmetric converging shock simulations. Furthermore, modifica- tion of the theory has been performed, and good agreement has been obtained when compared to both the Euler equations and previous experiments. 5.2 Future directions Although the current work fulfills the general objectives, there are more features that can be improved. TopredictthetransitionconditionfromRRtoIR,numericalresultshaveshowngood agreement with the sonic criterion. However, there is still some discrepancy between previous experimental data and the theory. To the best of the author’s knowledge, there are few experimental studies on the transition condition in the interaction of two blastwaves. Inmost experimental work, viscous effects cannotbeneglected since asolid objectwasinvolved during theinteraction. Thus, hopefullymoreexperimental workwill be done on the shock-shock interaction to truly understand the transition condition. In addition, although this study only focuses on converging shock configurations, the parallel schemes can be easily extended to solve more generalized shock setups by changing the boundary conditions. In the future, parallelization on three-dimensional GSD model can be investigated. 81 Regarding the modification of the GSD theory, similitude curves are required as a guideline to improve the accuracy of the shock expansion process. Under the current study, those curves are generated from individual sources (analytical solutions, experi- mental data and numerical data). For a single problem, multiple curves may be needed depending on the levels of interactions. Generating these curves will take extra compu- tation resources, which is notdesired for the design process. Thus, algebraic expressions that can describe most of the flow situations for different levels of interactions can im- prove the current model dramatically. To develop these expressions, both the modern computational techniques (e.g. non-linear regression) and the physical understanding of shock-shock interactions are required. Moreover, once these general expressions are es- tablished, it is possible that a pure analytical model can be developed without using the numerical theory, which will significantly improve the computational efficiency without losing much accuracy. Lastbutnotleast, theultimategoalistodevelop softwares thatcanserve forvarious types of optimization problems. A flow chart is shown in Figure 5.1 to illustrate the framework of this software. It can be seen that given the input, which can be the initial placement,theenergy,thenumberorthetimingoftheshockwaves, transitionconditions will be predicted by the analytical solver and the transition curve. Then, based on the transition condition, initial conditions for the modified GSD are set up. After the first run, a solution is then compared to an objective function to decide if it is optimized or not. If not, then a new input will be created by changing the design variables, which depend on the problem setup. This process will be repeated until the solution meets the objective. Depending on different optimization algorithms, the optimal solution may be different. Thus, a detailed study on choosing a proper optimization algorithm that can provide a globally optimal solution within an acceptable time is necessary. 82 Input Analytical Prediction Transition Curve Transition Condition Initial Condition for GSD Modified GSD Optimized? N Y Optimial Solution Figure 5.1 Optimization framework. 83 Bibliography [1] Aki, T., Higashino, F.: A numerical study on implosion of polygonally interacting shocks and consecutive explosion in a box. In: Current topics in shock waves; Proc. of the Int. Symp. on Shock Waves and Shock Tubes, 17th, Bethlehem, PA, July 17-21, (A91-40576 17-34), 167–172, AIP New York (1990) [2] Anderson, J. D.: Modern compressible flow: with historical perspective. Vol. 12. New York: McGraw-Hill (1990) [3] Apazidis, N.andLesser, M.B.: Ongenerationandconvergence ofpolygonal-shaped shock waves. Journal of Fluid Mechanics 309, 301–319 (1996) [4] Bach,G.G.,andLee,J.H.S.: Ananalyticalsolutionforblastwaves. AIAAjournal 8(2), 271–275 (1970) [5] Balasubramanian, K., Eliasson, V.: Numerical investigations of the porosity effect on the shock focusing process. Shock Waves 23(6), 583–594 (2013) [6] Banks, J.W., Aslam, T.D.: Richardson extrapolation for linearly degenerate dis- continuities. Journal of Scientific Computing 57(1), 1–18 (2013) [7] Bazhenova, T.V., Gvozdeva, L.G. and Nettleton, M.A.: Unsteady interactions of shock waves. Progress in Aerospace Sciences 21, 249–331 (1984) [8] Best, J. P.: A generalisation of the theory of geometrical shock dynamics. Shock Waves 1(4), 251–273(1991) [9] Betelu, S., Aronson, D.: Focusing of noncircular self-similar shock waves. Physical review letters 87(7), 074501 (2001) [10] Bethe, H. A., Fuchs, K., Hirschfelder, J. O., Magee, J. L., Neumann, J . V. : Blast wave No. LA-2000, Los Alamos Scientific Laboratory. (1958) [11] Ben-Dor,G.,Takayama,K.andKawauchi, T.: ThetransitionfromregulartoMach reflexion and from Mach to regular reflexion in truly non-stationary flows. Journal of Fluid Mechanics 100(01), 147–160 (1980) 84 [12] Ben-Dor, G. and Takayama, K.: Streak camera photography with curved slits for the precise determination of shock wave transition phenomena. Canadian Aeronau- tics and Space Journal 27, 128–134 (1981) [13] Ben-Dor, Gabi.: Shock wave reflection phenomena. Berlin: Springer (2007) [14] Blackmore, John T. : Ernst Mach; his work, life, and influence. University of California Press, Los Angeles (1972) [15] Book, D., L¨ ohner, R.: Simulation and theory of the quatrefoil instability of a converging cylindrical shock. In: Current topics in shock waves; Proc. of the Int. Symp. on Shock Waves and Shock Tubes, 17th, Bethlehem, PA, July 17-21, (A91- 40576 17-34), 149–154. AIP New York (1990) [16] Brode, H.L.: Quick estimates of peak overpressure from two simultaneous blast waves. Tech.rep., Tech.Rep.DNA4503T,DefenseNuclearAgency, AberdeenProv- ing Ground, MD (1977) [17] Brown, D.L., Henshaw, W.D., Quinlan, D.J.: Overture: An object-oriented frame- workforsolving partialdifferentialequations onoverlapping grids. Object Oriented Methods for Interoperable Scientific and Engineering Computing, SIAM, 245–255 (1999) [18] Bryson, A.E. and Gross, R.W.F.: Diffraction of strong shocks by cones, cylinders, and spheres. Journal of Fluid Mechanics 10(1), 1–16 (1961) [19] Chesshire, G., Henshaw, W.: Composite overlapping meshes for the solution of partial differential equations. J. Comput. Phys. 90(1), 1–64 (1990) [20] Demmig, F., Hemmsoth, H.H.: Model computation of converging cylindrical shock waves - initial configurations, propagation, and reflection. In: Current topics in shock waves; Proc. of the Int. Symp. on Shock Waves and Shock Tubes, 17th, Bethlehem, PA, July 17-21, (A91-40576 17-34), 155–160. AIP New York (1990) [21] Dimotakis, P.E., Samtaney, R.: Planar shock cylindrical focusing by a perfect-gas lens. Phys. Fluids 18(3), 031705 (2006) [22] Eliasson, V., Apazidis, N., Tillmark, N., Lesser, M.: Focusing of strong shocks in an annular shock tube. Shock Waves 15, 205–217 (2006) [23] Eliasson V., Apazidis, N., Tillmark, N.: Shaping converging shock waves by means of obstacles. Journal of Visualization 9, 240 (2006) 85 [24] Eliasson, V., Apazidis, N., Tillmark, N.: Controlling the form of strong converging shocks by means of disturbances. Shock Waves 17, 29–42 (2007) [25] Ghosh, D., Constantinescu, E.M. and Brown, J.: Efficient implementation of non- linearcompact schemes onmassively parallelplatforms. SIAM JournalonScientific Computing 37(3), C354–C383 (2015) [26] Glass, I.: Shock Waves and Man. The University Toronto Press, Toronto (1974) [27] Goldstine, H. H., Neumann, J. V.: Blast wave calculation. Communications on Pure and Applied Mathematics 8(2), 327-353 (1955) [28] Guderley, G.: Starke kugelige und zylindrische Verdichtungsst¨ oße in der N¨ ahe des Kugelmittelpunktes bzw. der Zylinderachse. Luftfahrt Forsch. 19, 302–312 (1942) [29] Hakkaki-Fard, A.: Study on the sonic point in unsteady shock reflections via nu- merical flowfield analysis. McGill University (Canada) (2012) [30] Hakkaki-Fard, A. and Timofeev, E.: On numerical techniques for determination of the sonic point in unsteady inviscid shock reflections. International Journal of Aerospace Innovations 4(1), 41–52 (2012) [31] Heilig, W.H.: Diffraction of a shock wave by a cylinder. The Physics of Fluids 12(5), 154–157 (1969) [32] Henderson, L. F. and Lozzi, A.: Experiments on transition of Mach reflection. Journal of Fluid Mechanics 68(01), 139–155 (1975) [33] Henshaw, W. D., Smyth, N. F., and Schwendeman, D. W.: Numerical shock propa- gationusing geometrical shock dynamics. Journal ofFluid Mechanics171, 519–545 (1986) [34] Henshaw, W.D., Schwendeman, D.W.: Parallel computation of three-dimensional flows using overlapping grids with adaptive mesh refinement. Journal of Computa- tional Physics 227(16), 7469–7502 (2008) [35] Higashino, F., Henderson, L.F. and Shimizu, F.: Experiments on the interaction of a pair of cylindrical weak blast waves in air. Shock Waves 1(4), 275–284 (1991) [36] Hikida, S., Needham, C.E.: Low amplitude multiple burst (lamb) model. Tech. rep., S-cubed Final Report, S-CUBED-R-81-5067 (1981) [37] Hornung, H. G. and Kychakoff, G.: Transition from regular to Mach reflexion of shock waves in relaxing gases. Shock Tube and Shock Wave Research, 296–302, (1978) 86 [38] Hornung, H. G., Oertel, H. and Sandeman, R. J.: Transition to Mach reflection of shock waves in steady andpsuedo-steady flow withand withoutrelaxation. Journal of Fluid Mechanics 90, 541–547 (1979) [39] Hosseini, S.H.R., Takayama, K.: Implosion of a spherical shock wave reflected from a spherical wall. Journal of Fluid Mechanics 530, 223–239 (2005) [40] Iosilevskii, G., and Weihs, D.: Speed limits on swimming of fishes and cetaceans. Journal of The Royal Society Interface 5(20), 329–338 (2008) [41] Itoh, S., Okazaki, N. and Itaya, M.: On the transition between regular and Mach reflection in truly non-stationary flows. Journal of Fluid Mechanics 108, 383–400 (1981) [42] Jiang, Z., Takayama, K., Moosad, K.P.B., Onodera, O., Sun, M.: Numerical and experimental study of a micro-blast wave generated by pulsed-laser beam focusing. Shock waves 8(6), 337–349 (1998) [43] Kandula, M., Freeman, R.: On the interaction and coalescence of spherical blast waves. Shock waves 18(1), 21–33 (2008) [44] Keefer, J.H., Reisler, R.E.: Simultaneous and non-simultaneous multiple detona- tions. In: Shock waves and shock tubes; Proc. of the 14th Int. Symp., New South Wales, Australia, 543–552 (1984) [45] Kjellander, M.K., Tillmark, N.T., Apazidis, N.: Experimental determination of self-similarity constant for converging cylindrical shocks. Physics of Fluids 23(11), 116103 (2011) [46] Kjellander, M.K., Tillmark, N.T., Apazidis, N.: Energy concentration by spherical converging shocks generated in a shock tube. Physics of Fluids 24(12), 126103 (2012) [47] Kleine, H., Timofeev, E., Hakkaki-Fard, A. and Skews, B.: The influence of Reynolds number on the triple point trajectories at shock reflection off cylindri- cal surfaces. Journal of Fluid Mechanics 740, 47–60 (2014) [48] Klein, R. I., Budil, K. S., Perry, T. S., Bach, D. R.: The interaction of supernova remnants with interstellar clouds: experiments on the NOVA laser. The Astrophys- ical Journal 583(1), 245 (2003) [49] Knystautas,R.,Lee,B.,Lee,J.: Diagnosticexperiments onconvergingdetonations. The Physics of Fluids 12(5), 165–168 (1969) 87 [50] Li, C.K., Zylstra, A.B., Frenje, J.A., S´ eguin, F.H., Sinenian, N., Petrasso, R.D., Amendt, P.A., Bionta, R., Friedrich, S., Collins, G.W. and Dewald, E.: Observa- tion of strong electromagnetic fields around laser-entrance holes of ignition-scale hohlraums in inertial-confinement fusion experiments at the National Ignition Fa- cility. New Journal of Physics 15(2), 025040 (2013) [51] Liverts, M. and Apazidis, N.: Limiting temperatures of spherical shock wave im- plosion. Physical review letters 116(1), 014501 (2016) [52] Lock, Gary David, and J. M. Dewey.: An experimental investigation of the sonic criterion for transition from regular to Mach reflection of weak shock waves. Ex- periments in Fluids 7(5), 289–292 (1989) [53] Matsuo, M., Ebihara, K., Ohya, Y.: Spectroscopic study ofcylindrically converging shock waves. Journal of applied physics 58(7), 2487–2491 (1985) [54] Mulley, A.G.: Shock-wave lithotripsy The New England Journal of Medicine 314, 845–847 (1986) [55] Needham,C.E.: BlastWaves. ShockWaveandHighPressurePhenomena.Springer- Verlag Berlin Heidelberg (2010) [56] Neemeh, R.A., Ahmad, Z.: Stability and collapsing mechanism of strong and weak converging cylindrical shock waves subjected to external perturbation. In: Shock waves and shock tubes; Proc. of the 15th Int. Symp., Berkeley, CA, 423–430. Stan- ford Univ. Press (1986) [57] Noumir, Y., Le Guilcher, A., Lardjane, N., Monneau, R. and Sarrazin, A.: A fast- marching like algorithm for geometrical shock dynamics. Journal of Computational Physics 284, 206–229 (2015) [58] Ofengeim, D. Kh, and Drikakis, D.: Simulation of blast wave propagation over a cylinder. Shock Waves 7(5), 305–317 (1997) [59] Oshima, K.: Blast waves produced by exploding wire. Aeronautical Research In- stitute, University of Tokyo q‹’x*z v1J26(9), 137 (1960) [60] Perry, R.W., Kantrowitz, A.: The production and stability of converging shock waves. Journal of Applied Physics 22, 878–886 (1951) [61] Polizzi, E. and Sameh, A.H.: A parallel hybrid banded system solver: the SPIKE algorithm. Parallel Computing, 32(2), 177–194 (2006) 88 [62] Polizzi, E.andSameh, A.: SPIKE:Aparallelenvironment forsolvingbandedlinear systems. Computers & Fluids, 36(1), 113–120 (2007) [63] Qiu, S., Eliasson, V.: Interaction and coalescence of multiple simultaneous and non-simultaneous blast waves. Shock Waves 26(3), 287–297 (2015) [64] Qiu, S., Liu, K. and Eliasson, V.: Parallel implementation of geometrical shock dynamics for two dimensional converging shock waves. Computer Physics Commu- nications 207, 186–192 (2016) [65] Ram, O., Geva, M. and Sadot, O.: High spatial and temporal resolution study of shock wave reflection over a coupled convex–concave cylindrical surface. Journal of Fluid Mechanics 768, 219–239 (2015) [66] Roig, R., Glass, I.: Spectroscopic study of combustion-driven implosions. Phys. Fluids 20(10), 1651–1656 (1977) [67] Roy, C.J.: Review of discretization error estimators in scientific computing. AIAA Paper 2010-0126 (2010) [68] Saito, T. and Glass, I.I.: Temperature measurements atan implosion focus. In Pro- ceedingsoftheRoyalSocietyofLondonA:Mathematical,PhysicalandEngineering Sciences 384(1786), 217–231, The Royal Society (1982) [69] Sakurai, A.: On the propagation and structure of the blast wave, I. Journal of the Physical Society of Japan 8(5), 662–669 (1953) [70] Sakurai, A.: On the propagation and structure of the blast wave, II. Journal of the Physical Society of Japan 9(2), 256–266 (1954) [71] Sakurai, A.: Blast Wave Theory. Modern Developments in Fluid Dynamics, edited by M. Holt, Academic Press, New York, 309–375 (1965) [72] Schwendeman, D.W. and Whitham, G.B.: On converging shock waves. In Proceed- ings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 413(1845), 297–311, The Royal Society (1987) [73] Schwendeman, D. W.: A numerical scheme for shock propagation in three dimen- sions. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 416(1850), 179–198, The Royal Society (1988) [74] Schwendeman, D.W.: Numericalshockpropagationinnon-uniformmedia. Journal of Fluid Mechanics 188, 383–410 (1988) 89 [75] Schwendeman, D. W.: A new numerical method for shock wave propagation based on geometrical shock dynamics. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 441(1912), 331–341, The Royal Society (1993) [76] Schwendeman, D.W.: On converging shock waves ofspherical and polyhedral form. Journal of Fluid Mechanics 454, 365–386 (2002) [77] Sedov, L.I.: Propagation of intense (strong) blast waves (in Russian). Prikl. Mat. Mek. 10, 241 (1946) [78] Shao-Lin, L.: Cylindrical shock waves produced by instantaneous energy release. Journal of Applied Physics 25(1), 54–57 (1954) [79] Stone, H.S.: An efficient parallel algorithm for the solution of a tridiagonal linear system of equations. Journal of the ACM (JACM) 20(1), 27–38 (1973) [80] Stowe,NicholasA.: FocusingEnergyUnderwaterThroughOptimizationofaSpher- ical Source Array. Diss. The University of Michigan (2013) [81] Starkenberg, J.K., Benjamin, K.J.: Predicting coalescence of blast waves from se- quentially exploding ammunition stacks. Tech. rep., Army Research Lab Report ARL-TR-645 (1994) [82] Takayama,K.,Kleine,H.,Gr¨ onig,H.: Anexperimentalinvestigationofthestability of converging cylindrical shock waves in air. Experiments in Fluids 5(5), 315–322 (1987) [83] Takayama, K. and Sasaki, M.: Effects of radius of curvature and initial angle on theshock transitionoverconcave andconvex walls. ReportsoftheInstitute ofHigh Speed Mechanics 46, 1–30, Tohoku University (1983) [84] Takayama, K., Onodera, O., Hoshizawa, Y.: Experiments on the stability of con- verging cylindrical shock waves. Journal of Theoretical and Applied Mechanics 32, 305–329 (1984) [85] Taylor,G.: Theformationofablastwavebyaveryintenseexplosion.I.Theoretical discussion. Proceedings ofthe Royal Society ofLondon. Series A Mathematical and Physical Sciences, 159–174 (1950) [86] Timofeev, E., Skews, B.W., Voinovich, P.A. and Takayama, K.: The influence of unsteadiness and three–dimensionality on regular–to–Mach reflection transitions: a high–resolution study. In Proceedings of the 22nd International Symposium on Shock Waves (ed. GJ Ball, R. Hillier & GT Roberts), 1231–1236 (1999) 90 [87] von Neumann, J.: Oblique reflection of shocks. Explosives Research Report 12, Bureau of Ordnance (1943) [88] Watanabe, M., Onodera, O., Takayama, K.: Shock wave focusing in a vertical annular shock tube. Journal of Theoretical and Applied Mechanics 32, 99–104 (1995) [89] Wang,C.,Qiu,S.,andEliasson,V.: Investigationofshockwavefocusinginwaterin a logarithmic spiral duct, Part 1: Weak coupling. Ocean Engineering 102, 174–184 (2014) [90] Wang, H.H.: A parallel method for tridiagonal equations. ACM Transactions on Mathematical Software (TOMS) 7(2), 170–183 (1981) [91] Whitham, G. B.: A new approach to problems of shock dynamics Part I Two- dimensional problems. Journal of Fluid Mechanics 2(02), 145–171(1957) [92] Whitham, G.B.: Linear and nonlinear waves. J. Wiley, New York (1974) [93] Wu, J., Neemeh, R., Ostrowski, P.: Experiments on the stability of converging cylindrical shock waves. AIAA Journal 19(3), 257–258 (1981) [94] Yeghiayan, R.P., Lee, W.N., Walsh, J.P.: Blast and thermal effects of multiple nuclear burst exposure of aircraft in a base-escape mode. Tech. rep., No. KA-TR- 146. Kaman Avidyne Burlington, MA (1977) 91

Abstract (if available)

Abstract Field of shock focusing has been studied for a long time, and its application has been expanded to a lot of areas such as industrial and civil engineering, biomedical, anti-terrorist and aerospace industry. Due to the nonlinear nature of shock waves and their interactions with nonuniform or moving media, solid or porous surfaces or other shock waves, predictions of shock dynamics and shock focusing events are far from trivial. Our current study aims to deepen the understanding of shock wave focusing processes by studying the interaction of multiple cylindrical shock fronts. The study has been performed using a combination of numerical, theoretical and analytical tools. For the numerical part, simulations of multiple point energy sources were performed, and the resulting shock interaction and coalescence behavior were explored. Three to ten energy sources were placed concentrically around the target, and conditions at the target area were monitored and compared to those obtained using a single source. For each simulation, the energy summed over all sources was kept constant, while the radial distances between target and energy sources and the energy source initiation times were varied. Each energy source was modeled as a point source explosion. The resulting shock wave propagation and shock front coalescence were solved using the inviscid Euler equations of gas dynamics on overlapping grids employing a finite difference scheme. Results show that multiple energy sources can be beneficial for creating extreme conditions at the intended target area

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

An experimental study of shock wave attenuation

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

PDF

Shock wave response of in situ iron-based metallic glass matrix composites

PDF

Converging shocks in water and material effects

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Passive rolling and flapping dynamics of a heaving Λ ﬂyer

PDF

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Accuracy and feasibility of combustion studies under engine relevant conditions

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

$On the dynamic fracture behavior of polymeric materials subjected to extreme conditions$

PDF

On the dynamic fracture behavior of polymeric materials subjected to extreme conditions

PDF

Development and applications of a body-force propulsor model for high-fidelity CFD

PDF

The development of a hydraulic-control wave-maker (HCW) for the study of non-linear surface waves

PDF

Shock sensitivity of energetic material and nanometric damage mechanisms in silica glass

PDF

Mechanical behavior and deformation of nanotwinned metallic alloys

PDF

Document understanding based design support (DocUDS) for augmenting engineering design

PDF

Integration of energy-efficient infrastructures and policies in smart grid

PDF

Molecular dynamics simulation study of initial protein unfolding induced by the photo-responsive surfactants, azoTAB

Asset Metadata

Creator Qiu, Shi (author)

Core Title Numerical study of focusing effects generated by the coalescence of multiple shock waves

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Publication Date 01/28/2018

Defense Date 08/08/2017

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag geometrical shock dynamics,numerical study,OAI-PMH Harvest,shock focusing

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Eliasson, Veronica (committee chair), Bermejo-Moreno, Ivan (committee member), Nakano, Aiichiro (committee member)

Creator Email shiq@usc.edu,undeadstone@yahoo.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-419511

Unique identifier UC11265265

Identifier etd-QiuShi-5654.pdf (filename),usctheses-c40-419511 (legacy record id)

Legacy Identifier etd-QiuShi-5654.pdf

Dmrecord 419511

Document Type Dissertation

Rights Qiu, Shi

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA