STABILITY of CONVERGING SHOCK WAVE by KENNETH SAU-KIN FONG B.Sc, 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 t h i s t h e s i s as conforming to the r e q u i r e d 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 B r i t i s h 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 Physics The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 April 4. 1978 ii ABSTRACT The theory stability and by of experiment. FORTRAN l a n g u a g e shock waves and ideal gas. limit of s t a b i l i t y The code. set up These Good were and computer and the procedure converging waves from found model in waves. an of waves in were d e f i n e d and the obtained using speed shapes o f framing was shock into the a waves camera pictures. experimentally observed calculation, method electrothermal made t o p r o p a g a t e of the between t h e is a valid written in t h e m o d e l , an e x p e r i m e n t were high with positions shock shock waves was the v a r i a t i o n s measured agreement was velocity of shock waves code was converging to t e s t shock analysed of wavefront for s t a b i l i t y imploding order to generate wedged c h a n n e l , front of was A computer cylindrical criteria In tube. waves f o r the c a l c u l a t i o n plane the converging indicating to analyse the that the stability TABLE OF CONTENTS Chapter 1 I n t r o d u c t i o n 1.1 M o t i v a t i o n For S t a b i l i t y 1.2 1 Investigation 1 O u t l i n e Of T h e s i s 3 Chapter 2 D i s c u s s i o n Of P e r t u r b a t i o n And S t a b i l i t y 5 2.1 Volume P e r t u r b a t i o n 7 2.2 R a d i a l P e r t u r b a t i o n 9 2.3 10 P e r t u r b a t i o n Angle 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 S i m u l a t i o n 26 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 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 P o i n t s 36 4.6 Streamlines 36 4.7 Time I n t e r v a l 37 4.8 L i m i t a t i o n s Of The Code 38 Chapter 5 R e s u l t Of Numerical Computation 5.1 Plane Shock Propagating Into Wedge-like 5.2 Symmetrical 41 Channel C i r c u l a r Converging Shock 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 Chapter 6 Experiment ...41 45 ..47 57 iv 6.1 Shock Tube And Vacuum System 58 6.2 Power Source And T r i g g e r C i r c u i t 61 6.3 Method Of Measurement 61 6.4 R e s u l t 64 Chapter 7 Conclusion 69 Bibliography 71 Appendix A Repeated Linear Interpolation 73 Appendix B Runge-Kutta I n t e g r a t i o n Method 76 Appendix C Remark On Accuracy 79 Appendix D Code For Wavefront C a l c u l a t i o n 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 vi Figure(B.1) 77 1 CHAPTER 1 INTRODUCTION 1.1 M o t i v a t i o n For S t a b i l i t y Hydrodynamic s t a b i l i t y plays of an of imploding shock and heat fusion. cannot be Classical applied to Rayleigh-Taylor the fronts in B u l t e r / 4 / on i n f i n i t e l y question. mass transport The s t a b i l i t y strong converging Guderley/11/ s i m i l i a r i t y stability shock f r o n t s or a b l a t i o n heat f r o n t s such as those produced by l a s e r due to across waves important r o l e i n the d i s c u s s i o n of the f e a s i b i l i t y laser theory Investigation shock method c o n s i d e r s only theory o f waves using infinitesimally small p e r t u r b a t i o n s , b e s i d e s , important d i s s i p a t i v e mechanisms were This not considered. theory sufficient to determine the s t a b i l i t y by f u s i o n , s i n c e the imploding laser are n e i t h e r i n f i n i t e l y strong nor expected to be i n f i n i t e s i m a l l y It i s well perturbations. Lighthill/16/, by Huni(1968)/l/ cylindrical therefore conditions can as not required wave produced by l a s e r s the perturbations be small. known that plane shock f r o n t s are s t a b l e to This has Freeman/10/ Lapworth/14/. detonations is been derived theoretically and Chester/6/ and Experimental results by experimentally on concentric from Perry and Krantrowitz(1959)/18/, Ahlborn and and Knystantaus imploding and detonation Lee(1970)/13/ show that f r o n t s are s t a b l e to f i n i t e 2 sized perturbations. observed sized as with an In image-converter perturbations i n the time p r o g r e s s e d . many a n g u l a r front as a i t c o l l a p s e d towards behave l i k e collapse, 2) . What detonation strong this 1) . Whether is leads the extent This that stability as wavefronts will calculate the or to the the stable; words and centers of no longer method d i f f e r s fundamemtally i n which evolution of of a to such an restored. questions is conditions from A fast, various described. The and channel is able and to plane cylindrical perturbations. scale for discussed, Mach numbers o f standard both of language small be positions different size defined. the a wedged from is is p o s i t i o n s and waves w i t h what above two collapse into and wavefront fusion FORTRAN wavefront with cylindrical the The stability in a front a p e r t u r b a t i o n f o r which shape can fronts coded diminish strongly overdriven entire calculate shock waves p r o p a g a t i n g the of i n other t o answer of front smooth experimentally. the shock finite i t s center. size its original procedure is that questions: modify meaning as camera beginning perfectly a p p l i c a b l e to l a s e r procedure converging i n the shock waves a r e lost, and the numerical to two thesis tries numerically Huni detonation waves a r e maximum i s not perturbation and shock waves when c l o s e t o converging stability which became Ahlborn framing imploding What was regions Converging particular The stability analysis, perturbations are 3 predicted. Here the e v o l u t i o n of l a r g e p e r t u r b a t i o n s are calculated. The c o n d i t i o n s f o r s t a b i l i t y series then obtained of such c a l c u l a t i o n s of v a r i o u s c y l i n d r i c a l with d i f f e r e n t from a wavefronts perturbations. An experiment of a plane shock wave wedged channel i s then performed result is obtained propagating into a and compared to the numerical by the numerical method to t e s t the v a l i d i t y of the code. 1.2 O u t l i n e Of T h e s i s This thesis consists investigation of two main parts: theorectical of the s t a b i l i t y problem with numerical r e s u l t s and an experiment to t e s t p a r t s of the numerical r e s u l t . Chapter 2 c o n t a i n s the c r i t e r i a of hydrodynamic as a p p l i c a b l e to l a s e r stability fusion. In Chapter 3 the equation of motion of d i f f e r e n t types of wavefronts i s d e r i v e d using the method developed by Whitham. In Chapter 4 wavefronts are these then calculations yield equations solved the development a l s o l i m i t s f o r the s t a b i l i t y The experiment with for the calculation numerical method. of The of perturbed wavefronts and i s given i n t h i s chapter. used to t e s t the accuracy of the computer simulation i s described Chapter further 6 works. then i n Chapter draws 5. the c o n c l u s i o n w i t h s u g g e s t i o n f o r 5 CHAPTER 2 DISCUSSION OF PERTURBATION AND Most p e r t u r b a t i o n s on converging dampened out shock waves w i l l e n t i r e l y as shown by B u l t e r / 4 / and They concluded higher STABILITY that infinitesimally small harmonic modes i n a strong converging grow i n s t e a d of decay. symmetry of due spread be Whitham/23/. disturbances of shock f r o n t w i l l that the original not be r e s t o r e d but instead out and modify the e n t i r e front. 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 m o d i f i c a t i o n tends to cause f r o n t to achieve time. implies the shock f r o n t w i l l the p e r t u r b a t i o n w i l l However This not some s o r t of d i f f e r e n t symmetry This stabilizing mechanism at arises a the later from the dependence of the f r o n t v e l o c i t y on the changes of area as front progresses calculated in waves. a The function of a method This developed of the dependence The front the l o c a l s t r e n g t h of a wave with i n c r e a s e s more r a p i d l y due velocity of decreasing width. to pressure Consider center of a a radius build up into converging then at t h i s p e r t u r b a t i o n , where the wave i s concave, the f r o n t v e l o c i t y wavefront due is flowing p e r t u r b a t i o n which causes a wavefront to l a g behind, the be curvature a smaller the f r o n t from the compression of the gas channel can by Whitham/21/ for shock the s t r e n g t h of the wave and curvature behind time. rate of change of the wave. of in the i n c r e a s e s more r a p i d l y than the r e s t of to i t s higher curvature. This d i f f e r e n c e the in 6 front velocity increase w i l l which l a g behind. slowed down. i s reduced and along speed up the p a r t s of the f r o n t s S i m i l a r l y the p a r t s which lead The d i f f e r e n c e in curvatures along be the f r o n t at the same time the p e r t u r b a t i o n i s spread out the f r o n t . Stability i n the usual sense means that a p e r t u r b a t i o n i n a wavefront w i l l be damped out and its will original symmetry. shocks are unstable the wavefront w i l l T h i s i m p l i e s that most according restore converging to the above d e f i n i t i o n . many purposes, i n p a r t i c u l a r where one expects the For existence of f i n i t e s i z e d p e r t u r b a t i o n s as f o r i n s t a n c e i n l a s e r f u s i o n , the definition restrictive, necessary of stability because a symmetrical for high density d e s i r a b l e that the converging spherical at radii compressed mass. parts radii is obviously smaller heating and than R continue sphere wavefront greater have is than unstable. is not to to decrease as to not stays more or less the r a d i u s R of the final into several Furthermore, s t a b i l i t y required since the to "absolute" functions of R. One may increase for r<R. geometries of the p e r t u r b a t i o n We by at necessary when the can thus for p e r t u r b a t i o n s which R and "partial" p e r t u r b a t i o n s which decrease as f u n c t i o n s of r a d i u s r f o r but so absolutely have a l r e a d y been achieved be be However, i t i s A wavefront that breaks up compression stability not compression. r a d i u s of the shock f r o n t i s equal classify does can also define for r>R, the 7 1/. I t s volume 2/. I t s maximum p e r t u r b a t i o n 3/. I t s d e v i a t i o n from 2.1 Volume We AV.(area AA i n 2-dimension) amplitude similarity Ar. t o the original shape. Perturbation define the dimensionless extent of parameter the PERT volume , perturbation which is the as ratio the of the v difference i n volume unperturbed reference unperturbated Absolute before volume r e f by r e where zero. R. Notice that for laser shock does pellet i s not However, p a r t i a l compressed obtain term the in area is Eq(2.1) the an at a r a d i u s For two should perturbation volume be a shock dimensional replaced equation. zero stability converging the infinitesimal the a l a r g e r t h a n some comparable for t o be approaches v purposes, to necessary fusion. PERT absolutely stable since compressed stability radius scheme i n l a s e r volume be PERT for radius fusion t o be required Partial PERT^ d e c r e a s e s to by i.e. would that have enclosed an (2.1) only not f r o n t and r e f time, reaches perturbed volume £, stability of the the front, V AV/V function v requires = v enclosed f r o n t to reference PERT decreasing AV to the fuel size. final compression geometry, with This an area the A to definition 8 will be used later in Chapter 5 for c y l i n d r i c a l wavefront calculations. - A X I S OF — SYMMETRY V i \ PERTURBED AREA AA r : L O C A L RADIUS OF CURVATURE Fig. (2.1) P e r t u r b a t i o n Notice only the that parameters i n t h e above d e f i n i t i o n t h e volume o f t h e p e r t u r b a t i o n of the perturbation i s considered. g e o m e t r y o f t h e p e r t u r b a t i o n nor t h e d i s t r i b u t i o n Neither of front 9 velocity along the consideration. which With has s p l i t also be c a l l e d of volume guarantee would reference y alone be been alone, taken a into perturbation decreases. The is insufficient therefore greatly requirement fusion since impeded i f the to energy the wavefront was parts. Perturbation supplement ratio i f PERT has i n t o two n o n - c o h e r e n t p a r t s c a n t o b r e a k up i n t o d i f f e r e n t dimensionless the definition the wavefront stable front a smooth c o l l a p s e i n l a s e r 2.2 R a d i a l To this stability absorption permitted perturbed t h e above parameter "radial o f t h e maximum unperturbed PERT "Radial r definition, define perturbation"(PERT^), perturbation f r o n t with we amplitude Ar namely to the radius R = Ar/R stability" another (2.2) then r e q u i r e s PERT t o be a decreasing r function of R, and i f PERT goes t o z e r o before R does, then r we c l a s s i f y radial i t as " a b s o l u t e stability" i f PERT radial stability", decreases for radius and "partial r greater than r some R. For area the laser stability radius of f u s i o n p u r p o s e , one would and p a r t i a l ignition. radial r e q u i r e both stability, at least partial up to 10 2.3 P e r t u r b a t i o n So far dimensionless not Angle we give a complete we only restrict a general parameter an many more front Such which 1, w h i c h parameters the PERT^, may be do of the required. only parameter this one is extra Y; the spanned c a n be used by to d e f i n e of times circle. give a general t h e Mach number density can be front a Thus collapes o f the should But the d i s t r i b u t i o n shock f r o n t i s n o t apparent derivations. idea which vary is to the centre. front, For which field behind the o f t h e mach number considered the from a l s o be s p e c i f i e d from a snap shot not as t h e o f t h e shape o f and t e m p e r a t u r e o f t h e f l o w obtained, information y, PERT^, and completeness, this a i n t o a complete as front. still d e s c r i p t i o n o f t h e shape i n d i c a t e s t h e number perturbation front two 1 = 180°/y collapses, will pressure, by shape i s one h a l f o f t h e a n g l e c a n be f i t t e d thus the i s required, Furthermore, number perturbation of parameters o f t h e shape perturbation. The perturbation t o t w o - d i m e n s i o n a l g e o m e t r y , and i f needed. angle, azimuthal description ourselves idea is perturbation the a For a comprehensive the p e r t u r b a t i o n , If described p a r a m e t e r s PERT^ and PERT^, b u t t h e s e perturbation. of have along along the the the o f t h e f r o n t , and in the following 11 CHAPTER 3 EQUATION OF MOTION OF WAVEFRONT T h i s chapter d e s c r i b e s the c o n s e r v a t i o n equation across a wavefront and a l s o the d i f f e r e n t i a l c o n s e r v a t i o n equations f o r continuous equations gas flow. The of mass, momentum and energy behind are a s e t of h y p e r b o l i c p a r t i a l combined differential and rewritten c o n d i t i o n s across a conservation front by independently. ordinary are obtained These Chisnell/7/, are from f r o n t and the area equation, the A equation of complicated the one of are The jump the integral combined to o b t a i n a This Chester/6/ By t h i s method differential equations, i n c h a r a c t e r i s t i c form. equations. developed the f r o n t , which differential equation of motion d e s c r i b i n g the f r o n t . first conservation method was and Whitham/21/ obtains a non-linear the Mach number M o f the front. flow behind In deriving the f r o n t this i s neglected through the assumption that the jump c o n d i t i o n s can be a p p l i e d to the C+ c h a r a c t e r i s t i c o f the flow, and that the area A i s a f u n c t i o n of the l o c a l Mach number M of assumptions would not be v a l i d the propagation the front. i n e x p l o s i o n problems i n which towards the region of d i s t u r b a n c e must prime consideration. extent i n s p h e r i c a l and c y l i n d r i c a l converging indicated by Bulter's These However, comparison i t works be of to a very good shock waves as of the r e s u l t 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 d e s c r i b e s the area A 12 as a f u n c t i o n of the Mach number M , y e t i t i s the change in area that causes a corresponding change i n 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/, model Chisnell for calculating area. the front /21/ developed A waves i n g e n e r a l . a quasi-one similar a v e l o c i t y as a f u n c t i o n of Lee & Lee /6/ a p p l i e d t h i s method f o r the detonation. for /6/ and Whitham procedure The d i f f e r e n t i a l case of can be a p p l i e d to heat conservation equation dimensional gas flow with no heat source or sink i s : Momentum: Energy: |H. + |E + 11£. = 3t or p 9r U -|| u||+ a2(§| • uf|) = 0 where P,u,p,a i s the d e n s i t y , gas v e l o c i t y pressure and l o c a l (3.2) 0 i n the sound v e l o c i t y of the gas. (3.3) lab frame, 13 Manipulating these equations a one gets 2 x the (3.1) + according to the expression: (3.3) tap characteristic (3.3) (3.2) x forms of the conservation equation: d £ du dt * ""dt ± where the a +ve + a 2 pu sign 3A. A 3r = (3.4) Q applies to the positive applies to the negative c h a r a c t e r i s t i c : du/dt=u+a and the -ve sign c h a r a c t e r i s t i c : du/dt=u-a changing v a r i a b l e s along the 1 3A A 3r and substituting 1 3A 3t A 3t 3r (3.5) d i f f e r e n t i a l equations into Space-time characteristics has a s e t of ordinary a p u dA ^J+a A (3.6) _ = 0 are to be a p p l i e d to regions which do diagram is (3.5) characteristics: 2 + 1 u±a (3.4), one have any d i s c o n t i n u i t y , i . e . A 13A A 3t along the . , , dp ± apdu These equations characteristics given of in i n regions behind the wave F i g . (3.1). and not the wavefront. the The C+, C- schematic 14 diagram of is in Fig.(3.2). shown converging flow the p h y s i c a l front immediately Whitham / 2 1 / behind the In the next density will of a converging has t h e C+ c h a r a c t e r i s t i c approximation. and situation be e x p r e s s e d as a wave that for shown c a n be a p p l i e d front section shock with very the p r e s s u r e , function to the of good velocity the front velocity. 3.1.2 Equations The Of M o t i o n integral Across conservation equations Mass: P Momentum: Energy: Where measured per measured gives and the + p l v l = P V x 2 = and specific 2 + P are 2 across the f r o n t a r e : (3.7) 2 2 2 (3.8) V (3.9) + the d e n s i t y , the gas v e l o c i t y frame, the p r e s s u r e , W . h P = ^ 2 p,v,p, i n Watt/m l v l hj + % v + mass. 2 refer resp. p i n the f r o n t unit The F r o n t is For energy the input chemical energy the o f and enthalpy at the f r o n t , reactions, i n p u t Q=W/pv. t o t h e c o n d i t i o n s ahead and one often The s u b s c r i p t s behind the 1 front Fig. (3.1) Space-time diagram of converging shock 16 PRESSURE DENSITY VELOCITY Fig. (3.2) Schematic diagram of converging The general rearranged solution of equations shock wave (3.7) , (3.8), (3.9) can be as /!/: _ 2 = J U l 2 v 1 + p 1 «?i n (3.10a) 3 £i = 1 - nB Pl (3.10b) n =2 " l l g +1 g g (3.10c) M 2 J * = 1 ± t 1 2g^ + G - ( g 2-l) W (8 -l)(g -g Mj)^P v h 1 2 1 1 1 1 17 2(g +DCg -g ) fg 2 e = 1 2 M 1 (3.10e) (g -D(g -g Mp^ 1 2 1 (3.10f) 2 a ^1 l iP7 =g (3.10g) In general where e i s always small. root in (3.10d) has the p o s i t i v e shocks and W=0 for Supersonic heat waves f o r M ]<1. 3.1.3 Area The C+ waves rule characteristic terms negative are found together sign The square with M ^>1 a n d heat waves. for f o r M ^>1 a n d s u b s o n i c heat Rule area front. the sign 0<3<2. i s obtained equation b y t h e jump The s u b s t i t u t e d o f t h e Mach number This c a n be e a s i l y done gas ahead case Mi=M of M i s then of front the expressed and the initially. of the front, i nthe across of reference i s stationary speed p, u and p conditions equation i n t h e frame the front i s t h e mach by r e p l a c i n g the only i n area A. i n which the In this 18 and t h e gas v e l o c i t y b e h i n d U= V the f r o n t i s r 2 - i ^ V = (3.11a) P P ii \I^MJ+n8/ 7?+,R a (3.11b) } a,P - L 1 2 (1-nS) (3.11c) 8 for t h e c a s e o f shock , 8=2 2na rr~ 1 1 u «= 8 1 1 l l ^ 2 n " 8 - p p M . V i p (3.12a) M 8 (1 _ ( 3 ' 1 2 b ) (3.12c) 2n) 1 one g e t s Whitham's a r e a rule by s u b s t i t u t i n g (3.12) into (3.6) X (M)MdM . + ^=0 (M^-g / ) A (3.13) + 2 where X (M) s gl i ^ ) ( - 1 + i^T^ " (g -l)M +2 2 and a 2 = — £ _ _ 2g M -(g -l) i! 2 2 X (M) is a dependence of X (M) g gas. slowly g varying function of M. on M i s shown i n F i g . ( 3 . 3 ) f o r an The ideal 20 For i n f i n t e l y strong shocks i n i d e a l gas, t h i s reduces to dM , dA + — s M A _ n — (3.14) = 0 n » 1+ — • f c ^ ^ - 4.4361 for y = 5/3 s Y where » 5.0743 for y = 7/5 i n t h i s case M a A t h e r e f o r e the r u l e (3.15) s gives M a r^ _ 1 / ' s^ n f o r c y l i n d r i c a l shock ( 2 /n ) M a r these r e l a t i o n s solution may be compared the Guderley (3.16b) similarity yields CJ ^dA (M^-g / ) A X ( M ) M 2 and to shock the case of a Chapman-Jouguet d e t o n a t i o n , where 3=1 area r u l e where f o r spherical where Mar.exp(l/n). For the s (3.16a) X d M . gl CJ i t s graph i s shown i n F i g . ( 3 . 4 ) . (3.17) 21 22 For infintiy strong detonation n where n_ T dM . dA -A Cj1l + wave one finds A = (3.18) 0 = 3(1+ -) = 4.80 for Y = 5/3 - 5.14 for Y = 7/5 which g i v e s M a r^ ^ CJ^ for cylindrical CJ detonation (3.19a) M a r*~ C.P for spherical CJ detonation (3.19b) -1 n 2/n Notice the shock f r o n t and It close f o r Chapman-Jouguet detonation should strengthening be is the i n c r e a s e i n converges. converging initial The without for only moderate. number is is front. weak But quite fronts, the f o r f r o n t s with rapid as the M>5, area a consequence of t h i s f a c t that strong shows shock waves the as strengthening a function shock s t r e n g t h s , i n an i d e a l gas result be that for shock wave can be s t a b l e i n some aspect. i s obtained n o t i c e t h a t the can noted Mach It F i g . (3.5) converging resemblance between the equations applied by n u m e r i c a l l y i n f i n i t e l y strong of area of Y =7/5 cylindrical for various and 5/3. i n t e g r a t i n g Eq.(3.13), shock to shocks of i n i t i a l of approximation s t r e n g t h M=5 and i n c u r r i n g an e r r o r of more than 5% for the worst (M=°°) higher case Fig. (3.5) Area strengthening o f shock wave 24 in the final strengthening s t r e n g t h of the shock, of CJ d e t o n a t i o n waves. Fig.(3.6) shows the 25 26 CHAPTER NUMERICAL 4.1 I n t r o d u c t i o n And G e n e r a l Equations(3.13) the Mach number and system cylindrical converging of along of equation each of mesh the local can and be the t h e n be t r e a t e d however front of solid uniformly the mesh represent across the points are The of waves is flow as front positions. actual particle that the The r a y s c a n tubes which no the problem of d e s c r i b i n g reduced element c o n f i n e d The l o c a l at time introduced the f r o n t r e q u i r e walls a i n space numerically are successive I n t h i s way of evaluated considers i n normal d i r e c t i o n . of a wavefront walls. Rays o f these as t h e s o l i d convergence be a The numbers o f t h e f r o n t . plotted. can c r o s s . propagation with can on the f r o n t p o s i t i o n s a t s u c c e s s i v e p a t h s as t h e jump c o n d i t i o n s to difficult. w a v e f r o n t method Mach trajectories leave respect perturbation very motion Rays c l o s e t o t h e f r o n t would the is wavefront The point orthogonal particles single o f motion o f the f r o n t i s i n t e g r a t e d intervals flow a wave with An a n a l y t i c s o l u t i o n t o an p o i n t s which are d i s t r i b u t e d front. functions as This numerical mesh the Considerations (3.17) a r e n o n l i n e a r such the numerically. set SIMULATION M o f the f r o n t . asymmetric propagation 4 to decribing within Mach number flow o f each the tubes r a y tube 27 i s however a f u n c t i o n of the ray-tube given by Eg.(3.13) and (3.17). Using l o c a l Mach number at l a t e r times successive the method positions i s similar of located can t h i s concept, be calculated the new and f r o n t can be p l o t t e d . along the are not used and the This to the method o f c h a r a c t e r i s t i c s except the c h a r a c t e r i s t i c equations points area and the r e l a t i o n i s only that those wavefront are c o n s i d e r e d . The domain of dependence f o r the i n t e g r a t i o n i s determined by the velocity the of propagation of the disturbance along wavefront. Furthermore, i n c o n t r a s t to the usual moving method, fixing the their mesh can angular separations be to a prescribed center. to vary and new r a d i i are to be obtained interpolation. then point p o i n t i n the code are p a r t l y p r e s c r i b e d by The r a d i i are allowed by mesh P h y s i c a l q u a n t i t i e s at the mesh p o i n t s evaluated by using Eq.(3.11a), (3.11b) and (3.11c). The advantage of t h i s semi-moving mesh p o i n t method i s two-fold: (i) The p o s s i b i l i t y o f overcrowding of the mesh p o i n t s due too l a r g e a d i s c r e t i z a t i o n ( i i ) The grid size to i n time step i s removed. i s determined mainly by the r a d i u s , thus simplifying error control. However, Firstly there are also several i t needs more i n t e r p o l a t i o n s than disadvantages. i f the mesh p o i n t s 28 were to be determined by the t r a j e c t o r i e s calculated points, and readily available. disadvantages, secondly the of the streamlines previously are not But the advantages outweight a l l these e s p e c i a l l y i n the case where a shock-shock 1 is p r e s e n t , t r a d i n g o f f a c e r t a i n degree of accuracy of the f i x e d flow tube approach with the s t a b i l i t y of the v a r y i n g flow tube approach. Briefly, one c a l c u l a t i o n step c o n s i s t s of the f o l l o w i n g calculations: 1. E v a l u a t i o n o f new p o s i t i o n s o f mesh p o i n t s by the wavefront method. 2. C a l c u l a t i o n o f Mach number at new mesh p o i n t s by Whitham's area r u l e . 3. Choosing new mesh p o i n t s . 4. C a l c u l a t i o n of normals and Mach numbers at the new mesh points. A flow chart The the following process integration boundary of the o p e r a t i o n of of points i s shown i n F i g . ( 4 . 1 ) . s i x sections give a d e t a i l discretization, the and time wavefront. r e s p e c t i v e l y , while the choice integration Special streamlines of are given time A\ s;hock-shock i s an abrupt change of t r a v e l l i n g sideways across a shock f r o n t . 1 description of and area treatments of i n 4.5 and 4.6 intervals physical used is variables 29 f START j READ SHAPE READ I.A , M • CIRCLE r ^HAT .SHAPE •WEDGE- initialize X.Y.M.T initialize X.Y.M.T (total number of t i m e steps)] J= J• 1 j= 1 XN = X*M adt YN=Y* MyOdt WHITHAM AREA RULE _Y£5_ C STOP MN :X,XN,Y,YN,M J: NEW NORMALS T : X.Y ) / PLOT 7 YES / REPEATED LINEAR INTERPOLATION X : XN.YN Y : XN.YN M : XN, Y N,MN \ N O MOD(J,20)^ 0 YES NO STOP Fig. (4.1) Flow c h a r t o f numerical procedure 30 discussed i n 4.7. S p e c i a l treatments of 4.8 then d i s c u s s e s the l i m i t a t i o n s of the code to p h y s i c a l problems. 4.2 G r i d S i z e C o n t r o l In order to make the code a p p l i c a b l e to a wide v a r i e t y of problems, dimensionless parameters are used wherever p o s s i b l e . For a s e t of mesh p o i n t s on the converging ( r£ , 6 k ) for k = 1 and where where k and n are f r o n t g i v e n by K : (4.1) 6,. = kA6 non-negative i n t e g e r s and where i e i s a constant.(Fig.(4.2)) k i s the l a b e l of a mesh p o i n t and n i s the l a b e l of step. dividing cartesian time The r" are obtained by the a c t u a l p h y s i c a l dimensions of the wavefront by a normalizing let a f a c t o r K such that r j i s one u n i t of l e n g t h . t n + 1 - t n =At n coordinate and l e t < > j£ denote < M r £ , V V ' I n the p o s i t i o n of a mesh p o i n t i s d e f i n e d by R " , where R£ = ( x£ , Y £ ) (4.2a) x£ - r£cosG (4.2b) k (4.2c) Fig. (4.2) Coordinate system for Grid Size Control 32 between mesh p o i n t s R £ ^, Assuming the wavefront and t o D e segment of a c i r c l e with the three p o i n t s on a i t s circumference. wavefront at and R£ The d i r e c t i o n is then the of propagation o f the normal N £ of the c o n s t r u c t e d circle.(Fig.(4.2)) At the next time step t=t ,, n+1 r described initially *r - *r 'yr 1 ( 1 R£ at would the be front at which R^ + , was where 1> X^ *k = x£ + M|W cosN£ 1 (4.3a) n +1 = Y k + M k a A t n 8 l n N k (4.3b) where M " i s the Mach number at the mesh p o i n t R " a i s the dimensionless speed the actual speed of of sound, sound obtained by by the p r o p o r t i o n a l i t y dividing constant K Transforming back i n t o p o l a r c o o r d i n a t e s r f i f 1 1 =V(xfV = T a r f V f + ttfV Vxf > ; X i ( 4 ( 4 - * 4 a ) 4 b ) 33 the new evaluated The mesh points by repeated at the prescribed angles are then linear interpolation new d i r e c t i o n of propagation at each new mesh p o i n t i s the normal of the wavefront at the mesh p o i n t . The normals are approximated by as segments o f arcs considering Three neighbouring which passes of circles the wavefront connecting the composed of mesh points. mesh p o i n t s can d e s c r i b e an a r c of a c i r c l e through the three p o i n t s . The normal of the wavefront at the middle mesh p o i n t i s then approximated by the normal of the c i r c l e thus formed. one o f the neighbouring At the boundary mesh p o i n t s i s absent. points, T h i s has to be t r e a t e d d i f f e r e n t l y as o u t l i n e d i n the f o l l o w i n g s e c t i o n . 4.3 Changes Of Mach Number The trajectory of R£ to is the p a r t i c l e s which are c l o s e to the shock f r o n t . the boundaries of a ray-tube through which not and cross. the t r a j e c t o r y of This defines particles do The area of the ray-tube between mesh p o i n t s k k+1 f o r 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 i n area of the ray-tube from t A A ^ ^ + 1 to t - A ^ n + 1 i s then (4.7) Due to the changes i n area, the gas i n the tube behind the f r o n t w i l l experience a pressure change which i n turn causes a change i n the Mach number o f the f r o n t . the The Mach number at center o f the tube i s taken to be the average of the Mach numbers at the mesh p o i n t s M£ - (M£ + M£)/2 +1 (4.8) The Mach number at the middle of the ray tube at the next time step can be obtained by i n t e g r a t i n g M Eq.(3.13). For t h i s i n t e g r a t i o n i n t e g r a t i o n method i s used. numerically using a second order Runge-Kutta The Mach number at the new mesh p o i n t i s then evaluated by l i n e a r interpolation. 1 A d e t a i l e d d e s c r i p t i o n of the method i s g i v e n i n Appendix B. 35 4.4 Initial The start necessary of 1) . The the 2) . The of going later the a Two of the could can i n the shape o f the a one perturbed by shock giving from region initial uniform Two shock front calculation is second strengths the initial maximum shape of can be size of displacement unperturbed number o f onto a c y l i n d e r . assumes the M^.' front. initial sine R^. A perturbation to s p e c i f y the is the tried. The front. initial fitted are This results. used wavefront If at the is giving i s a plane channel. number, namely t h e be i s taken numbers first s h a p e s were u s e d . wavefront. wavefront Mach wavefront azimuthal unperturbed be front cylindrical are by t o be Ar i s the perturbed perturbation the the which shock The different onto points, front, specified semicircular perturbation. perturbation initiate of used. parameters i s the to detailed information to experimental the 1 of wedged-shape s i z e s and superimposed required mesh the various imploding various the of f r o n t , the and compared an conditions wave f r o n t , s p e c i f i e d strength shapes are into the absence strength, general is shape o f the strength initial p o s i t i o n of initial In —- c a l c u l a t i o n are initial initial in Conditions wavefront. times such Two general I f a smooth d e v i a t i o n a n t i c i p a t e d , the curve an i s replaced from perturbation superimposed abrupt a on the deviation, the by a circular 36 arc of an 4.5 Boundary At the appropriate Points the boundary neigbouring Eq.(4.6) for to treating change the can be makes o f new the it its the can as impossible Mach numbers and case direction either may as p a r a l l e l then be the to the calculated of solid be, by a b s e n c e o f one force of to use the new the shock propagation. By w a l l s or direction as of the axis propagation boundary. The c o n s i d e r i n g the new symmetry problem. Streamlines As the mesh the w a v e f r o n t addition can in determined Mach number 4.6 calculation boundaries symmetry, as the points or N, A l s o b o u n d a r y c o n d i t i o n s may of of p o i n t s k=l mesh the directions. front radius. be of traced are a p o i n t s are not readily partially available. fixed, streamlines However, " 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 i n the program. by of the streamline 37 4.7 Time I n t e r v a l The very choice o f the time important truncation will for the interval stability, reduce points the error of the d i f f e r e n c e equations. time interval will Too l o n g program time and a time i f i s only the In neighbours. the over each o t h e r . before such separation a numbers which constant time where order the is interval this will the separations o f magnitude interval compensate also c a n be used the t h e domain o f time. the r e d u c t i o n will than to cross one can depends and wedge remain But f o r a its oscillate. use on t h e the computation. the between mesh p o i n t s f o r a long In t o two o f i t s n e a r e s t points in a unacceptably unacceptable: therefore neighbouring has t o be changed for time interval during the short t h e mesh p o i n t s g e o m e t r y , where t h e wave i s c o n v e r g i n g time too i s too long, cause interval discretizing the s o l u t i o n occurs varying interval time includes neighbours other condition are and However, The maximum t i m e between discretization and a l s o t h e t r u n c a t i o n coupled Therefore extreme case time A short interval interval d e p e n d e n c e o f a mesh p o i n t nearest is make t h e c o m p u t a t i o n a l a mesh p o i n t neighbours, step the e r r o r i n c u r r e d i n mesh the and each the e r r o r o f the c a l c u l a t i o n . obviously long. for Mach A geometry, i n t h e same cylindrical towards the c e n t e r , the i n every time s t e p i n separations i n order to between t h e mesh points. In conjunction with the use of dimensionless length 38 s c a l e , K, the dimensionless time paratmeter D i s d e f i n e d D=aAt (4.9) where At i s the r e a l time and a in 4.2. R k a n d interval i s the dimensionless We R speed of sound d e f i n e d A R £ to be the s e p a r a t i o n between define as previously mesh points k 1 + k A R It = | R was k" R k l" ( 4 + found that in order - > 1 0 for the c a l c u l a t i o n to be s t a b l e , the f o l l o w i n g i n e q u a l i t y must be obeyed: —-> 4.8 20 L i m i t a t i o n s Of The The code uses the Whitham's area rule wavefront. achieved by T h i s time saving flexibility to calculate reduction the in together with propagation computer the of the time is only the p h y s i c a l q u a n t i t i e s at calculating of the code. to method the entire flow Besides numerical subjecting i t s e l f methods, the code the field. i s however not without some s a c r i f i c e s usual l i m i t a t i o n s of and wavefront considering than (4.11) Code Tremendous wavefront rather subjected for a l l n,k is i n the to the also the l i m i t a t i o n s imposed on the wavefront method the Whitham's area rule. Numerical s i m u l a t i o n cannot g i v e an o v e r a l l p i c t u r e 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 particular q u a n t i t y develops as a f u n c t i o n system. Also, since all numerical results applicable to a p a r t i c u l a r criteria often can and be the initial to simulation shape of the at the by only are but different this point, and asymmetrical s i m i l a r i t y rules. to the front. considered and front. for exemplifies perturbation the Similarity individually code a the are condition. Whitham's area r u l e i s obtained by terminates C+ The a p p l y i n g the characteristics flow f i e l d i s assumed to have which behind the no effect jump front on the T h i s c r i t e r i o n i s true with a converging shock wave where the changes convergence as the it be can also flow f i e l d front. at But the applied the f r o n t has To mesh down the a certain extent, effect on the the s a i d to a d i v e r g e n t shock appreciably shock decays point to area steady shocks where only a s l i g h t same cannot be shock slows semi-moving problems i f the are p r i m a r i l y due to plane and the p r e s s u r e d r i v i n g The front f r o n t propagates. behind the wave where the since the initialized, s t r e n g t h d i s t r i b u t i o n of the wavefront c o n d i t i o n s across a f r o n t i s not be r e s t of a wider c l a s s of problems, calculated cannot be g e n e r a l i z e d The from Our initial must initial applied conditions. because the also obtained r e s u l t s must be initial data of the how method as it expands rapidly. a l s o p r e s e n t s some p e r t u r b a t i o n s tend to separate from p a r t s of the wavefronts to form separate f r o n t s . the main The code, 40 as the mesh p o i n t s a r e consideration Such in to a the off code and terminate the the since for front will easy be both task information front field behind the front at by any an of is plotted last i s not shock time the the wavefront multiple abnormality in the the break seriously method itself. calculated, flow behind quantities it the not gas only the using or is and o f the calculation. shock code termination and non-uniform, instant curves. more the wavefront front into plot. however the p h y s i c a l a p p l i c a b l e to calculations. is converging to c a l c u l a t e cause in this n o n - i s e n t r o p i c and provided i s not code the thus closed However, b e f o r e wavefront limitation simple d e t e c t i o n system visible o f the a non-uniform flow behind t h e code usually the flow error calculated are by the p o l a r angles, take which are calculation. flexibility restricted Since cause last in their off condition will the processes The an only wavefronts breaking occurs, fixed Thus shock-heat 41 CHAPTER 5 RESULT OF NUMERICAL COMPUTATION A FORTRAN code of the above numerical method was and the programme was UBC A computing =5/3 The computer used calculation r e q u i r e s about time i s approximately wavefront Y compiled using the *FTN Compiler of Computing Center. typical calculations were p l o t t e d supported by the written of shock directly UBC 5 i s the IBM 300 Results waves using 370/168. time steps and sec. a in the of the the an i d e a l gas of plotting Computing Center. subroutine The output p l o t c o n s i s t s of the wavefront at d i f f e r e n t radial a f u n c t i o n of the unperturbed r a d i u s perturbation as and the graph of volume p e r t u r b a t i o n unperturbed 5.1 times, as a initial done of of a the Channel plane shock wave with d i f f e r e n t Mach number propagating i n t o a wedged i n order to t e s t the code by comparison shows calculations. the wavefront plots channel with of two typical The number of mesh p o i n t s i n both cases Mach number was 20 degrees. The wavefronts were p l o t t e d once every 20 The were experiment. 60, the i n i t i a l steps. of radius. calculations F i g . (5.1) graph function Plane Shock Propagating Into Wedge-like The the M=8.00 and the wedge angle was was time d i m e n s i o n l e s s time increments D per time step were held constant at 0.0001 f o r F i g . ( 5 . l . A ) and 0.0002 for 42 DIMENSIONLESS F i g . (5.1) Wavefront plot i n t o a wedged channel. (a) D=0.0001 (b) D=0.0002 LENGTH SCALE of a plane shock wave propagating 43 Fig.(5.1.B). of a In shock reflection, front was 1 Fig.(5.2) points the with Mach stablized reached observed and value. the upper This time however when compared mesh p o i n t s . seen wall t o be in a which wall, i s the part three-shock formed. the period maximum For the as t i p of Mach d e v e l o p s and a g a i n p r o p a g a t e s around 2% o f t h e as a much l a r g e r t h e v a l u e o f t h e Mach n u m b e r s . finally caused final along the front. i s too i n and t h i s long of the c a n be amplitude o s c i l l a t i o n At the 352th the t r a j e c t o r i e s has from interval sets 20 originating t o t h e Mach n u m b e r s and t h e s e p a r a t i o n s instability i n an t h e Mach stem stem the d i m e n s i o n l e s s time mesh D=0.0002, oscillates of o s c i l l a t i o n overshoot a second Numerical Mach t h e Mach numbers a t d i f f e r e n t After i n Fig.(5.2.B) instability stem, of the shock-shock. with the t h e upper a Mach a t t h e t i p o f t h e Mach stem fashion, wall a illustrates number steps case near the passage underdamped time each time in step, the o f t h e mesh points A Mach reflection i s differed from the regular shock 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 from 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 the wall i s a new s h o c k , c a l l e d a M a c h s t e m , w h i c h i s n o r m a l to the 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 shock-shock in shock dynamics, where there i s an abrupt change o f shock 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 the region behind t h e two fronts, originating from the shock-shock. T h i s t h i r d shock f r o n t i s c r e a t e d by disturbance originating from the w a l l t r a v e l l i n g i n t o the shocked r e g i o n . l A comparison of the r e s u l t with i n f i n i t e strong 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 ) . 2 shock theory rn _ (A) D= 0 . 0 0 0 1 CC LU CO JO o cr cn —I 0.0 NO 20.0 OF 40.0 TIME STEP 80 (XlO a (B) rn. D= 0.0002 on LU CD •a IE cr a 0.0 10.0 NO (5.2) (a) (b) OF —I— 20.0 TIME STEP CX10 30.0 r 40 ) Mach numbers at d i f f e r e n t mesh p o i n t s D=0.0001 D=0.0002 45 to c r o s s over, which was detected by the program in in a t e r m i n a t i o n of the program. The l a s t wavefront p l o t Fig.(5.1.B) shows t h i s c o n d i t i o n . exists in the mesh points, A which shock-shock also can be i d e n t i f i e d as an abrupt change i n Mach number and normal mesh points. The first observed at the 280th time travelled and r e s u l t i n g direction along the r e f l e c t i o n of the shock-shock i s step, when the shock-shock had a l l the way along the f r o n t and h i t the other w a l l . 5.2 Symmetrical C i r c u l a r Converging Shock F i g . (5.3) shows the wavefront p l o t f o r the geometry symmetrical converging M=8.00. The cylinder with initial number of mesh p o i n t s i s 90. Mach number The initial time increment D i s 0.0003, but i t i s made to vary during computation in separations between computing time From the p l o t remains the compensate mesh required for the points. to 5 complete i n F i g . (5.3) i t can be seen that cylindrical as the decreases i n seconds of 240 time s t e p s . the wavefront i t c o l l a p s e s towards the The Mach number o f the f r o n t as a f u n c t i o n of the can be used to e v a l u a t e the computational e r r o r of the code by comparing from to were perfectly center. radius order of a direct the obtained dependence with those numerical integration i n f i n i t e l y strong shock approximation. r e s u l t of such obtained of the Area Rule and the Fig(5.4) shows the comparison. Exact agreement i s obtained between the wavefront method 46 CM DIMENSIONLESS F i g . (5.3) Wavefront f r o n t with D=0.0003. and the direct shock these of two the area methods and approximation i s l e s s than 2%. i n d i c a t e that i f there i s no d i s c o n t i n u i t y i n the discretization UNIT of a c y l i n d r i c a l converging integration d i s c r e p a n c y between strong plot LENGTH of the wavefront shock rule, while the the infinitely These r e s u l t s the wavefront, and the approximation of the area change have caused n e g l i g i b l e computational An e s t i m a t i o n o f the computational e r r o r caused error. by the r e l o c a t i o n o f the mesh p o i n t s a f t e r each time step i s achieved by displacing from i t s wavefront. the natural o r i g i n o f the p o l a r c o o r d i n a t e s 0.2 u n i t s positsion, i.e. at the center of the T h i s displacement created an asymmetrical point d i s t r i b u t i o n grid i n an otherwise symmetrical geometry, hence 47 comparison with previous calculation the m a g n i t u d e o f t h e e r r o r c a u s e d points. since the The r e s u l t The perturbation onto a front circular similar front curvature shock approximation. circular circular is the wavefront. an a b r u p t Mach wavefront change number at infinitely a the center linear shock center o f the p e r t u r b a t i o n where t h e Mach number case, the superimposed along the second given case Eg.(3.16a) an of for o f the arc the infinitely case, of a part of a smaller two c a s e s there front. o f the p e r t u r b a t i o n The again follows approximation(Eq.(3.16a)), of t o t h e edge t h e Mach number of the were i n v e s t i g a t e d b a s e d different on t h e and from t h e perturbation, i s t h e same as t h e u n p e r t u r b e d f o r the three is as f u n c t i o n o f t h e In t h e t h i r d transition of s t a b i l i t y kinds b u t t h e Mach number i n the curvature is perturbations The by that error. diferent displacment Thus u n l i k e t h e f i r s t strong limit using negligible t h e Mach number case, i s replaced there The with mesh i s not p l o t t e d In t h e f i r s t but i s rather radius of the indicated three throughout. i n shape t o t h e f i r s t strong the wavefront, constant of Perturbation f r o n t s with i s a sinsuidal i s not constant, local has c a u s e d were i n v e s t i g a t e d . amplitude being 0.01% a g a i n Shock Wave W i t h e v o l u t i o n o f shock of p e r t u r b a t i o n s by r e l o c a t i o n than r e l o c a t i o n o f t h e mesh p o i n t s Converging an i n d i c a t i o n o f o f the comparison which the d i f f e r e n c e are l e s s 5.3 C i r c u l a r can g i v e front. cases following of three 48 o.s ~i 0.6 A WAVEFRONT CALCULATION • NUMERICAL INTEGRATION OF AREA RULE • INFINITELY STRONG SHOCK APPROXIMATION 1 0.7 RADIUS 1 OF 0.8 FRONT r— 0.9 1.0 F i g . (5.4) Comparison of Mach number increase between W a v e f r o n t melt'hod w i t h I n f i n i t e l y S t r o n g Shock A p p r o x i m a t i o n . 49 criteria for s t a b i l i t y : 1) . Shape s t a b i l i t y . ( W a v e f r o n t does not break 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 3) . P a r t i a l The area s t a b i l t y . ( A A / A second and d i r e c t l y from the condition the decreases.) decreases.) third wavefront up.) criteria can calculations, must be observed and i n t e r p r e t e d p l o t s because the code i s not designed to be while obtained the first from the wavefront handle broken up wavefronts. Careful inspection of time d i f f e r e n t cases of wavefronts instead that The p e r t u r b a t i o n angles were perturbation depending either revealed of in the three a l l cases, of d i s a p p e a r i n g , the p e r t u r b a t i o n s spreaded out along the f r o n t s . the evolutions amplitudes on the i n i t i a l were conditions, decreased. PERT, wavefront the PERT would r 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 o f the f i r s t case i s shown i n F i g . (5.5). 90 and and while However, A i n c r e a s e or decrease. A increased initial time The number of mesh p o i n t s interval D was 0.0003. was The p e r t u r b a t i o n had an i n i t i a l amplitude of 0.1 and substended an angle of 120 time steps were c a l c u l a t e d and the y 20 degrees. = wavefront was p l o t t e d once every 20 time s t e p s . the wavefront p l o t , i t can spreaded At out and decreased the 120th time s t e p , small that the the wavefront be seen that From the p e r t u r b a t i o n i n amplitude as time progressed. perturbation amplitude became roughly c i r c u l a r was so i n shape, 50 2.25 2.5 F i g . (5.5) Wavefront plot of a t y p i c a l p e r t u r b a t i o n of 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 w i t h c o n s t a n t shock s t r e n g t h ) indicating In that order perturbations superimposed the initial to appreciate of onto a number), the Ar changed, were by open condition for stability are a region of circular Ar traced stability and area i n the on shape of stability with region left constant was of Fig.(5.6). stability satisfies for perturbation perturbation stability in Fig.(5.6). which stable. case(sinusoidal angle triangles partial limit wavefront partial and the first perturbation of y combinations marked the p e r t u r b a t i o n was the and Together a l l three Mach amplitude found the for broken line Similarly partial the radial t h e n one criteria. finds 51 10 20 30 40 PERTURBATION ANGLE 7 50 60 IN DEGREE Fig. (5.6) 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 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) F i g . (5.7) shows the wavefront p l o t of a the second case. and the The i n i t i a l plot revealed that of at the Inspection become circular. s t a b l e because both PERT of radius. Also A the of the 160th time step a Mach stem perturbation Mach stem grew i n l e n g t h as time progressed, not of p e r t u r b a t i o n amplitude was 0.3 the p e r t u r b a t i o n angle was 20 degrees. appeared at the center did perturbation This and PERT r existence (y=0). This and the wavefront p e r t u r b a t i o n i s however decreased of the as functions Mach stem had 52 Fig. (5.7) W a v e f r o n t plot o f a t y p i c a l p e r t u r b a t i o n o f 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 w i t h non-uniform shock s t r e n g t h ) prevented the perturbation stability f o r the second from b r e a k i n g Notice second even though t h a t shape, the s t a b i l i t y that the Obviously area initial that for limits Mach they instability For they are d i f f e r e n t , in first and of due number d i s t r i b u t i o n s are plotted may have t h e same case and a m p l i t u d e - w i s e , due t o t h e h i g h e r yet is perturbations the p e r t u r b a t i o n o f the second perturbation, The l i m i t o f case o f p e r t u r b a t i o n Fig. (5.8). cases, off. slighly more to were is the fact different. more Mach number prone initial to stable at the shape b e c a u s e o f t h e same. the perturbation of the third case(circular 53 LU PERTURBATION ANGLE F i g . (5.8) L i m i t o f s t a b i l i t y 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 perturbation wavefront stability plot is i s shown F i g . (5.9) had perturbation the superimposed plotted an wavefront, part the wavefront, on therefore kind a the typical limit of perturbation in amplitude o f 0.3 and From t h e o b s e r v a t i o n seemed l i k e l y Observations this and perturbation i t strength) front), The o f the p e r t u r b a t i o n shape u n s t a b l e . parameters Fig.(5.10). o f t h e second shock Fig.(5.9) a n g l e o f 36 d e g r e e s . time a s m a l l parent non-uniform in initial IN DEGREE for perturbations onto a c i r c u l a r shown in / would that a t some broke away the p e r t u r b a t i o n on different of perturbation of later from i s termed perturbation indicates that they 54 2.5 Fig. (5.9) 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 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) the are more susceptable and second to shape i n s t a b i l i t y than the first kind. The l i m i t s of s t a b i l i t y perturbations are common region i n the different f o r the three d i f f e r e n t cases of a l l plotted figure perturbations. i n Fig.(5.11). that is Since stable the i n i t i a l There are a to the three c o n d i t i o n s of the three p e r t u r b a t i o n s cover a wide v a r i e t y of p e r t u r b a t i o n s , we can say that, the s t a b l e region in F i g . (5.11) i s stable to most p e r t u r b a t i o n s . Several conclusions can be drawn from i n s p e c t i n g the 55 UJ PERTURBATION ANGLE / IN DEGREE F i g . (5.10) L i m i t o f 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 o f 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 w i t h n o n - u n i f o r m shock s t r e n g t h ) limits of o f s t a b i l i t y f o r the three perturbations. related, and t h e y distribution Radial depend of Mach p e r t u r b a t i o n with radially and than that of perturbations since the a with than along small on the angle the i s more perturbation due t o a r e a perturbation perturbation with A unstable a small angle, convergence angle. angles initial front. for a large perturbation front small a cases s t a b i l i t y are c l o s e l y significantly number area-wise for and a r e a a large perturbation perturbation angle, strengthening very c r i t e r i a of the three the i s less However, will tend to 56 LU PERTURBATION ANGLE F i g . (5.11) Limits perturbat ions break off therefore only be earlier, unfortunate computing wavefront time is i s actually of shape that the l i m i t from judgements discontinuous stability causing determined subjective handle of the are wavefronts breaking the and before up. IN DEGREE three cases instability. of It i s f o r shape s t a b i l i t y c a n observations reguired, required / since in which t h e code furthermore some cannot too one c a n a c t u a l l y much see t h e 57 CHAPTER 6 EXPERIMENT An and experiment letting result. a plane channel with direction produced cross several in the an from validity be at formed the front betweeen t h e n u m e r i c a l numerical will plot discontinuity, be as amplified, oscillations shaped wave and is reproducibitiy is Thirdly, since a circular dividing a plane towards each shock and other(Perry/18/), will the shock bending it is sensitive result. numerical show up a the not at the or as symmetrical generate ,difficult to is a This wavefront to of channel, from difficult imploding geometry because very Secondly, especially front shock experimental this along wavefronts. relatively the a error and was wavefront. give the a rectangular converging the will shock chamber the at with converging in changes i n camera p h o t o s f o r The imploding by wedged i n Mach number o f tube the produced The result. Firstly of a converging framing e l e c t r o t h e r m a l shock will unphysically imploding into change computer reasons. simulation wavefront the instead of a c i r c u l a r comparison the constructed shock wave was A two-dimensional discontinuity any and were o b t a i n e d with shock-shock to t e s t was a rectangular cross section. section. used apparatus shock wave p r o p a g a t e shock comparison the A converging of propagation reflected At designed, measurements were p e r f o r m e d computer was was achieve. produced separated clear and by section how the 58 turbulence produced artificially 6.1 induced Shock Tube And An Fig.(6.1b). three along The of the rectangular one inch thick center p i e c e has The shock inches, a f t e r which one tube electrodes are filling i s Argon and the range of 0.5 MKS The digital mounted tube a channel of one inch to 10 t o r r . The which Two is i s measured by a the pressure pressure controlled by reproducible acetylene visible. can to 0.2%. T h i s automatic gas In be added i n order is in transducer. and measured voltmeter. pressure The pressure i s attenuated valve end. pressure capacitative displays tungsten filling the d i g i t a l presettable cross for a length of 8 which compares a p r e s e t reading with a from inch width square section. initial filling one uniform the voltmeter, in p l a t e s sandwiched comparator produces i s shown i n o p p o s i t e l y i n the other of the transducer An automatic gas the s i d e of the channel i s bent 20 degrees I n t r u . Inc. Type 221M output shock lucite is inward to from a wedge-shaped end gas affect shock tube i s c o n s t r u c t e d i t to act as the shock tube of section. will a c r o s s - s e c t i o n a l diagram of i t i s shown of together. bending perturbation. diagram The pieces the Vacuum System isometric F i g . ( 6 . 1 a ) , while by by a directly. a digital the output filling of system i n the shock tube which i s addition a small trace of to make the shock f r o n t more Fig. (6.1a) Isometric diagram o f t h e shock tube Fig. (6.1b) Cross sectional diagram o f t h e shock tube 61 6.2 Power Source And Tr igger C i r c u i t The power source line constructed inductors. voltage automatic up an unbalanced lumped transmission of 5 m i c r o f a r a d c a p a c i t o r s and 5 The of is capacitor to charging 30 KV can a preset by a s a t u r a b l e transformer and an control. be charged microhenry to A N a t i o n a l E l e c t r o n i c s high v o l t a g e s w i t c h i n g i g n i t r o n model NL-1037 i s used i n the output of the c a p a c i t o r bank. along The d i s c h a r g e c u r r e n t then a RG-8/U c o a x i a l c a b l e and d i s c h a r g e s i n the shock tube through tungsten electrodes. The t y p i c a l d i s c h a r g e c u r r e n t i s a f l a t p u l s e which ranges from 2KA to 37KA f o r of 80 microseconds. trigger capable The i g n i t r o n a duration i s t r i g g e r e d by a k r y t r o n u n i t which i s i n t u r n c o n t r o l l e d by a d e l a y of g e n e r a t i n g delays of up to 100 schematic in travels generator millisecond. 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 A i s given Fig.(6.2) 6.3 Method Of Measurement Space-time measurments of the f r o n t p o s i t i o n are obtained by an TRW a image-convertor 5-frame microsecond midway camera model 1-D fitted with p l u g - i n u n i t 26-B with delay between frames i n the region. between piece l u c i t e framing the The camera i s front focused onto a plane p l a t e and the backplate of the 3- plate. The framing camera i s t r i g g e r e d by a +300v pulse produced 62 MANUAL TRIGGER •10V DELAY GENERATOR Jl • 60V • 300V ITRW 2 6 B 5-FRAME PLUG-IN UNIT [TRW 4 6 A i DELAY GENERATOR T R W 1D IMAGE CONVERTOR CAMERA KRYTRON TRIGGER UNIT Ll •5KV PRESETTABLE HIGH VOLTAGE SUPPLY CAPACITOR BANK IGNITRON SWITCH SHOCK F i g . (6.2) circuitry Schematic diagram of the TUBE electronic control 63 from a TWR delay 46A d e l a y generator which p r o v i d e s a v a r i a b l e from 0.5 microsecond to 100 microsecond. time Exposure times of 50,100,200 or 500 nanoseconds per frame can be The used. i n t e r - f r a m e d e l a y and the exposure time are independently a d j u s t a b l e on the camera. The TWR 46A i s i n turn triggered by the d e l a y generator of the power source, so that the spacetime measurement is synchronized by the trigger of the of the d i s c h a r g e were taken using P o l a r o i d land discharge.(Fig.(6.2)) Pictures film Type 107 measurements, for the quick observations. For d i s c h a r g e s are recorded by Kodak Royal RX120 negative and developed by A c u f i n e f i l m developer a film speed of ASA The Pan giving 3000. speed of the unperturbed f r o n t and the Mach stem obtained from the p o s i t i o n s of the wavefront the accurate l a r g e s c a l e enlargements as determined was on of the framing camera n e g a t i v e s . 64 6.4 R e s u l t Various filling pressure and discharge tried and i t was d i s c o v e r e d that f i l l i n g torr gave the most visible shock v o l t a g e s were pressure of front 3 to 8 with a d i s c h a r g e v o l t a g e o f 10 KV, e q u i v a l e n t to a d i s c h a r g e c u r r e n t o f 20 KA. The framing camera photos of the d i s c h a r g e showed a plane shock f r o n t followed by a heat uniform speed of approximately section with a front propagating t y p i c a l p i c t u r e i s shown i n F i g . ( 6 . 3 ) . down the 3km/s. The i n i t i a l A filling pressure was 5 t o r r argon and the charging v o l t a g e was 10 The of d i s c h a r g e c u r r e n t was measured to be 20KA. The s e r i e s photographs shows c l e a r l y the development of the shock heat wave. Inter-frame delays are 1 microsecond exposure times are 500 nanoseconds. wave i n the uniform seen as The heat The speed of the shock The plane shock f r o n t can be f r o n t i s a b r i g h t luminuous cloud f o l l o w i n g the shock I t s boundary i s not w e l l d e f i n e d . As the shock f r o n t e n t e r s the nonuniform r e g i o n , appears, At l a t e r spread to a a t h i n luminuous region and i s very s h a r p l y d e f i n e d . front. stem and and the region i s 2.98+0.15km/s corresponding Mach number o f M=9.37+0.47. be KV. time indicating at unperturbed the front. front. intersection The Mach e x i s t e n c e of Mach r e f l e c t i o n . i n t e r v a l s , the Mach stem i s observed along the unperturbed seen the a of Mach to grow and A shock-shock can a l s o the Mach stem stem i s observed and the to have a 65 Fig. ( 6 . 3 ) Framing camera picture of shock front produced from a discharge current o f 2 0 K A in 5 torr argon. Inter-frame delay: lysec. 66 higher velocity pictures the than the unperturbed velocity is Due the M=10.77+1.00. the to because of higher observed to be b r i g h t e r than front. 3.42+0.32km/s, higher Mach the the equivalent temperature number, the In to attained Mach stem can be unperturbed front i n the photograph. The ratio of the velocity v e l o c i t y o f the unperturbed numbers were measured r e s u l t s are p l o t t e d 5% discrepancy considers compared i n Fig(6.4). This the the with experiment discrepancy can the calculation the i o n i z a t i o n at and behind shock the although than the unperturbed the shock front the and region. the front, and i f one caused by also Mach stem has a higher speed a higher speed following i n the r e g i o n where the the Mach This means perturbation the stem due to higher gas d e n s i t y achieved by the Mach stem that although there i s a s t a b i l i z i n g i t caused in the flow shock front. This m u l t i p l e step wave compression is very scheme. front, behind the f r o n t perturbes the subsonic heat f r o n t , which i s the d r i v i n g for the tube. mechanism at work to r e s t o r e the symmetry o f the shock the a wavefront reduced lowering of g shock was not perturbed than the r e g i o n where has formed The f r o n t , the subsonic heat f r o n t has Mach theory. be boundary e f f e c t of the w a l l s of the shock that initial There i s approximately in Notice Mach stem and the front for different and between calculation. of force undesirable for a However, i f one 67 F i g . (6.4) C o m p a r i s o n o f Mach number change due t o r e f l e c t i o n between infinitely strong shock approximation, wavefront c a l c u l a t i o n s and e x p e r i m e n t . 68 turns of on this result the p o s i t i o n t h e shape front the o f t h e p r e c e d i n g shock through t o the shock velocity. pressure at inhomogenuity reaches the front front, i s only head of in the front front. is slightly the heat itself leading effect the flow subsonic, are c a r r i e d T h i s means t h a t front, While i . e . the v e l o c i t y the t h e shock shock information perturbations little front. and t h e h e a t f r o n t which C+ c h a r a c t e r i s t i c s , front one can c o n c l u d e t h a t o f t h e s u b s o n i c h e a t wave have between t h e shock velocity around the from t h e h e a t characterized by h i g h e r t h a n t h e shock the inhomogenuity front will caused i n the by be d i s p e r s e d t o a u n i f o r m Mach the as i t number on 69 CHAPTER 7 CONCLUSION The wavefront calculations propagations of wavefronts are present. wavefronts initial They to the also show the numerical numbers between are v i s i b l e . and one can numerical computer The expect 9.0, performed of r e s u l t of indicates with initial obtained nelgects ionization the even effects of d i s c o n t i n u i t y i n front has an adverse better the c i r c u l a r agreement converging Hence between the f r o n t with the situation. further code of or the the range i n which the A l s o the presence similar results work, improvements f o r the wavefront can The mesh.points can be non-uniform, more region. with Furthermore, s e v e r a l c e n t e r s of convergence, be made to the 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 r e l o c a t i o n of the mesh p o i n t s . perturbation of comparison the r e s u l t of the numerical c a l c u l a t i o n . real physical For of the stability Such good agreement was the d i r e c t i o n of propagation on the i n t o a wedged channel T h i s was 8.0 boundary e f f e c t . effect idea r e s u l t and'the experimental though i n the numerical method one and that perturbation. an agreement of w i t h i n 5%. shocks good when v a r i o u s kinds of p e r t u r b a t i o n s plane shock waves propagating Mach a p e r t u r b a t i o n s depends to a l a r g e extent on c o n d i t i o n s of between gave so that s e p a r a t i o n s between emphasis the wavefront the breaking on the can have up of wavefront for unstable perturbations can be observed. 71 BIBLIOGRAPHY 1. Ahlborn,B. and J.P. Huni. 1969. S t a b i l i t y Measurements of C o n c e n t r i c Detonations. 1191-1192 and Space-Time AIAA J o u r n a l 1_, 2. Belokon,V.A. et a l . 1964. Entrance of a strong shock wave i n t o a wedge-like c a v i t y . Sov. Phys. JETP 2A, 33-40 3. Brueckner,K.A. et a l . 1974. Hydrodynamic stability l a s e r - d r i v e n plasma. Phys. F l u i d 1_, 1554-1559 of a 4. B u t l e r , D . S . 1954. Converging spherical and cylindrical shocks. Report No. 54/54 Armament Research and Development Establishment, M i n s t r y of Supply. Fort H a l s t e a d , Kent. 5. B u t l e r , D . S . 1955. Symposium on b l a s t and shock waves. Armament Research and Development Establishment, M i n s t r y of Supply. F o r t H a l s t e a d , Kent. 6. Chester,W. 1954. The quasi-cylindrical P h i l . Mag.(7)45, 1293-1301 shock tube. 7. C h i s n e l l , 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. I n t e r s c i e n c e P u b l i s h e r s , New York. 9. Dorn W.S. and D.D McCracken. 1972. Numerical Methods with F o r t r a n IV case s t u d i e s . 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 . F l u i d Mech. 6, 469-480 11. Guderley,G. 1942. Starke Kugelige und Zylindrishe Verdichtungstoesse i n der Naehe des K u g e l m i t t e l p u n k t e s 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 d e c r e a s i n g density. Sov. Phys. JETP 31, 747-749 13. Knystautas,R and J.H. Lee. Spark I g n i t i o n of Detonation Waves. AIAA J o u r n a l 5, 1209-1211 Converging 14. Lapworth,K.C. 1959. An Experimental i n v e s t i g a t i o n of the s t a b i l i t y o f plane shock waves. J . F l u i d Mech. 6_, 469480 15. Lee,H.L.and B.H.K. Lee. 1965. C y l i n d r i c a l Waves. Phys. F l u i d 8, 2148-2152 Imploding Shock 72 16. L i g h t h i l l , M . J . 1949. Proc. Roy. Soc. A198, 454-470 17. Murakami,H. e t a l . 1977. S t a b i l i t y o f D e c e l e r a t i n g Shock Wave. Trans. Japan Soc. Aero. Space S c i . 2Q_, 100-107 18. Perry,R.W.and A.Kantrowitz. 1951. The s t a b i l i t y o f converging shock waves. 878-886 production and J . Appl. Phys.22, 19. Russel,D.A. 1966. Shock Wave Strengthening convergence. J . F l u i d . Mech. 2_7, 305-314 20. Swan,G.W. and G.R. Fowles. 1975. Phys. F l u i d 18, 28-35 Shock Wave by area Stability. 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 through regions of non-uniform area J . F l u i d . Mech. 4, 337-360 23. Whitham,G.B. 1973. L i n e a r and Sons, New York. shock or waves flow. and Nonlinear Waves. John Wiley 73 APPENDIX A REPEATED LINEAR In p r i n c i p a l , n-1 degree polynomial therefore of given if one of x T h i s method on will give saving. One the repeated linear (x procedure using the value. The first two y from (x _ k 1 # of a n d one these w value order h requires 3 linear and r k and use the the value the given value for of y. large of y with three interpolations, the last (x ,y ), k points, n. x great is points e i n t e r p o l a t i o n s g i v e the k of f i ta interpolations Take e can " s t r a d d l e " the v a l u e interpolation. value can consuming second ^k+l'Yk + l* sets to f i n d approximation y _ ^ , U ,y ) k o f x, time one interpolation whose v a l u e s order given interpolated of very plane entire polynomial (x^,y^) a second ^ k - l ' ^ k - l * ' k'yk' n-1 is time two i n a x-y the to some v a l u e s however three points The through the Using n points wishes" to o b t a i n by y corresponding value INTERPOLATION k one the first using the interpolated values ( X f c + l ' V k + l ) resp. (A.2) (A.3) 74 In figure(A.l) the v a l u e s y^}, interpolated three y5^ given of points are L l and L2 the two represent k - l k + 1 the lines of both l i n e a r We functions o f x. interpolated once ^ k + l'y"k+l^ as p o i n t s . x shown by l i n e above interpolation. L3 Note t h a t more using in fig.(A.l). P i s y£ and three points ( x k 1 k 1 n d 7 The y- k a linear ^k' ^ directly (A.4) + k polynomial . fy . )# P on L3 interpolation. , + interpolated a by k-l\ y ^ j ^ i s no l o n g e r Note .that quadratic y W l~\-l Yk are (x^ I'^k^l^ The p o i n t i s given Pk+1 (2) and The r e s u l t o f t h e i n t e r p o l a t i o n i s x i s the r e s u l t o f t h i s f i n a l coordinate y£*j ' that a n d function passes ^k+l^k+l o f x, i t i s a through 1 . the Fig. (A.l) Repeated Linear Interpolation 76 APPENDIX B RUNGE-KUTTA The Runge-Kutta using 1) forward method integrated. permits steps. thepast Only integration 2) I t i s a second and backward I t does n o t r e q u i r e of I N T E G R A T I O N METHOD order I t sadvantages history thevalue integration method include of the function to be of thefunction at the point i srequired. an easy change in the step size of the integration. Suppose point x=x a differential y - y (x ,y ) m m line y this . m + 1 ,y<lj), , slope obtained, then theslope integration r L2 a t ( x m+1 step, ,,V^1) m+1 J 1 slope of thelast 2 ' through once Y^t\ integration (B step ^ s c a n be e v a l u a t e d (1) aT -^!.^ the - (B.3) J forward d at where \ = f ( x ,y ) ( x ..-x ) + y m+1 m'^m m+1 m m the B y' p a s s i n g ( 1 i s ' ' L l c a n be drawn w i t h and < x 3 y (B.l) m ' ' V df a and a s o l u t i o n i s given, m then equation - 4) i st h e average o f 77 the forward and the backward steps. (2) (B.5) dx (2) (B.6) 78 Fig. (B.l) Runge-Kutta Integration 79 APPENDIX C REMARK ON A numerical position of a approximation errors major and in simulation nonlinear D equation, the can denotes solution the wavefront best termed an solution. three The into quantities, solution from numerical the c a n be c l a s s i f i e d exact the s o l u t i o n actual be inacessible I f we c o n s i d e r X denotes N denotes calculating problem numerical categories. N. for t o an o t h e r w i s e the ACCURACY two D, X of the nonlinear a " p e r f e c t " computer and solution of the problem. Then t h e e r r o r i s D-N= (D-X) + (X-N) The setting being order e r r o r has two p a r t s , D-X up t h e n o n l i n e a r p r o b l e m w i t h used. terms difference This i s caused i n any n o n l i n e a r equations. discretization off i s the i n s t a b i l i t y o f time of intermediate The numerical by t h e equations other interval the p a r t i c u l a r part, size value due t o t h e f i n i t e and t h e f i n i t e computer. X-N, This t h a t we must use a g r i d the algorithm i n c o n v e r t i n g them fact speed by truncation of higher is and mesh p o i n t s and values. caused and time i s due interval memory due into to roundingto the of f i n i t e size of 80 APPENDIX D CODE FOR WAVEFRONT CALCULATION The following calculation. to Center, 'ALGRAF , 'ALSIZE', C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C FORTRAN and 'FTNCMD' so code for of detonation the subroutine subroutine Computing 1 the For c a l c u l a t i o n modifications MTS@ Command is are 'FUNC is fronts, minor are required. supported the shockfront plotting 'ALSCAL', 'ALDASH' AND by The the UBC subroutines 'PLOTND'. THIS PROGRAM EMPLOYES WHITHAM'S AREA RULE TO CALCULATE THE TIME EVOLUTION OF A SHOCK FRONT. IT COMBINES THE CALCULATION OF WEDGED AND CYLINDRICAL FRONTS. THE BOUNDARY CONDITIONS ARE TREATED AS AXIS OF SYMMETRY TO CALCULATE THE VALUE OF THE BOUNDARY POINTS. THE NEW MACH NUMBERS ARE CALCULATED BY RUNGE-KUTTA INTGERATION AND THE NEW MESH POINTS ARE CALCULATED BY REPEATED LINEAR INTERPOLATION. R=X,Y,M,T CARTESIAN COORDINATES, STRENGTH, AND NORMAL OF MESH POINTS. RN=XN,YN,MN,TN CARTESIAN COORDINATES, STRENGTH AND NORMAL OF MESH POINTS AFTER ONE TIME STEP. RADIUS,ANGLE POLAR COORDINATES OF MESH POINTS. PANGLE ANGLE SUBSTAINED BY SHOCK FRONT IN RADIANS. TANG ANGLE OF WAVEFRONT IN WEDGED GEOMETRY. DEL,L S I Z E AND L # OF PERTURBATION. GAMMA RATIO OF S P E C I F I C HEATS. XC CENTER OF POLAR COORDINATES. D TIME INTERVAL. J # OF TIME STEPS. II CONTAINS INFORMATION OF THE SHAPES OF THE WAVEFRONT. XLENG,YLENG X,Y VALUES OF RADIAL PERTURBATION. XAREA,YAREA X,Y VALUES PERTURBATION. OF AREA 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/ N'/,FIT,FRONT COMMON /LA/R /LB/RN /LE/D /LG/GAMMA /LPI/PI / L J / 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 1 LOGICAL UNIT 2 I S DEFAULTED TO 'REF.DATA , WHICH CONTAINS THE DATA OF THE REFERENCE CIRCLE. LOGICAL UNIT 5 I S DEFAULTED TO *SOURCE*. IT RECEIVES THE NECESSARY INPUT PARAMETERS. LOGICAL UNIT 6 IS ASSIGNED TO A SCRATCH F I L E " - F I L E " , THIS CONTAINS A L L THE PRINTED OUTPUT. LOGOCAL UNIT 7 IS DEFAULTED TO *SINK*. I T PROMPTS FOR THE NECESSARY INPUT PARAMETERS. THE PLOTFILE IS DEFAULTED TO "-PLOT#". 1 CALL FTNCMD('DEFAULT 2=REF.DATA ,18) CALL FTNCMD('ASSIGN 6=-FILE',14) CALL FTNCMD('DEFAULT 7=*SINK*',16) 1 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 READ(5,FMT) N NN=N-1 PI=3.14159265358 WRITE(7,102) N') 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 ' I I ' 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'' , OR ' N''') READ(5,881) FRONT FORMAT(A4) 1 1 1 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 IF (FRONT.EQ.NO) CALL CIRCLE(X,Y,T,M,N) (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 400 C C C C GOTO 401 CALL WEDGE(X,Y,T,M,N) THE SIZE OF THE OUTPUT PLOT IS 5" X 5" FOR WEDGED FRONT. CALL ALSIZE(5.0,5.0) CALL ALSCAL(0.0,1.0,0.0,1.0) C C C C C C C FOR WEDGED FRONT, 'FIT'='NO' SO THAT A BEST-FIT CIRCLE WILL NOT BE EVALUATED. 401 FIT=NO CONTINUE THE INITIAL SHAPE IS PLOTTED. CALL DATAPL(X,Y,N,0) WRITE(7,100) C C C C C C C C C C C C C C C C C C C C C C C C C C THE # OF TIME STEPS 'IT' IS REQUESTED. 100 FORMAT ('1','ENTER NUMBER OF TIME STEPS') READ (5,FMT) IT "GEOM" IS CALLED TO INITIALIZE THE VALUES OF T. CALL GEOM(X,Y,T,N) IF(II.EQ.1) READ(2,15)D,AREARE,RADREF FOR CYLINDRICAL WAVEFRONTS, THE AREA 'AREARE' AND THE RADIUS 'RADREF' OF THE REFERENCE CIRCLE IS READ. THE VALUES 'AREARE' AND 'RADREF' WERE COMPUTED USING THIS SAME PROGRAM WITH NO PERTURBATION AND STORED IN A STORAGE FILE. 'D', THE TIME INTERVAL IS ALSO READ. 15 2 FORMAT(1X,F10.8,2X,F10.7,2X,F10.8) DO 2 1=1,N RADIUS(I)=1. THE CALCULATION IS STARTED IN THE FOLLOWING DO-LOOP WHICH CALCULATES THE NEW CHARACTERISTIC POINTS, THE NEW SHOCK STRENGTH AT THE CHARACTERISTIC POINTS AND THEN THE NEW MESH POINTS. DO 20 J=1,IT 84 C C C C FOR CYLINDRICAL WAVEFRONTS, 'AREARE' AND 'RADREF' IS READ TO CALCULATE THE SIZE OF THE PERTURBATION. IF C C C C (II.EQ.l) READ(2,15)D,AREARE,RADREF THE NEW POSITION OF THE MESH POINTS IS CALCULATED. 25 C C C C C C C C DO 25 1=1,N XN (I)=X(I)+M(I)*DCOS(T(I))*D YN(I)=Y(I)+M(I)*DSIN(T(I) ) *D AFTER THE NEW POSITIONS ARE OBTAINED, "SHOCK" IS CALLED TO CALCULATE THE NEW MACH NUMBERS. THE RELOCATED MESH POINTS ARE THEN CALCULATED BY "MESH". FINALLY THE NEW DIRECTION OF PROPAGATION ARE CALCULATED BY "GEOM". 300 C C C C C C C C C C C C C C C C C 131 20 C C C CALL SHOCK(N) CALL MESH (N,5,998) CALL GEOM(X,Y,T,N) FOR CYLINDRICAL WAVEFRONTS, THE FUNCTION "AREA" CALCULATES THE AREA ENCLOSED BY THE WAVEFRONT AND THE AREA IS COMPARED TO THAT OF THE REFERENCE FRONT. 'YAREA' IS THE AREA PERTURBATION. XAREA(J)=RADREF YAREA(J)=AREA(N)-AREARE YAREA(J)=YAREA(J)/AREARE IF (II.EQ.2) GOTO 131 THE LENGTH IS COMPARED TO THAT OF THE REFERENCE FRONT. 'YLENG' IS THE RADIAL PERTURBATION. CALL MAXM(RADIUS,K,N) XLENG(J)=RADREF YLENG(J)=RADIUS(K)-RADREF YLENG(J)=YLENG(J)/RADREF THE WAVEFRONT IS PLOTTED WHEN THE TIME STEPS ARE MULTIPLES OF 20. IF (MOD(J,20).EQ.O) CALL DATAPL(X,Y,N,0) CONTINUE GOTO 997 AFTER THE COMPLETION OF THE REQUIRED NUMBER OF STEPS, THE RADIAL PERTURB., 85 C C C THE AREA PERTURB. VS RADIUS CURVE IS CURVE IS PLOTTED. 998 997 505 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C CONTINUE CALL DATAPL(X,Y,N,0) 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) CALL PLOTND STOP END THIS SUBROUTINE CALCULATES THE LOCATIONS OF THE NEW MESH POINTS GIVEN THE CONDITION THAT THEIR ANGULAR SEPARATION REMAIN THE SAME. SECOND ORDER REPEATED LINEAR INTERPOLATION IS USED TO CALCULATE THE LOCATIONS OF THE MESH POINTS, AND LINEAR INTERPOLATION IS USED TO CALCULATE THE VALUE OF THE MACH NUMBER AT THE NEW MESH POINT. 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 / L J / 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) XN,YN INPUT CARTESIAN COORDINATES OF POSITION RN,THN CALCULATED POLAR COORDINATES FROM INPUT MN INPUT STRENGTH OF MACH NUMBER X,Y REDISTRIBUTED CARTESIAN COORDINATES OF POSITION TH PRESCRIBED POLAR ANGLE FOR INTERPOLATION R INTERPOLATED RADIAL DISTANCE M REDISTRIBUTED STRENGTH OF MACH NUMBER XC CENTER OF POLAR COORDINATE S,SN DISTANCE BETWEEN MESH POINT FOR INTERPOLATION OF MACH NUMBER MFIR,MLAS MACH NUMBER OF THE BOUDARY POINTS. THE POLAR COORDINATES OF THE INPUT MESH POINTS IS EVALUATED. 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 C C C C C C C C THE FOLLOWING DO-LOOP CHECKS TO SEE IF ANY ABNORMAL CONDITION EXISTS. IF SUCH CONDITION EXISTSM THEN THE COMPUTATION IS STOPPED , THE POINT AT WHICH THE ABNORMALITY EXIST IS PRINTED, AND THE WAVEFRONT AND THE PERTURBATION CURVES ARE PLOTTED BY THE MULTIPLE RETURN 998. 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('0 ,'*****ABNORMAL TERMINATION AT J=',I3, AND AT 1=',13/' THE', 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 1 1 1 C C C C C C C C C C C C C C C C C R ARE OBTAINED BY REPEATED LINEAR INTERPOLATION OF RN FROM THE ANGLES THN AND TH. THE SUBROUTINE "RELINP" PERFORMS THE INTERPOLATION. AFTER R IS OBTAINED, M IS CALCULATED USING LINEAR INTERPOLATION OF THE DISTANCE BETWEEN R AND RN. THE DISTANCES ARE S AND SN. CALL RELINP(THN,RN,TH,R,N) THE NEW MESH POINTS ARE TO BE REPRESENTED IN CARTESIAN COORDINATES. 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 ) - ( Y N ( D - Y N ( 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 INSTEAD THE VALUES ARE OBTAINED DIRECTLY FROM "SHOCK" THROUGH COMMON BLOCK #LBP#. C C C M(1)=MFIR M(N)=MLAS RETURN END C C C C C C C C C C C 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". 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 C C C C C Al A2 MA OLD AREA BETWEEN MESH POINTS I AND I + l . NEW AREA BETWEEN MESH POINTS I AND I+l STRENGTH AT MID POINTS OF MESH POINT I I+l. 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 C C C C C C C C C THE THE THE AND FOR THE INTERIOR POINTS, THE NEW MACH NUMBER 'MA' BETWEEN TWO MESH POINTS ARE EVALUATED. IF (I.EQ.l.OR.I.EQ.NN) GOTO 111 CALL RUNGE(MA,Al,A2) MN(I)=MA GOTO 10 FOR THE BOUNDARY POINTS, THE INITIAL VALUE FOR INTEGRATION IS THE MACH NUMBER AT THE MESH POINT SINCE ITS VALUE IS NOT TO BE INTERPOLATED. 111 IF (I.EQ.l) MA=M(1) GOTO 112 88 C C C C C THE MACH NUMBER AT N-1 IS TO BE INTERPOLATED SO SPECIAL PROCEDURE IS TAKEN. 112 C C C C C SPECIAL VARIABLE NAMES MFIR' AND 'MLAS' ARE GIVEN TO THE MACH NUMBERS AT THE BOUNDARY POINTS. 1 10 C C C C C C C C C C CALL RUNGE(MA,Al,A2) MN(NN)=MA MA=M(N) CALL RUNGE(MA,Al,A2) I F ( I . E Q . l ) MFIR=MA IF (I.EQ.NN) MLAS=MA CONTINUE RETURN END THE SUBROUTINE "FUNC" CONTAINS THE DERVATIVE OF MACH NUMBER AS FUNCTION OF AREA TO BE USED IN RUNGE-KUTTA INTEGRATION PROCEDURE. THIS IS THE WHITHAM'S AREA RULE OF SHOCK WAVE. MS DEPENDENT VARIABLE,MACH NUMBER DF DERIVATIVE A INDEPENDENT VARIABLE,AREA SUBROUTINE FUNC(MS,DF,A) I M P L I C I T REAL*8(A-H,M,O-Z) COMMON /LG/GAMMA C C C C C C 'GAMMA' IS THE IDEAL GAS CONSTANT I T IS AVAIBLE IN THE COMMON BLOCK #LG# AND IS EQUAL TO 5/3. YL SIGMA SQUARED 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 C C C C C C C C THIS SUBROUTINE IS THE RUNGE-KUTTA INTEGRATION PROCEDURE. THIS PROCEDURE USES MODIFIED EULER'S METHOD TO ESTIMATE THE MACH NUMBER AT A NEW AREA GIVEN THE DERVATIVE OF THE MACH NUMBER AS A FUNCTION OF AREA IN THE SUBROUTINE ' F U N C . Y FUNCTION TO BE INTEGRATED. XI,X2 LIMITS OF INTEGRATION. 89 SUBROUTINE RUNGE(Y,XI,X2) IMPLICIT REAL*8(A-H,M,0-Z) C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C DF1 IS THE SLOPE AT XI. CALL FUNC(Y,DFl,X1) THE 1ST APPROXIMATION ' Y l IS OBTAINED. 1 Y1=DF1*(X2-X1)+Y1 THE SLOPE 'DF2' AT Y l IS OBTAINED. CALL FUNC(Y2,DF2,X2) THE COMPLETE INTEGRATION IS DONE BY THE AVERAGE SLOPE. Y=(DF1+DF2)/2*(X2-X1)+Y RETURN END THIS SUBROUTINE CALCULATES THE ANGLES OF THE MESH POINTS USING THE INFORMATION OBTAINED. THREE MESH POINTS ARE USED TO CALCULATE THE DIRECTION OF PROPAGATION OF THE MIDDLE MESH POINT. THEY ARE TREATED TO BE SEGMENTS OF A CIRCLE. AFTER THE CENTER OF THE CIRCLES ARE LOCATED, THEN THE NORMALS OF THE CIRCLES CAN BE FOUND. SUBROUTINE GEOM(X,Y,T,N) IMPLICIT REAL*8 (A-H,0-Z) COMMON /LPI/PI /LS/II /LANGEL/PANGLE REAL*8 X(N),Y(N),T(N) T(1)=0. /LXC/XC/LNN/NN THE FOLLOWING ARE PROCEDURES TO EVALUATE THE CENTERS OF CIRCLES FORMED BY THE MESH POINTS. 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 IF THE CURVATURE OF THE CIRCLE IS TOO SMALL, THEN THE THREE MESH POINTS ARE TREATED AS 90 C C C POINTS ON A STRAIGHT LINE. THE NORMAL IS THEN EVALUATED FROM THE ST. LINE FORMULA. IF IF C C C C C (II.EQ.2.AND.DABS(V).LE.1D-3) GOTO 10 (II .EQ.LAND.DABS (V) .LE.lD-9) GOTO 10 'A' IS THE X-COORDINATE OF THE CENTER. 'B' IS THE Y-COORDINATE OF THE CENTER. 'TN' IS THE NORMAL OF THE CIRCLE. 20 C C C C C C C C C C C A=(P*U-Q*TT)/V B=(P*S-Q*R)/(-V) TN=DATAN2(B-Y(I),A-X(I)) GOTO 2 IF THE CURVATURE IS TOO SMALL, THEN THE NORMAL IS TAKEN TO BE PERPENDICULAR TO A LINE JOINING THE THREE MESH POINTS. 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) THE NORMALS ARE CHECKED TO MAKE SURE THAT THE CHANGES IN DIRECTION DO NOT EXCEED 90 DEGREE, MAKING SURE THAT THE WAVEFRONT REMAINS CONVERGING. 2 C C C C C C C C C C C C C C C IF (DABS(TN-T(I-l)).GE.1.5) T(I)=TN 1 CONTINUE TN=TN-PI THE BOUNDARY CONDITION THAT THE FLOW MUST BE PARELLEL TO THE WALL IS IMPOSED UPON THE BOUNDARY POINT. T(N)=-PANGLE RETURN END THIS SUBROUTINE TAKES IN THE DATA STORED ARRAYS X AND Y WHICH ARE IN DOUBLE PRECISION AND CONVERT THE VALUES INTO SINGLE PRECISION, THEN PLOT THEM WITH A SYMBOL J AT EACH DATA POINT. SUBROUTINE DATAPL(X,Y,N,J) REAL*8 X(N),Y(N) REAL A(100) ,B (100) THE FOLLOWING DO-LOOP DOES THE 91 C C CONVERSION INTO SINGLE PRECISION. 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 C C C C THIS SUBROUTINE GIVES THE I N I T I A L SHAPE OF THE WAVEFRONT AS A SEMI-CIRCLE. THE PERTURBATION IS SMOOTH AND SUPERIMPOSED ON THE CIRCULAR WAVEFRONT IN SINUSOIDAL FORM. SUBROUTINE CIRCLE(X,Y,T,SM,N) I M P L I C I T REAL*8 (A-H,0-Z) COMMON / L P I / P I /LP/TH /LXC/XC /LL/DELTA,L COMMON /LRA/R /LG/GAMMA /LANGEL/ANGEL REAL*8 T H ( 1 0 0 ) , R ( 1 0 0 ) REAL*8 XN(100),YN(100),MN(100),TN(100) REAL*8 X ( N ) , Y ( N ) , T ( N ) , S M ( N ) INTEGER FMT(1)/'*'/,YES/ Y'/fNO/'N'/rFORM COMMON /LB/XN,YN,MN,TN WRITE(7,1) 1 FORMAT('0','ENTER L,DELTA, AND I N I T I A L ' , . ' STRENGTH OF SHOCK') 1 C C C C C C THE THE THE ARE READ AZIMUTHAL NUMBER ' L , PERTURBATION AMPLITUDE 'DELTA' AND I N I T I A L SHOCK STRENGTH 'SMI' READ. 1 (5,FMT) L,DELTA,SMI C C C C 'MM' IS THE LARGEST MESH POINT NUMBER THAT IS MODIFIED BY THE PERTURBATION. MM=N*2/L ANGEL=PI C C C C THE POLAR COORDINATES 'R' AND T H ' A CIRCLE IS FIRST CALCULATED. 1 OF DO 30 1=1,N TH(I)=PI*(N-I)/(N-1) R(D=1. C C C C THE MACH NUMBERS AT THE MESH POINTS ARE INITIALIZED. 30 SM(I)=SMI CONTINUE 'XC IS THE CENTER OF CONVERGENCE. XC=R(1)+DELTA THE PERTURBATION IS INSERTED INTO THE I N I T I A L 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 'ANSWER ''Y'' OR ''N '') READ (5,101) FORM FORMAT(A4) I F (FORM.EQ.NO) GOTO 102 STRENGTH?'/ 1 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 THE PERTURBATION. AT THE CENTER OF 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 I N I T I A L SHAPE IS WRITTEN INTO OUTPUT UNITS '6' AND '7'. WRITE(7,200) CURV,(SM(I),1=1,10) 93 200 102 2 C C C C C FORMAT(1X,11F10.5) WRITE(6,2) N, L,DELTA,SMI FORMAT('0' CIRCLE'/'NO. OF POINTS=' , 13,2X ,'L=' ,12/ 'INITIAL PERTURBATION=',F7.5/ 'INITIAL SHOCK STRENGTH=',F7.3) RETURN END THIS SUBROUTINE CALCULATES THE INITIAL SHAPE OF A WEDGED CHANNEL WITH INPUT WEDGE ANGLE AND FRONT ANGLE. 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('0 ,'ENTER ANGLE OF WEDGE, ANGLE OF FRONT, AND', .' INITIAL STRENGTH OF SHOCK') 1 C C C C THE WEDGE ANGLE, THE FRONT ANGLE, AND THE MACH NUMBER IS READ. READ(5,FMT) ANGLE,TANG,SMI PANGLE=ANGLE*P1/180 C C C C 'A' IS THE X-COORDINATE MESH POINT. OF THE NTH A=DTAN(TANG*PI/180)/(DTAN(TANG*PI/180)+DTAN(PANGLE)) IF (ANGLE.EQ.0.0) A=l.-DCOTAN(TANG*PI/180) C C C XC = 1. C C C THE CENTER OF CONVERGENCE IS SET TO BE EQUAL TO 1. THE BOUNDARY POINTS ARE FIRST CALCULATED. 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 C THE REST OF THE MESH POINTS ARE CALCULATED. 20 DO 20 1=2,NN X(I)=X(N)*(I-1)/(N-1) Y(I)=Y(N)*(I-1)/(N-1) 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 THE I N I T I A L CONDITIONS ARE PRINTED IN 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 L E F T . 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 / L P I / P I REAL*8 T ( 1 0 0 ) , T H ( 1 0 0 ) , U ( 1 0 0 ) , P I REAL*4 A,B 'A' IS I N I T I A L I Z E D TO BE ZERO. A=0 . DO 10 1 = 1,N U(I)=TH(I)-PI B=DABS(T(I)-U ( I ) ) I F (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 I N I T I A L MAX VALUE TO EQUAL TO ZERO. . 'G' IS SET TO G=0.0 DO 10 1=1,N IF ( A ( I ) . L T . G ) GOTO 10 IF 'A' IS BIGGER THAN IS THE NEW MAX. 'G , THEN 1 'A' 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 R A D I U S ( 1 0 0 ) , A N G L E ( 1 0 0 ) 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 15 C C C C C D015 1 = 1,N X(I)=XC+R(I)*DCOS(TH (I)) Y(I)=R(I)*DSIN(TH(I)) CALL GEOM(X,Y,TN,N) THE PERTURBATION IS INSERTED INTO THE CIRCULAR WAVEFRONT. THE PERTURBATION IS IN THE SHAPE OF AN CIRCULAR ARC. 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 C C THE 'X' Y' COORDINATES ARE BEING MODIFIED TO TAKE THE PERTURBATION. 1 10 C C C C C C C DO 10 1=2,MMM X (I)=SLEN-SLEN*DCOS(ANG*(1-1)/(MMM)) Y(I)=SLEN*DSIN(ANG*(I-1)/(MMM)) CONTINUE ADDITION INFORMATION ABOUT THE DISTRIBUTION OF MACH NUMBERS ARE ASKED. IF THE ANSWER FORM EQUALS TO 'NO', THEN THE MACH NUMBER IS UNIFORM AND IS NOT BEING MODIFIED, OTHERWISE THEY ARE TO BE MODIFIED BY THE FOLLOWING PROCEDURE. 1 10 0 101 C C C C C C C C C C C 103 C C C 1 WRITE(7,100) FORMAT('0','DO YOU WANT A NON-UNIFORM STRENGTH?'/ 'ANSWER ''Y'' OR ''N''') READ (5,101) FORM FORMAT(A4) IF (FORM.EQ.NO) GOTO 102 IF A NON-UNIFORM SHOCK STRENGTH IS WANTED, THEN THE STRENGTH AT THE CENTER OF THE PERTURBATION IS EVALUATED USING THE STRONG SHOCK APPROXIMATION AND THE REST OF THE PERTURBATION IS LINEARLY INTERPOLATED. CURV=SLEN GN=l+2/GAMMA+DSQRT(2*GAMMA/(GAMMA-1)) SM(1)=SMI*CURV**(-1/GN) MODIFICATION OF THE MESH POINT MACH NUMBERS. DO 103 1=2,MM SM(I)=SM(1) PRINTOUT ON UNIT '6' INDICATES SUCESSFULNESS OF MACH NUMBER MODIFICATION. 98 C C C THE I N I T I A L CONDITIONS ARE PRINTED IN THE OUTPUT UNIT ' 6 . WRITE(7,200) CURV,(SM(I),1=1,10) FORMAT(lX 11F10.5) WRITE(6,2) N,L,DELTA,SMI F O R M A T ( 0 , CIRCLE'/'NO. OF POINTS= ,13 , 2X,'L=' ,12/ 'INITIAL PERTURBATION=*,F7.5/ 'INITIAL SHOCK STRENGTH=' ,F7 . 3) RETURN END 1 200 102 2 f 1 1 1 1
- 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
Page Metadata
Item Metadata
Title | Stability of converging shock wave |
Creator |
Fong, Kenneth Sau-Kin |
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 |
Aggregated Source Repository | 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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0094383/manifest