STABILITY of CONVERGING SHOCK WAVE by KENNETH SAU-KIN FONG B . S c , University of B r i t i s h Columbia, 1975 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (DEPARTMENT OF PHYSICS) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA March, 1978 © Kenneth Sau-kin Fong, 1978 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of Brit ish Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of P h y s i c s The University of Brit ish Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 A p r i l 4. 1978 i i ABSTRACT The s t a b i l i t y of converging waves was analysed with theory and by experiment. A computer code was w r i t t e n i n FORTRAN language f o r the c a l c u l a t i o n of wavefront p o s i t i o n s of plane shock waves and c y l i n d r i c a l converging shock waves i n i d e a l gas. The c r i t e r i a for s t a b i l i t y were d e f i n e d and the l i m i t of s t a b i l i t y of imploding shock waves was obtained using the code. In order to t e s t the model, an experiment was set up to generate shock waves i n an e l e c t r o t h e r m a l shock tube. These shock waves were made to propagate i n t o a wedged channel, and the v a r i a t i o n s of the shapes of the waves f r o n t were measured from high speed framing camera p i c t u r e s . Good agreement was found between the e x p e r i m e n t a l l y observed v e l o c i t y and the model c a l c u l a t i o n , i n d i c a t i n g that the computer procedure i s a v a l i d method to analyse the s t a b i l i t y of converging waves. TABLE OF CONTENTS Chapter 1 Introduction 1 1.1 Motivation For S t a b i l i t y Investigation 1 1.2 Outline Of Thesis 3 Chapter 2 Discussion Of Perturbation And S t a b i l i t y 5 2.1 Volume Perturbation 7 2.2 Radial Perturbation 9 2.3 Perturbation Angle 10 Chapter 3 Equation Of Motion Of Wavefront 11 3.1 The C.C.W. model Of An Imploding Wave 12 3.1.1 Equations Of Motion Behind The Front 12 3.1.2 Equations Of Motion Across The Front 14 3.1.3 Area Rule 17 Chapter 4 Numerical Simulation 26 4.1 Introduction And General Considerations 26 4.2 Grid Size Control 30 4.3 Changes Of Mach Number 33 4.4 I n i t i a l Conditions 35 4.5 Boundary Points 36 4.6 Streamlines 36 4.7 Time Interval 37 4.8 Limitations Of The Code 38 Chapter 5 Result Of Numerical Computation 41 5.1 Plane Shock Propagating Into Wedge-like Channel ...41 5.2 Symmetrical Circular Converging Shock 45 5.3 Circular Converging Shock Wave With Perturbation ..47 Chapter 6 Experiment 57 i v 6.1 Shock Tube And Vacuum System 58 6.2 Power Source And Trigger C i r c u i t 61 6.3 Method Of Measurement 61 6.4 Result 64 Chapter 7 Conclusion 69 Bibliography 71 Appendix A Repeated Linear Interpolation 73 Appendix B Runge-Kutta Integration Method 76 Appendix C Remark On Accuracy 79 Appendix D Code For Wavefront Calculation 80 TABLE OF FIGURES Figure(2.1) 8 Figure(3.1) 14 Figure(3.2) 14 Figure(3.3) 18 Figure(3.4) 20 Figure(3.5) 22 Figure(3.6) 24 Figure(4.1) 28 Figure(4.2) 30 Figure(5.1) 41 Figure(5.2) 43 .Figure(5.3) 45 Figure(5.4) . 47 Figure(5.5) 49 Figure(5.6) 50 Figure(5.7) 51 Figure(5.8) 52 Figure(5.9) 53 Figure(5.10) 53 Figure(5.11) 54 Figure(6.1a) 58 Figure(6.1b) 58 Figure(6.2) 6 1 Figure(6.3) 6 4 Figure(6.4) 6 6 Figure(A.l) 7 4 Figure(B.1) v i 77 1 CHAPTER 1 INTRODUCTION 1.1 Motivation For S t a b i l i t y Investigation Hydrodynamic s t a b i l i t y of imploding shock and heat waves plays an important role in the discussion of the f e a s i b i l i t y of laser fusion. C l a s s i c a l Rayleigh-Taylor s t a b i l i t y theory cannot be applied to shock fronts or ablation heat fronts such as those produced by laser due to mass transport across the fronts in question. The s t a b i l i t y theory of Bulter/4/ on i n f i n i t e l y strong converging shock waves using Guderley/11/ s i m i l i a r i t y method considers only i n f i n i t e s i m a l l y small perturbations, besides, important d i s s i p a t i v e mechanisms were not considered. This theory i s therefore not s u f f i c i e n t to determine the s t a b i l i t y conditions as required by laser fusion, since the imploding wave produced by lasers are neither i n f i n i t e l y strong nor can the perturbations be expected to be i n f i n i t e s i m a l l y small. It i s well known that plane shock fronts are stable to perturbations. This has been derived t h e o r e t i c a l l y by L i g h t h i l l / 1 6 / , Freeman/10/ and Chester/6/ and experimentally by Lapworth/14/. Experimental results on concentric detonations from Perry and Krantrowitz(1959)/18/, Ahlborn and Huni(1968)/l/ and Knystantaus and Lee(1970)/13/ show that c y l i n d r i c a l imploding detonation fronts are stable to f i n i t e 2 s i z e d p e r t u r b a t i o n s . In p a r t i c u l a r Ahlborn and Huni observed with an image-converter framing camera that f i n i t e s i z e d p e r t u r b a t i o n s i n the imploding d e t o n a t i o n f r o n t d i m i n i s h as time progressed. What was i n the beginning a f r o n t with many angular regions became a p e r f e c t l y smooth c y l i n d r i c a l f r o n t as i t c o l l a p s e d towards i t s c e n t e r . Converging d e t o n a t i o n waves are s t r o n g l y o v e r d r i v e n and behave l i k e s trong shock waves when c l o s e to the c e n t e r s of c o l l a p s e , t h i s leads to two q u e s t i o n s : 1) . Whether converging shock waves are s t a b l e ; and 2) . What i s the maximum s i z e of a p e r t u r b a t i o n f o r which s t a b i l i t y i s not l o s t , or i n other words what s i z e of a p e r t u r b a t i o n w i l l modify the e n t i r e wavefront to such an extent that i t s o r i g i n a l shape can no longer be r e s t o r e d . This t h e s i s t r i e s to answer the above two q u e s t i o n s both n u m e r i c a l l y and e x p e r i m e n t a l l y . The c o n d i t i o n s f o r s t a b i l i t y as a p p l i c a b l e to l a s e r f u s i o n i s d i s c u s s e d , from which the meaning of s t a b i l i t y i s d e f i n e d . A fast, numerical procedure to c a l c u l a t e the p o s i t i o n s of v a r i o u s wavefronts as the f r o n t s c o l l a p s e i s d e s c r i b e d . The procedure i s coded i n FORTRAN language and i s able to c a l c u l a t e the wavefront p o s i t i o n s and Mach numbers of plane shock waves propagating i n t o a wedged channel and c y l i n d r i c a l converging shock waves with d i f f e r e n t p e r t u r b a t i o n s . The method d i f f e r s fundamemtally from standard s t a b i l i t y a n a l y s i s , i n which the e v o l u t i o n of small s c a l e p e r t u r b a t i o n s are 3 predicted. Here the evolution of large perturbations are calculated. The conditions for s t a b i l i t y i s then obtained from a series of such calculations of various c y l i n d r i c a l wavefronts with d i f f e r e n t perturbations. An experiment of a plane shock wave propagating into a wedged channel i s then performed and compared to the numerical result obtained by the numerical method to test the v a l i d i t y of the code. 1.2 Outline Of Thesis This thesis consists of two main parts: theorectical investigation of the s t a b i l i t y problem with numerical results and an experiment to test parts of the numerical r e s u l t . Chapter 2 contains the c r i t e r i a of hydrodynamic s t a b i l i t y as applicable to laser fusion. In Chapter 3 the equation of motion of d i f f e r e n t types of wavefronts i s derived using the method developed by Whitham. In Chapter 4 these equations for the c a l c u l a t i o n of wavefronts are then solved with numerical method. The calculations y i e l d the development of perturbed wavefronts and also l i m i t s for the s t a b i l i t y i s given in th i s chapter. The experiment used to test the accuracy of the computer s i m u l a t i o n i s d e s c r i b e d i n Chapter 5. Chapter 6 then draws the c o n c l u s i o n with suggestion f o r f u r t h e r works. 5 CHAPTER 2 DISCUSSION OF PERTURBATION AND STABILITY Most perturbations on converging shock waves w i l l not be dampened out e n t i r e l y as shown by Bulter/4/ and Whitham/23/. They concluded that i n f i n i t e s i m a l l y small disturbances of higher harmonic modes in a strong converging shock front w i l l grow instead of decay. This implies that the o r i g i n a l symmetry of the shock front w i l l not be restored but instead the perturbation w i l l spread out and modify the entire front. However due to a s t a b i l i z i n g mechanism which tends to smooth the converging wavefront, t h i s modification tends to cause the front to achieve some sort of d i f f e r e n t symmetry at a later time. This s t a b i l i z i n g mechanism arises from the dependence of the front v e l o c i t y on the changes of area as the front progresses in time. This dependence can be calculated in a method developed by Whitham/21/ for shock waves. The rate of change of the front v e l o c i t y i s a function of the strength of the wave and the l o c a l curvature of the wave. The strength of a wave with a smaller radius of curvature increases more rapidly due to pressure build up behind the front from the compression of the gas flowing into a channel of decreasing width. Consider a converging perturbation which causes a wavefront to lag behind, then at the center of this perturbation, where the wave i s concave, the front v e l o c i t y increases more rapidly than the rest of the wavefront due to i t s higher curvature. This difference in 6 front v e l o c i t y increase w i l l speed up the parts of the fronts which lag behind. S i m i l a r l y the parts which lead w i l l be slowed down. The difference in curvatures along the front is reduced and at the same time the perturbation i s spread out along the front. S t a b i l i t y in the usual sense means that a perturbation in a wavefront w i l l be damped out and the wavefront w i l l restore i t s o r i g i n a l symmetry. This implies that most converging shocks are unstable according to the above d e f i n i t i o n . For many purposes, in particular where one expects the existence of f i n i t e sized perturbations as for instance in laser fusion, the d e f i n i t i o n of s t a b i l i t y does not have to be so r e s t r i c t i v e , because a symmetrical sphere i s not absolutely necessary for high density compression. However, i t i s desirable that the converging wavefront stays more or less spherical at r a d i i greater than the radius R of the f i n a l compressed mass. A wavefront that breaks up into several parts i s obviously unstable. Furthermore, s t a b i l i t y at r a d i i smaller than R i s not required since the necessary heating and compression have already been achieved when the radius of the shock front is equal to R. One can thus c l a s s i f y s t a b i l i t y to be "absolute" for perturbations which continue to decrease as functions of R and " p a r t i a l " for perturbations which decrease as functions of radius r for r>R, but may increase for r<R. We can also define the geometries of the perturbation by 7 1/. I t s volume AV.(area AA i n 2-dimension) 2/. I t s maximum p e r t u r b a t i o n amplitude Ar. 3/. I t s d e v i a t i o n from s i m i l a r i t y to the o r i g i n a l shape. 2.1 Volume P e r t u r b a t i o n We d e f i n e the extent of the volume p e r t u r b a t i o n as the dimensionless parameter PERT , which i s the r a t i o of the v d i f f e r e n c e i n volume AV enclosed by the perturbed f r o n t and an unperturbed r e f e r e n c e f r o n t to the volume enclosed by the unperturbated r e f e r e n c e f r o n t , V r e £ , i . e . PERT v = A V /V r e f (2.1) Absolute volume s t a b i l i t y would r e q u i r e d PERT to be a d e c r e a s i n g f u n c t i o n of time, where PERT v approaches zero before v r e f reaches zero. P a r t i a l volume s t a b i l i t y r e q u i r e s o n l y that PERT^ decreases f o r r a d i u s l a r g e r than some R. N o t i c e that f o r l a s e r f u s i o n purposes, a converging shock does not have to be a b s o l u t e l y s t a b l e s i n c e the f u e l p e l l e t i s not to be compressed to an i n f i n i t e s i m a l s i z e . However, p a r t i a l s t a b i l i t y at a r a d i u s comparable to the f i n a l compressed r a d i u s i s necessary f o r the shock compression scheme i n l a s e r f u s i o n . For two dimensional geometry, the volume term i n Eq(2.1) should be r e p l a c e d with an area A to o b t a i n the area p e r t u r b a t i o n equation. T h i s d e f i n i t i o n 8 w i l l be used l a t e r i n Chapter 5 f o r c y l i n d r i c a l wavefront c a l c u l a t i o n s . - A X I S OF V i — SYMMETRY \ PERTURBED AREA AA r :LOCAL RADIUS OF CURVATURE F i g . (2.1) P e r t u r b a t i o n parameters N o t i c e that i n the above d e f i n i t i o n of the p e r t u r b a t i o n only the volume of the p e r t u r b a t i o n i s c o n s i d e r e d . Neither the geometry of the p e r t u r b a t i o n nor the d i s t r i b u t i o n of f r o n t 9 v e l o c i t y along the perturbed f r o n t has been taken i n t o c o n s i d e r a t i o n . With t h i s d e f i n i t i o n alone, a p e r t u r b a t i o n which has s p l i t the wavefront i n t o two non-coherent p a r t s can a l s o be c a l l e d s t a b l e i f PERT y decreases. The requirement of volume s t a b i l i t y alone i s t h e r e f o r e i n s u f f i c i e n t to guarantee a smooth c o l l a p s e i n l a s e r f u s i o n s i n c e the energy a b s o r p t i o n would be g r e a t l y impeded i f the wavefront was permitted to break up i n t o d i f f e r e n t p a r t s . 2.2 R a d i a l P e r t u r b a t i o n To supplement the above d e f i n i t i o n , we d e f i n e another di m e n s i o n l e s s parameter " r a d i a l p e r t u r b a t i o n " ( P E R T ^ ) , namely the r a t i o of the maximum p e r t u r b a t i o n amplitude Ar to the re f e r e n c e unperturbed f r o n t with r a d i u s R PERT = Ar/R (2.2) r " R a d i a l s t a b i l i t y " then r e q u i r e s PERT to be a d e c r e a s i n g r f u n c t i o n of R, and i f PERT goes to zero before R does, then r we c l a s s i f y i t as "absolute r a d i a l s t a b i l i t y " , and " p a r t i a l r a d i a l s t a b i l i t y " i f PERT decreases f o r r a d i u s r g r e a t e r than r some R. For l a s e r f u s i o n purpose, one would r e q u i r e both p a r t i a l area s t a b i l i t y and p a r t i a l r a d i a l s t a b i l i t y , at l e a s t up to the r a d i u s of i g n i t i o n . 10 2.3 P e r t u r b a t i o n Angle So f a r we have d e s c r i b e d a p e r t u r b a t i o n by two dimensionless parameters PERT^ and PERT^, but these s t i l l do not g i v e a complete d e s c r i p t i o n of the shape of the p e r t u r b a t i o n . For a comprehensive d e s c r i p t i o n of the shape of the p e r t u r b a t i o n , many more parameters may be r e q u i r e d . I f we r e s t r i c t o u r s e l v e s to two-dimensional geometry, and i f only a gen e r a l idea of the shape i s r e q u i r e d , only one e x t r a parameter i s needed. Such a parameter i s Y; the p e r t u r b a t i o n angle, which i s one h a l f of the angle spanned by the p e r t u r b a t i o n . Furthermore, t h i s can be used to d e f i n e an azimuthal number 1, which i n d i c a t e s the number of times a p e r t u r b a t i o n can be f i t t e d i n t o a complete c i r c l e . Thus thus 1 = 180°/y The parameters PERT^, PERT^, and y, which vary as the f r o n t c o l l a p s e s , w i l l g i v e a gen e r a l idea of the shape of the p e r t u r b a t i o n as the f r o n t c o l l a p e s to the c e n t r e . For completeness, the Mach number of the f r o n t , from which the pr e s s u r e , d e n s i t y and temperature of the flow f i e l d behind the f r o n t can be obta i n e d , should a l s o be s p e c i f i e d along the f r o n t . But the d i s t r i b u t i o n of the mach number along the shock f r o n t i s not apparent from a snap shot of the f r o n t , and t h i s i n f o r m a t i o n i s not considered i n the f o l l o w i n g d e r i v a t i o n s . 11 CHAPTER 3 EQUATION OF MOTION OF WAVEFRONT This chapter describes the conservation equation across a wavefront and also the d i f f e r e n t i a l conservation equations for continuous gas flow. The d i f f e r e n t i a l conservation equations of mass, momentum and energy behind the front, which are a set of hyperbolic p a r t i a l d i f f e r e n t i a l equations, are combined and rewritten in c h a r a c t e r i s t i c form. The jump conditions across a front are obtained from the integral conservation equations. These are combined to obtain a equation of motion describing the front. This method was f i r s t developed by Chisnell/7/, Chester/6/ and Whitham/21/ independently. By this method one obtains a non-linear ordinary d i f f e r e n t i a l equation of the Mach number M of the front and the area A of the front. In deriving this equation, the complicated flow behind the front i s neglected through the assumption that the jump conditions can be applied to the C+ c h a r a c t e r i s t i c of the flow, and that the area A i s a function of the l o c a l Mach number M of the front. These assumptions would not be v a l i d in explosion problems in which the propagation towards the region of disturbance must be of prime consideration. However, i t works to a very good extent in spherical and c y l i n d r i c a l converging shock waves as indicated by Bulter's comparison of the result obtained by using the CCW model and the method of s i m i l a r i t y . It should be noted that the equation describes the area A 12 as a function of the Mach number M , yet i t i s the change in area that causes a corresponding change in pressure behind the front, which then causes the speed of the wavefront and the Mach number to change. 3.1 The C.C.W. model Of An Imploding Wave 3.1.1 Equations Of Motion Behind The Front Chester /I/, C h i s n e l l /6/ and Whitham /21/ developed a model for calculating the front v e l o c i t y as a function of area. Lee & Lee /6/ applied this method for the case of detonation. A similar procedure can be applied to heat waves in general. The d i f f e r e n t i a l conservation equation for a quasi-one dimensional gas flow with no heat source or sink i s : Momentum: |H. + U|E + 11£. = 0 (3.2) 3t or p 9r Energy: -|| + u||- a 2 ( § | • uf|) = 0 (3.3) where P,u,p,a i s the density, gas ve l o c i t y in the lab frame, pressure and l o c a l sound v e l o c i t y of the gas. 13 Manipulating these equations according to the expression: a 2 x (3.1) + (3 .3 ) t a p x (3 .2 ) (3.3) one gets the c h a r a c t e r i s t i c forms of the conservation equation: d £ ± a du + a 2 p u 3A. = Q dt * ""dt A 3r (3.4) where the +ve sign applies to the pos i t i v e c h a r a c t e r i s t i c : du/dt=u+a and the -ve sign applies to the negative c h a r a c t e r i s t i c : du/dt=u-a changing variables along the c h a r a c t e r i s t i c s 1 3A 1 3A 3t 1 3 A 1 A 3r A 3t 3r A 3t u±a (3.5) and substituting (3.5) into (3.4), one has a set of ordinary d i f f e r e n t i a l equations along the c h a r a c t e r i s t i c s : . , , a 2pu dA _ dp ± apdu + J^+a A = 0 (3.6) These equations are to be applied to regions which do not have any discontinuity, i . e . in regions behind the wavefront. A Space-time diagram of the wave and the C+, C-c h a r a c t e r i s t i c s is given in F i g . (3.1). The schematic 1 4 diagram of the p h y s i c a l s i t u a t i o n of a converging shock wave i s shown i n F i g . ( 3 . 2 ) . Whitham / 2 1 / has shown that f o r converging f r o n t the C+ c h a r a c t e r i s t i c can be a p p l i e d to the flow immediately behind the f r o n t with very good approximation. In the next s e c t i o n the p r e s s u r e , v e l o c i t y and d e n s i t y w i l l be expressed as a f u n c t i o n of the f r o n t v e l o c i t y . 3 . 1 . 2 Equations Of Motion Across The Front The i n t e g r a l c o n s e r v a t i o n equations across the f r o n t are: Mass: P l v x = P 2 V 2 ( 3 . 7 ) Momentum: p l + p l v l = P 2 + P 2 V 2 ( 3 . 8 ) Energy: hj + %v 2 + = ^ 2 + ( 3 . 9 ) Where p,v,p, and h are the d e n s i t y , the gas v e l o c i t y measured i n the f r o n t frame, the p r e s s u r e , and the enthalpy per u n i t mass. W i s the input energy at the f r o n t , measured i n Watt/m . For chemical r e a c t i o n s , one o f t e n g i v e s the s p e c i f i c energy input Q=W/pv. The s u b s c r i p t s 1 and 2 r e f e r to the c o n d i t i o n s ahead of and behind the f r o n t resp. F i g . (3.1) Space-time diagram of converging shock 16 P R E S S U R E D E N S I T Y V E L O C I T Y F i g . (3.2) Schematic diagram of converging shock wave The general solution of equations (3.7) , (3.8), (3.9) can be rearranged as /!/: _ 2 = J U 1 + 1 n 3 (3.10a) v l p2 «?i £i = 1 - nB (3.10b) Pl n = g2 " g l M l (3.10c) g 2 + 1 J 2 g ^ ( g 2 - l ) W * = 1 ± t 1 + G - (8 1-l)(g 2-g 1Mj)^P 1v 1h 1 17 e = 2(g 2 +DCg 1 -g 2 ) M fg 1 (g 1 -D(g 2 -g 1Mp^ (3.10e) (3. 1 0 f ) 2 ^1 al = giP7 (3.10g) where e i s a l w a y s s m a l l . In g e n e r a l 0<3<2. The s q u a r e r o o t i n (3.10d) has t h e p o s i t i v e s i g n t o g e t h e r w i t h M ^ >1 and W=0 f o r s h o c k s and t h e n e g a t i v e s i g n f o r h e a t waves. S u p e r s o n i c h e a t waves a r e f o u n d f o r M ^ >1 and s u b s o n i c h e a t waves f o r M ]<1. 3.1.3 A r e a R u l e The a r e a r u l e i s o b t a i n e d by r e p l a c i n g p, u and p i n t h e C+ c h a r a c t e r i s t i c e q u a t i o n by t h e jump c o n d i t i o n s a c r o s s t h e f r o n t . The s u b s t i t u t e d e q u a t i o n i s t h e n e x p r e s s e d o n l y i n te r m s o f t h e Mach number M o f t h e f r o n t and t h e a r e a A. T h i s c a n be e a s i l y done i n t h e frame o f r e f e r e n c e i n w h i c h t h e gas ahead o f t h e f r o n t i s s t a t i o n a r y i n i t i a l l y . In t h i s c a s e Mi=M i s t h e mach speed o f t h e f r o n t , 18 and the gas v e l o c i t y behind the f r o n t i s U = V r V2 = - i ^ (3.11a) } i \I^MJ+n8/ P a P i 7 ? + , R (3 .11b) a,P2 - L - (1-nS) (3.11c) 81 f o r the case of shock , 8=2 2na u «= r r ~ (3.12a) 81 M1 81 M1 p - p l l ^ 2 n " ( 3 ' 1 2 b ) p . V i ( 1 _ 2 n ) (3.12c) 81 one gets Whitham's area r u l e by s u b s t i t u t i n g (3.12) i n t o (3.6) X (M)MdM . + ^=0 (3.13) (M^-g2/gl) + A where Xs(M) - i ^ ) ( 1 + i ^ T ^ " (g -l)M2+2 and a 2 = — £ _ _ 2g 2M i !-(g 2-l) X g(M) i s a s l o w l y v a r y i n g f u n c t i o n of M. The dependence of X g(M) on M i s shown i n Fig.(3.3) f o r an i d e a l gas. 20 For i n f i n t e l y strong shocks in ideal gas, this reduces to dM , dA _ n — + — = 0 s M A (3.14) where n » 1+ — s Y in t h i s case •f c ^ ^ - 4.4361 for y = 5/3 » 5.0743 for y = 7/5 M a A s (3.15) therefore the rule gives M a r ^ _ 1 / ' n s ^ f o r c y l i n d r i c a l shock (3.16a) ( 2 /n ) M a r s f o r s p h e r i c a l shock (3.16b) these relations may be compared to the Guderley s i m i l a r i t y solution where Mar.exp(l/n). For the case of a Chapman-Jouguet detonation, where 3=1 the area rule y i e l d s X C J ( M ) M d M ^ d A . (M^-g 2/ g l) A (3.17) where X CJ and i t s graph is shown in Fig.(3.4). 21 22 For i n f i n t i y strong d e t o n a t i o n wave one f i n d s dM . dA A n C j 1 l + -A = 0 (3.18) where n_T = 3(1+ -) = 4.80 for Y = 5/3 - 5.14 for Y = 7/5 which gives M a r^ - 1^ nCJ^ for c y l i n d r i c a l (3.19a) CJ detonation M a r*~2/nC.P for spherical (3.19b) CJ detonation Notice the close resemblance between the equations for shock front and for Chapman-Jouguet detonation front. It should be noted that for weak fronts, the strengthening i s only moderate. But for fronts with M>5, the increase in Mach number i s quite rapid as the area converges. It is a consequence of this fact that strong converging shock wave can be stable in some aspect. Fi g . (3.5) shows the strengthening of c y l i n d r i c a l converging shock waves as a function of area for various i n i t i a l shock strengths, in an ideal gas of Y =7/5 and 5/3. The result i s obtained by numerically integrating Eq.(3.13), notice that the i n f i n i t e l y strong shock approximation (M=°°) can be applied to shocks of i n i t i a l strength M=5 and higher without incurring an error of more than 5% for the worst case F i g . (3.5) Area st r e n g t h e n i n g of shock wave 24 in the f i n a l strength of the shock, strengthening of CJ detonation waves. Fig.(3.6) shows the 25 26 CHAPTER 4 NUMERICAL SIMULATION 4.1 I n t r o d u c t i o n And General C o n s i d e r a t i o n s Equations(3.13) and (3.17) are n o n l i n e a r with r e s p e c t to the Mach number M of the f r o n t . An a n a l y t i c s o l u t i o n to an asymmetric system such as a s i n g l e p e r t u r b a t i o n on a c y l i n d r i c a l converging wave i s very d i f f i c u l t . The propagation of the wavefront can however be evaluated n u m e r i c a l l y . T h i s numerical wavefront method c o n s i d e r s a set of mesh p o i n t s which are d i s t r i b u t e d u n i f o r m l y i n space along the f r o n t . The motion of the mesh p o i n t s are f u n c t i o n s of the l o c a l Mach numbers of the f r o n t . The equation of motion of the f r o n t i s i n t e g r a t e d n u m e r i c a l l y at each mesh p o i n t and the f r o n t p o s i t i o n s at s u c c e s s i v e time i n t e r v a l s can be p l o t t e d . Rays are introduced as orthogonal t r a j e c t o r i e s of these s u c c e s s i v e f r o n t p o s i t i o n s . Rays c l o s e to the f r o n t would represent the a c t u a l p a r t i c l e paths as the jump c o n d i t i o n s across the f r o n t r e q u i r e that the flow leave the f r o n t i n normal d i r e c t i o n . The rays can then be t r e a t e d as the s o l i d w a l l s of flow tubes which no p a r t i c l e s can c r o s s . In t h i s way the problem of d e s c r i b i n g the convergence of waves i s reduced to d e c r i b i n g the propagation of a wavefront element c o n f i n e d w i t h i n flow tubes with s o l i d w a l l s . The l o c a l Mach number of each ray tube 27 is however a function of the ray-tube area and the r e l a t i o n i s given by Eg.(3.13) and (3.17). Using t h i s concept, the new l o c a l Mach number at later times can be calculated and the successive positions of the front can be plotted. This method i s similar to the method of c h a r a c t e r i s t i c s except that the c h a r a c t e r i s t i c equations are not used and only those points located along the wavefront are considered. The domain of dependence for the integration i s determined by the ve l o c i t y of propagation of the disturbance along the wavefront. Furthermore, in contrast to the usual moving mesh point method, the mesh point in the code are partly prescribed by f i x i n g their angular separations to a prescribed center. The r a d i i are allowed to vary and new r a d i i are to be obtained by in t e r p o l a t i o n . Physical quantities at the mesh points can then be evaluated by using Eq.(3.11a), (3.11b) and (3.11c). The advantage of this semi-moving mesh point method i s two-fold: (i) The p o s s i b i l i t y of overcrowding of the mesh points due to too large a d i s c r e t i z a t i o n in time step i s removed. ( i i ) The grid size i s determined mainly by the radius, thus simplifying error control. However, there are also several disadvantages. F i r s t l y i t needs more interpolations than i f the mesh points 28 were to be determined by the t r a j e c t o r i e s of the previously calculated points, and secondly the streamlines are not readily available. But the advantages outweight a l l these disadvantages, es p e c i a l l y in the case where a shock-shock 1 i s present, trading off a certain degree of accuracy of the fixed flow tube approach with the s t a b i l i t y of the varying flow tube approach. B r i e f l y , one calc u l a t i o n step consists of the following cal c u l a t i o n s : 1. Evaluation of new positions of mesh points by the wavefront method. 2. Calculation of Mach number at new mesh points by Whitham's area rule. 3. Choosing new mesh points. 4. Calculation of normals and Mach numbers at the new mesh points. A flow chart of the operation i s shown in Fig.(4.1). The following six sections give a d e t a i l description of the process of d i s c r e t i z a t i o n , time integration and area integration of the wavefront. Special treatments of boundary points and streamlines are given in 4.5 and 4.6 respectively, while the choice of time intervals used i s 1 A\ s;hock-shock i s an abrupt change of physical variables t r a v e l l i n g sideways across a shock front. 29 READ I.Ar, M initialize X.Y.M.T C STOP ) • CIRCLE _Y£5_ / PLOT 7 f START j READ SHAPE ^HAT .SHAPE •WEDGE-(total number of time steps)] j= 1 XN = X*M adt YN=Y* MyOdt J: NEW NORMALS T : X.Y initialize X.Y.M.T Y E S / \ N O MOD(J,20)^ 0 YES J = J • 1 WHITHAM AREA RULE MN :X,XN,Y,YN,M REPEATED LINEAR INTERPOLATION X : XN.YN Y : XN.YN M: XN, Y N,MN N O S T O P F i g . (4.1) Flow chart of numerical procedure 30 discussed in 4.7. Special treatments of 4.8 then discusses the l i m i t a t i o n s of the code to physical problems. 4.2 Grid Size Control In order to make the code applicable to a wide variety of problems, dimensionless parameters are used wherever possible. For a set of mesh points on the converging front given by ( r£ , 6 k ) for k = 1 K : (4.1) and where 6,. = kA6 where k and n are non-negative integers and where i e i s a constant.(Fig.(4.2)) k i s the label of a mesh point and n i s the label of a time step. The r" are obtained by dividing the actual physical dimensions of the wavefront by a normalizing factor K such that r j i s one unit of length. l e t t n + 1 - t n = A t n and l e t <j>£ denote < M r £ , V V ' I n cartesian coordinate the position of a mesh point i s defined by R", where R£ = ( x£ , Y £ ) (4.2a) x£ - r£cosGk (4.2b) (4.2c) F i g . (4.2) Coordinate system for G r i d S i z e C o n t r o l 32 Assuming the wavefront between mesh points R £ ^, and and t o D e a segment of a c i r c l e with the three points on i t s circumference. The d i r e c t i o n of propagation of the wavefront at R£ i s then the normal N£ of the constructed ci r c l e . ( F i g . ( 4 . 2 ) ) At the next time step t=t ,, the front which was r n+1 described i n i t i a l l y at R£ would be at R ^ + , where *r 1- (*r 1 'yr 1 > X ^ 1 = x£ + M|WncosN£ (4.3a) *k + 1 = Yk + M k a A t n 8 l n N k (4.3b) where M " i s the Mach number at the mesh point R " a i s the dimensionless speed of sound, obtained by dividing the actual speed of sound by the proportionality constant K Transforming back into polar coordinates r f 1 =V(xf V + ttfV ( 4 - 4 a ) i f 1 = T a r f V f Vxf X> ; ( 4 * 4 b ) i 33 the new mesh points at the prescribed angles are then evaluated by repeated linear interpolation The new d i r e c t i o n of propagation at each new mesh point i s the normal of the wavefront at the mesh point. The normals are approximated by considering the wavefront as composed of segments of arcs of c i r c l e s connecting the mesh points. Three neighbouring mesh points can describe an arc of a c i r c l e which passes through the three points. The normal of the wavefront at the middle mesh point i s then approximated by the normal of the c i r c l e thus formed. At the boundary points, one of the neighbouring mesh points i s absent. This has to be treated d i f f e r e n t l y as outlined in the following section. 4.3 Changes Of Mach Number The trajectory of R£ to i s the trajectory of p a r t i c l e s which are close to the shock front. This defines the boundaries of a ray-tube through which the p a r t i c l e s do not cross. The area of the ray-tube between mesh points k and k+1 for c y l i n d r i c a l or wedge geometry i s given by the norm (4.6) 34 The changes in area of the ray-tube from t to t n + 1 i s then A A ^ ^ + 1 - A ^ (4.7) Due to the changes in area, the gas in the tube behind the front w i l l experience a pressure change which in turn causes a change in the Mach number of the front. The Mach number at the center of the tube i s taken to be the average of the Mach numbers at the mesh points M£ - (M£ + 1 + M£)/2 (4.8) The Mach number at the middle of the ray tube at the next time step can be obtained by integrating M numerically using Eq.(3.13). For this integration a second order Runge-Kutta integration method i s used. The Mach number at the new mesh point i s then evaluated by linear interpolation. 1 A detailed description of the method i s given in Appendix B. 35 4.4 I n i t i a l C o n d i t i o n s — -The necessary i n i t i a l c o n d i t i o n s r e q u i r e d to i n i t i a t e the s t a r t of the c a l c u l a t i o n are 1) . The i n i t i a l p o s i t i o n of the mesh p o i n t s , which i s the i n i t i a l shape of the wave f r o n t , s p e c i f i e d by g i v i n g R^. 2) . The i n i t i a l s t r e n g t h of the f r o n t , s p e c i f i e d by g i v i n g M^ .' In the absence of d e t a i l e d i n f o r m a t i o n of the i n i t i a l s t r e n g t h of the f r o n t , the shock f r o n t i s taken to be uniform i n s t r e n g t h , and v a r i o u s Mach numbers are t r i e d . Two g e n e r a l shapes are used. The f i r s t i s a plane shock f r o n t going i n t o a wedged-shape channel. T h i s c a l c u l a t i o n i s l a t e r compared to experimental r e s u l t s . The second shape i s an imploding s e m i c i r c u l a r f r o n t . A p e r t u r b a t i o n of v a r i o u s s i z e s and d i f f e r e n t i n i t i a l shock s t r e n g t h s can be superimposed onto the c y l i n d r i c a l f r o n t . Two parameters are used to s p e c i f y the i n i t i a l s i z e of the p e r t u r b a t i o n . Ar i s the i n i t i a l maximum displacement of the perturbed wavefront from the unperturbed wavefront. 1 i s the azimuthal number, namely the number of times such a p e r t u r b a t i o n can be f i t t e d onto a c y l i n d e r . Two general p e r t u r b a t i o n shapes were used. I f a smooth d e v i a t i o n from the unperturbed wavefront i s a n t i c i p a t e d , the p e r t u r b a t i o n could be i n the shape of a s i n e curve superimposed on the wavefront. I f one assumes an abrupt d e v i a t i o n , the wavefront at the perturbed r e g i o n i s r e p l a c e d by a c i r c u l a r 36 arc of an a p p r o p r i a t e r a d i u s . 4.5 Boundary P o i n t s At the boundary p o i n t s k=l or N, the absence of one of the neigbouring mesh p o i n t s makes i t impossible to use Eq.(4.6) f o r the c a l c u l a t i o n of new Mach numbers and the new d i r e c t i o n s . A l s o boundary c o n d i t i o n s may f o r c e the shock f r o n t to change i n i t s d i r e c t i o n of p r o p a g a t i o n . By t r e a t i n g the boundaries as e i t h e r s o l i d w a l l s or as the a x i s of symmetry, as the case may be, the d i r e c t i o n of propagation can be determined as p a r a l l e l to the boundary. The new Mach number can then be c a l c u l a t e d by c o n s i d e r i n g the symmetry of the problem. 4.6 Streamlines As the mesh p o i n t s are p a r t i a l l y f i x e d , s t r e a m l i n e s of the wavefront are not r e a d i l y a v a i l a b l e . However, by the a d d i t i o n of a " f l o a t i n g mesh p o i n t " , a p a r t i c u l a r s t r e a m l i n e can be t r a c e d i n the program. 37 4.7 Time I n t e r v a l The choice of the time i n t e r v a l f o r each time step i s very important f o r the s t a b i l i t y , the d i s c r e t i z a t i o n and t r u n c a t i o n e r r o r of the c a l c u l a t i o n . A shor t time i n t e r v a l w i l l o b v i o u s l y reduce the e r r o r i n c u r r e d i n d i s c r e t i z i n g the mesh p o i n t s and the time i n t e r v a l and a l s o the t r u n c a t i o n e r r o r of the d i f f e r e n c e equations. However, too short a time i n t e r v a l w i l l make the computational time unacceptably l o n g . Too long a time i n t e r v a l i s a l s o unacceptable: In the program a mesh p o i n t i s only coupled to two of i t s nearest neighbours, and i f the i n t e r v a l i s too long, the domain of dependence of a mesh p o i n t i n c l u d e s neighbours other than i t s nearest neighbours. Therefore the s o l u t i o n w i l l o s c i l l a t e . In the extreme case t h i s w i l l cause the mesh p o i n t s to cr o s s over each o t h e r . The maximum time i n t e r v a l one can use before such a c o n d i t i o n occurs t h e r e f o r e depends on the s e p a r a t i o n between the neighbouring p o i n t s and the Mach numbers which are v a r y i n g d u r i n g the computation. A constant time i n t e r v a l can be used i n the wedge geometry, where the s e p a r a t i o n s between mesh p o i n t s remain i n the same order of magnitude f o r a long time. But f o r a c y l i n d r i c a l geometry, where the wave i s converging towards the c e n t e r , the time i n t e r v a l has to be changed i n every time step i n order to compensate f o r the r e d u c t i o n i n s e p a r a t i o n s between the mesh p o i n t s . In c o n j u n c t i o n with the use of dimensionless l e n g t h 38 scale, K, the dimensionless time paratmeter D i s defined as D=aAt (4.9) where At i s the real time in t e r v a l and a i s the dimensionless speed of sound defined previously in 4.2. We define AR£ to be the separation between mesh points Rk a n d R k + 1 A R k = | R k " R k + l " ( 4- 1 0> It was found that in order for the c a l c u l a t i o n to be stable, the following inequality must be obeyed: —-> 20 for a l l n,k (4.11) 4.8 Limitations Of The Code The code uses the wavefront method together with the Whitham's area rule to calculate the propagation of the wavefront. Tremendous reduction in computer time is achieved by considering only the physical quantities at the wavefront rather than calculating the entire flow f i e l d . This time saving i s however not without some s a c r i f i c e s in the f l e x i b i l i t y of the code. Besides subjecting i t s e l f to the usual l i m i t a t i o n s of numerical methods, the code i s also subjected to the l i m i t a t i o n s imposed on the wavefront method and the Whitham's area rule. Numerical simulation cannot give an overa l l picture of a 39 physical problem and i t i s often v e r y ' d i f f i c u l t to trace how a particular quantity develops as a function of the rest of the system. Also, since a l l data must be i n i t i a l i z e d , the numerical results obtained from the simulation are only applicable to a p a r t i c u l a r i n i t i a l condition. S i m i l a r i t y c r i t e r i a can be applied to a wider class of problems, but often results must be calculated i n d i v i d u a l l y for d i f f e r e n t i n i t i a l conditions. Our code exemplifies this point, because the i n i t i a l strength d i s t r i b u t i o n of the wavefront and also the i n i t i a l shape of the perturbation are asymmetrical and cannot be generalized by s i m i l a r i t y rules. The Whitham's area rule i s obtained by applying the jump conditions across a front to the C+ c h a r a c t e r i s t i c s which terminates at the front. The flow f i e l d behind the front is not considered and i s assumed to have no e f f e c t on the front. This c r i t e r i o n is true with a converging shock wave where the changes at the front are primarily due to area convergence as the front propagates. To a certain extent, i t can also be applied to plane and steady shocks where the flow f i e l d behind the front has only a s l i g h t effect on the front. But the same cannot be said to a divergent shock wave where the shock slows down appreciably as i t expands since the pressure driving the shock decays rapidly. The semi-moving mesh point method also presents some problems i f the perturbations tend to separate from the main parts of the wavefronts to form separate fronts. The code, 40 as the mesh p o i n t s are f i x e d i n t h e i r p o l a r angles, take i n t o c o n s i d e r a t i o n o n l y wavefronts which are simple c l o s e d curves. Such a breaking o f f c o n d i t i o n w i l l thus cause an abnormality i n the code and cause the e r r o r d e t e c t i o n system i n the code to terminate the c a l c u l a t i o n . However, before t e r m i n a t i o n o c c u r s , the l a s t c a l c u l a t e d wavefront i s p l o t t e d and the break o f f processes are u s u a l l y v i s i b l e i n t h i s l a s t p l o t . The f l e x i b i l i t y of the code i s however more s e r i o u s l y r e s t r i c t e d by the l i m i t a t i o n of the wavefront method i t s e l f . Since the flow f i e l d behind the f r o n t i s not c a l c u l a t e d , and s i n c e f o r a non-uniform converging shock the flow behind the f r o n t w i l l be both n o n - i s e n t r o p i c and non-uniform, i t i s not an easy task to c a l c u l a t e the p h y s i c a l q u a n t i t i e s of the gas flow behind the f r o n t at any time i n s t a n t using o n l y the i n f o r m a t i o n provided by the wavefront c a l c u l a t i o n . Thus the code i s not a p p l i c a b l e to m u l t i p l e shock or shock-heat f r o n t c a l c u l a t i o n s . 41 CHAPTER 5 RESULT OF NUMERICAL COMPUTATION A FORTRAN code of the above numerical method was written and the programme was compiled using the *FTN Compiler of the UBC Computing Center. The computer used i s the IBM 370/168. A t y p i c a l c a l c u l a t i o n requires about 300 time steps and the computing time i s approximately 5 sec. Results of the wavefront calculations of shock waves in an ideal gas of Y =5/3 were plotted d i r e c t l y using a p l o t t i n g subroutine supported by the UBC Computing Center. The output plot consists of the wavefront at d i f f e r e n t times, the graph of ra d i a l perturbation as a function of the unperturbed radius and the graph of volume perturbation as a function of the unperturbed radius. 5.1 Plane Shock Propagating Into Wedge-like Channel The calculations of a plane shock wave with d i f f e r e n t i n i t i a l Mach number propagating into a wedged channel were done in order to test the code by comparison with experiment. Fig. (5.1) shows the wavefront plots of two t y p i c a l c a l c u l a t i o n s . The number of mesh points in both cases was 60, the i n i t i a l Mach number was M=8.00 and the wedge angle was 20 degrees. The wavefronts were plotted once every 20 time steps. The dimensionless time increments D per time step were held constant at 0.0001 for Fig.(5.l.A) and 0.0002 for 42 DIMENSIONLESS L E N G T H S C A L E Fig. (5.1) Wavefront plot into a wedged channel. (a) D=0.0001 (b) D=0.0002 of a plane shock wave propagating 43 F i g . ( 5 . 1 . B ) . I n e a c h c a s e a Mach s t e m , w h i c h i s t h e p a r t o f a s h o c k f r o n t n e a r a w a l l i n a t h r e e - s h o c k Mach r e f l e c t i o n , 1 was o b s e r v e d t o be f o r m e d . F i g . ( 5 . 2 ) i l l u s t r a t e s t h e Mach numbers a t d i f f e r e n t mesh p o i n t s w i t h t h e p a s s a g e o f t h e s h o c k - s h o c k . F o r D=0.0002, t h e Mach number a t t h e t i p o f t h e Mach s t e m o s c i l l a t e s i n an underdamped f a s h i o n , w i t h t h e p e r i o d o f o s c i l l a t i o n a r o u n d 20 t i m e s t e p s and t h e maximum o v e r s h o o t as 2% o f t h e f i n a l s t a b l i z e d v a l u e . A f t e r t h e t i p o f t h e Mach stem h a s r e a c h e d t h e u p p e r w a l l , a s e c o n d Mach stem o r i g i n a t i n g f r o m t h e u p p e r w a l l d e v e l o p s and a g a i n p r o p a g a t e s a l o n g t h e f r o n t . T h i s t i m e however t h e d i m e n s i o n l e s s t i m e i n t e r v a l i s t o o l o n g when com p a r e d t o t h e Mach numbers and t h e s e p a r a t i o n s o f t h e mesh p o i n t s . N u m e r i c a l i n s t a b i l i t y s e t s i n and t h i s c a n be s e e n i n F i g . ( 5 . 2 . B ) as a much l a r g e r a m p l i t u d e o s c i l l a t i o n i n t h e v a l u e o f t h e Mach numbers. A t t h e 3 5 2 t h t i m e s t e p , t h e i n s t a b i l i t y f i n a l l y c a u s e d t h e t r a j e c t o r i e s o f t h e mesh p o i n t s lA Mach r e f l e c t i o n i s d i f f e r e d f r o m t h e r e g u l a r s h o c k r e f l e c t i o n i n w h i c h t h e r e f l e c t i o n o c c u r s away f r o m a. w a l l . From t h e p o i n t o f r e f l e c t i o n , a l s o c a l l e d t h e t r i p l e p o i n t , t o t h e w a l l i s a new s h o c k , c a l l e d a Mach s t e m , w h i c h i s n o r m a l t o t h e w a l l . The t r i p l e p o i n t i s c a l l e d a s h o c k - s h o c k i n s h o c k d y n a m i c s , where t h e r e i s an a b r u p t c h a n g e o f s h o c k s t r e n g t h f r o m t h e i n c i d e n t s h o c k f r o n t t o t h e Mach s t e m . A s e c o n d , much w e a k e r r e f l e c t e d s h o c k i s f o r m e d i n t h e r e g i o n b e h i n d t h e two f r o n t s , o r i g i n a t i n g f r o m t h e s h o c k - s h o c k . T h i s t h i r d s h o c k f r o n t i s c r e a t e d by d i s t u r b a n c e o r i g i n a t i n g f r o m t h e w a l l t r a v e l l i n g i n t o t h e s h o c k e d r e g i o n . 2 A c o m p a r i s o n o f t h e r e s u l t w i t h i n f i n i t e s t r o n g s h o c k t h e o r y and w i t h e x p e r i m e n t i s p l o t t e d i n F i g . ( 6 . 4 ) . rn _ CC LU CO JO o cr cn ( A ) D= 0 . 0 0 0 1 0.0 NO —I 20.0 OF T I M E 40.0 S T E P (XlO 80 a rn. on LU CD •a IE cr a ( B ) D = 0 . 0 0 0 2 0.0 NO 10.0 OF T I M E — I — 20.0 S T E P CX10 30.0 r ) 40 (5.2) Mach numbers at d i f f e r e n t mesh points (a) D=0.0001 (b) D=0.0002 45 to cross over, which was detected by the program and resulting in a termination of the program. The l a s t wavefront plot in Fig.(5.1.B) shows this condition. A shock-shock also exists in the mesh points, which can be i d e n t i f i e d as an abrupt change in Mach number and normal d i r e c t i o n along the mesh points. The f i r s t r e f l e c t i o n of the shock-shock i s observed at the 280th time step, when the shock-shock had travelled a l l the way along the front and h i t the other wall. 5.2 Symmetrical Circular Converging Shock Fig. (5.3) shows the wavefront plot for the geometry of a symmetrical converging cylinder with i n i t i a l Mach number M=8.00. The number of mesh points i s 90. The i n i t i a l time increment D i s 0.0003, but i t i s made to vary during the computation in order to compensate for the decreases in separations between the mesh points. 5 seconds of computing time were required to complete 240 time steps. From the plot in Fig. (5.3) i t can be seen that the wavefront remains p e r f e c t l y c y l i n d r i c a l as i t collapses towards the center. The Mach number of the front as a function of the radius can be used to evaluate the computational error of the code by comparing the obtained dependence with those obtained from d i r e c t numerical integration of the Area Rule and the i n f i n i t e l y strong shock approximation. Fig(5.4) shows the result of such comparison. Exact agreement i s obtained between the wavefront method 4 6 CM D I M E N S I O N L E S S L E N G T H UNIT F i g . (5.3) Wavefront plot of a c y l i n d r i c a l converging shock front with D=0.0003. and the d i r e c t integration of the area rule, while the discrepancy between these two methods and the i n f i n i t e l y strong shock approximation i s less than 2%. These results indicate that i f there i s no discontinuity in the wavefront, the d i s c r e t i z a t i o n of the wavefront and the approximation of the area change have caused n e g l i g i b l e computational error. An estimation of the computational error caused by the relocation of the mesh points after each time step i s achieved by displacing the o r i g i n of the polar coordinates 0.2 units from i t s natural positsion, i . e . at the center of the wavefront. This displacement created an asymmetrical grid point d i s t r i b u t i o n in an otherwise symmetrical geometry, hence 47 comparison with p r e v i o u s c a l c u l a t i o n can g i v e an i n d i c a t i o n of the magnitude of the e r r o r caused by r e l o c a t i o n of the mesh p o i n t s . The r e s u l t of the comparison which i s not p l o t t e d s i n c e the d i f f e r e n c e are l e s s than 0.01% again i n d i c a t e d that the r e l o c a t i o n of the mesh p o i n t s has caused n e g l i g i b l e e r r o r . 5.3 C i r c u l a r Converging Shock Wave With P e r t u r b a t i o n The e v o l u t i o n of shock f r o n t s with three d i f e r e n t kinds of p e r t u r b a t i o n s were i n v e s t i g a t e d . In the f i r s t case, the p e r t u r b a t i o n amplitude i s a s i n s u i d a l displacment superimposed onto a c i r c u l a r wavefront, with the Mach number along the f r o n t being constant throughout. The second case i s s i m i l a r i n shape to the f i r s t case, but the Mach number of the f r o n t i s not consta n t , but i s rather given as f u n c t i o n of the l o c a l r a d i u s of cu r v a t u r e using Eg.(3.16a) f o r i n f i n i t e l y s t rong shock approximation. In the t h i r d case, a p a r t of the c i r c u l a r wavefront i s repl a c e d by an arc of a smaller c i r c u l a r wavefront. Thus u n l i k e the f i r s t two cases there i s an abrupt change i n the curv a t u r e of the f r o n t . The Mach number at the center of the p e r t u r b a t i o n again f o l l o w s the i n f i n i t e l y s trong shock approximation(Eq.(3.16a)), and there i s a l i n e a r t r a n s i t i o n of the Mach number from the center of the p e r t u r b a t i o n to the edge of the p e r t u r b a t i o n , where the Mach number i s the same as the unperturbed f r o n t . The l i m i t of s t a b i l i t y f o r the three d i f f e r e n t cases of p e r t u r b a t i o n s were i n v e s t i g a t e d based on the f o l l o w i n g three 48 A WAVEFRONT CALCULATION • NUMERICAL INTEGRATION OF AREA RULE • INFINITELY STRONG SHOCK APPROXIMATION o.s ~ i 1 1 r— 0.6 0.7 0.8 0.9 R A D I U S OF F R O N T 1.0 F i g . (5.4) Comparison of Mach number i n c r e a s e between Wavefront melt'hod with I n f i n i t e l y Strong Shock Approximation. 49 c r i t e r i a for s t a b i l i t y : 1) . Shape stability.(Wavefront does not break up.) 2) . P a r t i a l r a d i a l s t a b i l i t y . ( A r / r decreases.) 3) . P a r t i a l area stabilty.(AA/A decreases.) The second and the th i r d c r i t e r i a can be obtained d i r e c t l y from the wavefront calculations, while the f i r s t condition must be observed and interpreted from the wavefront plots because the code i s not designed to handle broken up wavefronts. Careful inspection of time evolutions of the three d i f f e r e n t cases of wavefronts revealed that in a l l cases, instead of disappearing, the perturbations spreaded out along the fronts. The perturbation angles were increased while the perturbation amplitudes were decreased. However, depending on the i n i t i a l conditions, PERT, and PERT would A r either increase or decrease. A wavefront plot of a t y p i c a l perturbation of the f i r s t case i s shown in Fi g . (5.5). The number of mesh points was 90 and the i n i t i a l time interval D was 0.0003. The perturbation had an i n i t i a l amplitude of 0.1 and substended an angle of y =20 degrees. 120 time steps were calculated and the wavefront was plotted once every 20 time steps. From the wavefront p l o t , i t can be seen that the perturbation spreaded out and decreased in amplitude as time progressed. At the 120th time step, the perturbation amplitude was so small that the wavefront became roughly c i r c u l a r in shape, 50 2.25 2.5 F i g . (5.5) Wavefront p l o t of a t y p i c a l p e r t u r b a t i o n of the f i r s t case ( s i n s u s o i d a l p e r t u r b a t i o n with constant shock strength) i n d i c a t i n g that the i n i t i a l p e r t u r b a t i o n was s t a b l e . In order to a p p r e c i a t e the l i m i t of s t a b i l i t y f o r p e r t u r b a t i o n s of the f i r s t c a s e ( s i n u s o i d a l p e r t u r b a t i o n superimposed onto a c i r c u l a r wavefront with constant Mach number), the p e r t u r b a t i o n angle and p e r t u r b a t i o n amplitude Ar were changed, p a r t i a l area s t a b i l i t y was found f o r combinations of y and Ar i n the r e g i o n l e f t of the broken l i n e marked by open t r i a n g l e s on F i g . ( 5 . 6 ) . S i m i l a r l y the c o n d i t i o n f o r p a r t i a l shape s t a b i l i t y and p a r t i a l r a d i a l s t a b i l i t y are t r a c e d i n F i g . ( 5 . 6 ) . Together then one f i n d s a region of s t a b i l i t y which s a t i s f i e s a l l three c r i t e r i a . 51 10 20 30 40 50 60 PERTURBATION ANGLE 7 IN DEGREE F i g . (5.6) Limit of s t a b i l i t y for perturbations of the f i r s t case (sinsusoidal perturbation with constant shock strength) Fig. (5.7) shows the wavefront plot of a perturbation of the second case. The i n i t i a l perturbation amplitude was 0.3 and the perturbation angle was 20 degrees. Inspection of the plot revealed that at the 160th time step a Mach stem appeared at the center of the perturbation (y=0). This Mach stem grew in length as time progressed, and the wavefront did not become c i r c u l a r . This perturbation is however stable because both PERTA and PERTr decreased as functions of radius. Also the existence of the Mach stem had 52 F i g . (5.7) Wavefront p l o t of a t y p i c a l p e r t u r b a t i o n of the second case ( s i n s u s o i d a l p e r t u r b a t i o n with non-uniform shock strength) prevented the p e r t u r b a t i o n from breaking o f f . The l i m i t o f s t a b i l i t y f o r the second case of p e r t u r b a t i o n i s p l o t t e d i n F i g . (5.8). N o t i c e that f o r p e r t u r b a t i o n s of f i r s t and second cases, even though that they may have the same i n i t i a l shape, the s t a b i l i t y l i m i t s are d i f f e r e n t , due to the f a c t that the i n i t i a l Mach number d i s t r i b u t i o n s were d i f f e r e n t . Obviously the p e r t u r b a t i o n of the second case i s more s t a b l e area and amplitude-wise, due to the higher Mach number at the p e r t u r b a t i o n , yet they are s l i g h l y more prone to shape i n s t a b i l i t y because of the same. For the p e r t u r b a t i o n of the t h i r d c a s e ( c i r c u l a r 53 L U PERTURBATION ANGLE / IN DEGREE F i g . (5.8) L i m i t of s t a b i l i t y f o r p e r t u r b a t i o n s of the second case ( s i n s u s o i d a l p e r t u r b a t i o n with non-uniform shock strength) p e r t u r b a t i o n superimposed onto a c i r c u l a r f r o n t ) , a t y p i c a l wavefront p l o t i s shown i n Fig.(5.9) and the l i m i t of s t a b i l i t y i s shown i n F i g . ( 5 . 1 0 ) . The p e r t u r b a t i o n i n F i g . (5.9) had an i n i t i a l p e r t u r b a t i o n amplitude of 0.3 and p e r t u r b a t i o n angle of 36 degrees. From the o b s e r v a t i o n of the p l o t t e d wavefront, i t seemed l i k e l y t h a t at some l a t e r time a small p a r t of the p e r t u r b a t i o n would broke away from the parent wavefront, t h e r e f o r e the p e r t u r b a t i o n i s termed shape u n s t a b l e . Observations on d i f f e r e n t p e r t u r b a t i o n parameters on t h i s kind of p e r t u r b a t i o n i n d i c a t e s that they 54 2.5 F i g . (5.9) Wavefront plot of a t y p i c a l perturbation of the th i r d case (circular perturbation with non-uniform shock strength) are more susceptable to shape i n s t a b i l i t y than the f i r s t and second kind. The l i m i t s of s t a b i l i t y for the three d i f f e r e n t cases of perturbations are a l l plotted in Fig.(5.11). There are a common region in the figure that is stable to the three d i f f e r e n t perturbations. Since the i n i t i a l conditions of the three perturbations cover a wide variety of perturbations, we can say that, the stable region in F i g . (5.11) i s stable to most perturbations. Several conclusions can be drawn from inspecting the 55 UJ PERTURBATION ANGLE / IN D E G R E E F i g . (5.10) L i m i t of s t a b i l i t y f o r p e r t u r b a t i o n s of the t h i r d case ( c i r c u l a r p e r t u r b a t i o n with non-uniform shock strength) l i m i t s of s t a b i l i t y f o r the three c r i t e r i a of the three cases of p e r t u r b a t i o n s . R a d i a l and area s t a b i l i t y are c l o s e l y r e l a t e d , and they depend very s i g n i f i c a n t l y on the i n i t i a l d i s t r i b u t i o n of Mach number along the f r o n t . A p e r t u r b a t i o n with a l a r g e p e r t u r b a t i o n angle i s more unstable r a d i a l l y and area-wise than a p e r t u r b a t i o n with a small p e r t u r b a t i o n angle, s i n c e f o r a l a r g e p e r t u r b a t i o n angle, the stren g t h e n i n g of the f r o n t due to area convergence i s l e s s than that f o r a small p e r t u r b a t i o n angle. However, p e r t u r b a t i o n s with small p e r t u r b a t i o n angles w i l l tend to 56 L U PERTURBATION ANGLE / IN DEGREE F i g . (5.11) L i m i t s of s t a b i l i t y of the three cases of pe r t u r b a t ions break o f f e a r l i e r , causing shape i n s t a b i l i t y . I t i s t h e r e f o r e unfortunate that the l i m i t f o r shape s t a b i l i t y can only be determined from the o b s e r v a t i o n s i n which some s u b j e c t i v e judgements are r e g u i r e d , s i n c e the code cannot handle d i s c o n t i n u o u s wavefronts and furthermore too much computing time i s r e q u i r e d before one can a c t u a l l y see the wavefront i s a c t u a l l y breaking up. 57 CHAPTER 6 EXPERIMENT An experiment was designed, the apparatus was c o n s t r u c t e d and measurements were performed to t e s t the v a l i d i t y of the computer r e s u l t . A converging shock wave was produced by l e t t i n g a plane shock wave propagate i n t o a converging wedged channel with a r e c t a n g u l a r c r o s s s e c t i o n . The changes i n d i r e c t i o n of propagation and the change i n Mach number of the r e f l e c t e d shock were obtained from framing camera photos f o r comparison with the computer r e s u l t . The shock was produced i n an e l e c t r o t h e r m a l shock tube with a r e c t a n g u l a r c r o s s s e c t i o n . A two-dimensional converging shock geometry was used i n s t e a d of a c i r c u l a r imploding chamber because of s e v e r a l reasons. F i r s t l y i n the converging channel, a shock-shock w i l l be formed at the wavefront. T h i s d i s c o n t i n u i t y at the f r o n t w i l l g i v e a very s e n s i t i v e comparison betweeen the numerical and the experimental r e s u l t . At any numerical d i s c o n t i n u i t y , e r r o r from numerical s i m u l a t i o n w i l l be a m p l i f i e d , and t h i s w i l l show up at the wavefront p l o t as o s c i l l a t i o n s along the wavefront or as u n p h y s i c a l l y shaped wavefronts. Secondly, a symmetrical imploding wave i s r e l a t i v e l y d i f f i c u l t to generate and r e p r o d u c i b i t i y i s e s p e c i a l l y , d i f f i c u l t to achieve. T h i r d l y , s i n c e a c i r c u l a r imploding shock i s produced by d i v i d i n g a plane shock f r o n t and bending the separated s e c t i o n towards each o t h e r ( P e r r y / 1 8 / ) , i t i s not c l e a r how the 58 turbulence produced by the bending w i l l a f f e c t the a r t i f i c i a l l y induced perturbation. 6.1 Shock Tube And Vacuum System An isometric diagram of the shock tube i s shown in Fig.(6.1a), while a cross-sectional diagram of i t is shown in Fig.(6.1b). The rectangular shock tube is constructed from three pieces of one inch thick l u c i t e plates sandwiched together. The center piece has a channel of one inch width along i t to act as the shock tube of one inch square cross section. The shock tube i s uniform for a length of 8 inches, after which one side of the channel is bent 20 degrees inward to from a wedge-shaped end section. Two tungsten electrodes are mounted oppositely in the other end. The f i l l i n g gas i s Argon and the i n i t i a l f i l l i n g pressure i s in the range of 0.5 to 10 t o r r . The pressure i s measured by a MKS Intru. Inc. Type 221M capacitative pressure transducer. The output of the transducer i s attenuated and measured by a d i g i t a l voltmeter, which displays the pressure d i r e c t l y . An automatic gas f i l l i n g valve i s controlled by a d i g i t a l comparator which compares a preset reading with the output of the d i g i t a l voltmeter. This automatic gas f i l l i n g system produces a presettable pressure in the shock tube which i s reproducible to 0.2%. In addition a small trace of acetylene can be added in order to make the shock front more v i s i b l e . F i g . (6.1a) Iso m e t r i c diagram of the shock tube F i g . (6.1b) Cross s e c t i o n a l diagram of the shock tube 61 6.2 Power Source And Tr igger C i r c u i t The power source i s an unbalanced lumped transmission l i n e constructed of 5 microfarad capacitors and 5 microhenry inductors. The capacitor can be charged to a preset voltage of up to 30 KV by a saturable transformer and an automatic charging control. A National Electronics high voltage switching ignitron model NL-1037 i s used in the output of the capacitor bank. The discharge current then travels along a RG-8/U coaxial cable and discharges in the shock tube through tungsten electrodes. The t y p i c a l discharge current is a f l a t pulse which ranges from 2KA to 37KA for a duration of 80 microseconds. The ignitron is triggered by a krytron trigger unit which i s in turn controlled by a delay generator capable of generating delays of up to 100 millisecond. A schematic diagram of the electronic control c i r c u i t r y i s given in Fig.(6.2) 6.3 Method Of Measurement Space-time measurments of the front position are obtained by an TRW image-convertor framing camera model 1-D f i t t e d with a 5-frame plug-in unit 26-B with delay between frames in the microsecond region. The camera i s focused onto a plane midway between the front plate and the backplate of the 3-piece l u c i t e plate. The framing camera i s triggered by a +300v pulse produced 62 M A N U A L TRIGGER • 1 0 V DELAY GENERATOR Jl • 6 0 V P R E S E T T A B L E HIGH VOLTAGE S U P P L Y CAPACITOR BANK [TRW 4 6 A i D E L A Y GENERATOR KRYTRON TRIGGER UNIT • 5 K V IGNITRON S W I T C H • 3 0 0 V ITRW 2 6 B 5 - F R A M E P L U G - I N UNIT TRW 1D I M A G E C O N V E R T O R C A M E R A L l SHOCK TUBE F i g . (6.2) Schematic diagram of the e l e c t r o n i c c o n t r o l c i r c u i t r y 6 3 from a TWR 46A delay generator which provides a variable time delay from 0.5 microsecond to 100 microsecond. Exposure times of 50,100,200 or 500 nanoseconds per frame can be used. The inter-frame delay and the exposure time are independently adjustable on the camera. The TWR 46A is in turn triggered by the delay generator of the power source, so that the space-time measurement i s synchronized by the trigger of the discharge.(Fig.(6.2)) Pictures of the discharge were taken using Polaroid land fil m Type 107 for quick observations. For accurate measurements, the discharges are recorded by Kodak Royal Pan RX120 negative and developed by Acufine f i l m developer giving a f i l m speed of ASA 3000. The speed of the unperturbed front and the Mach stem was obtained from the positions of the wavefront as determined on the large scale enlargements of the framing camera negatives. 64 6.4 Result Various f i l l i n g pressure and discharge voltages were trie d and i t was discovered that f i l l i n g pressure of 3 to 8 torr gave the most v i s i b l e shock front with a discharge voltage of 10 KV, equivalent to a discharge current of 20 KA. The framing camera photos of the discharge showed a plane shock front followed by a heat front propagating down the uniform section with a speed of approximately 3km/s. A ty p i c a l picture is shown in Fig.(6.3). The i n i t i a l f i l l i n g pressure was 5 torr argon and the charging voltage was 10 KV. The discharge current was measured to be 20KA. The series of photographs shows c l e a r l y the development of the shock and heat wave. Inter-frame delays are 1 microsecond and the exposure times are 500 nanoseconds. The speed of the shock wave in the uniform region i s 2.98+0.15km/s corresponding to a Mach number of M=9.37+0.47. The plane shock front can be seen as a thin luminuous region and i s very sharply defined. The heat front i s a bright luminuous cloud following the shock front. Its boundary i s not well defined. As the shock front enters the nonuniform region, a Mach stem appears, indicating the existence of Mach r e f l e c t i o n . At later time i n t e r v a l s , the Mach stem is observed to grow and spread along the unperturbed front. A shock-shock can also be seen at the intersection of the Mach stem and the unperturbed front. The Mach stem is observed to have a 6 5 Fig. ( 6 . 3 ) Framing camera picture of shock front produced from a discharge current o f 2 0 K A in 5 torr argon. I n t e r - f r a m e d e l a y : l y s e c . 66 higher v e l o c i t y than the unperturbed front. In the pictures the v e l o c i t y i s 3.42+0.32km/s, equivalent to M=10.77+1.00. Due to the higher temperature attained because of the higher Mach number, the Mach stem can be observed to be brighter than the unperturbed front in the photograph. The r a t i o of the v e l o c i t y of the Mach stem and the ve l o c i t y of the unperturbed front for diff e r e n t i n i t i a l Mach numbers were measured and compared with theory. The results are plotted in Fig(6.4). There i s approximately a 5% discrepancy between the experiment and the wavefront c a l c u l a t i o n . This discrepancy can be reduced i f one considers in the cal c u l a t i o n the lowering of g caused by ionization at and behind the shock front, and also the boundary e f f e c t of the walls of the shock tube. Notice that although the Mach stem has a higher speed than the unperturbed front, the subsonic heat front following the shock front has a higher speed in the region where the shock was not perturbed than the region where the Mach stem has formed due to higher gas density achieved by the Mach stem region. This means that although there i s a s t a b i l i z i n g mechanism at work to restore the symmetry of the shock front, the perturbation i t caused in the flow behind the front perturbes the subsonic heat front, which i s the driving force for the shock front. This i s very undesirable for a multiple step wave compression scheme. However, i f one 67 F i g . (6.4) Comparison of Mach number change due to r e f l e c t i o n between i n f i n i t e l y strong shock approximation, wavefront c a l c u l a t i o n s and experiment. 68 turns t h i s r e s u l t around one can conclude that p e r t u r b a t i o n s of the p o s i t i o n of the subsonic heat wave have l i t t l e e f f e c t on the shape of the preceding shock f r o n t . While the flow between the shock f r o n t and the heat f r o n t i s subsonic, the v e l o c i t y through which i n f o r m a t i o n are c a r r i e d from the heat f r o n t to the shock f r o n t , i . e . the v e l o c i t y c h a r a c t e r i z e d by the C+ c h a r a c t e r i s t i c s , i s only s l i g h t l y higher than the shock f r o n t v e l o c i t y . T h i s means that the inhomogenuity i n the pressure at the head of the heat f r o n t caused by the inhomogenuity i n the f r o n t i t s e l f w i l l be d i s p e r s e d as i t reaches the shock f r o n t , l e a d i n g to a uniform Mach number on the shock f r o n t . 69 CHAPTER 7 CONCLUSION The wavefront calculations gave a good idea of the propagations of wavefronts when various kinds of perturbations are present. They also show that the s t a b i l i t y of wavefronts to perturbations depends to a large extent on the i n i t i a l conditions of the perturbation. The comparison between the numerical res u l t and'the experimental result of plane shock waves propagating into a wedged channel indicates an agreement of within 5%. This was performed with i n i t i a l Mach numbers between 8.0 and 9.0, the range in which the shocks are v i s i b l e . Such good agreement was obtained even though in the numerical method one nelgects ionization effects and boundary e f f e c t . Also the presence of discontinuity in the d i r e c t i o n of propagation of the front has an adverse effect on the res u l t of the numerical c a l c u l a t i o n . Hence one can expect similar or better agreement between the numerical results of the c i r c u l a r converging front with the real physical s i t u a t i o n . For further work, improvements can be made to the computer code for the wavefront c a l c u l a t i o n , p a r t i c u l a r l y in the relocation of the mesh points. The separations between mesh.points can be non-uniform, with more emphasis on the perturbation region. Furthermore, the wavefront can have several centers of convergence, so that the breaking up of wavefront fo r unstable p e r t u r b a t i o n s can be observed. 71 BIBLIOGRAPHY 1. Ahlborn,B. and J.P. Huni. 1969. S t a b i l i t y and Space-Time Measurements of Concentric Detonations. AIAA Journal 1_, 1191-1192 2. Belokon,V.A. et a l . 1964. Entrance of a strong shock wave into a wedge-like cavity. Sov. Phys. JETP 2A, 33-40 3. Brueckner,K.A. et a l . 1974. Hydrodynamic s t a b i l i t y of a laser-driven plasma. Phys. Fl u i d 1_, 1554-1559 4. Butler,D.S. 1954. Converging spherical and c y l i n d r i c a l shocks. Report No. 54/54 Armament Research and Development Establishment, Minstry of Supply. Fort Halstead, Kent. 5. Butler,D.S. 1955. Symposium on blast and shock waves. Armament Research and Development Establishment, Minstry of Supply. Fort Halstead, Kent. 6. Chester,W. 1954. The q u a s i - c y l i n d r i c a l shock tube. P h i l . Mag.(7)45, 1293-1301 7. Chisnell,R.F. 1955. The normal motion of a shock wave through a nonuniform one-dimensional medium. Proc. Roy. Soc. A232, 350-370 8. Courants,R. and K.O. F r i e d r i c h . 1948. Supersonic flow and shock waves. Interscience Publishers, New York. 9. Dorn W.S. and D.D McCracken. 1972. Numerical Methods with Fortran IV case studies. John Wiley and Sons, New York. 10. Freeman,N.C. 1957. On the s t a b i l i t y of plane shock wave. J. Fluid Mech. 6, 469-480 11. Guderley,G. 1942. Starke Kugelige und Zylindrishe Verdichtungstoesse in der Naehe des Kugelmittelpunktes bzw der Zyl inder ache . Luf t f ahr t f or schung 302-312 12. Gurevich,B. and A.A. Rumyanstev. 1970. Propagation of shock waves in a medium of decreasing density. Sov. Phys. JETP 31, 747-749 13. Knystautas,R and J.H. Lee. Spark Ignition of Converging Detonation Waves. AIAA Journal 5, 1209-1211 14. Lapworth,K.C. 1959. An Experimental investigation of the s t a b i l i t y of plane shock waves. J. Fluid Mech. 6_, 469-480 15. Lee,H.L.and B.H.K. Lee. 1965. C y l i n d r i c a l Imploding Shock Waves. Phys. Fluid 8, 2148-2152 72 16. L i g h t h i l l , M . J . 1949. Proc. Roy. Soc. A198, 454-470 17. Murakami,H. et a l . 1977. S t a b i l i t y of Decelerating Shock Wave. Trans. Japan Soc. Aero. Space S c i . 2Q_, 100-107 18. Perry,R.W.and A.Kantrowitz. 1951. The production and s t a b i l i t y of converging shock waves. J. Appl. Phys.22, 878-886 19. Russel,D.A. 1966. Shock Wave Strengthening by area convergence. J. F l u i d . Mech. 2_7, 305-314 20. Swan,G.W. and G.R. Fowles. 1975. Shock Wave S t a b i l i t y . Phys. Fl u i d 18, 28-35 21. Whitham,G.B. 1957. A new approach to problems of shock dynamics, Part I. Two dimensional problems. J. F l u i d . Mech. 2, 146-171 22. Whitham,G.B. 1958. On the propagation of shock waves through regions of non-uniform area or flow. J. F l u i d . Mech. 4, 337-360 23. Whitham,G.B. 1973. Linear and Nonlinear Waves. John Wiley and Sons, New York. 73 APPENDIX A REPEATED LINEAR INTERPOLATION In p r i n c i p a l , given n p o i n t s i n a x-y plane one can f i t a n-1 degree polynomial through the e n t i r e s e t s of p o i n t s , t h e r e f o r e i f one wishes" to o b t a i n by i n t e r p o l a t i o n the value of y corresponding to some values of x, one can use the g i v e n value of x on the n-1 polynomial to f i n d the value of y. This method however i s very time consuming for l a r g e n. Using three p o i n t s (x^,y^) whose values " s t r a d d l e " the value x w i l l g i v e a second order approximation value of y with great time s a v i n g . One of these second order i n t e r p o l a t i o n s i s the repeated l i n e a r i n t e r p o l a t i o n . Take three p o i n t s ^ k - l ' ^ k - l * ' ( xk'yk' a n d ^k+l'Yk + l * w h e r e The procedure r e q u i r e s 3 l i n e a r i n t e r p o l a t i o n s , the f i r s t two using the given value and the l a s t one using the i n t e r p o l a t e d v a l u e . The f i r s t two i n t e r p o l a t i o n s g i v e the i n t e r p o l a t e d values of y from (x k_ 1 # y k_ ^ , U k , y k ) and ( x k , y k ) , ( X f c + l ' V k + l ) resp. (A.2) (A.3) 74 In f i g u r e ( A . l ) the three g i v e n p o i n t s are the two i n t e r p o l a t e d values y^}, y 5 ^ of L l and L2 represent k - l k + 1 the l i n e s of i n t e r p o l a t i o n . Note that y £ * j and are both l i n e a r f u n c t i o n s of x. We i n t e r p o l a t e d once more using (x^ I'^k^l^ a n d ^ xk + l'y"k+l^ as p o i n t s . The r e s u l t of the i n t e r p o l a t i o n i s shown by l i n e L3 i n f i g . ( A . l ) . The p o i n t P on L3 d i r e c t l y above x i s the r e s u l t of t h i s f i n a l i n t e r p o l a t i o n . The y-coo r d i n a t e P i s y£ and i s given by (2) Pk+1 y k - l \ , + (A.4) Yk W+l~\-l k k Note .that y^ j^is no longer a l i n e a r f u n c t i o n of x, i t i s a q u a d r a t i c i n t e r p o l a t e d polynomial that passes through the three p o i n t s ( x k . 1 f y k . 1 ) # ^ k ' 7 ^ ' a n d ^ k + l ^ k + l 1 . F i g . (A.l) Repeated Linear I n t e r p o l a t i o n 76 APPENDIX B RUNGE-KUTTA INTEGRATION METHOD The R u n g e - K u t t a method i s a s e c o n d o r d e r i n t e g r a t i o n method u s i n g f o r w a r d and b a c k w a r d s t e p s . I t s a d v a n t a g e s i n c l u d e 1) I t d o e s n o t r e q u i r e t h e p a s t h i s t o r y o f t h e f u n c t i o n t o be i n t e g r a t e d . O n l y t h e v a l u e o f t h e f u n c t i o n a t t h e p o i n t o f i n t e g r a t i o n i s r e q u i r e d . 2) I t p e r m i t s an e a s y c h a n g e i n t h e s t e p s i z e o f t h e i n t e g r a t i o n . S u p p o s e a d i f f e r e n t i a l e q u a t i o n and a s o l u t i o n y a t p o i n t x=x i s g i v e n , m y - y m ( B . l ) d f ' ' V . ' , B - 2 ' t h e n a l i n e L l c a n be drawn w i t h s l o p e y' p a s s i n g t h r o u g h ( x m , y m ) and < x m + 1 , y < l j ) , where y ( 1 \ = f ( x ,y ) ( x ..-x ) + y (B.3) 3 m+1 m'^m m+1 m Jm t h i s i s t h e f o r w a r d i n t e g r a t i o n s t e p , o n c e Y^t\ ^s o b t a i n e d , t h e n t h e s l o p e L2 a t (x ,,V^1) c a n be e v a l u a t e d ' r m+1 Jm+1 d (1) a T 1 - ^ ! . ^ ( B-4 ) t h e s l o p e o f t h e l a s t i n t e g r a t i o n s t e p i s t h e a v e r a g e o f 77 t h e f o r w a r d and t h e b a c k w a r d s t e p s . (2) dx (B.5) (2) (B.6) 78 F i g . (B.l) Runge-Kutta I n t e g r a t i o n 79 APPENDIX C REMARK ON ACCURACY A numerical s i m u l a t i o n f o r c a l c u l a t i n g the wavefront p o s i t i o n of a no n l i n e a r problem can be best termed an approximation to an otherwise i n a c e s s i b l e s o l u t i o n . The e r r o r s i n the numerical s o l u t i o n can be c l a s s i f i e d i n t o two major c a t e g o r i e s . I f we consid e r three q u a n t i t i e s , D, X and N. D denotes the exact s o l u t i o n of the no n l i n e a r equation, X denotes the s o l u t i o n from a " p e r f e c t " computer and N denotes the a c t u a l numerical s o l u t i o n of the problem. Then the e r r o r i s D-N= (D-X) + (X-N) The e r r o r has two p a r t s , D-X i s the i n s t a b i l i t y caused by s e t t i n g up the n o n l i n e a r problem with the p a r t i c u l a r a l g o r i t h m being used. T h i s i s caused by the t r u n c a t i o n of higher order terms i n any n o n l i n e a r equations i n c o n v e r t i n g them i n t o d i f f e r e n c e equations. The other p a r t , X-N, i s due to d i s c r e t i z a t i o n of time i n t e r v a l and mesh p o i n t s and rounding-o f f of intermediate numerical v a l u e s . T h i s i s due to the f a c t t h a t we must use a g r i d s i z e and time i n t e r v a l of f i n i t e value due to the f i n i t e speed and the f i n i t e memory s i z e of the computer. 80 APPENDIX D CODE FOR WAVEFRONT CALCULATION The f o l l o w i n g i s the FORTRAN code for shockfront c a l c u l a t i o n . For c a l c u l a t i o n of d e t o n a t i o n f r o n t s , minor m o d i f i c a t i o n s to the subroutine 'FUNC are r e q u i r e d . The MTS@ Command subroutine 'FTNCMD' i s supported by the UBC Computing Center, and so are the p l o t t i n g subroutines 'ALGRAF 1, 'ALSIZE', 'ALSCAL', 'ALDASH' AND 'PLOTND'. C C C THIS PROGRAM EMPLOYES WHITHAM'S AREA RULE TO C CALCULATE THE TIME EVOLUTION OF A SHOCK FRONT. IT C COMBINES THE CALCULATION OF WEDGED AND CYLINDRICAL C FRONTS. C C THE BOUNDARY CONDITIONS ARE TREATED AS AXIS OF C SYMMETRY TO CALCULATE THE VALUE OF THE BOUNDARY C POINTS. C C THE NEW MACH NUMBERS ARE CALCULATED BY RUNGE-KUTTA C INTGERATION AND THE NEW MESH POINTS ARE CALCULATED C BY REPEATED LINEAR INTERPOLATION. C C R=X,Y,M,T CARTESIAN COORDINATES, STRENGTH, C AND NORMAL OF MESH POINTS. C RN=XN,YN,MN,TN CARTESIAN COORDINATES, C STRENGTH AND NORMAL OF MESH POINTS C AFTER ONE TIME STEP. C RADIUS,ANGLE POLAR COORDINATES OF MESH C POINTS. C PANGLE ANGLE SUBSTAINED BY SHOCK FRONT IN C RADIANS. C TANG ANGLE OF WAVEFRONT IN WEDGED GEOMETRY. C DEL,L SIZE AND L # OF PERTURBATION. C GAMMA RATIO OF SPECIFIC HEATS. C XC CENTER OF POLAR COORDINATES. C D TIME INTERVAL. C J # OF TIME STEPS. C II CONTAINS INFORMATION OF THE SHAPES OF C THE WAVEFRONT. C XLENG,YLENG X,Y VALUES OF RADIAL C PERTURBATION. XAREA,YAREA X,Y VALUES OF AREA PERTURBATION. IMPLICIT REAL*8 (A-H,0-Z) REAL*4 A(200) ,B(200) REAL*4 XLENG(251),YLENG(251),XAREA(251),YAREA(251) DIMENSION X(100) ,Y(100) ,T(100) ,ANGLE(100) DIMENSION XN(100),YN(100),TN(100),RADIUS(100) DIMENSION R(100,4) ,RN(100,4) ,S(100) REAL*8 M(100) ,MN(100) INTEGER*4 YES/'Y'//NO/1N'/,FIT,FRONT COMMON /LA/R /LB/RN /LE/D /LG/GAMMA /LPI/PI /LJ/J /LDE/DE COMMON /LRA/RADIUS /LXC/XC /LP/ANGLE /LL/DEL,L/LU/U COMMON /LS/II /LANGEL/PANGLE /LTANG/TANG /LNN/NN EQUIVALENCE (R(l) ,X) , (R(101) ,Y) ,(R(201) ,M) , (E(301) ,T) EQUIVALENCE (RN(1) ,XN) , (RN(101) ,YN) ,(RN(201) ,MN) EQUIVALENCE (RN(301),TN) INTEGER FMT(1)/'* 1/ LOGICAL UNIT 2 IS DEFAULTED TO 'REF.DATA1 , WHICH CONTAINS THE DATA OF THE REFERENCE CIRCLE. LOGICAL UNIT 5 IS DEFAULTED TO *SOURCE*. IT RECEIVES THE NECESSARY INPUT PARAMETERS. LOGICAL UNIT 6 IS ASSIGNED TO A SCRATCH FILE "-FILE", THIS CONTAINS ALL THE PRINTED OUTPUT. LOGOCAL UNIT 7 IS DEFAULTED TO *SINK*. IT PROMPTS FOR THE NECESSARY INPUT PARAMETERS. THE PLOTFILE IS DEFAULTED TO "-PLOT#". CALL FTNCMD('DEFAULT 2=REF.DATA 1,18) CALL FTNCMD('ASSIGN 6=-FILE',14) CALL FTNCMD('DEFAULT 7=*SINK*',16) THE IDEAL GAS CONSTANT 'GAMMA' IS EQUAL TO 5/3. GAMMA=5.D00/3.D00 FORMAT(IX,16A1) WRITE(7,101) THE # OF MESH POINTS IS REQUESTED. 01 FORMAT('0','ENTER N') READ(5,FMT) N NN=N-1 PI=3.14159265358 WRITE(7,102) THE INITIAL TIME INTERVAL 'DD' IS REQUESTED. FORMAT('0','ENTER TIME INTERVAL'/ 'FOR CIRCLE,SUGGEST USE l.D-3' /'FOR WEDGE,SUGGEST USE l.D-4') READ(5,FMT)DD D=DD THE ORIGINAL SHAPE OF THE WAVEFRONT IS READ IN USING EITHER THE SUBROUTINES "WEDGE" OR "CIRCLE". THE VARIABLE 'II' DETERMINES THE SHAPE. IF THE INITIAL SHAPE IS A WEDGE,A '2' IS ENTERED, IF THE INITIAL SHAPE IS A CIRCLE, A '1' IS ENTERED. WRITE (7,5) FORMAT('0','CHOOSE SHAPE REQUIRED'/ 'ENTER FOR CIRCLE', ' AND '' 2'' FOR WEDGE') READ(5,6) II FORMAT(II) IF THE INITIAL SHAPE IS A WEDGE, FURTHER INFORMATION ABOUT THE SHAPE OF THE PERTURBATION IS NOT NEEDED, SO THE PROGRAM GOES DIRECTLY TO MACH STRENGTH INITIALIZATION AND PLOTTING OF THE INITIAL SHAPE. IF (II.NE.l) GOTO 400 WRITE(7,880) FORMAT('0','DO YOU WANT ABRUPT CHANGE IN FRONT?'/ 'ENTER ''Y''1,1 OR '1N''') READ(5,881) FRONT FORMAT(A4) THE TWO DIFFERENT TYPES OF CIRCULAR CONVERGING FRONT WITH PERTURBATION IS DETERMINED. AN ANSWER 'Y' WILL GIVE A CIRCULAR PERTURBATION WITH SUBROUTINE 'ABRUPT' WHILE 'N' WILL GIVE A SINUSOIDAL PERTURBATION WITH SUBROUTINE 'CIRCLE'. IF (FRONT.EQ.NO) CALL CIRCLE(X,Y,T,M,N) IF (FRONT.EQ.YES) CALL ABRUPT(X,Y,T,M,N) THE SIZE OF THE OUTPUT PLOT IS 10" X 5" FOR CYLINDRICAL FRONTS. CALL ALSIZE(10.0,5.0) CALL ALSCAL(0.0,2.5,0.0,1.25) 83 GOTO 401 400 CALL WEDGE(X,Y,T,M,N) C C THE SIZE OF THE OUTPUT PLOT IS C 5" X 5" FOR WEDGED FRONT. C CALL ALSIZE(5.0,5.0) CALL ALSCAL(0.0,1.0,0.0,1.0) C C FOR WEDGED FRONT, 'FIT'='NO' SO THAT A C BEST-FIT CIRCLE WILL NOT BE EVALUATED. C FIT=NO 401 CONTINUE C C THE INITIAL SHAPE IS PLOTTED. C CALL DATAPL(X,Y,N,0) WRITE(7,100) C C THE # OF TIME STEPS 'IT' IS REQUESTED. C 100 FORMAT ('1','ENTER NUMBER OF TIME STEPS') READ (5,FMT) IT C C "GEOM" IS CALLED TO INITIALIZE THE C VALUES OF T. C CALL GEOM(X,Y,T,N) IF(II.EQ.1) READ(2,15)D,AREARE,RADREF C C FOR CYLINDRICAL WAVEFRONTS, C THE AREA 'AREARE' AND THE RADIUS C 'RADREF' OF THE REFERENCE CIRCLE C IS READ. C THE VALUES 'AREARE' AND 'RADREF' C WERE COMPUTED USING THIS SAME PROGRAM C WITH NO PERTURBATION AND STORED IN C A STORAGE FILE. C 'D', THE TIME INTERVAL IS ALSO READ. C 15 FORMAT(1X,F10.8,2X,F10.7,2X,F10.8) DO 2 1=1,N 2 RADIUS(I)=1. C C THE CALCULATION IS STARTED IN THE C FOLLOWING DO-LOOP WHICH CALCULATES C THE NEW CHARACTERISTIC POINTS, THE C NEW SHOCK STRENGTH AT THE CHARACTER-C ISTIC POINTS AND THEN THE NEW MESH C POINTS. C DO 20 J=1,IT 84 C FOR CYLINDRICAL WAVEFRONTS, 'AREARE' C AND 'RADREF' IS READ TO CALCULATE C THE SIZE OF THE PERTURBATION. C IF (II.EQ.l) READ(2,15)D,AREARE,RADREF C C THE NEW POSITION OF THE MESH POINTS C IS CALCULATED. C DO 25 1=1,N XN (I)=X(I)+M(I)*DCOS(T(I))*D 25 YN(I)=Y(I)+M(I)*DSIN(T(I) ) *D C C AFTER THE NEW POSITIONS ARE OBTAINED, C "SHOCK" IS CALLED TO CALCULATE THE C NEW MACH NUMBERS. THE RELOCATED MESH C POINTS ARE THEN CALCULATED BY "MESH". C FINALLY THE NEW DIRECTION OF PROPAGATION C ARE CALCULATED BY "GEOM". C CALL SHOCK(N) 300 CALL MESH (N,5,998) CALL GEOM(X,Y,T,N) C C FOR CYLINDRICAL WAVEFRONTS, C THE FUNCTION "AREA" CALCULATES THE C AREA ENCLOSED BY THE WAVEFRONT AND C THE AREA IS COMPARED TO THAT OF THE C REFERENCE FRONT. 'YAREA' IS THE AREA C PERTURBATION. C XAREA(J)=RADREF YAREA(J)=AREA(N)-AREARE YAREA(J)=YAREA(J)/AREARE IF (II.EQ.2) GOTO 131 C C THE LENGTH IS COMPARED TO THAT OF THE C REFERENCE FRONT. 'YLENG' IS THE RADIAL C PERTURBATION. C CALL MAXM(RADIUS,K,N) XLENG(J)=RADREF YLENG(J)=RADIUS(K)-RADREF YLENG(J)=YLENG(J)/RADREF C C THE WAVEFRONT IS PLOTTED WHEN THE TIME C STEPS ARE MULTIPLES OF 20. C 131 IF (MOD(J,20).EQ.O) CALL DATAPL(X,Y,N,0) 20 CONTINUE GOTO 997 C C AFTER THE COMPLETION OF THE REQUIRED C NUMBER OF STEPS, THE RADIAL PERTURB., 85 C THE AREA PERTURB. VS RADIUS CURVE IS C CURVE IS PLOTTED. C 998 CONTINUE CALL DATAPL(X,Y,N,0) 997 CONTINUE CALL ALSCAL(0.0,0.0,0.0,0.0) CALL ALGRAF(XLENG,YLENG,IT,0) CALL ALSCAL(0.0,0.0,0.0,0.0) CALL ALGRAF(XAREA,YAREA,IT,0) 505 CALL PLOTND STOP END C C C THIS SUBROUTINE CALCULATES THE LOCATIONS OF C THE NEW MESH POINTS GIVEN THE CONDITION THAT C THEIR ANGULAR SEPARATION REMAIN THE SAME. C SECOND ORDER REPEATED LINEAR INTERPOLATION C IS USED TO CALCULATE THE LOCATIONS OF THE MESH C POINTS, AND LINEAR INTERPOLATION IS USED TO C CALCULATE THE VALUE OF THE MACH NUMBER AT THE C NEW MESH POINT. C SUBROUTINE MESH(N,*) IMPLICIT REAL*8 (A-H,M,0-Z) COMMON /LA/X,Y,M,T /LB/XN,YN,MN,TN /LP/TH /LXC/XC COMMON /LPI/PI COMMON /LJ/J /LDE/DE /LRA/R /LNN/NN COMMON /LBP/MFIR,MLAS DIMENSION S (100) , SN (100) DIMENSION XN(100),YN(100),MN(100),TN(100) DIMENSION X(100),Y(100),M(100),T(100) DIMENSION TH(100),THN(100),R(100),RN(100) C C C XN,YN INPUT CARTESIAN COORDINATES OF POSITION C RN,THN CALCULATED POLAR COORDINATES FROM INPUT C MN INPUT STRENGTH OF MACH NUMBER C X,Y REDISTRIBUTED CARTESIAN COORDINATES OF C POSITION C TH PRESCRIBED POLAR ANGLE FOR INTERPOLATION C R INTERPOLATED RADIAL DISTANCE C M REDISTRIBUTED STRENGTH OF MACH NUMBER C XC CENTER OF POLAR COORDINATE C S,SN DISTANCE BETWEEN MESH POINT FOR C INTERPOLATION OF MACH NUMBER C MFIR,MLAS MACH NUMBER OF THE BOUDARY POINTS. C C THE POLAR COORDINATES OF THE INPUT C MESH POINTS IS EVALUATED. C DO 5 1=1,N THN(I)=DATAN2(YN(I),XN(I)-XC) 86 5 RN(I)=DSQRT((XN(I)-XC)**2+YN(I)**2) C C THE FOLLOWING DO-LOOP CHECKS TO SEE C IF ANY ABNORMAL CONDITION EXISTS. C IF SUCH CONDITION EXISTSM THEN THE C COMPUTATION IS STOPPED , THE POINT C AT WHICH THE ABNORMALITY EXIST IS C PRINTED, AND THE WAVEFRONT AND THE C PERTURBATION CURVES ARE PLOTTED BY C THE MULTIPLE RETURN 998. C DO 6 1=2,N IF (THN (I) .LT.THN(1-1)) GOTO 6 WRITE(6,601)J,I,THN(I-1),THN(I),THN(I+1), RN(1-1) ,RN (I) ,RN(1 + 1) WRITE(7,601)J,I,THN(I-1),THN(I),THN(I+1), RN(I-l) ,RN(I) ,RN(I+1) RETURN1 601 FORMAT('01,'*****ABNORMAL TERMINATION AT J=',I3, 1 AND AT 1=',13/' THE', 1 FOLLOWINGS ARE VAULES OF', ' THN(I-l),THN(I),AND THN(I+1):' ,3F14.5/' THE FOLLOWINGS ARE VALUES OF' ' RN(I-l) ,RN(I) ,RN(I + 1) ' : ,3F14.10) 6 CONTINUE C C C R ARE OBTAINED BY REPEATED LINEAR C INTERPOLATION OF RN FROM THE ANGLES C THN AND TH. C THE SUBROUTINE "RELINP" PERFORMS C THE INTERPOLATION. C AFTER R IS OBTAINED, M IS CALCULATED C USING LINEAR INTERPOLATION OF THE C DISTANCE BETWEEN R AND RN. THE C DISTANCES ARE S AND SN. C C CALL RELINP(THN,RN,TH,R,N) C C THE NEW MESH POINTS ARE TO BE C REPRESENTED IN CARTESIAN COORDINATES. C DO 10 1=1,N X(I)=XC+R(I)*DCOS(TH(I)) 10 Y(I)=R(I)*DSIN(TH (I) ) DO 20 1=2,NN S(I)=DSQRT((X(I)-(XN(I)-XN(1-1))/2)**2+ (Y(I)-(YN ( D-YN(I-l) )/2)**2) SN(I)=DSQRT((X(I)-(XN(I)-XN(1+1))/2)**2+ (Y(I)-(YN(I)-YN(I + 1) )/2)**2) 20 M(I) = (S(I)*MN(I)+SN(I)*MN(1-1))/(S (I)+SN(I)) C C THE BOUNDARY POINTS ARE NOT INTERPOLATED 87 C C C M(1)=MFIR M(N)=MLAS RETURN END C C C C C C C C C C C SUBROUTINE SHOCK(N) IMPLICIT REAL*8 (A-H,M,0-Z) DIMENSION XN(100) ,YN(100) ,MN(100) , TN(100) ,SM (100) DIMENSION X(100),Y(100),M(100),T(100) COMMON /LA/X,Y,M,T /LB/XN,YN,MN,TN /LG/GAMMA /LU/U COMMON /LNN/NN /LBP/MFIR,MLAS U=1.0 C C C Al C A2 C MA C C C DO 10 1=1,NN A2=DSQRT((XN(I+1)-XN(I))**2+(YN(I+l)-YN(I))**2) A1=DSQRT((X(I+1)-X(I))**2+(Y(I+l)-Y(I))**2) MA=(M(I)+M(I+1))/2 C C FOR THE INTERIOR POINTS, THE NEW MACH C NUMBER 'MA' BETWEEN TWO MESH POINTS C ARE EVALUATED. C IF (I.EQ.l.OR.I.EQ.NN) GOTO 111 CALL RUNGE(MA,Al,A2) MN(I)=MA GOTO 10 C C FOR THE BOUNDARY POINTS, THE INITIAL C VALUE FOR INTEGRATION IS THE MACH C NUMBER AT THE MESH POINT SINCE ITS C VALUE IS NOT TO BE INTERPOLATED. C 111 IF (I.EQ.l) MA=M(1) GOTO 112 INSTEAD THE VALUES ARE OBTAINED DIRECTLY FROM "SHOCK" THROUGH COMMON BLOCK #LBP#. THIS SUBROUTINE IS TO CALCULATE THE NEW SHOCK STRENGTH GIVEN THE AREA RATIO AND THE ORGINAL STRENGTH BY WHITHAM'S AREA RULE. THE INTEGRATION IS DONE BY RUNGE-KUTTA METHOD CALLED FROM SUBROUTINE "RUNGE". THE OLD AREA BETWEEN MESH POINTS I AND I+l. THE NEW AREA BETWEEN MESH POINTS I AND I+l THE STRENGTH AT MID POINTS OF MESH POINT I AND I+l. 88 C C THE MACH NUMBER AT N-1 IS TO BE C INTERPOLATED SO SPECIAL PROCEDURE C IS TAKEN. C CALL RUNGE(MA,Al,A2) MN(NN)=MA MA=M(N) 1 1 2 CALL RUNGE(MA,Al,A2) C C SPECIAL VARIABLE NAMES 1MFIR' AND C 'MLAS' ARE GIVEN TO THE MACH NUMBERS C AT THE BOUNDARY POINTS. C IF (I.EQ.l) MFIR=MA IF (I.EQ.NN) MLAS=MA 1 0 CONTINUE RETURN END C C C THE SUBROUTINE "FUNC" CONTAINS THE C DERVATIVE OF MACH NUMBER AS FUNCTION OF AREA C TO BE USED IN RUNGE-KUTTA INTEGRATION PROCEDURE. C THIS IS THE WHITHAM'S AREA RULE OF SHOCK WAVE. C MS DEPENDENT VARIABLE,MACH NUMBER C DF DERIVATIVE C A INDEPENDENT VARIABLE,AREA C SUBROUTINE FUNC(MS,DF,A) IMPLICIT REAL*8(A-H,M,O-Z) COMMON /LG/GAMMA C C 'GAMMA' IS THE IDEAL GAS CONSTANT C IT IS AVAIBLE IN THE COMMON BLOCK #LG# C AND IS EQUAL TO 5/3. C YL SIGMA SQUARED C YL=((GAMMA-1)*MS**2+2)/(2*GAMMA*MS**2-GAMMA+l) GAM=MS/(MS**2-1)*(1+2/(GAMMA+1)*(1-YL)/DSQRT(YL)) *(1+2*DSQRT(YL)+1/MS**2) DF=-1/(GAM*A) RETURN END C C THIS SUBROUTINE IS THE C RUNGE-KUTTA INTEGRATION PROCEDURE. C THIS PROCEDURE USES MODIFIED EULER'S METHOD C TO ESTIMATE THE MACH NUMBER AT A NEW AREA C GIVEN THE DERVATIVE OF THE MACH NUMBER AS A C FUNCTION OF AREA IN THE SUBROUTINE 'FUNC. C Y FUNCTION TO BE INTEGRATED. C XI,X2 LIMITS OF INTEGRATION. C 89 SUBROUTINE RUNGE(Y,XI,X2) IMPLICIT REAL*8(A-H,M,0-Z) C C DF1 IS THE SLOPE AT XI. C CALL FUNC(Y,DFl,X1) C C THE 1ST APPROXIMATION ' Y l 1 IS OBTAINED. C Y1=DF1*(X2-X1)+Y1 C C THE SLOPE 'DF2' AT Yl IS OBTAINED. C CALL FUNC(Y2,DF2,X2) C C THE COMPLETE INTEGRATION IS DONE BY C THE AVERAGE SLOPE. C Y=(DF1+DF2)/2*(X2-X1)+Y RETURN END C C C THIS SUBROUTINE CALCULATES THE ANGLES OF THE C MESH POINTS USING THE INFORMATION OBTAINED. C THREE MESH POINTS ARE USED TO CALCULATE THE C DIRECTION OF PROPAGATION OF THE MIDDLE MESH C POINT. THEY ARE TREATED TO BE SEGMENTS OF A C CIRCLE. AFTER THE CENTER OF THE CIRCLES ARE C LOCATED, THEN THE NORMALS OF THE CIRCLES CAN C BE FOUND. C SUBROUTINE GEOM(X,Y,T,N) IMPLICIT REAL*8 (A-H,0-Z) COMMON /LPI/PI /LS/II /LANGEL/PANGLE /LXC/XC/LNN/NN REAL*8 X(N),Y(N),T(N) T(1)=0. C C THE FOLLOWING ARE PROCEDURES TO C EVALUATE THE CENTERS OF CIRCLES C FORMED BY THE MESH POINTS. C DO 1 I=2,NN P=X(I)**2+Y(I)**2-X(1-1)**2-Y(1-1)**2 Q=X(I)**2+Y(I)**2-X(1+1)**2-Y(1+1)**2 R=X (I)-X(1-1) S=X(I)-X(1+1) TT=Y(I)-Y(1-1) U=Y(I)-Y(I + 1) V=(R*U-TT*S)*2 C C IF THE CURVATURE OF THE CIRCLE IS TOO SMALL, C THEN THE THREE MESH POINTS ARE TREATED AS 90 C POINTS ON A STRAIGHT LINE. THE NORMAL IS C THEN EVALUATED FROM THE ST. LINE FORMULA. C IF (II.EQ.2.AND.DABS(V).LE.1D-3) GOTO 10 IF (II .EQ.LAND.DABS (V) .LE.lD-9) GOTO 10 C C 'A' IS THE X-COORDINATE OF THE CENTER. C 'B' IS THE Y-COORDINATE OF THE CENTER. C 'TN' IS THE NORMAL OF THE CIRCLE. C A=(P*U-Q*TT)/V B=(P*S-Q*R)/(-V) 20 TN=DATAN2(B-Y(I),A-X(I)) GOTO 2 C C IF THE CURVATURE IS TOO SMALL, THEN THE C NORMAL IS TAKEN TO BE PERPENDICULAR TO C A LINE JOINING THE THREE MESH POINTS. C 10 Al=DSQRT((X(I-l)-X(I))**2+(Y(1-1)-Y(I))**2) A2=DSQRT((X(I+1)-X(I))**2+(Y(I+l)-Y(I))**2) TN=(-DATAN((X(I+l)-X(I))/(Y(I+l)-Y(I)))*Al -DATAN((X(I)-X(I-l) )/(Y(I)-Y(I-l) ) ) *A2) /(A1+A2) C C THE NORMALS ARE CHECKED TO MAKE SURE C THAT THE CHANGES IN DIRECTION DO NOT C EXCEED 90 DEGREE, MAKING SURE THAT C THE WAVEFRONT REMAINS CONVERGING. C 2 IF (DABS(TN-T(I-l)).GE.1.5) TN=TN-PI T(I)=TN 1 CONTINUE C C THE BOUNDARY CONDITION THAT THE FLOW MUST C BE PARELLEL TO THE WALL IS IMPOSED UPON C THE BOUNDARY POINT. C T(N)=-PANGLE RETURN END C C C THIS SUBROUTINE TAKES IN THE DATA STORED C ARRAYS X AND Y WHICH ARE IN DOUBLE PRECISION C AND CONVERT THE VALUES INTO SINGLE C PRECISION, THEN PLOT THEM WITH A SYMBOL J C AT EACH DATA POINT. C SUBROUTINE DATAPL(X,Y,N,J) REAL*8 X(N),Y(N) REAL A(100) ,B (100) C C THE FOLLOWING DO-LOOP DOES THE 91 C CONVERSION INTO SINGLE PRECISION. C DO 10 1=1,N A(I)=X(I) 10 B(I)=Y(I) CALL ALDASH(0.0,0.0,0.0,0 . 0 ) CALL ALGRAF(A,B,-N,J) RETURN END C C C THIS SUBROUTINE GIVES THE INITIAL SHAPE OF THE C WAVEFRONT AS A SEMI-CIRCLE. THE PERTURBATION IS C SMOOTH AND SUPERIMPOSED ON THE CIRCULAR WAVEFRONT C IN SINUSOIDAL FORM. C SUBROUTINE CIRCLE(X,Y,T,SM,N) IMPLICIT REAL*8 (A-H,0-Z) COMMON /LPI/PI /LP/TH /LXC/XC /LL/DELTA,L COMMON /LRA/R /LG/GAMMA /LANGEL/ANGEL REAL*8 TH(100),R(100) REAL*8 XN(100),YN(100),MN(100),TN(100) REAL*8 X(N),Y(N),T(N),SM(N) INTEGER FMT(1)/'*'/,YES/ 1Y'/fNO/'N'/rFORM COMMON /LB/XN,YN,MN,TN WRITE(7,1) 1 FORMAT('0','ENTER L,DELTA, AND INITIAL', . ' STRENGTH OF SHOCK') C C THE AZIMUTHAL NUMBER ' L 1 , C THE PERTURBATION AMPLITUDE 'DELTA' AND C THE INITIAL SHOCK STRENGTH 'SMI' C ARE READ. C READ (5,FMT) L,DELTA,SMI C C 'MM' IS THE LARGEST MESH POINT NUMBER C THAT IS MODIFIED BY THE PERTURBATION. C MM=N*2/L ANGEL=PI C C THE POLAR COORDINATES 'R' AND 1TH' OF C A CIRCLE IS FIRST CALCULATED. C DO 30 1=1,N TH(I)=PI*(N-I)/(N-1) R ( D=1. C C THE MACH NUMBERS AT THE MESH POINTS ARE C INITIALIZED. C SM(I)=SMI 30 CONTINUE 'XC IS THE CENTER OF CONVERGENCE. XC=R(1)+DELTA THE PERTURBATION IS INSERTED INTO THE INITIAL SHAPE OF THE SEMICIRCLE. DO 10 1=1,MM R(I)=DELTA*(DCOS(PI/2*(1-1)/(MM-1)))**2+R(I) 'X', 'Y' ARE THE CARTESIAN COORDINATES OF THE MESH POINTS. DO 25 1=1,N X (I)=XC+R(I)*DCOS(TH(I)) Y(I)=R(I) *DSIN (TH (I) ) CONTINUE THE NORMALS OF THE MESH POINTS ARE EVALUATED. CALL GEOM(X,Y,TN,N) WRITE(7,100) FORMAT('0','DO YOU WANT A NON-UNIFORM STRENGTH?'/ 'ANSWER ''Y'' OR ''N1'') READ (5,101) FORM FORMAT(A4) IF (FORM.EQ.NO) GOTO 102 THE RADIUS OF CURVATURE IS CALCULATED, AND THE MACH NUMBER AS FUNCTION OF CURVATURE IS GIVEN. ANG=DATAN(Y(2)/(X(2)-X(1) ) ) SLEN=DSQRT(Y(2)**2+(X(2)-X(1))**2)/2 'CURV IS THE CURVATURE AT THE CENTER OF THE PERTURBATION. CURV=SLEN/DCOS(ANG) GN=l+2/GAMMA+DSQRT(2*GAMMA/(GAMMA-1)) THE MACH NUMBERS OF THE MESH POINTS AFFECTED BY THE PERTURBATION ARE MODIFIED ACCORDINGLY. SM(1)=SMI*CURV**(-1/GN) DO 103 1=2,MM SM(I)=(SM(1)-SM(N))/MM*(MM-I+1.)+SM(N) THE INFORMATION ABOUT THE INITIAL SHAPE IS WRITTEN INTO OUTPUT UNITS '6' AND '7'. WRITE(7,200) CURV,(SM(I),1=1,10) 93 200 FORMAT(1X,11F10.5) 102 WRITE(6,2) N, L,DELTA,SMI 2 FORMAT('0' CIRCLE'/'NO. OF POINTS=' , 13,2X ,'L=' ,12/ 'INITIAL PERTURBATION=',F7.5/ 'INITIAL SHOCK STRENGTH=',F7.3) RETURN END C C THIS SUBROUTINE CALCULATES THE INITIAL SHAPE OF A C WEDGED CHANNEL WITH INPUT WEDGE ANGLE AND FRONT C ANGLE. C SUBROUTINE WEDGE(X,Y,T,SM,N) IMPLICIT REAL*8(A-H,0-Z) REAL*8 TH(100) REAL*8 X(N),Y(N),T(N),SM(N) INTEGER FMT(1)/'*'/ COMMON /LPI/PI /LP/TH /LXC/XC /LS/II /LANGEL/PANGLE COMMON /LTANG/TANG WRITE(7,1) 1 FORMAT('01,'ENTER ANGLE OF WEDGE, ANGLE OF FRONT, AND', .' INITIAL STRENGTH OF SHOCK') C C THE WEDGE ANGLE, THE FRONT ANGLE, AND THE C MACH NUMBER IS READ. C READ(5,FMT) ANGLE,TANG,SMI PANGLE=ANGLE*P1/180 C C 'A' IS THE X-COORDINATE OF THE NTH C MESH POINT. C A=DTAN(TANG*PI/180)/(DTAN(TANG*PI/180)+DTAN(PANGLE)) IF (ANGLE.EQ.0.0) A=l.-DCOTAN(TANG*PI/180) C C THE CENTER OF CONVERGENCE IS SET TO C BE EQUAL TO 1. XC = 1. C C THE BOUNDARY POINTS ARE FIRST CALCULATED. C X(N)=l.-A Y (N)=X(N)*DTAN(TANG*PI/180) IF (ANGLE.EQ.0.0) Y(N)=1. X (1)=0. Y(1)=0. NN=N-1 C C THE REST OF THE MESH POINTS ARE CALCULATED. C DO 20 1=2,NN X(I)=X(N)*(I-1)/(N-1) Y(I)=Y(N)*(I-1)/(N-1) 20 CONTINUE THE MACH NUMBER OF THE MESH POINTS ARE INITIALIZED. DO 10 1=1,N SM(I)=SMI TH(I)=DATAN2(Y(I),X(I)-XC) THE INITIAL CONDITIONS ARE PRINTED IN THE OUTPUT UNIT '6'. WRITE(6,2) N,ANGLE,TANG,SMI FORMAT('0','WEDGE'/'NO. OF POINTS=' ,I 3,2X/ 'ANGLE OF WEDGE=',F6.2/'ANGLE OF FRONT=',F6.2/ 'INITIAL SHOCK STRENGTH=',F7.3) RETURN END REPEATED LINEAR INTERPOLATION THE INPUTS ARE X AND Y, AND YN IS INTERPOLATED USING THE VALUE OF XN. YN(1) AND YN(N) RETAINED THE VALUE OF Y ( l ) AND Y(N) . SUBROUTINE RELINP(X,Y,XN,YN,N) IMPLICIT REAL*8(A-H,0-Z) DIMENSION X(N) ,Y(N) ,XN (N) ,YN(N) DIMENSION DX(100),DY(100) NN=N-1 DO 10 1=1,NN DX(I)=X(I + 1)-X(I) DY(I)=Y(I+1)-Y(I) DO 20 1=2,NN YKP IS THE FIRST APPROXIMATION ON THE RIGHT. YKN IS THE FIRST APPROXIMATION ON THE LEFT. YN IS THE SECOND APPROXIMATION EMPLOYING THE VALUES OF YKP AND YKN. YKP=Y(I)+(XN(I)-X(I))/DX(I)*DY(I) YKN=Y(1-1)+(XN(I)-X(1-1))/DX(1-1)*DY(1-1) YN(I)=YKN+(XN(I)-X(1-1))/(X(I+l)-X(1-1))*(YKP-YKN) THE BOUNDARY POINTS ARE NOT INTERPOLATED. YN(1)=Y (1) YN(N)=Y(N) RETURN END THE SUBROUTINE "MAXANG" LOCATES THE MAXIMUM ANGLE A AT MESH POINT # K FROM THE SET OF ANGLES T WITH DIMENSION N. SUBROUTINE MAXANG(A,K,T,N) COMMON /LP/TH /LPI/PI REAL*8 T ( 1 0 0 ) , T H ( 1 0 0 ) , U ( 1 0 0 ) , P I REAL*4 A,B 'A' IS INITIALIZED TO BE ZERO. A=0 . DO 10 1 = 1,N U(I)=TH(I)-PI B=DABS(T(I)-U (I)) IF (B.LE.A) GOTO 10 ANY 'B' BIGGER THAN 'A' IS TO BE THE NEW 'A' AND ITS MESH POINT NUMBER IS MEMORIZED. A=B K=I CONTINUE RETURN END THE SUBROUTINE "MAXM" LOCATES THE MAXIMUM OF THE SET A OF DIMENSION N # K. SUBROUTINE MAXM(A,K,N) REAL*8 A(N) THE INITIAL MAX VALUE 'G' IS SET TO TO EQUAL TO ZERO. . G=0.0 DO 10 1=1,N IF (A(I).LT.G) GOTO 10 IF 'A' IS BIGGER THAN 'G 1, THEN 'A' IS THE NEW MAX. G=A(I) K=I CONTINUE RETURN END THE AREA OF THE CIRCULAR WAVEFRONT IS CALCULATED BY EVALUATING THE RIEMANN SUM OF THE AREA FORMULA FOR A SEGMENT OF A CIRCLE. REAL FUNCTION AREA*8(N) COMMON /LRA/RADIUS /LP/ANGLE REAL*8 RADIUS(100),ANGLE(100) AREA=0. DO 10 1 = 2,N AREA=AREA+(RADIUS(I)**2)*(ANGLE(1-1)-ANGLE(I))/2 RETURN END THIS SUBROUTINE WILL GIVE AN INITIAL SHAPE OF A CIRCULAR CONVERGING SHOCK FRONT WITH AN ABRUPT CHANGE IN THE FRONT AT THE PERTURBATION. SUBROUTINE ABRUPT(X,Y,T,SM,N) IMPLICIT REAL*8 (A-H,0-Z) COMMON /LPI/PI /LP/TH /LXC/XC /LL/DELTA,L COMMON /LRA/R /LG/GAMMA /LB/XN,YN,MN,TN /LANGEL/ANGEL REAL*8 TH(100),R(100) REAL*8 XN(100),YN(100),MN(100),TN(100) REAL*8 X(N),Y(N),T(N),SM(N) INTEGER FMT(1)/'*'/,YES/'Y'/,NO/' N ' /, FORM WRITE(7,1) 1 FORMAT('0','ENTER L,DELTA, AND', . ' INITIAL STRENGTH OF SHOCK') THE AZIMUTHAL NUMBER 'L', THE PERTURBATION AMPLITUDE 'DELTA' AND THE INTIAL MACH STRENGTH IS READ. READ (5,FMT) L,DELTA,SMI THE SEMICIRCLE IS DIVIDED INTO N MESH POINTS AND L # REPRESENTS HOW MANY SUCH PERTURBATIONS CAN BE FITTED IN A COMPLETE CIRCLE. •MM* IS THE HIGHEST NUMBER OF THE MESH POINTS BEING AFFECTED BY THE PERTURBATION. MM=N*2/L ANGEL=PI 'TH' AND 'R' ARE THE POLAR COORDINATES OF A CIRCLE WHICH IS TO BE MODIFIED. DO 30 1=1,N TH(I)=PI*(N-I)/(N-1) R(I)=1. SM(I)=SMI 30 CONTINUE THE CENTER OF THE POLAR COORDINATE IS LOCATED. 'XC IS THE CENTER OF CONVERGENCE. XC=R(1)+DELTA 'X' AND 'Y' ARE THE CARTESIAN COORDINATES OF THE CIRCLE. 97 C D015 1 = 1,N X(I)=XC+R(I)*DCOS(TH (I)) 15 Y(I)=R(I)*DSIN(TH(I)) CALL GEOM(X,Y,TN,N) C C THE PERTURBATION IS INSERTED INTO THE CIRCULAR C WAVEFRONT. THE PERTURBATION IS IN THE SHAPE C OF AN CIRCULAR ARC. C X(1)=X(1)-DELTA SLEN= (X (MM) **2+Y (MM) **2-X (1) **2 ) / ( 2* (X (MM) -X (1) ) ) ANG=DARSIN(Y(MM)/SLEN) MMM=MM-1 C C THE 'X' 1Y' COORDINATES ARE BEING MODIFIED TO C TAKE THE PERTURBATION. C DO 10 1=2,MMM X (I)=SLEN-SLEN*DCOS(ANG*(1-1)/(MMM)) Y(I)=SLEN*DSIN(ANG*(I-1)/(MMM)) 10 CONTINUE C C ADDITION INFORMATION ABOUT THE DISTRIBUTION C OF MACH NUMBERS ARE ASKED. IF THE ANSWER 1 FORM1 C EQUALS TO 'NO', THEN THE MACH NUMBER IS UNIFORM C AND IS NOT BEING MODIFIED, OTHERWISE THEY ARE C TO BE MODIFIED BY THE FOLLOWING PROCEDURE. C WRITE(7,100) 10 0 FORMAT('0','DO YOU WANT A NON-UNIFORM STRENGTH?'/ 'ANSWER ''Y'' OR ''N''') READ (5,101) FORM 101 FORMAT(A4) IF (FORM.EQ.NO) GOTO 102 C C IF A NON-UNIFORM SHOCK STRENGTH IS WANTED, THEN C THE STRENGTH AT THE CENTER OF THE PERTURBATION C IS EVALUATED USING THE STRONG SHOCK APPROXIMATION C AND THE REST OF THE PERTURBATION IS LINEARLY C INTERPOLATED. C C CURV=SLEN GN=l+2/GAMMA+DSQRT(2*GAMMA/(GAMMA-1)) SM(1)=SMI*CURV**(-1/GN) C C MODIFICATION OF THE MESH POINT MACH NUMBERS. C DO 103 1=2,MM 103 SM(I)=SM(1) C C PRINTOUT ON UNIT '6' INDICATES SUCESSFULNESS C OF MACH NUMBER MODIFICATION. 98 C C THE INITIAL CONDITIONS ARE PRINTED IN THE C OUTPUT UNIT '6 1. WRITE(7,200) CURV,(SM(I),1=1,10) 200 FORMAT(lX f11F10.5) 102 WRITE(6,2) N,L,DELTA,SMI 2 FORMAT( 10 1 , 1 CIRCLE'/'NO. OF POINTS= 1 ,13 , 2X,'L=' ,12/ 'INITIAL PERTURBATION=*,F7.5/ 'INITIAL SHOCK STRENGTH=' ,F7 . 3) RETURN END
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Stability of converging shock wave
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Stability of converging shock wave Fong, Kenneth Sau-Kin 1978
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Stability of converging shock wave |
Creator |
Fong, Kenneth Sau-Kin |
Publisher | University of British Columbia |
Date Issued | 1978 |
Description | The stability of converging waves was analysed with theory and by experiment. A computer code was written in FORTRAN language for the calculation of wavefront positions of plane shock waves and cylindrical converging shock waves in ideal gas. The criteria for stability were defined and the limit of stability of imploding shock waves was obtained using the code. In order to test the model, an experiment was set up to generate shock waves in an electrothermal shock tube. These shock waves were made to propagate into a wedged channel, and the variations of the shapes of the waves front were measured from high speed framing camera pictures. Good agreement was found between the experimentally observed velocity and the model calculation, indicating that the computer procedure is a valid method to analyse the stability of converging waves. |
Subject |
Shock waves |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-26 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0094383 |
URI | http://hdl.handle.net/2429/21073 |
Degree |
Master of Science - MSc |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1978_A6_7 F65.pdf [ 5.69MB ]
- Metadata
- JSON: 831-1.0094383.json
- JSON-LD: 831-1.0094383-ld.json
- RDF/XML (Pretty): 831-1.0094383-rdf.xml
- RDF/JSON: 831-1.0094383-rdf.json
- Turtle: 831-1.0094383-turtle.txt
- N-Triples: 831-1.0094383-rdf-ntriples.txt
- Original Record: 831-1.0094383-source.json
- Full Text
- 831-1.0094383-fulltext.txt
- Citation
- 831-1.0094383.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0094383/manifest