PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 5/08 K:IProj/Aoc&PrelelRC/DateDue,indd BLAST SIMULATION WITH SHOCK TUBE TESTING AND COMPUTATIONAL FLUID DYNAMICS ANALYSIS By Kai Long A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Mechanical Engineering 2008 ABSTRACT BLAST SIMULATION WITH SHOCK TUBE TESTING AND COMPUTATIONAL FLUID DYNAMICS ANALYSIS By Kai Long This research was focused on using a shock tube to simulate blast waves. To begin with, it was necessary to understand the characteristics of the shock tube. Besides some basic testing, a computational fluid dynamics (CFD) model based on the commercial code FLUENT was established for numerical simulations. The numerical results agreed favorably with the testing results and those from a one-dimensional analysis. The CFD model was then used to explore more advanced uses and possible upgrading of the shock tube. Venting holes in multiple plates were found to be feasible for mitigating shock waves. The effects of water and sand in blast waves were also investigated in this research. Acknowledgments I would like to begin by extending a sincere thanks to my advisor, Dr. Dahsin Liu for his profound insight, constant encouragement, and generous support in all aspects of my professional development. My appreciation also goes to my committee members, Dr. Mei Zhuang and Dr. Alfred Loos, for their valuable advice on my research and career. I would like to thank Dr. Guojin Li for providing the testing facility and the experimental data, and his generous help on this project. I am immensely grateful to this great department for providing resources to make studying and working here an unforgettable experienCe, to the faculty for nurturing me in the field of mechanical engineering, and to my fellow class mates and lab mates for sharing their knowledge and countless helpful discussion session. Lastly, I would like to thank my parents and my brother, without whose support, I would never have been here. iii TABLE OF CONTENTS LIST OF TABLES ................................................................................... vii LIST OF FIGURES ............................................................................... viii CHAPTER 1 INTRODUCTION .................................................................. 1 1 .1 . Background ................................................................................ 1 1.2. Statement of Problem ................................................................... 8 1.3. Organization of Thesis ................................................................. 11 CHAPTER 2 SHOCK TUBE BASED FACILITY FOR BLAST TESTING INTRODUCTION ..................................................................................... 12 2.1. The MSU Shock Tube ................................................................. 12 2.2. Diaphragms .............................................................................. 13 2.3. Piston-assisted Shock Tube ......................................................... 14 2.4. Sensors for Pressure Measurement .............................................. 15 CHAPTER 3 THEORY AND CFD ANALYSIS TOOLS ................................. 16 3.1. Shock Tube Theory .................................................................... 16 3.2. 2D/3D CFD Code FLUENT ........................................................... 23 3.3. 1D CFD Code JIANG ................................................................... 31 3.4. Post Processing ........................................................................ 32 3.5. Summary ................................................................................... 33 CHAPTER 4 CASE STUDIES AND JUSTIFICATIONS ................................. 34 4.1. The MSU Shock Tube .................................................................. 34 4.2. The MSU Piston-assisted Shock Tube ........................................... 47 4.3. The HEK shock tunnel and the T3 shock tunnel.................................56 4.4. Summary ................................................................................... 63 CHAPTER 5 UPGRADING THE MSU SHOCK TUBE .................................. 64 5.1. Driving Gas/Application Gas Combination Studies..............................64 5.2. Geometry Modification Studies 69 5.3. Summary ................................................................................... 72 CHAPTER 6 SIMULATION OF REAL BLAST WAVES ................................. 73 61 Real Blast Waves Simulations” 73 6.2. Shock Wave Venting on Plate381 6.3. Summary ................................................................................... 96 CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS .......................... 97 7.1 . Conclusions ................................................................................. 97 7.2. Recommendations........................................................................99 APPENDIX A CASE STUDIES ............................................................... 100 iv APPENDIX B PROGRAMES110 APPENDIX C SHOCK TUBE FACILITIES..................................................141 BIBLIOGRAPHY ................................................................................. 142 LIST OF TABLES Table 3.1 Compatibility between 20 element types and mesh schemes ........... 25 Table 3.2 Compatibility between 2D element types and mesh schemes ........... 26 Table 4.1 Mesh parameters .................................................................... 39 Table 4.2 Operating conditions of the free piston shock tube ......................... 48 Table 4.3 Geometry parameters of MSU, HEK and T3 shock tubes ................ 57 Table 4.4 Operating conditions of HEK and T3 ............................................ 58 Table A.1 Shock tubes comparison ......................................................... 141 vi LIST OF FIGURES (Images in this thesis are presented in color) Figure 2.1 The shock tube based testing facility .......................................... 13 Figure 3.1 Distribution of pressure and temperature in shock tube at initial time ................................................................................................... 17 Figure 3.2 Distribution of pressure and temperature in shock tube after diaphragm rupture ................................................................................. 18 Figure 3.3 Distribution of pressure and temperature in shock tube after shock wave reflection ..................................................................................... 21 Figure 3.4 Approaches in segregated solver and coupled solver .................... 29 Figure 3.5 CFD analysis using FLUENT .................................................... 32 Figure 4.1 Initial conditions ..................................................................... 35 Figure 4.2 Mesh for the simple shock tube ................................................. 36 Figure 4.3 x-t Characteristics from the CFD simulation .................................. 36 Figure 4.4 Simulated and measured P5 ..................................................... 37 Figure 4.5 Grid dependence of P5 ............................................................. 40 Figure 4.6 Numerical schemes dependence................................................42 Figure 4.7 Viscous effect ........................................................................... 43 Figure 4.8 Performance curves ................................................................... 44 Figure 4.9 Three-chamber model46 Figure 4.10 Three-chamber model vs. two-chambers model .............................. 46 Figure 4.11 Mesh for the free piston shock tube .......................................... 48 Figure 4.12 Simulated P549 Figure 4.13 Simulated P651 vii Figure 4.14 Simulated piston motions (condition 1) ............................................. 52 Figure 4.15 Influence from sensor (condition 1).. ................................................ 53 Figure 4.16 Mesh grids considering the sensor ........................................... 54 Figure 4.17 Sensor’s influence on P6 simulation ......................................... 55 Figure 4.18 Schematic diaphragm for HEK and T3 ...................................... 56 Figure 4.19 Modeling of the diaphragm rupture process ............................... 57 Figure 4.20 Simulation results for HEK ...................................................... 60 Figure 4.21 Simulated P6 for T3 .............................................................. 61 Figure 5.1 Influence of Driving gas/Application gas combination on P5 ............ 65 Figure 5.2 Performance curves ............................................................... 66 Figure 5.3 Parabola end ......................................................................... 69 Figure 5.4 Annular shock tube ................................................................ 70 Figure 5.5 Simulated P5 ........................................................................ 71 Figure 6.1 Problem setup for shock waves under water and sand .................. 73 Figure 6.2 Total pressure time history at the blast tube outlet ........................ 75 Figure 6.3 Shock wave-sand .................................................................. 76 Figure 6.4 Shock wave-sand (t=0.98 ms) ................................................... 78 Figure 6.5 Shock wave-water (t=0.98 ms ................................................... 79 Figure 6.6 Case 1 ................................................................................. 81 Figure 6.7 Case 2-2 plates ..................................................................... 82 Figure 6.8 Simulated average pressures on plate 1 ..................................... 84 Figure 6.9 Simulated pressures on right end .............................................. 85 Figure 6.10 Contours of pressure on plate 1 .............................................. 87 viii Figure 6.11 Contours of pressure on plate 2 .............................................. 89 Figure 6.12 Contours of pressure on right end ............................................ 91 Figure 6.13 Average pressure on plate 1 ................................................... 94 Figure 6.14 Average pressure on plate 2 ................................................... 94 Figure 6.15 Pressure at the center of the right end ...................................... 95 Figure A.1 P5 .................................................................................... 100 Figure A.2 P6 .................................................................................... 101 Figure A.3 Problem setup ..................................................................... 102 Figure A.4 Simulation results of case 1 ................................................... 103 Figure A.5 Simulation results of case 2 ................................................... 105 Figure A.6 Pressure distribution along the X-axis ....................................... 106 Figure A.7 Problem setup ..................................................................... 107 Figure A.8 P6 .................................................................................... 108 ix CHAPTER 1 INTRODUCTION 1.1 Background A. Blast wave A blast wave is the pressure and flow resulting from the sudden release of a finite amount of energy in a small concentrated volume, which results a strong shock wave propagating outwards. Geoffrey Ingram Taylor developed a similarity analysis of the idealized point source explosion in his pioneering work in 1950 [1], which reduced the partial differential governing equations to ordinary differential equations using similarity assumptions. The so- called Sedov-Taylor solution gave the quantitative estimation of the outcome from the point explosion, such as the pressure distribution and the temperature distribution. Baker summarized the physical model of the blast wave, scaling laws, curves and tables of compiled blast parameters in his book [2]. The blast wave history at one point can be described by the phenomenological models. Baker summarized four phenomenological models in his book [2]. Brode proposed another two phenomenological models [3], which showed more accurate fit for the experimental data using more fitting parameters. Due to the advancements of computation capability in recent decades, the numerical simulations, Computational Fluid Dynamics (CFD) and Finite Element Analysis (FEA), can provide a cost-effective method to study the blast waves, in comparison to the expensive cost of the real blast experiments. The numerical simulations have been implemented in many researches on blast waves. Chapman [4] utilized the AUTODYN2D CFD code to investigate the blast waves generated in a point-explosion, which solves the one-dimentional spherical propagation. Pehrson [5] numerically investigated the interaction between blast wave and structures by incorporating the CONWEP blast model into the FEA code DYNA2D and DYNA3D, which included the angle of incidence of the blast wave on the surface of the structure. Ofengeim [6] performed viscous and inviscid simulations for the blast wave propagation over a cylinder, in which he used a planar shock wave instead of a real blast wave. Trelat [7] examined the response of a large plane surface to an explosion both experimentally and numerically, and the Hopkinson scaling law was employed for the validation of the approach on a large scale. Kambouchev [8] employed a 1D numerical model to compute the flow field corresponding to blast waves propagating in air and impacting on free-standing plates, which had different incident profiles (uniform planar, exponential and exact blast profiles). Recently, Kato [9] proposed a CFD scheme in non-conservative form, which accepts more flexible forms for additional terms such as chemical reactions and phase changes. Based on ANSYS/LS-DYNA, Xia [10] numerically investigated the propagation of blast wave in solid rock, and found that the presence of water in the rock damage the rock more than without water. B. Shock Tube The blast wave front is actually a strong shock wave. Besides the real explosion trials and the corresponding numerical simulations for the blast waves, the shock tube facility, which utilizes a difference in pressure to generate high-speed and high-enthalpy shock wave, is a possible method for simulating the real blast waves. Although the first shock tube was developed in 1861, it did become general uses for the blast wave studies during World War II [11]. Shock tubes developed into the applications of hypervelocity testing [12, 13, 14]. Hornung [15] discussed the important role of shock tube in the hypervelocity flow studies, which included the ground testing, the computation and flight testing, and showed the great advantage of shock tubes over the conventional nozzle expansion tubes in the ground testing, considering the simplicity of the structure and the strength requirements. Shock tubes were shown to be able to generate a shock speed of 5 km/s in the tube and then expand to 7 km/s. Shock tubes were also utilized in the high temperature chemical and physics research [16]. Recently, shock tubes were developed as an impulse generator to simulate the blast wave in the impact testing of composite materials [17] [18] [19]. A hand held shock tube system, called gene gun as a needle for delivering genes and drug system, was investigated both experimentally and numerically [20] [21] [22]. Besides the standard simple shock tube [23], several kinds of shock tubes have been developed to improve the shock tube performance. The so—called Free Piston Shock Tunnel was a milestone in the shock tube facility development history. It greatly increased the shock tube pressure performance, using a heavy piston as the isentropic compression device of the driver gas, and was recognized to be the only facility on ground to simulate the high enthalpy flow of re-entry conditions [24]. Several Free Piston Shock Tunnels are under operating around the world, such as the T5 [25], HIEST [26] and HYPULSE [27]. The piston motion and pressure performance was carefully examined using laser measurement in experiments [27]. The Planar laser-induced fluorescence (PLIF) technique was employed to measure the temperature field at the nozzle outlet of a Free Piston Driven Tunnel [28]. To explore the possible reason, causing the loss of the plateau pressure in Free Piston Shock Tunnels, the loss mechanism was discussed in [29].To improve the performance of the Free Piston Shock Tunnel, a tuned operation was introduced by carefully ensuring a relative high speed of the piston at diaphragm rupture and ensures a low speed of piston at the impacting of the end wall [30]. Most recently, the detonation shock tunnels [31] [32] [33], which used detonation in the driver section, were developed to achieve high pressure and high temperature in the driver section, and thus improved the shock tube’s performance. The operating and structure of the detonation shock tube were shown to be simpler than the Free Piston Shock Tunnel, except for the danger involved with the explosive mixture and the cost of an ignition source [34]. Some researchers explored the annular shock tube, in which an annular shock wave converged. The converging process was studied numerically in reference [35]. A noticeable pressure increase due to the shock wave converging was observed in the experiments [36] [37] [38]. In [39] [40], an overview of the shock waves, such as the physical description, experiments, theoretical analysis and numerical analysis, was described. The physical model and theoretical analysis for shock tubes were discussed in details in reference [11], which gives the analytical solution and design standards for the simple shock tube. In addition to the experimental method and analytical method, the CFD simulation is a promising alternative with the advantage of low cost. Mundt [41] numerically investigated the nozzle pressure of several free piston shock tunnels using a quasi-one—dimensional Lagrangian code L1D, which solved Euler equations. Petrie and Repar [42] performed two dimensional simulations for a simple shock tube with a diaphragm opening by using the finite-volume code U2DE to solve the Euler equations, in which the shock wave speed exceeded the speeds predicted by 1D theoretical solution. Goozee [43] investigated the vortex mechanism for flow contamination in a simple shock tube using the multi-block CFD code, MB_CNS. Sheng [44] investigated an updated Riemann solver to handle discontinues within a simple shock tube, with consideration of the specific heat ratio change across the shock wave. Liu [45] investigated the effect of an obstacle within the driver section of a simple shock tube using an implicit solver for the Reynolds averaged Navier-Stokes equations, MIFVS. Using a novel CFD scheme for the solution of the Euler equation and incorporating with non-ideal thermodynamics (liquid-vapor), Guardone [46] investigated the non-classical three-dimensional flow within a simple shock tube, with a consideration of incomplete diaphragm rupture, a 19% lower shock wave pressure was found than that from a complete diaphragm rupture. Kaneko [47] incorporated the possible reactions between the gas components to the compressible Navier-Stokes equations, and numerically investigated the thermal and chemical nonequilibrium in a simple shock tube. Mizoguchi [48] numerically investigated the diaphragm rupture in a Free Piston Shock Tunnel using a 1D model. In addition to the CFD approaches above, the commercial CFD code FLUENT was widely used by the shock wave researchers [20, 21, 49]. Using FLUENT, Liu [50] numerically investigated the gas-particle interaction within a simple shock tube. As part of the efforts to mitigate the damage from blast waves, some attempts to attenuate the shock wave were investigated numerically. Schwer [51] numerically investigated the mechanism of the shock wave attenuation using water mist, where a two-step reaction model was introduced at the reaction-front. Andreopoulos [49] numerically studied the shock wave attenuation using metallic grids. The numerical results showed that a 50% increase in solidity might attenuate the impinging shock by more than 200%. The numerical simulation for the detonation shock tubes still remained a challenging problem, due to the difficulty of the physical modeling of the deflagration-to-detonation transition [33, 52]. 1.2 Statement of Problem Blast threats from traditional explosive devices have attracted more and more concerns due to the extreme and often catastrophic nature of these kinds of threats. These blast threats include the accidental explosions such as the fuel-oxygen detonation in the oil and coal industries, and the intentional explosions such as roadside bombs in the battlefield. The need for protective equipment against these kinds of threats is important. The fiber-reinforced composites show a great potential in this kind of application considering the weight and capability against blast waves. As part of the efforts to mitigate the blast waves, a deep understanding of the blast waves is needed, which requires detailed theoretical, experimental and numerical studies. Shock tube facility, which is easier to measure, control and repeat in comparison with the real blast trial, is considered to be an ideal experimental facility to simulate the real blast. A double mode shock tube facility, standard simple shock tube and piston assisted shock tube, has been developed at Michigan State University as an impact testing facility [17] [18]. The strong shock wave generated by this facility is used to simulate the blast wave, and hence to impact the testing materials. A couple of experiments show this facility can generate a very high pressure impulse, which is more than 100 MPa, and a high velocity of piston above 100 m/s is achieved. All these severe conditions require careful operation. The detailed theoretical analysis and CFD simulation can predict the details of the flow phenomena, such as the generated pressure impulse, temperature and the piston velocity, under proper modeling, and thus to ensure the safety of the experiments. To simulate the real blast wave, the higher the pressure impulse generated in the shock tube, the better. Upon the careful examination of the simulation results, possible upgrading method to improve the pressure performance can also be achieved. On the other hand, the possible mitigation of the blast wave can be explored using the CFD analysis. Explosions can happen under water or in sands in the real applications, CFD analysis provides cost-effective tool to simulate the blast wave under these scenanos. Another important part of this research is the dynamic response of the composite materials under blast wave loading. The numerical modeling and simulation for the dynamic behavior requires strong coupling of the blast wave and the material. One coupling method is using the phenomenological models [2,3] to provide the blast wave input to the material modeling. The other method, which is more accurate, is to construct a proper CFD model to simulate the blast wave, and thus to use the CFD simulation output as the impulse input in the simulation for the dynamic behavior of the testing materials. To meet the project requirements above, the objectives of this project are: 5" 9° . Construct the CFD models for the shock tube facility implemented with the CFD codes. Justify the CFD models for pressure wave analysis. Identify the properties and physics of pressure waves Explore the possible upgrade of the MSU shock tube. Investigate the possible mitigation of pressure waves. 1.3 Organization of Thesis Chapter 1 Introduction Chapter 2 Shock Tube Based Testing Facility Chapter 3 Shock Tube Theory and CFD Analysis Tools Chapter 4 Case Studies for Justifications Chapter 5 Upgrading the Shock Tube Chapter 6 Mitigation of the Strong Shock Waves Chapter 7 Conclusions II CHAPTER 2 SHOCK TUBE BASED FACILITY FOR BLAST TESTING 2.1 The MSU Shock Tube A shock tube was developed at MSU as a pressure generator for impact testing [18]. Figure 2.1 shows the schematic diagram of the three- section shock tube and accessories for blast testing. The shock tube was made of a steel alloy containing chromium and manganese and had high strength and high temperature resistance. The circular shock tube has a constant cross-sectional area. The outer . diameter of the tube is 120mm while the inner diameter 80mm. The total length of the shock tube is 6.25m. The far-left section, which is also called the high-pressure chamber, stores a relatively high-pressure gas, and has a length of 2m. The right section, which is called the low-pressure chamber, is 4m long. It consists of two 2m parts joined by flanges. The third section, called the intermediate-pressure chamber, is located between the high-pressure chamber and the low-pressure chamber and has a length of 0.25 m. For high-pressure testing, a piston can be inserted to the low-pressure chamber to increase the pressure significantly. Besides, a fourth section can be added to the end of the shock tube. The fourth section, which is called the blast tube, is 0.1537m long and has a inner diameter of 0.0128 m. specime holder blast tube high-pressure low-pressure chamber chamber [:6 dia phragm Figure 2.1 The shock tube based testing facility [20] 2.2 Diaphragms Two diaphragms were used to isolate the intermediate-pressure chamber from the high-pressure chamber and the low-pressure chamber, one on each side of the intermediate-pressure chamber, as shown in Figure 2.1. In the operation of the shock tube, the intermediate-pressure chamber is filled with a pressure approximately equal to the average of the high pressure and the low pressure. For example, if the pressure in the high-pressure chamber is set for 1.72MPa and that in the low-pressure chamber is 0.27MPa, the pressure in the intermediate-pressure chamber will be set as 1MPa. When all chambers reach designated pressure levels, the gas in the intermediate- pressure chamber is vented out. The two diaphragms separating the intermediate-pressure chamber from the high-pressure chamber and the low- pressure chamber are machined with grooves with well-calculated depth to allow the first diaphragm to burst at the pressure level equals to the high I3 pressure, and the second diaphragm to burst at the pressure level equal to the difference between the high pressure and the low pressure, respectively. Accordingly, an instantaneous shock wave can be generated. The use of the intermediate-pressure chambers, i.e. the double diaphragms, instead of single diaphragm, is aimed at a higher pressure ratio as mentioned in reference [46]. The diaphragms need to be carefully prepared. In this study, they were made of aluminum 6061 with dimensions of 150mm X 150mm and a thickness of 2mm. Two diagonal groves were machined with a carefully calculated depth to warrant the burst of the diaphragms under the designated pressures. Once the diaphragms burst, the high-pressure gas flows rapidly into the low- pressure chamber and moves with the low-pressure gas toward the right end of the shock tube. The diaphragms are designed for a complete burst to establish consistence for experiments. However, due to imperfection in machining, the diaphragm may not burst completely. As a result, the pressure rise caused by the shock wave and the speed of the shock wave will be smaller than that of a complete burst [48]. This makes the pressure rise less predictable by the analytical shock tube theory, which assumes a complete diaphragm burst. 2.3 Piston-assisted Shock Tube In order to operate the shock tube at higher pressure levels, a piston can be inserted into the low-pressure chamber. Besides, a blast tube was added to the end of the low-pressure chamber. For simulating a blast wave, an additional diaphragm is added to the front of the blast tube for creating a blunt pressure front. Figure 2.1 also shows the piston, the blast tube and the diaphragm. During the stage of filling gases, the initial pressure in the high- pressure chamber is designated as P4 while that in the low-pressure chamber P1. After the burst of the double diaphragms, the low pressure gas is compressed by the piston, and the pressure is increased from P1 to P5. The diaphragm for the blast tube bursts when P5 is higher than the rupture pressure, resulting in a blunt shock wave — the blast wave. 2.4 Sensors for Pressure Measurement Three pressure sensors, installed at the high-pressure chamber, intermediate-pressure chamber and low-pressure chamber, are required to monitor the corresponding pressure levels. To measure the pressure P5 generated in the shock tube without piston, i.e. the simple shock tube, and the compressed pressure P5 generated in the shock tube with a piston, i.e. piston- assisted shock tube, a pressure sensor should placed in the end of the low- pressure chamber. To measure the pressure of the blast wave P6, a pressure sensor should be placed at the outlet of the blast tube. Both P5 and P6 are important parameters for blast simulations. CHAPTER 3 THEORY AND CFD ANALYSIS TOOLS Shock tubes utilize difference in pressure to generate high-enthalpy and high-speed flows. This chapter outlines the theoretical analysis of the shock tube flow. In addition, two computational fluid dynamics (CFD) codes are used to investigate the waves in the MSU shock tube. They are one-dimensional (1D) code JIANG and two-dimensionaI/three-dimensional (2D/3D) code FLUENT. For simplicity, the three-section shock tube is simulated by two sections only. And the bursting process of the diaphragm is also neglected. 3.1 Shock Tube Theory The governing equations for the shock tube performance can be derived under several simplifying assumptions: (1) the flow is of one dimensional, (2) there is no viscosity or heat transfer involved, (3) the bursting process is neglected, (4) the expansion wave is isentropic, and (5) the gases are ideal gases and there is no loss involved. Three performance parameters can be derived from the idealized shock tube theory, which include the pressure performance%, the temperature performanceg—f and the Mach number of the shock wave M s . There are two initial parameters which determine the shock tube performance: the pressure ratio between the high pressure and low P pressure /P1 , and the sound speed ratio between the high pressure gas I6 a (the driving gas) and the low pressure (the application gas) /al . The initial conditions are shown in Figure 3.1. 4 1 I P4 Pressure P1 Position 4 l Tern re peratu T4=T1 Position ' Flgure 3.1 Dlstrlbutlon of pressure and temperature In shock tube at InItIaI tIme As shown in Figure 3.2, upon the rupture of the diaphragm, a right- travelling shock wave, followed by the driving - application gas interface, which is the interface between the driving gas and the application gas, and a set of expansion waves with a fan-shape are formed. The pressure and temperature I7 N .5 l we Expansion wave Contact face Shock wave "\ Pressure P2=P3 P1 Position 7 A T2 Temperature T4 T1 \ T3 Position ' Flgure 3.2 DIstrIbutIon of pressure and temperature . In shock tube after diaphragm rupture . . in the shock tube can be divided into four categories as shown In Figure 3.2. The shock wave forms and travels to the right: it increases the pressure in region 1 from P1 to P2. The expansion wave travels to the left and lowers the pressure from P4 to P3. It also lowers the temperature from T4 to T3. The driving-gas/application gas interface travels to the right and has a speed slower than the shock wave. The pressure across the driving gas/application gas interface is identical (P2=P3), while the temperature drops from T2 to T3. 18 The equations addressing the changes of pressure and temperature can be derived from two basic relations; one across the expansion wave and the other across the shock wave. The equations to follow are from reference [11]. The isentropic relation across the expansion wave is 2i P4 (1,, 74" P1 a3 (3.1) wherea is the speed of sound and 7 is the specific heat ratio. The relations across the normal shock wave are: P21=1+%_7—1(M52 —1) T = [27M —(7. —1)I- [(7. —1)M_f +2] (32> 21 (y, +1)2Ms2 where M s is the Mach number of the shock wave. Using these two basic relations, the dependence of the Mach number of the shock wave M s on the initial condition% can be obtained as follows: 2 ' — 7 -1 Ell = 1+ 7,1 (I‘ls2 _1) 1— 7,4 1a’l4 Ms —L 4 3!, +1 71, +1 M (3-3) S —— ‘ Using Equation 3.1, 3.2 and 3.3, the dependence of the pressure ratio across the normal shock wave (g?) on the initial pressure ratio(—:—:) can be obtained with Equation 3.4: a1P2, +1 (3.4) 6 3 5 2 5 : i—b 4— T r I Expansion waves Contact face Shock wave P5 / P2=P3 P6 T5 T2 T3 _—;r_§—-—-—-'/ ‘ ' Flgure 3.3 Distribution of pressure and temperature In shock tube after shock wave reflection As shown in Figure 3.3, the shock wave reflects at the right end and travels to the opposite end. The reflected shock wave raises the pressure from P2 to P5. The expansion wave reflects at the left end, and the reflect expansion wave decreases the pressure from P3 to P6. To derive the major performance (gs-J and (g), the basic relations across the reflected shock wave are shown in Equation 3.5: 2] _ (05+2)P2| —1 p _ 52 P2l+a _ (a + P52 ”’52 (3'5) 52 aPSZ+1 How the performance parameters (1;) and (If?) depends on the initial condition[-II—:%) can be obtained by combining Equation 3.4 and Equation 3.5 with Equation 3.6. P51 T51 P52‘P21 T52 'T21 (36) 22 3.2 2DI3D CFD code FLUENT FLUENT is a widely used CFD analysis package, which is designed for numerical simulations of various kinds of flows. It has been shown to be able to simulate the flow in shock tubes adequately [49] [20] [21] [50]. The package includes two major components: the geometry and mesh preprocessor GAMBIT 2.3.16 and the solver FLUENT 6.2.16. 3.2.1 The Geometry and Mesh Preprocessor GAMBIT 2.3.16 GAMBIT 2.3.16 [53] is a software package designed to build and mesh models for computational fluid dynamics (CFD) and other scientific applications. Using its Graphical User Interface (GUI), the geometry model and the followed mesh grids can be generated in GAMBIT. A. Geometry Generating a. Geometry Types: 1. Real Geometry - entities characterized by a direct definition of their geometry. For example, a vertex defined by its coordinates (0,0,0) 2. Virtual Geometry - entities characterized ONLY by an indirect definition, such as a reference to another entity. For example, a vertex is defined as the mid-point of an edge 3. Faceted Geometry - entities characterized by an indirect definition with respect to an underlying grid. For example, a vertex is defined as the corner of a mesh element 23 b. Geometrical Topology: - Vertex - Edge (2 or more vertices) . Face (3 or more edges) - Volume (4 or more faces) c. Geometry Generating Approaches: Bottom-up approach: generates low dimensional entities and builds on top of them higher dimensional entities. Top-bottom approach: generates upper dimensional entities and uses Boolean operation (unite, subtract and intersect) to define the other entities. Importing Geometries: the realistic geometries, which can be modeled in the CAD software, sometimes are very complicated to be generated from the simple geometries in the CFD mesh generator. Translation between CAD and CFD systems is a major bottleneck in mesh generation. GAMBIT provides a bridge to connect the CAD and CFD systems. GAMBIT can import geometry from virtually any CAD/CAE software in Parasolid, ACIS, STEP, IGES, and native CATIA V4N5 formats. Tolerant modeling and healing capabilities automatically provide connected solid geometry during import. Newly added CAD connections for SolidWorks, Pro/ENGINEER, UGS’NXT, Autodesk, and 24 Inventor use native surface, connectivity and topology information to translate geometry models. Element Quadrilateral Triangular Quadrilateral / Scheme Triangular Map X X Submap X Pave X X X Tri Primitive X Map Split X Submap Split X Wedge Primitive X Table 3.1 Compatibility between 20 element types and mesh schemes B. Grid Generating Upon the generation of the geometrical model, the mesh grids can be generated using GAMBIT. Table 3.1 and Table 3.2 shows the available 2D and 3D element types and the compatible mesh schemes in GAMBIT. The basic steps to generate mesh grids upon the generated geometry are summarized below: a. 1D - edges: Define the number of grid points and the distribution of grid points. b. 20 - faces: Define the number of elements (using the predefined edge mesh in the last step or define a uniform face element size). Choose the mesh scheme based on the element types and the constraints from the edge meshing. 25 Element Hexahedron HexahedronNVedge Tetrahedron Scheme /Hybrid Map X Submap X Tet Primitive X X Cooper X Stairstep X Tgrid X Hex Core X Table 3.2 Compatibility between 3D element types and mesh schemes c. 30 - volumes: Define the number of elements (using the predefined face mesh in last step or define a uniform volume element size). Choose the mesh scheme based on the element types and the constraints from the edge meshing. d. Special tools: Use block meshes, non-conformal meshes to generate the structured mesh to the best extent possible. Use boundary layers meshes to refine the boundary layer. e. Examine the mesh qualities: Use the ‘Examine The Mesh’ toolpad to examine the mesh qualities and make the necessary modifications on the gnds. f. Specify the zone types: Define the boundary conditions and the continuum types. 9. Save the jobs and export the examined mesh 26 3.2.2 Solver FLUENT 6.2.16 [54] The solver operates through an interactive, menu-driven, Graphical User Interface (GUI). In addition to the GUI, FLUENT consists of Text User Interface (TUI), which executes the customized text command. A. Governing Equations solved by FLUENT 6.2.16 FLUENT numerically solves the conservation of mass, momentum and energy along with the species transport equation and the equation of state. Neglecting the viscous effects, these governing equations are reduced to the Euler equations. a. Conservation of Mass: 28'?- + V - (1)-v): Sm (37) where Sm is the mass source. b. Conservation of Momentum: % ;)+ V - (oft): —Vp + V . CF p} + I? (3.8) where : is the stress tensor, ,0} and I; are the gravitational body force and external body force. 0. Conservation of Energy: at __a_(pE)+V.(§(pE+ p))= V-[keflVT—Zhjjj +G-3)]+Sh (39) 27 where kc}, is the effective conductivity, 7 ,- is the diffusion flux of species j. S, is the source term. d. Species Transport Equations: %( Y,)+V-(p§i§)=—V-7.+R,+S,( Where Y, is the local mass fraction of species i, R, is the not rate of 3.10) production of species i by chemical reaction and S, is the source. .7.- is the diffusion flux of the 1’" species. Equation of state for ideal gas: P = PRT (3.11) 8. Numerical Schemes In FLUENT Two numerical methods (pressure-based segregated solver and density-based coupled solver) are implemented in FLUENT. Using either method, the governing equations (Equations 3.7-3.11) are numerically solved. In both cases a control-volume-based technique is used that consists of: a. Division of the domain into discrete control volumes using a computational grid. b. Integration of the governing equations on the individual control volumes to construct algebraic equations for the discrete dependent variables such as velocities, pressure, temperature, and conserved scalar. 28 Update properties. 1 Solve momentum equations. 1 Solve Pressure-correction (continuity) equation. Update pressure, face mass flow rate I Solve energy. species, turbulence, and other scalar equafions l Converged? a Stop (a) Segregated solver Update properties. 1 Solve continuity, momentum, energy, and species equations simultaneously. l Solve turbulence and other scalar equations. V Converged? v Stop (b) Couples solver Figure 3.4 Approaches In segregated solver and coupled solver 29 c. Linearization of the discretized equations and solution of the resultant linear equation system to yield updated values of the dependent variables. The two numerical methods employ a similar discretization process (finite-volume), but the segregated solver solves the governing equations (Equations 3.3-3.7) sequentially, while the coupled solver solves the governing equations simultaneously, as shown in Figure 3.4. C. User Defined Functions in FLUENT Since the FLUENT solver is a general-purpose but cannot anticipate all needs, there is a need to develop and incorporate new models to the FLUENT solver. User-defined functions, or UDF, are C routines that can be dynamically loaded with the FLUENT solver to enhance the standard features of the code. In summary, UDF can: a. return a value. b. modify an argument. 0. return a value and modify an argument. (I. modify a FLUENT variable (not passed as an argument). e. write information to (or read information from) a case or data file. 30 3.3 1D CFD code JIANG $5— + V 1,017): 0 % ;)+V - 1,0175): -Vp + +p§ + 75 §;(pE)+V-(iz(pE+ p))= V'(keflVT—Zhj7j)+sh (3.12) Neglecting the viscous effect, the governing equations Equation 37-39 are reduced to the Euler Equations (Equation 3.12). a .. —’0 + V - (pv) = O at % ;)+v.(o;;)=_vp Sat-(pEHV . (;(pE+ p)).—. o (3.13) Neglecting the body force, the heat transfer and source in Equation 3.12, the simplified 1D Euler Equations (Equations 3.13) and the Equation of State (Equation 3.11) are numerically solved. A nonoscillatory, nonfree and dissipative (NND) numerical scheme [55, 56], which has no free parameters and has good stability characteristics and convergence, is introduced in this code. The computation cost is significantly reduced by using the NND scheme compared with some other high-resolution Total Variation Diminishing (TVD) scheme. 31 Details of the code are given in Appendix 83. 3.4 Post Processing EXCEL, MATLAB and TECPLOT were used as the post-processing tools for the data analysis and flow field visualization. k MATLAB) Solver (FLEUNT) F Pro-Processing / { -Transport Equations Define Mesh Mass . Modeling ‘” Generator Species mass fraction .Shtlflver Settings Goals (GAMBIT) Phase volume fraction umencal L Momentum Schemes Energy °Setting solution ~Eq.iations of state controls -Supporting physical models oRevise geometrical, I I numeric” 0' phi/Sic“ 'Material Properties Physical Models - ~Boundary - Turbulence NO conditions 0 Combustion , * ~Initial Condtions - Mdtiphase _ Ye oExamlne results , Moving Zones ~ (EXCEL: - Moving Mesh ' " TECPLOT, J Flgure 3.5 CFD analysis uslng FLUENT 32 3.5 Summary The governing equations (Equation. 3.3, 3.4 and 3.6) for the P5 T5 performance parameters, (E) (— T1)and(M3) were obtained through the theoretical analysis. The basic steps to perform the CFD analysis using FLUENT (Figure 3.5) are summarized below: 1. Define the modeling goals. 2. Create the model geometry and grid using GAMBIT. 3. Set up the solver and physical models. 4. Compute and monitor the solution. 5. Examine and save the results. 6. Consider revisions to the numerical or physical model parameters. 33 CHAPTER 4 CASE STUDIES AND JUSTIFICATIONS This chapter justifies the CFD simulation results for the CFD models, which describe the MSU shock tube, the HEK shock tunnel and the T3 shock tunnel, by comparing the simulation results to the measured data. The 2D and 3D model comparison, grid dependence, numerical schemes comparison and viscous effects are investigated for the MSU shock tube. The performance curves for the MSU shock tube are then generated by both analytical method and CFD simulations. The numerical models for the other two shock tunnels were hence generated and justified based on the study for the MSU shock tube. Finally several improvements for the numerical models are explored. 4.1 The MSU Shock Tube The MSU shock tube can be operated with or without a piston. The inner boundaries were set to be the simulation boundaries as shown in Figure 4.1 along with the initial conditions. The 2D and 3D mesh grids were generated with GAMBIT as shown in Figure 4.2. Some simplifying assumptions were made in the simulations. For example, the gases were assumed to be ideal gases. A single diaphragm, instead of two diaphragms, was used. A complete and instant rupture of diaphragm was assumed. And the gas leakage or loss was neglected. 34 Diaphragm High-pressure Low-pressure chamber chamber 80 P4=1.72 MPa (250 P1=0.27 MPa (40 mm Psi) Psi) 2m 4m Condition 1 Diaphragm High-pressure Low-pressure chamber chamber 80 P4=4.13 MPa (coo P1=0.55 MPa (80 mm Psi) Psi) 2m 4m Condltion 2 Figure 4.1 Initial conditions Under the initial condition 1, given in Figure 4.1, a shock wave with a speed of 487.8 m/s (1.42 Mach) was generated. In Figure 4.3, the x-t diaphragm was plotted to describe the shock wave characteristics in the shock tube. Contours of simulated pressure were plotted in the position-time plane. From Figure 4.3, the shock wave was initially generated around x=-1 m, where the diaphragm was also located. The shock wave traveled along the tube to the right end before reflecting back at x=3 m around t=8 ms. The driving gas - application gas interface, where the driving gas and application gas began mixing, is depicted by the dark line in Figure 4.3. The gas interface traveled slower than the incident shock wave. It interacted with the reflected shock wave at x=1.65m at around t=12.24 ms and was slowed down by the reflected shock wave. The reflected shock wave transmitted through the interface, and was strengthened by the gas interface. At the same 35 time, a reflected compression wave was generated at the interface, and traveled back to the right end. (a) 20 mesh (b) 30 mesh Figure 4.2 Mesh for the simple shock tube 35 0) O 7,6 Pressure(MPa) ”shirt, 1 .7 25 eds/7 1 .6 °°k . . i -5 "’6 Transmitted expansron 1 .4 a waves i .3 20 i .2 E l .1 F 4 E .9 .- 15 l.8 '- t; 25 10 [g 5 0 o X(m) Figure 4.3 x-t Characteristlcs from the CFD simulation On the other side, the expansion waves, which were generated at =-1 m initially, traveled to the left and then reflected at the left end (x=-3 m). The reflected expansion waves traveled to the right and interacted with the 36 2 P4=1.72MPa (250 psi) , P1=0.27MPa (40psi) 1.8 Air-Air 1.6 ~ A 1.4 0‘.“ 1.2 g 1 If 0_3 Measured 0.6 , —JIANG-1D 0.4 + ——FLUENT-2D 0.2 ; —FLUENT-3D 0 + + i i 0 10 20 30 Time (ms) (a) Condition 1 P4=4.1 3MPa(600psi), P1 =0.55MPa(80psi) 4'5 T Air-Air 4 -2 ’8 o. é —Measured “- —JIANG-1D —FLUENT-2D —FLUENT-SD 0 i + i 0 10 30 Time(ms) (a) Condition 2 Figure 4.4 Simulated and measured P5 reflected shock wave at x=0.91 m at approximately t=15 ms. This interaction increased the transmitted shock wave speed and the shock wave strength, as depicted in Figure 4.3 that the slope of the shock wave was reduced after the interaction with the expansion waves. Also shown in Figure 4.3, the pressure at the right end (x=3 m) was first increased by the incident shock wave and the reflected shock wave at time t =8 ms, and then was increased to the maximum pressure by the compression wave reflected from the gas interface at time t=15 ms. Consequently it was attenuated by the reflected expansion wave. 37 Figure 4.4 shows the standard pressure history of P5 at the end of low pressure chamber under condition 1 and condition 2 as shown in Figure 4.1. The 1D, 2D and 3D CFD simulated pressures were compared to the measured pressure. The overall agreement for the wave shape between the simulated pressure and measured pressure was fair. The simulation, however, seems to over - predict the peak pressure and the plateau pressure by about 0.2 MPa in condition 1 and by 0.5 MPa in condition 2. These discrepancies were due to the loss in the experiment, such as the tube sealing loss, the diaphragm opening loss and the heat loss, which were not included in the simulations. The arriving sequences of shock wave and expansion wave depicted in Figure 4.3 can also be found in Figure 4.4. In Figure 4.4, P5, the pressure at the right end, is increased by the incident shock wave and the reflected shock wave, which refers to the peak pressure and the first plateau pressure. The pressure is further increased by the compression wave due to the driving - application gas interface, and becomes the second plateau before it is attenuated by the reflected expansion wave. The 2D and 30 CFD modelings using FLUENT were able to catch the pressure oscillation at the arriving of the shock wave. It was observed in the measured data. The interaction between the reflected shock wave and the gas interface was modeled by the 2D and 3D models, and a reflected compression wave from the interaction was predicted, which agreed with the experiment. In comparison, the 1D model predicted a reflected shock wave from the shock 38 wave-interface interaction. Overall, the 2D and SD models predicted a better pressure increase rate than the 1D model. Under condition 1, the pressure-time history of P5, which is the pressure at the end of the low pressure chamber, was used to evaluate the grid dependence, numerical schemes dependence and the viscous effects since P5 is the most concerned parameter in the shock tube study. Then the performance curves were generated by both the analytical method and CFD analysis. A. Grid Dependence The grid dependence was evaluated by comparing three different mesh densities for both the 2D model and 3D model, which are shown in Table 4.1. Number of cells Mesh size (m2) 800 6.23 e-4 4000 1.21 e-4 8000 6.05 e-4 a)2D mesh Number of Mesh size (m3) cells 2364 2.02 e-5 8382 3.6 e-6 16848 1.89 e-6 b)3D mesh Table 4.1 Mesh parameters 39 P4=1.72MPa (250 psi) , P1 =0.27MPa (40psi), 20 2 __ Air-Air 1.8 + 1.6 ~— 1.4 ~ 312 -- s 1 —— E 03 _, —Measured 0.6 _ —800 cells ()4 . —4000 cells 02 W —8000 cells 0 l i l l 0 10 20 30 Time (ms) (a) 2D grid dependence 2 P4=1.72ll/Pa(250psi), P1=0.27IVPa(40psi), 30 1.8 } Air-Air 1.6 1.4 - ’8 1.2 — a. 52; 1 _ ——-Measured a. 0.8 ~ 05 _ —2364 cells 0.4 - —8382 cells 0.2 - —16848 cells 0 +— . i i . r 0 5 10 15 20 25 30 Time (ms) (b) SD grid dependence Figure 4.5 Grld dependence of P5 In Figure 4.5, the grid dependence in 2D and 3D models is examined. For the 2D model, 800 grids were not able to predict the peak pressure of P5 and it gives a smooth incline at the shock front. The second plateau pressure duration was less predicted by the 800 grids. The influence of grids decreased 40 from 4000 cells to 8000 cells. The 3D model showed less dependence on the grids and the influence of grids decreased from 8382 cells to 16838 cells. The 4000 grids in the 2D model, which had the cell size of 1.21 10'4 m2, and the 16838 grids in the 3D model, which had the cell size of 3.6 10'6 m3, were considered to be suitable for this case study in terms of the grid dependence and the computation load. B. Numerical Scheme Dependence Several numerical schemes, which included the coupled-implicit, coupled-explicit and segregated-implicit schemes, were evaluated for both the 2D model and the 30 model. For the 2D model shown in Figure 4.6(a), the coupled-explicit schemes predicted more accurate duration of the first plateau pressure and pressure decay than the coupled-implicit linearization. This meant that the coupled-explicit scheme predicted a better expansion wave speed. The coupled-implicit scheme was able to predict the peak pressure at the shock front as observed in the experiments. The segregated-implicit scheme was less effective in catching the steep shock front than the other two schemes. For the 3D model shown in Figure 4.6b, less numerical scheme dependence was exhibited. The coupled-explicit and coupled-implicit schemes predicted similar results except for the peak pressure at the shock front. The segregated-implicit scheme did not predict an obvious second plateau. The 41 coupled-implicit scheme was considered to be more suitable in this case study in terms of its ability to catch the peak pressure. 2 P4=1.72IVPa (250 p50 , P1=0.27IVPa (40psi), 20, 4000 cells Air-Air 1.5 — To? a. E 1 - to —— Measured 0- . — FLUENT-Coupled-irrpllcit 0-5 “ -—FLUEI\lT-Coupled-expliclt I -——- FLUBdT—Segregated-irrplicit JIANG 0 l I l 0 10 Time (ms) 20 30 (a) 2D numerical schemes comparison P4=1.72NPa(250psi), PI=0.27NPa(40psi), 3D, 8382 cells 2 _ Air-Air 1.8 - 1.6 - ’ 1.4 - g 1.2 - E 1 - —Measured a“ 0.8 — —FLUB\IT—coupled-irrplicit 0.6 ~ —FLUEI\IT-coupled-explicit 0'4 ‘ — FLUBdT-segregated-irrplicit 0.2 JIANG 0 7" | l I | V I ' T ‘ V 'I 0 5 10 15 20 25 30 Time (ms) (b) 3D numerical schemes comparison Flgure 4.6 Numerical schemes dependence C. Viscous Effects For 2D model, the viscous calculation, based on both laminar model and Spalart Allmaras turbulence model, and the inviscid calculation were 42 performed to evaluate the viscous effect. The simulation result depicted in Figure 4.7 showed that they predicted almost identical results. It was concluded that the viscous effect did not play an important role for the P5 simulation because of the high speed of the shock wave, which was similar to the result in reference [54]. 2 P4=1.72lVPa (250 psi) , P1=0.27NPa (40psi), 20, 4000 cells Air-Air 1.5 + E g 1 s —- Measured g —— FLUENT hvisid 05 a- — FLUBIT Lam‘nar t ' _ . . . . ‘ — FLUExIT SA-Turbulent 0 l i 1* i i 1' o 5 10 15 20 25 30 Time (ms) Figure 4.7 Viscous effect D. Performance Calculations Using both analytical method (Equations (3.3), (3.4) and (3.6)) and CFD analysis, several performance curves were generated for the experimental studies. Figures 4.8(a), 4.8(b) and 4.9(c) show the dependence of P5max, T5max and the Mach number M5 on the initial operating conditions, respectively. The 2D coupled-explicit scheme gave an overall better agreement with the analytical solution than other schemes. The CFD models were justified by comparing the simulated performance data with the analytical solution. P5i’P1 TSmaxfII F’Smaxii'P1 ~ P4IP1 Air-Air 30 . A A 25 "’ A _/~"fl A aft/‘8 20 - g/a/ . 15 - / - /fi/ Analytical 10 _ fa" Ci FLUENT ZD-coupled-implicit q ’ El FLUENT 2D-coupled-explicit 5 _ {A FLUENT SID-coupled-explicit . 1? JIANG code ID 0 1 1 l 1 1 0 20 4|] 60 80 100 PMPI a) P5/P1 T5max/T1~P4IP1 Air-Air 3.5 . . . 3 - _ftt/’a’/: E . . A - - — Analytical (-2 FLUENT 2D-coupled-implicit - El FLUENT 2D-coupled-explicit A. FLUENT SD-coupled-explicit 4 V JIANG code 10 40 60 80 100 P4i’P1 b) Mach number of the shock wave Figure 4.8 Performance curves Mach number ~ PMPI Air-Air 2.5 . s ‘7 ” a I 1'7 fit A 2 f A A -i ffi/é/l Analytical FLUENT2D-coupled-implicit FLUENT ZD-coupled-explicit - FLUENT 3D-coupled—explicit JIANG code 10 m in \ [.33 RW' QPUL Mach number of the shock wave l l 0 20 40 60 80 100 P4IP1 c) T5/T1 A i— Figure 4.8 Performance curves (continue) E. Further Modification on the CFD Model In previous studies, the three-chamber simple shock tube (Figure 2.1) was simplified to a two-chamber model (Figure 4.1) by assuming a single diaphragm. The influence from the intermediate-pressure chamber was also studied by using a three-chamber CFD model. The intermediate-pressure chamber in Figure 2.1 was established and filled with a gas pressure averaging that of the high-pressure chamber and the low pressure chamber. As shown in Figure 4.9(a), a 44.8 mm long intermediate-pressure chamber with a pressure of about 1 MPa was introduced. Furthermore, the incomplete rupture of the diaphragm was studied by introducing a 65 mm hole on the diaphragms, which was measured from the actual diaphragm used in the experiment, as shown in Figure 4.9(b). The comparison between the two- lntermediate-pressu re chamber 44.8 mm long 0.9997 Mlza (145 psi) High-pressu re Low-pressu re chamber chamber P4=1.72 MPa (250 P1=0.27 MPa (40 30 mm Psi) Psi) (a) 3-chambers model Intermediate-pressure chamber 44.8 mm long 0.9997 MP4a (145 pSI) Hhigh-greessure I I I Law-[Erasure 0 am r c am er P4=1.72 MPa (250 P1=0.27 MPa (40 8° mm psi) I I 1350 Incomplete rupture hole of 65 mm (b) 3-chambers with incomplete rupture model Figure 4.9 Three-chamber model chamber model, three-chamber model and three-chamber model with incomplete burst is shown in Figure 4.10. 2 P4=1.72MPa (250 psi), P1=0.27M Pa 1.8 ] (40psi),2o 1.6 d Air-Air 1.4 ~ g 1.2 ~ 2 1 fl Measured E 0.8 ~ =2chambers 0-6 “ Schambers 0.4 — —3 c ha mbe rs /w incomplete 0-2 7 diaphragm rupture o . . 1 o 10 _ 20 30 Time (ms) Figure 4.10 Three-chamber model vs. two-chambers model 46 In Figure 4.10, the simulations from the three-chamber model and the three-chamber model with incomplete burst gave better pressure wave shapes than the two-chamber model. The three-chamber with incomplete burst predicted a lower plateau pressure than the three-chamber model with complete burst, which was closer to the measured data. These results suggested that an incomplete burst of diaphragms produced a shock wave with weaker shock wave strength, which was also observed in Ref. [46]. 4.2 The MSU Piston-assisted Shock Tube As shown in Figure 2.2, the shock tube can be operated as a piston- assisted shock tube by installing a piston in the low-pressure chamber. The operating and initial conditions are shown in Table. 4.1. Condition 1 Condition 2 High pressure 5.52 MPa (800 psi) air 19.3 MPa (2800 psi) air chamber P4 Low pressure 0.689 MPa (100 psi) air 1.03 MPa (150 psi) air chamber P1 Diaphragm burst 12 MPa 24.47 MPa pressure Shock tube initial 0.10135 MPa (14.7 psi) air 0.10135 MPa (14.7 psi) air ressure Piston mass weight 2 kg 2 kg Table 4.2 Operating conditions of the free piston shock tube A. Numerical Models Some simplifying assumptions were made in the numerical models. The gas was assumed to be ideal gas. A single diaphragm was used for simplicity. An ideal, complete and instant rupture of diaphragm was also considered while the gas leakage loss was not included in the models. 47 The 1D model was able to predict the compressive process but not able to predict P6 at the outlet of the blast tube. The 2D and 3D mesh grids were generated in GAMBIT and shown in Figure 4.11 based on the grids dependence study in section 4.1A. A total of 11,132 quadrilateral and triangular cells were generated for the 2D model and 54,715 hexahedral and tetrahedral elements were generated for the 3D model. A dynamic mesh technique was employed to simulate the moving zone caused by the piston > ”*1"—‘.‘ fit. . .2 (a) 20 mesh . lTll lit i i .7 i ‘ " ’l‘ * i—i 7L]. P“. —.~ v . ‘1’] 4’4 7 i,‘ t3 1.. list-"‘41 I l I ‘ ‘ A in" (b) 3D mesh Figure 4.11 Mesh for the free piston shock tube motion. The segregated-implicit solver was found to be the only numerical scheme which was compatible with the dynamic mesh. The piston motion was controlled by two subroutines, which calculated the force on the piston from 43 P4=5.52MPa (800 psD , P1=0.689MPa (100psi) 16 ’ Piston (2kg) 14 _. Air-Air MeasuredP5 ‘2 r JIANG P510 10 ~ . ,. —FLUENT P520 l 8 - —FLUENT P530 5 -- 4 .. 2 _ 0 '1' I I I 0.005 0.025 0.045 0.065 Time (s) (a)Condition 1 P4=19.3 NPa (2800psi), P1=1.03 IVPa (150psi). 20,Piston(2kg) 10° " Air-Air go . 80 ~ — Measured P5 70 « E 60 , —JIANG-1D g 50 ‘ ——FLUENT-2D a 40 « 30 s 20 - 10 . O I I I ———I 0 0.01 _ 0.02 0.03 0.04 Time (s) (b) Condition 2 Figure 4.12 Simulated P5 the flow field solutions and then solved the corresponding acceleration and velocity of the piston. The laminar model was used as the viscous model. B. Results and Discussions Upon the rupture of the first diaphragm at the intermediate-pressure chamber, P5 at the end of the low-pressure chamber rose sharply due to the 49 gas compression by the piston. Once it surpassed the rupture pressure, it burst the second diaphragm. P5 is the important parameter in this study. A detailed description of P5 can help examine the compression process and should be helpful for the design of the second diaphragm and the blast tube. In Figure 4.12, the simulated P5 by 1D, 2D and 3D models were compared to the measured data. For the maximum P5, which burst the second diaphragm, good agreements between the simulations and the experiments were achieved for both condition 1 and condition 2. The compression process was well predicted. The 1D model showed an overall better prediction on P5. 50 P4=5.52MPa (800 psi), P1=0.689MPa (100psi), g Piston(2kg) 20 Air-Air 15 Measured P6 (15 CL —FLUB\ITP62D E10 (0 —FLUENTP63D IL 5 r l l 0 l i i 0.02 06 4 . hmem) (a)Condition 1 P4=19.3 MPa (2800psi), P1=1.03 MPa (150psi), 2D, Pisotn(2kg) 100 Air-Air 80 — Measured A (i —Fluent20 E 60 E l E 40 i. 20 - 0 l I l fl 7 0 0.01 0.02 0.03 0.04 Time (s) (b) Condition 2 Figure 4.13 Simulated P6 Upon the rupture of the second diaphragm, the interaction between the very high pressure P5 at the low-pressure chamber and the low pressure at the blast tube generated a strong shock wave, which impacted the testing materials at the outlet of the blast tube. The pressure at the outlet of the blast tube, which is P6, was used to evaluate the shock wave. Thus P6 was the most important parameter in studies concerning shock wave based loading. 51 P465M Pa (80msi), P1=0.689M Pa (”$0.30 Piston(2kg) 200 — - - - Air Air —J|ANG-1 D 1 50 a — FLUBVT 20 .5 so _ Q 0 T I I :9: O 0.01 0.02 0.03 II '50 “ Time (s) -100 s -150 _. (a) Piston velocity P4=5.5IVPa (BOOpsi), P1=0.689MPa (100psi), 30 6 _ Piston(2kg) Air-Air Piston-position (m) — FLUBVT-ao 1 _ o I— I I I TI 0 0.01 0.02 0.03 0.04 0.05 Time (s) (b) Piston position Figure 4.14 Simulated piston motions (condition 1) In Figure 4.13, the CFD simulation results show excellent agreement with the measured data. The 3D model shows a better decay rate than the 2D model. These 2D and 3D models were justified to be able to provide suggestions for future operations of the MSU shock tube and the simulated P6 can be used as a pressure input in numerical studies for the dynamic response of materials under the shock wave. 52 First diaphragm Piston Second diaphragm High-pressure Low-pressure chamber chamber P4=5.5 MPa P1=0.689 MPa (800 psi) (100 psi) (a) Without sensor Sensor First diaphragm Piston Second diaphragm I High-pressure Low-pressure chamber chamber P4=5.5 MPa P1=0.689 MPa (800 psi) (100 psi) (b) With sensor Figure 4.15 Influence from sensor (condition 1) The piston within the low-pressure chamber can be accelerated to a high speed of around 200 m/s in the compression process. Thus the detailed description of the piston motion was critical to the operation safety of the shock tube. The CFD simulation provided a useful method to predict the piston motion. In Figure 4.14, the simulated piston velocity and piston position were compared based on 1D, 2D and 3D models. It shows that the 1D model predicted higher velocity than the 2D and 3D model. The discrepancy accumulated was noticeable in the rear part of the simulation. It was believed 53 to be caused by the viscous effects. The 1D model solved the inviscid Euler Equations which neglected both viscous effect and boundary layer effect, while the 2D and 30 model using FLUENT accounted for viscous effect. 0.1- . -. c U. a...“ WC... O. ’~€:~.¢ .5 naotok-v -- 9- a. S/rm‘etric Axis Blasttube Figure 4.16 Mesh grids considering the sensor C. Further Modification on the CFD Model To measure P6, a pressure sensor was placed at the outlet of the blast tube. The sensor has a diameter of 4 mm, while the blast tube had a diameter of 12.8 mm as shown in Figure 4.15. The sensor’s influence in previous studies was neglected initially. To consider the influence from the sensor, a new 2D CFD model was introduced, as shown in Figure 4.16. The simulation result was plotted in Figure 4.17. It showed a better pressure decay rate than the model 54 P4=5.5IVPa (BOOpsi), P1=0.689NPa (100psi), 2D, 16 _ Piston(2kg) Air-Air 14 e 12 .4 —Measured " 10 s g 8 With Sensor Geometry IE 6 _, —Wlthout Sensor Geometry 4 a 2 J O h‘ “‘ ““ , , 0.02 0.04 0.06 Tim e (s) Figure 4.17 Sensor’s Influence on P6 simulation without considering the sensor, although there was a 10% computation time increase due to the finer mesh near the sensor geometry. 55 4.3 The HEK shock tunnel and the T3 shock tunnel Two other shock tubes, HEK and T3, were also used for the justifications of the CFD codes. The HEK shock tunnel was a medium-sized piston-assisted shock tunnel constructed at the National Aerospace Laboratory (NAL) Kakuda Research Center. T3 was a medium-sized piston- assisted shock tunnel at Australia National University. These two shock tunnels had similar structure as shown in Figure 4.18. The major dimensions are listed in Table 4.2 while the operating and initial conditions are given in Table 4.3. 2"“I reservoir Compression tube Orifice I_ Nozzle Shock tube Piston P5 \L . ' Diaphragm P6 Figure 4.18 Schematic diaphragm for HEK and T3 Table 4.3 Geometry parameters of MSU, HEK and T3 shock tubes Length Diameter Piston Testing- . . Driving Compressron Shock Compressmn Shock gas tube tube tube tube IMSU 4m 0.154m 80mm 12.8mm kg Air-Air HEK 16m 6.5m 210mm 72mm 30kg Helium-Air T3 .8m 7.7m 150mm 38mm 90kg Helium-Air 56 Dump tank Table 4.4 Operating conditions of HEK and T3 HEK T3 High pressure 5.46 MPa air 3.7 MPa air chamber P4 Low pressure 0.113 MPa helium 0.0573 MPa helium chamber P1 Diaphragm burst 52 MPa 44.2 MPa pressure Shock tube initial 0.113 MPa air 0.0356 MPa air pressure Piston mass weight 30 kg 90 kg A. Diaphragm Rupture Model . I: dr dd Piston high-pressure low-pressu re blast tube chamber chamber dump tank Ruptur motion of the hole on the diaphragm.Rigid motion controlled by Eq. 4.1 Figure 4.19 Modeling of the diaphragm rupture process In section 4.1E, the influence of the diaphragm was studied by assuming an ideal, instant and incomplete rupture at the diaphragm. In reality, there was indeed a rupture process in the real experiments, instead of an instant rupture, which involved a time span of 300 #3 [48] as shown in Figure 4.19. The effect from this diaphragm rupture process was studied. Outa [48] 57 proposed a model for the diaphragm rupture process. The time dependent ruptured diameter of the diaphragm (d,) was described by Equation (4.1) 2 2 fl' t ’Zzl—cos— — (41) dd 2 top - where dd is the diameter of the final rupture hole on the diaphragm, d, is the diameter of the rupturing hole on the diaphragm, tis the rupture process time, and top is the opening time of the diaphragm. B. Numerical Models Some simplifying assumptions were made in the numerical studies. The gas was assumed to be ideal gas. A single diaphragm was used for simplicity. The gas leakage loss was not included in this model. In these two shock tunnels, the piston was designed to impact the piston stopper with a velocity of around 10 m/s after the diaphragm rupture, and then bounce back. To avoid the piston to impact the end wall, which would generate negative cell volumes in the computation domain, the piston was assumed to be stopped when the piston moved to 5.2 cm away from the tube end. The 2D mesh grids were generated in GAMBIT. A total of 14,992 quadrilateral and triangular cells were generated for the HEK tunnel and a total of 6,260 quadrilateral and triangular cells were generated for the T3 tunnel. A dynamic mesh technique was employed to simulate the moving zone because 1st of the piston motion. The segregated-implicit solver with -order accuracy 58 SILIILJ. scheme was used for both cases. The piston motion was controlled by two subroutines (User Defined Functions), one for each side of the piston, which calculated the force on the piston from the flow field solution and then solved the corresponding acceleration and velocity of the piston. The diaphragm rupture model, Equation 4.1, was integrated to the CFD model by another two subroutines. A Non-Conformal grids technique was employed for the grids around the diaphragm rupture zone, which allowed the use of block mesh generating to reduce the convergence difficulty in the dynamic mesh. The laminar model was used as the viscous model. C. Results and discussions For the HEK shock tunnel, the simulated P5 and P6 with consideration of the diaphragm rupture process were shown in Figure 4.20. As shown in Figure 4.20(a), the compression process was well predicted. At the time of diaphragm rupture, the simulated piston velocity was 135.49 m/s, which was close to the measured piston speed of 134.3 m/s. Two kinds of diaphragm rupture processes were assumed with the orifice diameter of 50 mm. One was an instant rupture represented by the green line and the other was with a 300 p: diaphragm rupture process, represented by the orange line in Figure 4.20a. In the case of instant diaphragm rupture, the mass flow was relatively high in the simulation because the diaphragm became fully open instaneously. Thus, the numerical result with instant rupture was not able to predict the gradual increase of pressure after the diaphragm rupture. The numerical 59 model with a 300 ,us rupture time was able to predict the gradual pressure rise 70“ 60 50 40 P5 (MPa) 201 10 -~ 70 60 after the 30* HEK Shock tube P4=5.46MPa, P1 =0.1 13MPa, 2D, Piston(30kg) Helium-Air — Measured —— 0us opening time with 1 50mm diamter -— 300us opening time with - 50mm diameter — Ous opening time with 20mm diamter 0 2 4 6 8 Time (ms) (a) Simulated P5 HEK Shock tube P4=5.46MPa, P1 =0.1 13MPa, 20, Piston(30kg) 1 Helium-. Measured i --—— 0us opening tiem with 50mm f. diameter ‘ ~ 300us opening time with 50mm 1 diameter -— Ous opening tiem with 20mm - 2‘ diameter \ —4 “MR—‘H‘MH 1 1.5 2 2.5 3 3.5 4 4.5 Time (ms) (b) Simulated P6 Figure 4.20 Simulation results for HEK diaphragm rupture accurately, as shown in Figure 4.20(a). The 60 diaphragm rupture process was shown to have a noticeable influence on the pressure history of P5. Similar conclusions can be found at reference [48]. T3 Shock tube 40 - P4=3.7MPa, P1=0.0573MPa, 2D, piston(90kg) Helium-Air 35 « ——- Simulation in Reference —— Measured in Reference — Ous opening time 300us opening time 30‘ 25* 20* P6 (MPa) 15* 10~ Time3(ms) Figure 4.21 Simulated P6 for T3 The holding time of the maximum P5 in simulation was shorter than the experiments, which was considered to be caused by lack of the information of the geometry at the orifice connecting the low-pressure chamber and the shock tube chamber. With the reduction of the diameter of the orifice from 50 mm to 20 mm, the holding time became closer to the measured data. This suggested that the cross-sectional area of the fluid flow might be smaller than the diameter of the orifice used, Le. 50 mm. And an incomplete rupture of the diaphragm might have occurred In Figure 4.20(b), the magnitude of simulated P6 shows a good agreement with the measured data. Due to the lack of the nozzle geometry 6] information and the pressure sensor positions, the holding time of the simulated nozzle pressure P6 was much smaller than the measured data. However, the simulated shock wave speed was found to be 3,970 m/s, which was close to the measured shock wave speed of around 4000 m/s. In the simulation of the T3 shock tunnel, similar phenomena, such as a good pressure magnitude and a much shorter holding time, were observed for I the nozzle pressure P6 as shown in Figure 4.21. The simulated shock wave speed was found to be 4,032 m/s, which was close to the measured 4,150 m/s and the simulated 4,110 m/s in reference [41]. To improve the simulation accuracy, details of the nozzle geometries, the orifice geometries and the sensor positions are required for further studies. 62 4.4 Summery The 1D CFD code JIANG was justified to be suitable for the simulation of the simple shock tube and the simulation of P5 in the piston-assisted shock tube. The 2D and 3D CFD models using FLUENT, however, were justified to be feasible for the simple shock tube simulation and prediction of both P5 and P6 in the piston-assisted shock tube simulations. Several conclusions are drawn from the numerical studies and the following recommendations are offered: 1. 3D model with coupled-implicit scheme was suggested for the simulations for P5 in the simple shock tube. 2. 2D model with coupled-explicit scheme was suggested for the air-air performance studies for the simple shock tube. 3. 3D model with segregated-implicit scheme was suggested for the free piston simulations. 4. Non-ideal diaphragm rupture has noticeable influence on P5. 63 CHAPTER 5 UPGRADING THE MSU SHOCK TUBE Using the justified CFD analysis methods described in Chapter 4, several possible upgrades of the MSU shock tube facility to improve the pressure performance were explored in this chapter, which included the combination studies of the driving gas and the application gas and the geometry modifications. 5.1 Driving Gas/Application Gas Combination Studies Recall from the theoretical analysis of the shock tube performance in section 3.1. The shock tube performance depends on two major parameters, one is the initial operating pressure P4/P1 and the other is the sound speed ratio between the driving gas and the application gas a4/a1. The influence of the first parameter was described in section 4.1D. The influence from the second parameter will be discussed in this section. Because the strength of the shock wave increased as the ratio of the sound speed increased, it was desirable to have a driving gas with low molecular weight and hence a low ratio of specific heat such as Hydrogen or Helium. On the other hand, it was desirable to have an application gas with high molecular weight such as Argon. So, a light driving gas and a heavy application gas combination was desirable in generating strong shock waves. Several combinations were investigated using CFD analysis methods implemented with FLUENT. A. Numerical Models 64 The 3D CFD model for the simple shock tube described in section 4.1 was utilized to study this problem. The species transport model, which was described in reference [54], was introduced to model the species transport phenomena between the driving gas and application gas. Besides the basic simplification assumptions made in section 4.1, the possible chemical reactions between the driving gas and application gas were not included in the model. Condition 1 in Figure 4.1 was used as the initial condition. B. Results and Discussions P4=1.72M3a (250 psi) . p1=0.27IVPa (40psi) 20 4 - Air-Air — He'He Air-Air - I-Ie-Air Air-He —— H2-I-I2 — I-Ie-H2 —— i-2-l-le —I-2-02 No Reaction 0 0.005 0.01 0.015 0.02 — HZ-Air No Reaction Time (s) Figure 5.1 Influence of Driving gas/Application gas combination on P5 Nine combinations of driving gas - application gas were studied. They included Air-Air, Helium-Helium, Helium-Air, Air-Helium, Hydrogen-Hydrogen, Helium-Hydrogen, Hydrogen-Helium, Hydrogen-Oxygen and Hydrogen-Air. Shown in Figure 5.1, the Air-Air, Helium-Helium and Hydrogen-Hydrogen combinations produced shock waves with same strength at 1.5 MPa. Also shown in Figure 5.1, the Hydrogen-Oxygen and Hydrogen-Air combination 65 generated a shock wave pressure of 3.5 MPa, which was more than two times that from Air-Air combination (1.5 MPa). Similarly, the Helium-Air combination generated P5maxiP1 P5maxiP1 ~ P4iP1 Helium-Air 100 - . . 80 r fi/e/Ei «Jr/“.3 i" so - /€ - 40 - Analytical 1 (I) 3D-coupled-implicit 20- £3 3D-compled-explicit i El ZD-compled-explicit 0 '— ' ' ' ' ‘ 0 20 40 60 80 100 P4iP1 (a) P5/P1 vs initial conditions T5maxi'T1 ~ P4iP1 Helium-Air El . . . 1:1 7 - El 11- D 5 . El M. E 5 _ Cl X (B 1% 4 i "" .--. Analytical 3‘ a 3‘ iii SD-coupled-implicit 2 _ 9 {1‘- BD-coupled-explicit I III 2D—coupled-explicit 1 1 1 4 1 0 20 40 50 80 100 P4iPl (b) T5/T1 vs initial conditions Figure 5.2 Performance curves 66 a Mach number ~ P4iP1 Helium-Air 4 . . , m D t: 35 D (gas-"'7‘; g 3 A i9 C) CD If?) 5 We g 8 0 ——-Analytical . ' . ‘E BD-coupled-implicit . 2 A 30-coupled—explicit ‘5 Cl ZD-coupled-explicit . ('5 2 # 20 40 60 80 100 P4iP1 (c) Mach number vs initial conditions Figure 5.2 Performance curves shock wave pressure of 2.7 MPa, which was almost two times the shock wave pressure generated in the Air-Air combination. It can be concluded from these studies that the lower the molecular weight ratio between the driving gas and the application gas, the stronger the shock wave generated. Although the Hydrogen-Oxygen and Hydrogen-Air combinations produced the highest overpressure, they also attenuated much faster than other combinations because the expansion waves traveled much faster in the light gas, such as Hydrogen, than in the heavy gas, such as Air. This attenuation was also considered to be a result from the boundary layer viscous effects which consumed the energy from the driving gas [16]. Because this research sought to achieve high pressure shock wave and high energy from the shock wave, the Helium-Air combination, which generated higher overpressure than the Air-Air combination and had less 67 decay rate than the Hydrogen-Air, was chosen as the driving gas - application gas in the future experimental studies. Although Hydrogen-Air combination could provide higher overpressure, their possible chemical reactions (combustion) [33] were considered less controllable and more hazardous. To avoid these concerns, the Helium-Air combination was investigated for future experimental studies. The performance curves were generated using the method described in section 4.1 D. Figure 5.2 shows the dependence of P5max, T5max and Mach number of the shock wave on the initial operating conditions. Comparing Figure 5.2 to Figure 4.8, the pressure performance was greatly increased by using Helium as the driving gas. For example, when the initial operating condition of P4/P1 was 10, the generated P5/P1 was 15 in Helium-Air case, which was more than two times of P5/P1 in the Air-Air combination case. The numerical schemes were also compared in Figure 5.2, and the 3D coupled-explicit scheme showed an overall better agreement with the analytical solutions. 68 5.2 Geometry Modification Studies Parabola 13 Monit red at focus “4'33 Parabola 2: Monitored at focus Figure 5.3 Parabola end Besides the different driving gas - application gas combinations described in section 5.1, two types of geometry modifications to improve the pressure performance in the simple shock tube were investigated. The basic idea was to explore the possible pressure performance improvements by converging pressure waves as discussed in reference [36]. One modification was modifying the end of the low pressure chamber to a parabola shape as shown in Figure 5.3. Two parabola shapes were investigated, whose cross profiles were defined by Eqn. 5.1: A. Parabola 1: X = —187.5y2 +6 B. Parabola 2: x = ‘6253’2 + 6 (5.1) 69 solid core outerwall stream line solid core 250 psi air end Figure 5.4 Annular shock tube The other modification was modifying the MAU simple shock tube to an annular shock tube as shown in Figure 5.4. A 6 m long inner body with a diameter of 60 mm was mounted coaxially in the interior of the outer tube which was defined in Figure 2.1, thus forming an annular channel to generate an annular shock wave. The annular shock wave converged at the right end of the annular shock tube. A. Numerical Models The 3D CFD models were generated using GAMBIT as shown in Figure 5.3 and Figure 5.4. Quadrilateral and triangular cells of 10,560, 42,163 and 17,362 were generated for the parabola 1, parabola 2 and the annular shock tube respectively. Same simplification assumptions were made as in section 4.1. Initial condition 1 in section 4.1, which was 1.72 MPa (250psi) in the high- pressure chamber and 0.27 MPa (40 psi) in the low-pressure chamber, was used for this case study. The coupled-implicit with 2"d-order scheme was chosen as the numerical scheme. The laminar model was used as the viscous model and Air-Air was used as the driving gas - application gas combination. 70 B. Results and Discussions P4=1.72lVPa (250 psi) , P1=0.27MPa (40psi), so 3 “ Air-Air 2.5 - A 2 # E E1 5 — § _ E 1 _ No parabola end Parabola 1 focus 0.5 a Parabola 2 focus 0 ‘ annular converge focus 0.008 0.0085 0.009 0.0095 0.01 Time(s) Flgure 5.5 Simulated P5 In Figure 5.5, P5 at the focus of the two parabola ended and P5 at the right end of the annular shock tube were compared to evaluate the shock wave strength. As shown in Figure 5.5, the shock front pressure at the normal simple shock tube was 1.5 MPa. The P5 was increased to 2.8 MPa using parabola 1 while to 2 MPa when using parabola 2. There was no noticeable pressure performance improvement in the annular shock tube design. This result indicated that modifying the end to a parabola shape was a possible method to improve the pressure performance. 7] 5.3 Summery Driving gas - Application gas combinations were studied using the simple shock tube model. The low molecular weight ratio between the driving gas and the application gas was desirable for improving the shock wave strength. The Helium-Air was chosen as the future driving gas - application gas combination based on the high pressure performance and the low attenuation rate. The 3D coupled-explicit scheme was also recommended for this case study. Two types of geometry modifications were investigated to improve the pressure performance. The parabola end was shown to be a possible method to improve the pressure performance. 72 CHAPTER 6 SIMULATION OF REAL BLAST WAVES Explosions can happen under water or in sands in the real applications. Experimental simulations studies can be conducted in the MSU free piston shock tube to simulate real explosions under these scenarios. As a preliminary study, CFD simulations for these experiments were described in this chapter. As an effort to mitigate the blast wave, several methodologies, such as the structure reinforcement, water blanket and venting, have been studied and developed recently. In this chapter, one possible method, using venting holes to mitigate shock wave, was investigated using the CFD analysis. 6.1 Real Blast Waves Simulations Application gas (0.1013 MPa) AI Compressed hig'i pressure driving gas (11.15 MPa) Sand or water (74 mm in length) Blast tube Low-pressure chamber Dump-tank Figure 6.1 Problem setup for shock waves under water and sand To simplify the simulation, the compression process of the MSU free piston shock tube was not included. The compression results from condition 2 73 in Table 4.1 were adopted as the initial condition of these simulations. As shown in Figure 6.1, the compressed driving gas with a pressure of 11.15 MPa was stored in the low-pressure chamber. The application gas (air) was stored in the first half of the blast tube with atmosphere pressure, followed by water or 56% sand in the rest of the blast tube. The sand particles for numerical simulation are 111 min diameter, 2500 kg/m3 in density and 0.001003 kg/m.s in viscosity. The driving gas and application gas produced the shock wave, which impacts the water or sand in the blast tube. The motion of the sand and water were monitored in simulations. A. Numerical Models Some simplifying assumptions were made. The gas was assumed to be ideal gas. A single diaphragm was used for simplicity. An ideal, complete and instant rupture of diaphragm was also considered while the gas leakage loss was not included in the model. The 2D mesh grids were generated in GAMBIT. A total of 25,111 quadrilateral cells were generated. The segregated-implicit solver was found to be the only numerical scheme which was compatible with the mixture multiphase model. B. Particle Drag Coefficient The particle drag coefficients in non-stationary flow are found to be considerably higher than those in a steady flow case [50]. To correctly model the transient fluid — particle interaction in the shock wave — sand simulation, an accurate drag coefficient correlation needs to be provided. By experimentally 74 investigating the particle motion in a shock tube, lgra and Takayama [57] proposed a drag coefficient correlation, which could safely be used for 200 s Re 3 101000. This correlation is shown in Equation 6.1. 10g10 CD = 7.8231 -5.8137logIO Re+ 1.4129(log10 Re)2 —0.1146(1ogl0 Re)3 (6.1) Where CD is the drag coefficient, Re is the Reynolds number. Coupled Equation 6.1 with the fluid-solid exchange coefficient [54] as shown in Equation 6.2, the drag coefficient can be incorporated into the CFD model using FLUENT. Vs _Vg 3C Dasag ,0 g 465 K38 = 4d 8 S (6.2) Where K is the fluid-solid exchange coefficient, a is the volume fraction, v is velocity and p is density. The subscription sis solid, 9 is gas. P5=1 1 .15MPa, P(blasttube)=0.1013MPa 25 ‘ -—-Air 20 g ——Sand A j —Water -— 3315 .2, $10 — 5 a 0 ' ' J. l i 0 0.0005 0.001 Time (s) Figure 6.2 Total pressure time history at the blast tube outlet 75 C. Results and Discussions -0.02 -0.04 a V 0.04 0.02 -0.02 -0.04 b) Simulated velocity contour at time=0.71 ms after diaphragm rupture Flgure 6.3 Shock wave-sand 76 In Figure 6.2, the simulated total pressure at the blast tube outlet was compared to the pressure history of all air case, which is from Figure 4.13. It shows that the shock wave attenuated in the sands from 15 MPa to about 8 MPa since part of the energy was transferred to the kinetic energy of sand particles. The shock wave — water interaction generated a compression pressure wave since the air-water interface speed is much lower than the sound speed in water. Since the water possesses a much higher density than air, the peak total pressure of the water jet, which is around 23 MPa is much higher than the 15 MPa from all air case. For the shock wave - sand case, the shock wave arrived at the blast tube outlet around the time of 0.71 ms as shown in Figure 6.2. In Figure 6.3, the details of the sand distribution and velocity contours of sand were plotted for this time point. Figure 6.3 (a) shows that the high pressure and high speed air flow burts into sand in the blast tube. The contact face was about 3 cm after the shock wave. Figure 6.3b shows that the sand particles were accelerated to around 80 m/s. In Figure 6.2, the pressure drops to a plateau pressure around 0.98 ms. Details of the flow field were plotted in Figure 6.4 for time around 0.98 ms. In Figure 6.4 (a), it shows that the air jet has penetrated the sand. As shown in Figure 6.4 (b), the sand particles front was accelerated to 150 m/s. The spherical pressure wave was shown in Figure 6.4 (c). 77 (a) Simulated vo ume fraction of sand at time=0.98 ms after diaphragm rupture 0.04 . Voioaw (mfs) .12 6.14x6.18 6.18 6.2 (b) Simulated velocity contours at time=0.98 ms after diaphragm rupture Pressure 1.0E+07 9.5E+06 8.4E+06 7.4E+08 6.9808 4.3E+06 3.7E+06 3.2E+06 2.2E+06 1.7E+06 6.2E+05 1.0E+05 "5.1' ' 6.12 5.14. A (c) Simulated F'ressure contours Figure 6.4 Shock wave-sand (t=0.98 ms) 78 Volume fraction :1 Bis 8.18 62 (a) Simulated volume fraction of water at time=0.98 ms after diaphragm rupture Velocity (mfs) 8.1 . 6.14 $16 6.18 6.2 (b) Simulated velocities contours at time=0.98 ms after diaphragm rupture Pressure 7.7E+06 T.3E+06 6.6E+06 6.2E+06 5.4E+06 5.1E+06 4.3E+06 3.5E+06 2.4E+06 1.2E+06 8.6E+IJS 1.0E+IJS (c) Simulated Pressure contours Figure 6.5 Shock wave-water (t=0.98 ms) For the shock wave — water case, the peak pressure at the blast tube outlet happens around 0.98 ms after the diaphragm rupture, as shown in 79 Figure 6.2. The details of the flow field around 0.98 ms were plotted in Figure 6.5. Figure 6.5 (a) shows that the water has been pushed out by the air jet at this time point. And the water droplet was accelerated to around 160 m/s as shown in Figure 6.5 (b). The corresponding pressure contours are shown in Figure 6.5 (c). 80 6.2 Shock Wave Venting on Plates 1mm 142.157 mm — E E 51‘ Hole numbers and geometry r 1 Thickness of first plate 1.5 mm——I defined below I 000 coo coo 1-hole 5-hole 9.|10Ie 16-hole Flgure 6.6 Case 1 12.8 mm L A venting method was proposed to lower the shock wave pressure on the composite plates, which was also under experimental studies. This section was to describe the CFD analysis for this venting method, which included two cases studies. For case series 1, as shown in Figure 6.6, the small blast tube in Figure 2.2 was used for the case studies. An incident shock wave with 4.1 MPa overpressure traveled from the left to the right of the shock tube. The plate with venting holes was placed near the right end. Four kinds of venting holes were studied, which were 1-hole, 5-hole, 9-hole and 16-hole. The total area of the holes was same in these four cases, which was 3.927 mmz, while the total 81 cross-sectional area of the plate was 128.7 mmz. The simulated pressures on the plate and after the plate were compared to see the venting effect. 12.8 mm f 142.157 mm ( Thickness of the plate 1.5 mrnal l \L mm 1.... liZBmmj—r; Hole geometry defined below Plae2 Plae1 (a) 5-holes with 0 degree twisted O O O O O O O 0 Plate 1 Plate 2 (c) 4-hole with 0 degree twisted o o o o o Pla e 1 Plate 2 (b) 5-holes with 45 degree twisted o 0 ° 0 o o O o Pla e 1 Plate 2 (d) 4-hole with 45 degree twisted Flgure 6.6 Case 2-2 plates For case series 2, as shown in Figure 6.7, two types of plates with different number of venting holes were investigated, one had five holes on each plate (Figure 6.7a) and the other had four holes on each plate (Figure 6.7b). In addition, two rotation angles (0 degree and 45 degree) between plate 1 and plate 2 were investigated. Initial conditions were same as Case 1. The pressure on plate 1, plate 2, and after the two plates were monitored to compare the effectiveness of the venting method. 82 A. Numerical Models Some simplifying assumptions were made for the case studies. The plate was assumed to be fixed. The dynamic behavior of the plate material, and the interaction between the fluid and the plate were not included in the model. A single diaphragm was used for simplicity. An ideal, complete and instant rupture of diaphragm was assumed. The gas leakage loss was neglected. The mesh grids, with 21,726, 27,049, 41,220 and 70,353 quadrilateral cells, respectively, were generated for the 1-hole, 5-hole, 9-hole and 16-hole cases in case series 1. The mesh grids, with 27,260, 27,956, 25,090 and 27956 quadrilateral cells, were generated for the four cases in case series 2. The coupled-implicit scheme with 2'"d order accuracy was used as the numerical scheme. Laminar model was chosen for the viscous effect. B. Results and Discussions For case series 1, the average pressures on plate 1, as defined in Figure 6.6, was plotted in Figure 6.8. The results showed that the maximum pressure on plate 1 was about 5% lower than the incident shock wave, which was around 4.1 MPa. The speed of the shock wave increased as the number of holes increased. The incident shock wave impacted plate 1 and was transmitted through the venting holes. To evaluate the effectiveness of plate 1 on the shock wave mitigation, the pressure after plate 1 was plotted in Figure 83 A Average pressure on plate 1 .5 4 1 —Noplate —1hole-1 late 3.5 1 p E 3 j ———5 holes-1 plate 2 2 5 9 holes-1 plate g 2 - ‘ -———16 holes-1 plate a 1.5 . 1 a 0.5 4 0 “I T I T T 0.12 0.14 0'16Timecli118 0.2 0.22 0.24 .0 —L Figure 6.8 Simulated average pressures on plate 1 6.9. In comparison with the incident shock wave pressure, the average pressure on the right and significantly dropped from 4.1 MPa to around 1.6 MPa after the shock wave passed plate 1. The number of venting holes did not affect the average pressure on right end, with a cross area of 128.7 mm2, as shown in Figure 6.9a. Because the center of the right end was the most possible place to fall, which had been observed in the experiments, the simulated pressure at the center was plotted in Figure 6.9b. The maximum pressure at the center of the right end was decreased to 1.6 MPa in comparison with the incident shock wave with 4.1 MPa pressure. In Figure 6.9, the first arriving shock at the center of the right end was 1.6 MPa, 0.58 MPa, 0.4 MPa, and 0.36 MPa for 1-hole, 5-holes, 9- holes and 16-holes, respectively. Thus, the shock wave pressure at the center of the right end was greatly influenced by the number of the venting holes. From these results, most of the energy was reflected by plate 1, and only a small portion of the energy vented through the venting holes and impacted the right end. This suggested that the pressure energy between plate 1 and the 84 right end can be intentionally distributed by a proper design of the venting holes with the aid of CFD simulations. For example, in order to lower the pressure at the center of right end, the venting flow path through Plate 1 may be designed not aligned with the object, which needs to be protected, to vent the high pressure air to the peripheral air. 4.5 1 Average pressure at right-end wall __ no plate 4 1 —1 hole-1 plate 3 5 1 —— 5 holes-1 plate ' 9 holes-1 plate E 3 1 —— 16 holes-1 plate 2 v2.5 1 2 I 1.5 . 1 a 0.5 « 0 l T T T I 0 0.1 0.2 11mm”) 0.3 0.4 0.5 (a) Average pressure on right end 4.5 - Pressure at the right-end mll center 4 7 — no plate 3.5 . -—-—-1 hole-1 plate A 3 A —— 5 holes-1 plate E 2 9 holes-1 plate V '5 7 —- 16 holes-1 plate E 2 A 1.5 7 x 1 - , ’ 0.5 - J1} ' 0 “ ##2##. T r r 1 0 0.1 0.2 0.3 0.4 0.5 Time (ms) (b) Pressure at the center of the right end Figure 6.9 Simulated pressures on right end 85 In comparison to a single plate with venting holes, two plates with venting holes were simulated in case series 2. Figure 6.10 to Figure 6.12 showed the pressure contours on plate 1, plate 2, and the right end of the tube, respectively, for the 5-hole plates with a 45-degree mismatch (Figure 6.7b). Figure 6.10, showed that the shock wave with a pressure around 4.1 MPa impacted plate 1 (Figure 6.10a), and then vented through the five holes (Figure 6.10c). In Figure 6.11, the shock wave impacted plate 2 through the venting holes on plate 1. At the same time, its pressure was attenuated to around 0.7 MPa (Figure 6.11a). Since the external four holes of plate 2 had a 45-degree rotation with respect to those of plate 1, the venting holes between the two plates were not aligned except at the center hole. The shock wave could not vent through plate 2 without being attenuated first. In Figure 6.12, the shock wave impacted the center of the right end first since the two center holes in both plates were aligned. The shock wave strength was attenuated further to about 0.37 MPa (Figure 6.12a). The center shock wave then spread and interacted with the shock waves from the four external holes on plate 2, which arrived later because of the misalignment between the four external holes with those of plate 1. It was believed that the interaction between the two plates and that between plate 2 and the right end diffused the incident flow and subsequently attenuated the shock wave. 86 Pressu re . ,, 4.1E+06 [i 4.02143E+06 f j 3.94286E+06 ‘ 3.86429E+06 . 3.78571E+06 3.70714E+06 3.62857E+06 3.55E+06 3.47143E+06 L 3.39286E+06 .‘ 331429E+06 3.23571E+06 3.15714E+06 3.07657E+06 3E+06 a) 1:161? e-4 s Pressu re 4.1 E+06 : 4.02143E+06 3.94 286E+06 1 3.86429E+06 3.78571 E+06 3.70714E+06 3.62857E+06 3.55 E+06 3.47143E+OB i- 3.39 288E+08 ’ ’ 3.31429E+06 ‘ 3.23571 E+06 3.1 5714E+06 3.07 857E406 ' 3E+06 b) t=1.667 e-4 s Flgure 6.10 Contours of pressure on plate 1 87 c) t=1.735 e-4 s d) t=1.784 e-4 s Pressure 4.1 E+06 4.02143E+06 334286506 3.86429E+06 3.78571 E+06 3.70714E+06 3.82857 E+06 3.55 E+08 3.47143E+06 3.39 285E+06 3.31 429E+06 3.23571 E+06 3.15714E+06 3.07857E+06 3E+06 l‘ 1% E E! I i i 1 Pressure .. 4.1E+06 . . 4.02143E+06 . 334286506 336429306 ' 3.78571E+06 3.70714E+08 3.62857E+06 3.55506 347143306 3.39286E+06 3.31429E+08 ', 3.23571E+06 '“ 115714506 3.07857E+06 3E+06 Figure 6.10 Contours of pressure on plate 1 (continued) Pressure 710000 666429 622857 579286 535714 492143 448571 405000 361429 317857 274286 5 . 230714 , 187143 143571 100000 'I flii 1i! Ii if E E i. a) t=1.634 e-4 s b) t=1.667 e-4 s Figure 6.11 Contours of pressure on plate 2 Pressure 710000 . . 666429 - i 622857 .-9 579286 ' 535714 492143 448571 405000 361429 317857 274286 230714 187143 143571 100000 c) t=1.702 e-4 s d) t=2.46 e-4 s Pressure 710000 666429 822857 579286 535714 492143 448571 405000 361429 317857 274286 230714 187143 143571 100000 . E! if IN. E I I g unqgjyfivnrwr Pressure 950000 : 939286 ‘ 928571 917857 907143 ‘ 896429 885714 875000 864286 853571 842857 832143 821429 810714 800000 -unhnmmmfl!fliflw Figure 6.11 Contours of pressure on plate 2 (continued) 90 uwnhmwwmiifl77 a) T=1.578 e-4 s Pressure 200000 192857 ' 185714 178571 1 171429 184286 157143 150000 142857 135714 128571 121429 114288 107143 100000 filliflv '3 z! a E J... l: f b) T=1.618 e-4 s Flgure 6.12 Contours of pressure on right end c) T=1.7846-4 s d) T=3.265 e-4 s Pressure 800000 796429 792857 789286 785714 782143 778571 775000 771429 767857 764286 760714 757143 753571 750000 -—hHfi-ufi[fiig4 Figure 6.12 Contours of pressure on rlght end (contlnued) In Figure 6.13, the average pressures on plate 1 of four cases, including two plates with five holes aligned and misaligned and two plates with four holes aligned and misaligned, were compared. There was no significant difference in the average pressure on plate 1 between the four setups. In Figure 6.14, the average pressures on plate 2 of the four cases were plotted. It showed the setups (both four holes and five holes) with a 0-degree rotation angle, i.e. no misalignment of holes between plate 1 and plate 2, had much lower pressure than the setups (both four holes and five holes) with a 45- degree rotation, Le. a 45-degree misalignment between plate 1 and plate 2, since the fluids vented much faster in the 0-degree rotation than in the 45- degree rotation. In comparison to the maximum pressure of 1.6 MPa after plate 1 in Figure 6.9a, the maximum pressure after plate 1 was 0.76 MPa for 0-degree rotation, which is lower than the pressure of 1 MPa for 45-degree rotation, as shown in Figure 6.14. This is because the flow was easier to be vented out from Plate 1 in the O-degree rotation configuration than in the 45- degree rotation configuration. 93 Average pressure on plate 1 — N0 plate 4-5 ‘ — 5 holes-2 plates-0 degeree 4 ‘ — 5 holes-2 plates-45 degree A 3'5 _ — 4 holes-2 plates-0 degree 0 _ g 3 —— 4 holes-2 plates- 45 degree "' 2.5 2 3 2 — 2 a 1.5 — 1 _ 0.5 - O 2 7 v .ii -7 7v. A. 7k 7.7. .i77 , .Wfi 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 Time (ms) Figure 6.13 Average pressure on plate 1 4'5 " Pressure on plate 2 — No plate 4 — 5 holes-2 plates-Odegree 3.5 —- 5 holes-2 plates-45 degree —— 4 holes-2 plates-0 degree —— 4 holes-2 plates-45 degree Presssure (MPa) .-‘ I0 OI N 01 (a) .I ._l 2 I 1 2 —L Figure 6.14 Average pressure on plate 2 94 Pressure at the right-end center — No plate 4'5 7 — 5 holes-2 plates-0 degree 4 4 —— 5 holes-2 plates-45 degre — 4 holes-2 plates-0 degree a. 3'5 “ — 4 holes-2 plates-45 degrei o. 3 - ‘5’ 2 5 g s at 2 ‘ a) 2 1.5 - IL 1 a 0.5 4 O 7 I T I I 0 0.1 0.2 0.3 0.4 0.5 Time (ms) Figure 6.15 Pressure at the center of the right end The pressure at the right-end of the four cases was plotted in Figure 6.15. The fluids vented slower through plate 2 in the 45-degree misaligned setup and hence more energy was diffused between plate 1 and plate 2. Accordingly, the maximum pressure in the 45-degree setup was 0.8 MPa, which was lower than the maximum pressure of 1 MPa in the 0-degree setup. Comparing Figure 6.15 to Figure 6.9, the pressure at the right and center dropped more by using two plates with venting holes than using one plate. The setup with 4 holes and a 45-degree rotation showed a higher capability to decrease the pressure at the center in comparison with the 5-hole setup because there was no center hole in the 4-hole case, since the high pressure air flow did not impact at the center of the right and directly in the 4 holes configuration. 95 6.3 Summary As a preliminary study for experimental simulations of real blast under water and sands, the CFD simulations predict an impact speed of 150 m/s for sand particles and 160 m/s for water droplet under an incident shock wave with a pressure of 15 MPa. In comparison, the air was accelerated to around 800m/s in the all air case. The venting holes showed the capability to distribute the energy among plates and hence to mitigate the shock wave strength. With a proper flow path design, the shock wave strength could be attenuated through the interactions between flows in different directions. The two-plate configuration, with four holes on each plate and with 45 degree rotation, shows the best capability to attenuate the shock wave pressure. 96 CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS 7.1 Conclusions Overall, the pressure shape in the simple shock tube was well predicted by both the 1D code JIANG and 2D/3D code FLUENT. FLUENT predicted a better pressure increase by the interaction between the reflected shock wave and the driving-application gas interface. Both these two CFD codes and analytical studies over-predicted the plateau pressure around 15%. It was believed to be caused by the negligence of the loss. The dependences of performance parameters (g?- 1T; and M5) on the initial conditions were generated by both the CFD and the analytical calculations. For the free piston shock tube, the P5 at the low-pressure chamber was predicted very well by both JIANG and FLUENT. An overall difference was less than 5% when compared with the measured data. The P6 at the outlet of the blast tube was well predicted by FLUENT, and an overall difference was around 10% comparing to the measured data. The piston within the compression tube was accelerated to a high speed of around 135 m/s under a diaphragm burst pressure of 12 Mpa, condition 1 in Table 4.1. The non-ideal diaphragm rupture showed a noticeable influence on P5, which led to a continuing pressure increase of P5 after the diaphragm rupture. For shock tube simulations using FLUENT, the following conclusions were drawn. (1) 3D model with coupled-implicit scheme was recommended for the simulations for P5 in the simple shock tube. (2) 2D model with coupled- explicit scheme was recommended for the air-air performance studies for the 97 simple shock tube. (3) 3D model with segregated-implicit scheme was recommended for the free piston simulations. In the simulations for modification of the shock tube, the driving- application gas combination of Helium-Air showed a noticeable increase in the pressure performance in comparison with the current used Air-Air combination. The geometry studies showed that there was a noticeable peak pressure increase at the end of the shock tube if a parabolic end section was used. The plateau pressure was shown to be increased a little by both the parabolic end and an annular end. As a preliminary study for experimental simulations of real blast under water and sands, the CFD simulations predicted an impact speed of 150 m/s for sand particles and 160 m/s for water droplet under an incident shock wave with a pressure of 15 MPa. In the simulations, the shock wave strength was greatly attenuated using venting-hole method. The two-plate configuration, with four holes on each plate and with 45 degree rotation, is shown to be best setup of current test, which is able to reduce the pressure at the center of the plate from 4 MPa to 0.8 MPa. 98 7.2 Recommendations As shown in the simulation result in Chapter 5, the Helium-Air combination is recommended to improve the pressure performance. Details of the flow field, such as the vortex structures, need to be investigated for the shock wave convergence. The loss mechanism in the shock tube, such as the heat transfer loss, boundary layer loss and gas leaking loss, needs to be modeled properly in the numerical simulations. The venting method shows a great potential to attenuate the shock wave, although the venting holes affect the material strength. A coupling between the CFD analysis for the flow field and the FEA analysis for the testing material under different venting designs is highly recommended for future numerical simulations. 99 APPENDIX A Cases Studies A.1 Study On The Length of The Orifice of HEK Shock Tunnel As shown in Figure 4.18, the geometry parameters of the orifice is unknown. The numerical studies on the diameter of the orifice were discussed in section 4.38. To investigate the effect from the length of the orifice, a numerical study was performed. Three orifice lengths were investigated in the CFD simulation. They were 0.1m, 1m and 2m. Same numerical model was employed as described in section 4.3. HEK Shock tube 70 — P4=5.46 MPa, P1 =0.1 13 MPa, Piston (30kg) Helium-Afr — Measured 60 — Onfice-0.1 m 1 —— Orifice-1 m ——- Orifice-2m Time (ms) Flgure A.1 P5 As shown in Figure A.1, the holding time of P5 increases as the length of the orifice increases, but the holding time is still much less than the measured P5. It shows that the holding time of P5 strongly depends on the geometry of the orifice. 100 HEK Shock tube P4=5.46 MPa, P1=0.113 MPa, Piston (30kg) Helium-Air — Measured -— Orifice-0.1 m —Orifice-1 m — Orifice-2m As shown in Figure A.2, the peak of P6 decreases as the length of the As described in Chapter 1, the water mist is considered to be a possible Time (ms) Flgure A.2 P6 orifice increases from 0.1m to 1m, and then does not change much when the length of the orifice increases from 1m to 2m. The holding time of P6 shows less dependence on the length of the orifice. A.2 Shock Wave Propagation In Water method to attenuate the blast wave [40]. Also, the explosions under water are one of the research interests in this project. For an understanding of the fundamentals of the shock wave propagation in water, two basic cases conducted in the MSU simple shock tube (Figure 2.1) were investigated using CFD analysis as shown in Figure A3. 101 / /‘ / 250psi Ail’ 40 psi All' 40psi Water / Monitor: right end wall Monitor: X=1-5 Monitor: x=3.5 (a) Case 1 250psi Air 40ml A“: «Inst Wata 40psi Air Moriitor: x=1.5 Monitor: X=3.5 Monitor: X=4.5 (b) Case 2 Flgure A.3 Problem setup In case 1 (Figure A.3 a), 2 m of air was stored at left of the tube to generate the incident shock wave while 4 m of water was stored at the right of the tube. A shock wave with an overpressure of 0.8 MPa, which was generated at the left end, travelled right and impacted the water. The incident shock wave and the pressure histories in the water and at the right end of the tube were monitored. In case 2 (Figure A.3b), 2 m of air was stored in the left part of the tube, which was used to generate the incident shock wave. 2 m of water was stored in the center followed by another 2 m of air in the right. A shock wave with the same strength as case 1 was generated at the left end of the tube, and travelled right. A third case with all air in the tube as described in section 4.1 was used as the comparison to case 1 and case 2. 1.8 - Air-Air-Water 16 J —lncident Shock in ' air at x=1 .5 — 11 water x=3.5 1.4 4 75‘ — Right and a 1.2 - E 1 2 a 0.8 1 a: 2 0.6 1 D. 0.4 4 0.2 0 . . . , T , , 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Time(s) (a) Air-Water case 1 .8 1 Air-Air-Air-Air 1.6 4 —x=1.5 —-x=3.5 1'4 ——-Right end KI ‘\ 1.2 a 1 0.8 1 0.6 fl \ 0.4 ~ 0.2 — O r I 0 0.005 Time(s) 0.01 0.015 . V \ Preesure(MPa) (b) Air-Air case Flgure A.4 Simulation results of case 1 A. Numerical models Some simplifying assumptions were made in numerical analysis. The 103 gas was assumed to be ideal gas. A single diaphragm was used for simplicity. An ideal, complete and instant rupture of the diaphragm was considered. The gas leakage loss was not included in the model. The steam, which was generated by the interaction between the incident shock wave and the water, was not included in the numerical simulations. The 3D mesh shown in Figure 4.2b, with 8,382 quadrilateral cells, was used for this case study. Since the flow was stratified, the Volume of Fluid (VOF) multiphase model [49] was employed. The Geo-Reconstruct scheme was used as the VOF formation scheme, considering its accuracy in the time- dependent VOF simulation. The segregated solver with 1St order accuracy was found to be the only scheme that was compatible with the VOF multiphase model. The laminar model was used to include the viscous effect. B. Results and discussions For case 1, the incident shock wave, and pressure history in the water and at the right end were plotted in Figure A.4a. The all-air case was plotted in Figure A.4b for comparison. This result showed that the water acted as a wall, as the incident shock wave reflected at the interface of air and water. The pressure at the interface was increased by the incident shock wave and the reflect shock wave from 0.276 MPa to about 1.4 MPa. The pressure within the water was the pressure at the interface. 104 1'8 1 Air-Air-Water-Air 1.6 — —x=1.5 Incident shock 1-4 1 ——-—inwate x=3.5 E 12 g —Air after water x=4.5 E v 1 a g 1 0.8 33’ a 0.6 — 0.4 ~ 0.2 J 0 r f 1 O 0.005 Time(s) 0.01 0.015 (a) Air-Water-Air case 1'2 — Air-Air-Air-Air —x=1.5 1 ‘ ——x=3.5 A —-x=4.5 g 0.8 ~ II? n V 06 5 \J 0.4 « 0.2 — o , , fl 0 0.005 Time(s) 0.01 0.015 (b) Air-Air case Flgure A.5 Simulatlon results of case 2 105 1 Pressure distribution ‘—A"'A"'Wat9r'A" — Air-Air-Air-Air P/—\ _L a) _L .b _L N 1 A J Pressure (M Pa) 0.6 I 0.4 0.2 o r e e . o 1 2 3 4 5 6 Position (m) Figure A.6 Pressure distribution along the X-axls For case 2, the incident shock wave, and the pressure history in the water and in the air after the water were plotted in Figure A.5a. The all-air case was plotted in Figure A.5b for comparison. Comparing the pressure histories at the right end of the tube for the air-air-water-air and all-air cases (green line in Fig 6.3), the result showed that the incident shock wave generated at the left did not penetrate the water to the opposite air. In Figure A.6, the pressure distribution along the axis at 0.015 ms was plotted for both the Air-Air-Water- Air case and the Air-Air-Air—Air case. A similar result can be found in Figure A.6. The pressure in the air after the water stayed at the initial pressure of 0.27 MPa, which was beyond 4 m in Figure A.6. The pressure within the water distributed linearly along the distance. It indicated, the incident shock wave did 106 not penetrate the water since the air-water interface move much slower than the sound of speed in water. A.3 Testing Cases for the Hydrogen-Oxygen Detonations l H 0 21% Owen, 79% Nitrogen Igrition po'nt (200(k) 0.0(D1537 kg/s 50% Oxygen, 50%Hydrogen Figure A? Problem setup This problem is a testing simulation for the Hydrogen-Oxygen detonation experimental studies in reference [33]. As shown in Figure A.7, a mixture of 15% oxygen, 25% hydrogen and 60% Nitrogen was pumped into the testing tunnel with the flow rate at 0.0001537 kg/s. The testing tunnel has a length of 5m and has a diameter of 25 mm. Ignition position are 250 mm from the inlet. 107 Measured 5 I 4.5 ' —-— Point 4-measured 4 —-— Point 7measured 3_5 ' + Pomt 10—measured 3 Pressure(m pa) [0 U1 Time(ms) (a) Measured Simulated —— Point 4-FLUENT-Finite — Point 7-FLUENT-Finite —— Point 10-FLUENT-Flnite Time(ms) (b) Simulated Figure A.8 P6 108 A mesh with 2500 quadrilateral cells was generated in GAMBIT. The inlet boundary was set to be a mass flow inlet boundary with 0.0001537 kg/s flow rate. The outlet was set as the pressure outlet boundary. A species transport model with finite rate reaction model [54] was employed in the solver FLUENT. The simulated pressure was compared to the measured data at three points. They are 2450 mm (point 4), 3250 mm (point 7) and 4500 mm (point 10) from the ignition point respectively as shown in Figure A.8. At point 4, the simulated pressure wave arriving time agrees well with the measured data. But the simulated wave speed after point 4, which is around 320 m/s, is much higher than the measured wave speed, which is around 200 m/s. And the simulated over pressure of 0.09 MPa is much lower than the measured over pressure of 4.5 MPa, although similar attenuation of the pressure wave was observed in the simulation. As mentioned in reference [52], the deflagration-to- detonation transition simulation still remained a challenging problem. A more accurate combustion model may be needed to be incorporated to the FLUENT solver other than the conventional finite rate combustion model. 109 APPENDIX B Programs 3.1 User Defined Functions (UDF) UDFs are routines written in C, that can be dynamically loaded with FLUENT solver to enhance the standard feature of the code. 8.1.1 UDF for Piston Motion Calculation /*********‘k*ir**************~k***********‘k'k***~k‘k*********** **** * UDF for calculating the motion of the left side of the piston and adjusting the mesh* ****‘k**‘k*************~k*********************************k* ***/ #include "udf.h" static real v_prev = 0.0; DEFINE_CG_MOTION(piston,dt,vel,omega,time,dtime) { Domain *domain= Get_Domain (1); int IDl=7; /* ID for the left side of the piston, from the boudary setting in FLUENT */ int ID2=6; /* ID for the right side of the piston, from the boudary setting in FLUENT */ Thread *t=Lookup_Thread(domain, IDl); Thread *tl=Lookup_Thread(domain, ID2); face_t f; real NV_VEC(A); real force, dv, forcel, forceO, y, yl; real X[ND_ND], X1[ND_ND]; /* reset velocities */ NV_S(Vel, =, 0.0); NV_S(omega, =, 0.0); f=l; v_prev = F_U(f,tl); if (!Data_Valid_P()) return; /* reset velocities */ NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); if (!Data_Valid_P()) HO return; /* compute pressure force on the left side of the piston by looping through all faces */ force = 0.0; begin_f_loop(f,t) { F_AREA(A,f,t); F_CENTROID(X,f,t); force += F_P(f,t) * NV_MAG(A); } end_f_loop(f,t) /* compute pressure force on the right side of the piston by looping through all faces */ forcel = 0.0; begin_f_loop(f,tl) { F_AREA(A,f,tl); F_CENTROID(X1,f,tl); forcel += F_P(f,tl) * NV_MAG(A); } end_f_loop(f,tl) /* compute change in velocity, i.e., dv = F * dt / mass velocity update using explicit Euler formula */ forceO = 0.0; force0=(force—forcel)/0.08*3.1415926*0.08*0.08/4.0; /*calculate the average pressure and then times area to get pressure force*/ dv = dtime * (force0)/2.0; v;prev += dv; Message ("time = %f, x_vel = %f, force = %f, forcel = %f, forceO = %f\n", time, v_prev, force, forcel,force0); /* set x—component of velocity */ vel[0] = v_preV; } /******it*********~k********i'****‘k************************* **** * UDF for calculating the motion of the right side of the piston and adjusting the mesh* 1” need assign, a new 'variable: v_prevl other than the v_prev used in left side ************~k*****************~k-k****‘k*******************‘k ***/ #include "udf.h" static real v_prevl = 0.0; DEFINE_CG_MOTIONipiston,dt,vel,omega,time,dtime) { Domain *domain= Get_Domain (1); int ID1=7; /* ID for the left side of the piston, found from the boudary setting in FLUENT */ int ID2=6; /* ID for the right side of the piston, found from the boudary setting in FLUENT */ Thread *t=Lookup_Thread(domain, IDl); Thread *tl=Lookup_Thread(domain, ID2); face_t f; real NV_VEC(A); real force, dv, forcel, forceO, y, yl; real x[ND_ND], xl[ND_ND]; /* reset velocities */ NV_S(vel, =, 0.0); NV_S(omega, =, 0.0); f=l; v_prev = F_U(f,tl); if (lData_Valid_P()) return; /* reset velocities */ NV_S(Vel, =, 0.0); NV_S(omega, =, 0.0); if (iData_Valid_P()) return; /* compute pressure force on the left side of the piston by looping through all faces */ force = 0.0; begin_f_loop(f,t) { F_AREA(A,f,t); F_CENTROID(X,f,t); force += F_P(f,t) * NV_MAG(A); H2 end_f_loop(f,t) /* compute pressure force on the right side of the piston by lOOping through all faces */ forcel = 0.0; begin_f_loop(f,tl) { F_AREA(A,f,tl); F_CENTROID(Xl,f,tl); forcel += F_P(f,tl) * NV_MAG(A); } end_f_loop(f,tl) /* compute change in velocity, i e., dv = F * dt / mass velocity update using explicit Euler formula */ forceO = 0.0; forceO=(force—forcel)/0.08*3.1415926*0.08*0.08/4.0; /*calculate the average pressure and then times area to get pressure force*/ dv = dtime * (forceO)/2.0; v_prev += dv; Message ("time = %f, x_vel = %f, force = %f, forcel = %f, forceO = %f\n", time, v_prev, force, forcel,force0); /* set x—component of velocity */ vel[0] = v_prevl; } 3.1.2 UDF for The Diaphragm Rupture Process /************************‘k******************************* ***'k * UDF for calculating the radius of the burst hole and adjusting the mesh * ***************************************~k***************** *‘k'k/ #include "udf.h" DEFINE_GEOM(Dlposition,domain,dt,position) { real LT, LD, LA, LTl,LT2, LRT, LTT, Dr2; LT2=CURRENT_TIME; LT=0.0003; /*rupture process duration*/ LD=0.048; /*diaphragm diameter 0.05 -— 0.002 to avoid the diaphragm contact the shocktube upper and bottom boundary in the end*/ H3 LRT=0.09429 /*+0.0000403; /*0.0943, time of rupture; 0.0000403 time required to have the initial 1mm gap in diaphragm*/ LTT=0.09429+0.0003; /*ending time of the burst process*/ if ((LT2 > LRT) && (LT2< LTT)) { LTl=LT2—0.09429; LA=(1-COS(3.1415926/2*(LTl/LT)*(LTl/LT)))*LD*LD; Dr2=(l-COS(3.1415926/2*(LTl/LT)*(LTl/LT)))*LD*LD; /* Dr“2 from Equation 4.1 */ position[l] =0.5*pow(Dr2,0.5); /* Dr/2 */ } } 8.1.3 UDF for Pressure Characteristics Data and Driving-Application Gas Interface Data /********************‘k**‘k******‘k*‘k********‘k***‘k********** **** * UDF for writing pressure-time data along the centerline for characteristics * ******‘kkir‘k‘k************‘k‘k‘kir'k‘kic'kir‘k‘k‘kkir‘k‘k****************** ***/ #include "udf.h" #include "stdio.h" FILE *fp; DEFINE_EXECUTE_AT_END(Chracteritics) { Domain *domain= Get_Domain (l); face_t f; real time,p; real X[ND_ND]; int. ID=5; /* II) for time centerline, defined. in boundary from FLUENT*/ Thread *t=Lookup_Thread(domain, ID); time=CURRENT_TIME*lOOO; fp = fopen("chracteristic.dat","a"); begin_f_loop(f, t) /* loops over faces in a face thread */ { F_CENTROID(X,f,t); H4 p=F_P(f,t)/1000000; /*Message (“time = %f, pressure 2 %f", time, p);*/ if (fabs(x[l]) FILE *fp; DEFINE_EXECUTE_AT_END(Chracteritics) { Domain *domain= Get_Domain (1); cell_t c,c0,c1; face_t f; Thread *t0,*tl; real time,te; real Xl[ND_ND],X2[ND_ND],X[ND_ND]; int ID=5; /* ID for the centerline of the rectangle domain*/ Thread *t=Lookup_Thread(domain, ID); time=CURRENT_TIME*lOOO; fp = fopen("chracteristic.dat","a"); begin_f_loop(f, t) /* loops over faces in the face thread, which is the centerline of the rectangle domain */ { C_CENTROID(X,f,t); c0 = F_C0(f,t); t0 = THREAD_TO(t); C_CENTROID(Xl,cO,tO); cl = F_Cl(f,t); tl = THREAD_Tl(t); H5 C_CENTROID(X2,Cl,tl); te = CLJWCO,t0); /* cells that are below the center line*/ if (fabs(x[l]) SOUTPUT 2>&l Redirect program output to a file in your home directory. fluent -g SFLUENTSOLVER —i SINPUT > $OUTPUT 2>&l An example for the journal file from the solid-gas interaction case: file read—case—data O603.cas quit solve/set/time-step le-7 solve/dual—time—iterate 20 15 quit solve/set/time-step H8 5e—7 solve/dual-time-iterate 400 15 quit solve/set/time—step le-6 solve/dual-time-iterate 500 15 quit solve/set/time—step 2e—6 solve/dual-time-iterate 500 15 quit solve/set/time—step 5e—6 solve/dual—time-iterate 3000 15 quit file/write-case—data/0529final.cas exit 8.3 JIANG Codes for One Dimensional Shock Tube Calculation Two JIANG codes written in Fortran were employed in this research. One is for the MSU simple MSU shock tube. The other is for the compression process for the MSU piston assisted shock tube. 8.3.1 JIANG Code for The Simple Shock Tube DIMENSION PP(21,2010),X(2010),PN(2100),X1(2010),Ul(3,2010), * X2(2010),T(2100),Tl(2101),XH(2101),TN(2100),PZ(2101),TZ(2 101), * V(3,4),FE(3,4),Vl(3,4),FEl(3,4),VA(3,4) COMMON/A/W(2010),P(2010),A(2010),U(3,2010),FP(3,2010), * FM(3,2010),GA1,GA2,N1,N2 OPEN(2,FILE='CS.DAT') READ(2,*)ALl,PR,TR,GA1,CM1,E1,AL2,PX,TX,E2,PS,GT,TK CLOSE(2) C ALl-—QDONG GUANCHANG PFr—QDONG YALI TR--QDONG WENDU C GAl-—QDONG BIREBI CMl--QDONG FENZILIANG El-—GZ ZHIJING H9 C AL2~—GONGZUO GUANCHANG PX--GZ QITI YALI TX--GZ QITI WENDU C YALI C 10 E2--GZ ZHIJING BM--HUOSAI ZHILIANG PS—-PEMO GT——JISUAN SHIJIAN Rl=8315./CM1 R2=286.7241379 GA2=1.4 FX=(AL1+AL2)*0.0005 Nl=ALl/FX+3.5 Fl=ALl/(Nl-3.) N2=AL2/FX+3.5 F2=AL2/(N2-3.) N2=N2+Nl NZ=(AL1+AL2)/FX+1 Sl=0.785398l634*El**2 S2=O.7853981634*E2**2 SB=SZ/Sl SS=O.5787037037 GM=O.4 GM1=SS*SB*(1.+O.5*(GA2—l.)*GM**2)**3 XA=ABS(GMl-GM) GM=GM1 IF(XA.GT.0.00000l) GOTO l KB=-2 HX=AL1 DO 10 I=1,4001 X(I)=FX*(I-l.) CONTINUE AR=SQRT(GA1*R1*TR) AX=SQRT(GA2*R2*TX) XM1=1.0+0.165*SQRT((PR—PX)/PX) V01=XM1*AX YM=SQRT((XMl**2+2./(GA2—l.))/(2.*GA2/(GA2- l.)*XMl**2-l.)) PY=PX*(2.*GA2*XM1**2—GA2+1.)/(GA2+1.) TY=TX*2.*(GA2—l.)*(1.+.5*(GA2- l.)*XM1**2)*(2.*GA2/(GA2-l.)* * XMl**2-l.)/((GA2+l.)*XM1)**2 AY=SQRT(GA2*R2*TY) V02=YM*AY VY=VOl-V02 AYR=AR-.5*(GAl—l.)*VY TYR=AYR**2/(GA1*R1) PZ1=PY*(TR/TYR)**(GAl/(GAl-l.)) XM2=1.0+O.17*SQRT((PR-PX)/PX) VOl=XM2*AX 120 YM=SQRT((XM2**2+2./(GA2-l.))/(2.*GA2/(GA2- l.)*XM2**2-l.)) PY=PX*(2.*GA2*XM2**2—GA2+1.)/(GA2+1.) TY=TX*2.*(GA2-1.)*(l.+.5*(GA2- l.)*XM2**2)*(2.*GA2/(GA2-l.)* 101 * XM2**2-l.)/((GA2+1.)*XM2)**2 AY=SQRT(GA2*R2*TY) V02=YM*AY VY=VOl-V02 AYR=AR-.5*(GAl-l.)*VY TYR=AYR**2/(GA1*R1) PZ2=PY*(TR/TYR)**(GAl/(GAl—l.)) XM3=XM1+(PR—PZl)/(PZZ—PZ1)*(XM2—XM1) XM1=XM2 PZl=PZ2 XM2=XM3 VOl=XM2*SQRT(GA2*R2*TX) YM=SQRT((XM2**2+2./(GA2-l.))/(2.*GA2/(GA2- l.)*XM2**2-l.)) PY=PX*(2.*GA2*XM2**2—GA2+1.)/(GA2+l.) TY=TX*2.*(GA2-l.)*(l.+.5*(GA2- l.)*XM2**2)*(2.*GA2/(GA2-l.)* 20 * XM2**2-l.)/((GA2+1.)*XM2)**2 AY=SQRT(GA2*R2*TY) V02=YM*AY VY=V01-V02 AYR=AR-.5*(GAl—1.)*VY TYR=AYR**2/(GA1*R1) PZZ=PY*(TR/TYR)**(GAl/(GAl—l.)) XXl=ABS((PR-PZ2)/PR) IF(XX1.GT.0.000003) GOTO 101 U(l,Nl)=PY/(R1*TYR) U(2,Nl)=U(l,Nl)*VY U(3,Nl)=U(l,Nl)*(AYR**2/GA1/(GAl—l.)+.5*VY**2) U(l,Nl+l)=PY/(R2*TY) U(2,N1+1)=U(1,N1+1)*VY U(3,Nl+1)=U(1,Nl+l)*(AY**2/GA2/(GA2-l.)+.5*VY**2) AR=SQRT(GA1*R1*TR) D=PR/(R1*TR) DO 20 I=1,N1—1 U(l,I)=D U(2,I)=0.0 U(3,I)=D*(AR**2/(GA1*(GA1-l.))) X1(I)=(I—3.)*F1 CONTINUE AR=SQRT(GA2*R2*TX) D=PX/(R2*TX) 121 21 6O 18 ll 12 1)+FM(J, * DO 21 I=N1+2,N2 U(l,I)=D U(2,I)=0.0 U(3,I)=D*(AR**2/(GA2*(GA2—l.))) X1(I)=AL1+(I-Nl-l.)*F2 CONTINUE GT1=GT/2100. M=TK/GT1 DO 60 I=1,210l T1(I)=(I-l.)*GT1 CALL FF PZ(1)=P(N1+1) TZ(1)=A(N1+1)**2/(GA2*R2) DO 30 Kl=1,2100 HT=0.0 LK=—2 FT=Fl/(A(l)+ABS(W(l)))*O.7O DO 11 I=2,N1 XA=Fl/(A(I)+ABS(W(I)))*O.7O IF(XA.LT.FT) FT=XA CONTINUE DO 12 I=Nl+l,N2-2 XA=F2/(A(I)+ABS(W(I)))*O.7O IF(XA.LT.FT) FT=XA CONTINUE IF(HT+FT.GE.GT1) THEN LK=2 FT=GTl—HT ENDIF HT=HT+FT DO 36 I=3,Nl-2 DO 74 J=l,3 Ul(J,I)=U(J,I)-FT/Fl*(FP(J,I)—FP(J,I- I+l)—FM(J,I)+ 0.5*(XL(FP(J,I—l),FP(J,I),FP(J,I+1))—XL(FP(J,I- I—l), FP(J,I))— I),FM(J,I+1),FM(J,I+2))+XL(FM(J,I-l), FM(J,I),FM(J,I+1)))) CONTINUE CONTINUE I=Nl—l DO 76 J=l,3 U1(J,I)=U(J,I)-FT/Fl*(FP(J,I)—FP(J,I- I+l)—FM(J,I)+ O.5*(XL(FP(J,I-l),FP(J,I),FP(J,I+1))—XL(FP(J,I- 122 2),FP(J,I—l), * FP(J,I) ) )) 76 CONTINUE I=Nl+2 DO 78 J=l,3 Ul(J,I)=U(J,I)-—FT/F1*(FP(J,I)—FP(J,I— 1)+FM(J,I+1)—FM(J,I)+ * O.5*(XL(FM(J,I—l),FM(J,I),FM(J,I+1))— XL(FM(J,I),FM(J,I+1), * FM(J,I+2)))) 78 CONTINUE C ________________________________________ CALL ELA(U(1,Nl—l),U(2,N1-l),U(3,N1- 1),V(1,1),V(2,l),V(3,1)) CALL ELA(U(1,N1),U(2,N1),U(3,Nl),V(1,2),V(2,2),V(3,2)) CALL ELA(U(1,N1+1),U(2,N1+1),U(3,Nl+l),V(l,3),V(2,3),V(3,3)) CALL ELA(U(1,N1+2),U(2,N1+2),U(3,Nl+2),V(l,4),V(2,4),V(3,4)) CALL LAF(V(1,1),V(2,1),V(3,l),GA1,FE(1,1),FE(2,1),FE(3,1)) CALL LAF(V(1,2),V(2,2),V(3,2),GAl,FE(1,2),FE(2,2),FE(3,2)) CALL LAF(V(1,3),V(2,3),V(3,3),GA2,FE(1,3),FE(2,3),FE(3,3)) CALL LAF(V(1,4),V(2,4),V(3,4),GA2,FE(1,4),FE(2,4),FE(3,4)) XL1=F1*O.5*(l./V(l,l)+1./V(l,2)) XL3=F2*O.5*(l./V(l,3)+l./V(l,4)) DO 110 I=l,3 Vl(I,l)=V(I,l)-FT/XL1*(FE(I,2)—FE(I,1)) 110 Vl(I,3)=V(I,3)-FT/XL3*(FE(I,4)—FE(I,3)) CALL YT(V(1,2),V(2,2),V(3,2),Vl(l,3),Vl(2,3),V1(3,3), * GAl,Rl,GA2,Vl(l,2),V1(2,2),Vl(3,2)) XL1=(F1+FT*O.5*(V(2,2)+V1(2,2)-V(2,l)—V1(2,1)))* * O.5*(l./Vl(l,1)+l./Vl(1,2)) DO 111 K=1,2 111 CALL LAF(V1(1,K),V1(2,K),Vl(3,K),GA1,FEl(l,K),FE1(2,K) * ,FEl(3,K)) DO 112 I=l,3 112 VA(I,2)=O.5*(V(I,2)+Vl(I,2)—FT/XL1*(FE1(I,2)— FE1(I,1))) CALL YT(V(1,3),V(2,3),V(3,3),VA(1,2),VA(2,2),VA(3,2), 123 * GA2,R2,GA1,VA(1,3),VA(2,3),VA(3,3)) CALL LAE(VA(1,2),VA(2,2),VA(3,2),U1(1,Nl),U1(2,N1),U1(3,N1)) 1)+FM(J, 2),FP(J, XL(FM(J, 80 38 82 'k ‘k * * CALL LAE(VA(1,3),VA(2,3),VA(3,3),U1(1,N1+1), Ul(2,Nl+l),Ul(3,N1+l)) HX=O.5*FT*(U(2,Nl)/U(1,N1)+Ul(2,Nl)/Ul(l,Nl))+HX DO 38 I=N1+3,N2—2 DO 80 J=l,3 Ul(J,I)=U(J,I)-FT/Fl*(FP(J,I)—FP(J,I- I+l)-FM(J,I)+ O.5*(XL(FP(J,I-1),FP(J,I),FP(J,I+1))-XL(FP(J,I- I-l). FP(J,I))- I),FM(J,I+1),FM(J,I+2))+XL(FM(J,I—l), FM(J,I),FM(J,I+1)))) CONTINUE CONTINUE DO 82 I=1,N1-1 X2(I)=X1(I) X2(Nl)=HX X2(Nl+l)=HX DO 84 I=Nl+2,N2 X2(I)=X1(I) C ________________________________________________________ 42 43 86 N3=HX/FX+3.5 Fl=HX/(N3-3.) N4=(AL1+AL2—HX)/FX+3.5 IF(N4.LE.7) THEN TS=(K1—1.)*GT1+HT WRITE(*,*)'GT should small or eq TT' WRITE(*,*)'GT 6:,Aaaéfi»ousou TT' WRITE(*,*)'TT=' WRITE(*,*)TS GOTO 23 ENDIF F2=(AL1+AL2—HX)/(N4-3.) N4=N4+N3 DO 42 I=1,N3 X1(I)=(I—3.)*F1 DO 43 I=N3+1,N4 X1(I)=(I-N3—l.)*F2+HX DO 86 J=l,3 U(J,3)=U1(J,3) U(J,N3)=U1(J,Nl) U(J,N3+1)=U1(J,N1+1) U(J,N4—2)=Ul(J,N2—2) 124 39 41 88 37 45 46 9O 44 NH=3 DO 37 I=4,N3—1 DO 39 K2=NH,N1 IF(X2(K2).GE.X1(I)) THEN JJ=K2-l NH=K2-1 GOTO 41 ENDIF CONTINUE DO 88 J=l,3 U(J,I)=Ul(J,JJ)+(X1(I)—X2(JJ))/(X2(JJ+1)—X2(JJ)) *(Ul(J,JJ+1)-U1(J,JJ)) CONTINUE CONTINUE NH=Nl+l DO 44 I=N3+2,N4—3 DO 45 K2=NH,N2-3 IF(X2(K2).GE.X1(I)) THEN JJ=K2-1 NH=K2—1 GOTO 46 ENDIF CONTINUE DO 90 J=l,3 U(J,I)=Ul(J,JJ)+(X1(I)—X2(JJ))/(X2(JJ+l)—X2(JJ)) *(Ul(J,JJ+l)—U1(J,JJ)) CONTINUE CONTINUE N1=N3 N2=N4 U(l,2)=U(l,4) U(2,2)=—U(2,4) U(3,2)=U(3,4) U(l,l)=U(l,5) U(2,l)=-U(2,5) U(3,1)=U(3,5) IF(KB.LT.O) THEN U(l,N2-l)=U(l,N2-3) U(2,N2-1)=-U(2,N2-3) U(3,N2—l)=U(3,N2—3) U(1,N2)=U(1,N2—4) U(2,N2)=—U(2,N2-4) U(3,N2)=U(3,N2—4) ELSE D=U(l,N2—2) Vn=U(2,N2-2)/D PJ=(U(3,N2—2)—.5*D*Vn**2)*(GA2-1.) 125 AJ=SQRT((U(3,N2—2)/D-.5*Vn**2)*GA2*(GA2-1.)) TJ=AJ**2/R2/GA2 AMJ=Vn/AJ POJ=PJ*(1.+.5*(GA2—l.)*AMJ**2)**(GA2/(GA2—l.)) TOJ=TJ*(1.+.5*(GA2-1.)*AMJ**2) P(N2-2)=POJ/(l.+.5*(GA2—l.)*GM**2)**(GA2/(GA2-1.)) TJ=TOJ/(1.+.5*(GA2—l.)*GM**2) A(N2-2)=SQRT(GA2*R2*TJ) W(N2-2)=GM*A(N2-2) U(l,N2-2)=P(N2—2)/(R2*TJ) U(2,N2-2)=U(l,N2-2)*W(N2-2) U(3,N2—2)=U(1,N2-2)*(A(N2-2)**2/GA2/(GA2- l.)+.5*W(N2-2)**2) 9 2 DO 92 J=l,3 U(J,N2-l)=U(J,N2—2) U(J,N2)=U(J,N2-2) CONTINUE ENDIF CALL FF IF(KB.LE.O.AND.P(NZ—Z).GE.PS) THEN KB=2 WRITE(*,*)' ----- membrane break ————— ' ENDIF C ________________________________________________________ IF(LK.LE.0) GOTO 18 XH(K1+1)=HX PZ(K1+1)=P(N1+1) TZ(K1+1)=A(N1+1)**2/(GA2*R2) IF(K1.GT.M) THEN MN=Kl-M PN(MN)=P(N2-2) TN(MN)=A(N2-2)**2/(GA2*R2) T(MN)=K1*GT1 ENDIF LC=k1/20. kkb=LC*20 if(kkb.eq.kl) then SS1=u(2,nl)/u(l,n1) SS2=P(N2-2)*.000001 SSB=A(N2-2)**2/GA2/R2 SS4=P(N1+1)*.OOOOOl SSS=A(N1+1)**2/GA2/R2 WRITE(*,lOO)Kl,HX,SS1,SSZ,SS3 WRITE(*,lOZ)ss4,ssS ENDIF NKzKl/lOO. 126 kkb=NK*lOO if(kkb.eq.kl) then PP(NK,1)=P(3) PP(NK,NZ) =P(N2—2) NH=3 DO 14 I=2,NZ-1 DO 15 K2=NH,N2-2 IF(Xl(K2).GE.X(I)) THEN JJ=K2-1 NH=K2—1 GOTO 16 ENDIF 15 CONTINUE 16 PP(NK,I)=P(JJ)+(X(I)—X1(JJ))/(Xl(JJ+l)— X1(JJ))*(P(JJ+l)-P(JJ)) 14 CONTINUE endif 30 CONTINUE NN=2100—M NN1=2101 WRITE(*,118)GM 100 FORMAT(3X,2HJ=,IS,3X,2HX=,G9.4,lX,lX,3HUj=,GlO.3,1X,3HPN=, * GlO.3,lX,3HTN=,GlO.4) 102 FORMAT(3X,3hpj=,glO.4,1X,3htj=,g10.4) 118 FORMAT(3X,3HGM=,G16.7) OPEN(5,FILE='PN.DAT') DO 29 I=1,NN WRITE(5,*)T(I) ,PN(I) ,TN(I) 29 CONTINUE CLOSE(5) OPEN(5,FILE='PJ.DAT') DO 89 I=1,NN1 WRITE(5,*)T1(I),PZ(I),TZ(I) 89 CONTINUE CLOSE(5) OPEN(5,FILE='PIS.DAT') DO 61 I=1,NN1 WRITE(5,*)T1(I),XH(I) 61 CONTINUE CLOSE(5) OPEN(5,FILE='P1.DAT') DO 28 I=1,NZ WRITE(5,*)X(I),PP(1,I),PP(2,I),PP(3,I) 28 CONTINUE CLOSE(5) OPEN(5,FILE='P2 .DAT') 127 DO 40 I=1,NZ WRITE(5,*)X(I),PP(4,I),PP(5,I),PP(6,I) 4O CONTINUE CLOSE(5) OPEN(5,FILE='P3.DAT') DO 50 I=1,NZ WRITE(5,*)X(I),PP(7,I),PP(8,I),PP(9,I) 50 CONTINUE CLOSE(5) OPEN(5,FILE='P4.DAT') DO 51 I=1,NZ WRITE(5,*)X(I),PP(10,I),PP(11,I),PP(12,I) 51 CONTINUE CLOSE(5) OPEN(5,FILE='P5.DAT') DO 52 I=1,NZ WRITE(5,*)X(I),PP(13,I),PP(14,I),PP(15,I) 52 CONTINUE CLOSE(5) OPEN(5,FILE='P6.DAT') DO 53 I=1,NZ WRITE(5,*)X(I),PP(16,I),PP(17,I),PP(18,I) 53 CONTINUE CLOSE(5) OPEN(5,FILE='P7.DAT') DO 54 I=1,NZ WRITE(5,*)X(I),PP(19,I),PP(20,I),PP(21,I) 54 CONTINUE CLOSE(5) 23 STOP END SUBROUTINE FF COMMON/A/W(2010),P(2010),A(2010),U(3,2010),FP(3,2010), * FM(3,2010),GA1,GA2,N1,N2 DO 10 I=1,NZ IF(I.LE.N2) THEN GA=GA1 ELSE GA=GA2 ENDIF D=U(1,I) W(I)=U(2,I)/D V=W(I) P(I)=(U(3,I)—.5*D*V**2)*(GA—l.) A(I)=SQRT((U(3,I)/D—.5*V**2)*GA*(GA—l.)) F1=U(2,I) F2=D*V**2+P(I) F3=(U(3,I)+P(I))*V X1=ABS(V+A(I)) X2=ABS(V-A(I)) El=D*ABS(V)-D/GA*(ABS(V)—.5*x1—.5*X2) E2=D*V*ABS(V)-D/GA*(V*ABS(V)-.5*(V+A(I))*Xl— .5*(V—A(I))*X2) E3=D*(GA-1.)/GA*.5*V**2*(ABS(V)—.5*x1— .5*X2)+.5*D*(A(I)**2/ * GA/(GA—l.)+.5*V**2)*(X1+X2)+.5*D*V*A(I)/GA*(Xl—X2) FP(1,I)=.5*(F1+E1) FM(1,I)=.5*(Fl-El) FP(2,I)=.5*(F2+E2) FM(2,I)=.5*(F2-E2) FP(3,I)=.5*(F3+E3) FM(3,I)=.5*(F3-E3) 10 CONTINUE RETURN END FUNCTION XL(X1,X2,X3) Y1=X2-Xl Y2=X3-X2 IF(Y1*Y2.LT.0.0) THEN FM=0.0 ELSE Zl=ABS(Y1) Z2=ABS(Y2) IF(Z1.LE.Z2) THEN FM=Y1 ELSE FM=Y2 ENDIF ENDIF RETURN END SUBROUTINE YT(U11,U12,U13,U21,U22,U23,GA1,R1,GA2,U1, 129 U2,U3) Dl=l./Ull Vl=U12 Pl=(Ul3-.5*Vl**2)*(GAl-l.)*Dl Tl=Pl/(R1*D1) D2=l./U21 V2=U22 P2=(U23-.5*V2**2)*(GA2—l.)*D2 CALL TP(P1,T1,P2,GA1,T) D=P2/(R1*T) Ul=l./D U2=V2 U3=R1*T/(GAl—l.)+.5*V2**2 RETURN END SUBROUTINE ELA(U1,U2,U3,V1,V2,V3) Vl=l./Ul V2=U2/U1 V3=U3/Ul RETURN END SUBROUTINE LAE(V1,V2,V3,U1,U2,U3) Ul=l./Vl U2=V2*Ul U3=V3*U1 RETURN END SUBROUTINE LAF(Vl,V2,V3,GA,F1,F2,F3) P=(GA—l.)/Vl*(V3-.5*V2**2) Fl=-V2 F2=P F3=P*V2 RETURN END SUBROUTINE TP(PX,TX,PY,GA,TY) 130 PG=PY/PX IF(PG.GT.1.) THEN TY=TX*PG*((GA-l.)*PG+GA+1.)/((GA+1.)*PG+GA-l.) ELSE TY=PG**((GA-l.)/GA)*TX ENDIF RETURN END 8.3.2 JIANG Code for The Compression Process within the MSU Piston Assisted Shock Tube DIMENSION PP(21,2010),X(2010),PN(2100),X1(20lO),Ul(3,2010), * X2(2010),T(2100),T1(2101),XH(2101),UH(2101),AH(2101),TN(2 100), * PJ(2101),TJ(2101) COMMON/A/W(2010),P(2010),A(2010),U(3,2010),FP(3,2010), * FM(3,2010),GA1,GA2,N1,N2 OPEN(2,FILE='CS.DAT') READ(2,*)ALl,PR,TR,GA1,CMI,E1,AL2,PX,TX,E2,PS,GT,TK,BM CLOSE(2) C ALl—-QDONG GUANCHANG PFr-QDONG YALI TR——QDONG WENDU C GAl-—QDONG BIREBI CMl--QDONG FENZILIANG Efl--GZ ZHIJING C AL2--GONGZUO GUANCHANG PX--GZ QITI YALI TX—-GZ QITI WENDU C E2--GZ ZHIJING BM-—HUOSAI ZHILIANG PS--PEMO YALI C GT--JISUAN SHIJIAN Rl=8315./CM1 R2=286.7241379 GA2=1.4 FX=(AL1+AL2)*0.00050 Nl=ALl/FX+3.5 Fl=ALl/(Nl-3.) N2=AL2/FX+3.5 F2=AL2/(N2—3.) N2=N2+N1 131 10 20 21 6O 18 NZ=(AL1+AL2)/FX+1 Sl=0.785398l634*El**2 SZ=O.7853981634*E2**2 SB=S2/Sl SS=O.5787037037 GM=O.4 GM1=SS*SB*(1.+O.5*(GA2-l.)*GM**2)**3 XA=ABS(GM1-GM) GM=GM1 IF(XA.GT.0.00000I) GOTO 1 KB=-2 HU=0.0 HX=AL1 DO 10 I=1,4001 X(I)=FX*(I-l.) CONTINUE AR=SQRT(GA1*R1*TR) D=PR/(R1*TR) DO 20 I=1,Nl U(l,I)=D U(2,I)=0.0 U(3,I)=D*(AR**2/(GA1*(GAl-l.))) X1(I)=(I-3.)*F1 CONTINUE AR=SQRT(GA2*R2*TX) D=PX/(R2*TX) DO 21 I=Nl+l,N2 U(l,I)=D U(2,I)=0.0 U(3,I)=D*(AR**2/(GA2*(GA2—l.))) X1(I)=AL1+(I—Nl-l.)*F2 CONTINUE HA=Sl*(PR-PX)/BM AH(1)=HA UH(l)=HU XH(l)=ALl GT1=GT/2100. M=TK/GT1 DO 60 I=l,2101 T1(I)=(I-l.)*GT1 CALL FF PJ(1)=P(N1+1) TJ(l)=A(Nl+l)**2/GA2/R2 DO 30 Kl=l,2100 HT=0.0 LK=-2 FT=F1/(A(l)+ABS(W(l)))*O.7O 132 11 12 DO 11 I=2,N1 XA=Fl/(A(I)+ABS(W(I)))*0.70 IF(XA.LT.FT) FT=XA CONTINUE DO 12 I=Nl+l,N2-2 XA=F2/(A(I)+ABS(W(I)))*O.70 IF(XA.LT.FT) FT=XA CONTINUE IF(HT+FT.GE.GT1) THEN LK=2 FT=GTl—HT ENDIF HT=HT+FT DO 36 I=3,N1—2 DO 74 J=l,3 Ul(J,I)=U(J,I)-FT/Fl*(FP(J,I)-FP(J,I- l)+FM(J,I+l)-FM(J,I)+ * O.5*(XL(FP(J,I-l),FP(J,I),FP(J,I+l))-XL(FP(J,I- )IFP JII_1) I 2 ( * XL(FM( * 74 36 FP(J,I))- J,I),FM(J,I+1),FM(J,I+2))+XL(FM(J,I-l), FM(J,I),FM(J,I+1)))) CONTINUE CONTINUE I=Nl-l DO 76 J=l,3 Ul(J,I)=U(J,I)—FT/Fl*(FP(J,I)-FP(J,I- l)+FM(J,I+l)-FM(J,I)+ * 0.5*(XL(FP(J,I-l),FP(J,I),FP(J,I+1))-XL(FP(J,I- 2),FP(J,I-l), * 76 FP(J,I)))) CONTINUE I=Nl+2 DO 78 J=l,3 Ul(J,I)=U(J,I)-FT/Fl*(FP(J,I)-FP(J,I- l)+FM(J,I+l)-FM(J,I)+ * O.5*(XL(FM(J,I—1),FM(J,I),FM(J,I+1))- XL(FM(J,I),FM(J,I+1), * C ______ FM(J,I+2)))) CONTINUE UB=HU+FT*HA BLl=FT*Fl*(4.*A(N1)+(GAl-l.)*(W(Nl)- UB))/(4.*Fl+(GAl-l.) * *FT*(W(Nl)-W(N1—l))+4.*FT*(A(Nl)—A(N1—1))) A3=A(Nl)+BLl/Fl*(A(N1—1)-A(Nl)) 133 UB))/(4. 1|: UB))/(4. UB))/(4. U3=W(Nl)+BLl/Fl*(W(N1-1)—W(Nl)) Cl=A3+.5*(GAl—l.)*U3 AQl=C1—O.5*(GAl-l.)*UB CALL AP(A(N1),AQl,GAl,P(Nl),PQl) BL2=FT*F2*(4.*A(N1+l)-(GA2-1.)*(W(N1+l)— *F2+(GA2 —l.)*FT*(W(N1+2)—W(N1+1))+4.*FT*(A(N1+l)—A(N1+2))) A2=A(Nl+1)+BL2/F2*(A(Nl+2)—A(N1+l)) U2=W(N1+l)+BL2/F2*(W(N1+2)—W(N1+l)) C2=A2-.5*(GA2—l.)*U2 AQZ=C2+.5*(GA2—l.)*UB CALL AP(A(N1+1),AQ2,GA2,P(N1+1),PQ2) BA=Sl*(PQl—PQ2)/BM UB=HU+FT*(HA+BA)*O.5 HX=HX+O.5*FT*(HU+UB) BL1=FT*F1*(4.*A(N1)+(GA1—l.)*(W(Nl)— *Fl+(GA1-l.) *FT*(W(Nl)-W(Nl-l))+4.*FT*(A(Nl)—A(Nl-l))) A3=A(Nl)+BLl/Fl*(A(Nl-l)—A(Nl)) U3=W(Nl)+BLl/Fl*(W(Nl-l)-W(Nl)) Cl=A3+.5*(GAl—l.)*U3 AQl=C1-0.5*(GAl—l.)*UB WQ1=UB CALL AP(A(N1),AQ1,GA1,P(N1),PQ1) BL2=FT*F2*(4.*A(N1+1)—(GA2-1.)*(W(Nl+l)— *F2+(GA2 -l.)*FT*(W(N1+2)—W(N1+l))+4.*FT*(A(Nl+l)-A(Nl+2))) A2=A(Nl+1)+BL2/F2*(A(Nl+2)—A(Nl+l)) U2=W(N1+l)+BL2/F2*(W(N1+2)—W(N1+1)) C2=A2—.5*(GA2—l.)*U2 AQ2=C2+.5*(GA2-l.)*UB WQ2=UB CALL AP(A(N1+1),AQ2,GA2,P(N1+1),PQ2) HU=UB HA=S1*(PQl—PQ2)/BM TQl=AQ1**2/(R1*GA1) TQ2=AQ2**2/(R2*GA2) Ul(l,Nl)=PQl/(R1*TQ1) U1(l,N1+l)=PQ2/(R2*TQ2) Ul(2,N1)=Ul(l,Nl)*WQ1 Ul(2,N1+l)=U1(l,N1+l)*WQZ Ul(3,Nl)=U1(1,Nl)*(AQ1**2/(GA1*(GA1— l.))+.5*WQl**2) Ul(3,N1+1)=U1(l,N1+1)*(AQ2**2/(GA2*(GA2- l.))+.5*WQ2**2) C _______ DO 38 I=N1+3,N2—2 134 1)+FM 2),FP XL(FM 80 38 82 J. J. J, DO 80 J=l,3 Ul(J,I)=U(J,I)—FT/Fl*(FP(J,I)—FP(J,I— I+l)-FM(J,I)+ O.5*(XL(FP(J,I-l),FP(J,I),FP(J,I+1))-XL(FP(J,I— I-ll, FP(J,I))- I),FM(J,I+1),FM(J,I+2))+XL(FM(J,I-l), FM(J,I),FM(J,I+1)))) CONTINUE CONTINUE DO 82 I=1,N1—1 X2(I)=X1(I) X2(N1)=HX X2(Nl+l)=HX DO 84 I=Nl+2,N2 X2(I)=X1(I) C ________________________________________________________ 42 43 86 N3=HX/FX+3.5 Fl=HX/(N3—3.) N4=(AL1+AL2-HX)/FX+3.5 IF(N4.LE.7) THEN TS=(K1-l.)*GTl+HT WRITE(*,*)'GT should small or eq TT' WRITE(*,*)'GT O:,AD;OU»ouEOU TT' WRITE(*,*)'TT=' WRITE(*,*)TS GOTO 23 ENDIF F2=(AL1+AL2—HX)/(N4-3.) N4=N4+N3 DO 42 I=1,N3 X1(I)=(I—3.)*Fl DO 43 I=N3+1,N4 Xl(I)=(I-N3-l.)*F2+HX DO 86 J=l,3 U(J,3)=U1(J,3) U(J,N3)=U1(J,Nl) U(J,N3+1)=U1(J,N1+l) U(J,N4-2)=U1(J,N2-2) NH=3 DO 37 I=4,N3-1 DO 39 K2=NH,N1 IF(X2(K2).GE.X1(I)) THEN JJ=K2—l NH=K2-1 GOTO 41 135 39 41 88 37 45 46 9O 44 ENDIF CONTINUE DO 88 J=l,3 U(J,I)=Ul(J,JJ)+(X1(I)—X2(JJ))/(X2(JJ+l)-X2(JJ)) *(Ul(J,JJ+l)—U1(J,JJ)) CONTINUE CONTINUE NH=N1+1 DO 44 I=N3+2,N4—3 DO 45 K2=NH,N2—3 IF(X2(K2).GE.X1(I)) THEN JJ=K2-l NH=K2-l GOTO 46 ENDIF CONTINUE DO 90 J=l,3 U(J,I)=U1(J,JJ)+(X1(I)-X2(JJ))/(X2(JJ+1)-X2(JJ)) *(Ul(J,JJ+l)—U1(J,JJ)) CONTINUE CONTINUE N1=N3 N2=N4 U(l,2)=U(l,4) U(2,2)=—U(2,4) U(3,2)=U(3,4) U(l,l)=U(l,5) U(2,1)=-U(2,5) U(3,1)=U(3,5) IF(KB.LT.O) THEN U(l,N2-1)=U(l,N2—3) U(2,N2—l)=—U(2,N2—3) U(3,N2-1)=U(3,N2—3) U(l,N2)=U(l,N2-4) U(2,N2)=-U(2,N2-4) U(3,N2)=U(3,N2—4) ELSE D=U(l,N2-2) V=U(2,N2—2)/D PQ=(U(3,N2-2)-.5*D*V**2)*(GA2-l.) AJ=SQRT((U(3,N2—2)/D—.5*V**2)*GA2*(GA2—1.)) TQ=AJ**2/R2/GA2 AMJ=V/AJ POJ=PQ*(1.+.5*(GA2-1.)*AMJ**2)**(GA2/(GA2-l.)) TOJ=TQ*(l.+.5*(GA2-1.)*AMJ**2) P(N2—2)=POJ/(l.+.5*(GA2—l.)*GM**2)**(GA2/(GA2—l.)) TQ=TOJ/(l.+.5*(GA2—1.)*GM**2) 136 A(N2—2)=SQRT(GA2*R2*TQ) W(N2-2)=GM*A(N2—2) U(l,N2—2)=P(N2~2)/(R2*TQ) U(2,N2—2)=U(l,N2-2)*W(N2—2) U(3,N2—2)=U(1,N2—2)*(A(N2—2)**2/GA2/(GA2- 1.)+.5*W(N2-2)**2) DO 92 J=l,3 U(J,N2-1)=U(J,N2-2) U(J,N2)=U(J,N2-2) 92 CONTINUE ENDIF CALL FF IF(KB.LE.O.AND.P(N2—2).GE.PS) THEN KB=2 WRITE(*,*)' ————— membrane break ————— ' ENDIF C ________________________________________________________ IF(LK.LE.O) GOTO 18 XH(K1+1)=HX UH(K1+1)=HU AH(K1+1)=HA PJ(K1+1)=P(N1+1) TJ(K1+1)=A(N1+1)**2/(GA2*R2) IF(K1.GT.M) THEN MN=Kl-M PN(MN)=P(N2-2) TN(MN)=A(N2—2)**2/(GA2*R2) T(MN)=K1*GTl ENDIF LC=kl/20. kkb=LC*20 if(kkb.eq.k1) then SS=P(N2—2)*.000001 SSl=A(N2—2)**2/R2/GA2 SS2=P(N1+1)*.000001 SS3=A(N1+1)**2/R2/GA2 WRITE(*,lOO)K1,HX,HU,HA,SS WRITE(*,lOZ)SSl,SSZ,SS3 ENDIF NK=K1/100. kkb=NK*100 if(kkb.eq.k1) then PP(NK,1)=P(3) PP(NK,NZ)=P(N2-2) NH=3 DO 14 I=2,NZ—1 137 DO 15 K2=NH,N2—2 IF(Xl(K2).GE.X(I)) THEN JJ=K2—l NH=K2—1 GOTO 16 ENDIF 15 CONTINUE 16 PP(NK,I)=P(JJ)+(X(I)-x1(JJ))/(X1(JJ+1)- X1(JJ))*(P(JJ+1)—P(JJ)) 14 CONTINUE endif 3O CONTINUE NN=2100—M NN1=2101 WRITE(*,118)GM 100 FORMAT(3X,2HJ=,I5,1X,2HX=,G9.4, * 2HU=,G10.3,2HA=,G11.4,lx,3HPN=,G10.3) 102 FORMAT(3X,3HTN=,G10.3,lX,3HPj=,GlO.3,lX,3HTj=,GlO.3) 118 FORMAT(3X,3HGM=,G16.6) OPEN(5,FILE='PN.DAT') DO 29 I=1,NN WRITE(5,*)T(I),PN(I),TN(I) 29 CONTINUE CLOSE(5) OPEN(5,FILE='PJ.DAT') DO 89 I=1,NN1 WRITE(5,*)T1(I),PJ(I),TJ(I) 89 CONTINUE CLOSE(5) OPEN(5,FILE='PIS.DAT') DO 61 I=1,NN1 WRITE(5,*)T1(I),XH(I),UH(I),AH(I) 61 CONTINUE CLOSE(5) OPEN(5,FILE='P1.DAT') DO 28 I=1,NZ WRITE(5,*)X(I),PP(l,I),PP(2,I),PP(3,I) 28 CONTINUE CLOSE(5) OPEN(5,FILE='P2.DAT') DO 40 I=1,NZ WRITE(5,*)X(I),PP(4,I),PP(5,I),PP(6,I) 40 CONTINUE CLOSE(5) OPEN(5,FILE='P3.DAT') DO 50 I=1,NZ 138 WRITE(5,*)X(I),PP(7,I),PP(8,I),PP(9,I) 5O CONTINUE CLOSE(5) OPEN(5,FILE='P4.DAT') DO 51 I=1,NZ WRITE(5,*)X(I),PP(10,I),PP(11,I),PP(12,I) 51 CONTINUE CLOSE(5) OPEN(5,FILE='P5.DAT') DO 52 I=1,NZ WRITE(5,*)X(I),PP(13,I),PP(14,I),PP(15,I) 52 CONTINUE CLOSE(5) OPEN(5,FILE='P6.DAT') DO 53 I=1,NZ WRITE(5,*)X(I),PP(16,I),PP(l7,I),PP(18,I) 53 CONTINUE CLOSE(5) OPEN(5,FILE='P7.DAT') DO 54 I=1,NZ WRITE(5,*)X(I),PP(l9,I),PP(20,I),PP(21,I) 54 CONTINUE CLOSE(5) 23 STOP END SUBROUTINE AP(AX,AY,GA,PX,PY) TG=(AY/AX)**2 IF(TG.LT.1.) THEN PY=PX*TG**(GA/(GA-l.)) ELSE A=GA-l. B=(GA+1.)*(TG-l.) C=l.-GA PY=PX*0.5*(B+SQRT(B**2-4.*A*C))/A ENDIF RETURN END SUBROUTINE FF COMMON/A/W(2010),P(2010),A(2010),U(3,2010),FP(3,2010), * FM(3,2010),GA1,GA2,N1,N2 DO 10 I=1,N2 139 IF(I.LE.N2) THEN GA=GA1 ELSE GA=GA2 ENDIF D=U(1,I) W(I)=U(2,I)/D V=W(I) P(I)=(U(3,I)-.5*D*V**2)*(GA—l.) A(I)=SQRT((U(3,I)/D—.5*V**2)*GA*(GA-l.)) Fl=U(2,I) ’ F2=D*V**2+P(I) F3=(U(3,I)+P(I))*V Xl=ABS(V+A(I)) X2=ABS(V-A(I)) El=D*ABS(V)-D/GA*(ABS(V)—.5*Xl—.5*X2) E2=D*V*ABS(V)-D/GA*(V*ABS(V)—.5*(V+A(I))*Xl- .5*(V-A(I))*X2) .5*X2)+. lO * E3=D*(GA-l.)/GA*.5*V**2*(ABS(V)—.5*Xl— 5*D*(A(I)**2/ GA/(GA-l.)+.S*V**2)*(Xl+X2)+.5*D*V*A(I)/GA*(X1-X2) FP(1,I)=.5*(F1+E1) FM(l,I)=.5*(Fl-El) FP(2,I)=.5*(F2+E2) FM(2,I)=.5*(F2—E2) FP(3,I)=.5*(F3+E3) FM(3,I)=.5*(F3-E3) CONTINUE RETURN END FUNCTION XL(X1,X2,X3) Yl=X2-X1 Y2=X3-X2 IF(Y1*Y2.LT.0.0) THEN FM=0.0 ELSE Zl=ABS(Yl) ZZ=ABS