1 Introduction2 Governing Equations3 Numerical Methods3.1 Shepherd's Algorithm3.2 Reynolds' Algorithm with Specified Shock Velocity3.3 Chapman-Jouguet Detonation Velocity3.3.1 Reynolds' Algorithm3.3.2 Modified Reynolds' Algorithm3.3.3 Minimizing Initial Velocity4 Validation5 Matlab Scripts5.1 Subfunction eq_state5.2 Subfunction FHFP_noshk5.3 Subfunction FHFP_shk5.4 Subfunction noshk_calc5.5 Subfunction PostShock_noshk5.6 Subfuntion PostShock_shk5.7 Subfunction shk_calc5.8 Subfunction soundspeed_eq5.9 Subfunction state6 Initial Velocity as a Function of Density Ratio6.1 Derivation6.2 CJ Point Analysis6.3 Derivatives of Pressure6.4 Thermodynamic Analysis6.5 Perfect Gas AnalysisNumerical Solution Methods for Shock and Detonation JumpConditionsS. Browne and J. E. ShepherdGALCIT Report FM2006.006July 2004AbstractPrograms that simulate constant volume explosions and ZND profiles are very important toolsfor making predictions on the outcome of expensive experiments. Each of these programs usesthe post-shock state as initial conditions for its computation. The post-shock state which is thesecond intersection of the Rayleigh line and the reactant Hugoniot can be found analytically for gasmixtures with constant properties, but must be solved iteratively for gas mixtures with temperaturedependent properties. Our software has traditionally iterated on a single variable, the density ratio,to find the post-shock state. Unfortunately, for very strong shock waves, the output of this singlevariable algorithm is incorrect. Reynolds’ STANJAN program also has the capability of calculatingthe post-shock state and is successful in a much wider range than our one variable method. Thisreport outlines the algorithms that STANJAN uses to calculate both the post-shock state and theChapman-Jouguet detonation velocity.1 IntroductionIn the field of detonation physics, simulations are invaluable tools because they allow scientists to predictproperty and species profiles before running expensive experiments. Simulation programs for constantvolume explosions and ZND detonations1require a system of equations to solve and a set of initialconditions. Generally the initial conditions are a post-shock state which must be calculated correctly fromthe initial state. The shock jump conditions can be solved analytically for gas mixtures with constantspecific heat. Unfortunately, this is only a small subset of all gas mixtures, and most combustible gasmixtures of interest do not have constant specific heat. In systems with variable specific heat, the shockjump conditions must be solved iteratively.The solutions to the jump conditions occur when the Rayleigh line and the reactant Hugoniot in-tersect. To date, we have solved this problem by using Shepherd’s algorithm [7], a Newton scheme, toiterate on the density ratio. As the shock speed increases, the upper intersection of the Rayleigh lineand reactant Hugoniot occurs when the Hugoniot is nearly vertical. Shepherd’s algorithm fails in theseextreme cases while Reynolds’ algorithm [5] does not.STANJAN, the original implementation of Reynold’s algorithms, was written in Fortran. We haveimplemented Reynolds’ algorithms with Matlab to reconstruct our software. Originally, Reynolds andShepherd used the Chemkin libraries [3] to calculate certain thermodynamic properties based on anexplicit reaction mechanism. Recently, Goodwin has developed a similar library called Cantera [2],which we have chosen to use in our Matlab software because it is open source.In this report, we first present the governing equations of the post-shock problem in Section 2. InSection 3.1, we give Shepherd’s algorithm while Sections 3.2 and 3.3 outline Reynolds’ algorithms forcalculating both the post-shock state and the Chapman-Jouguet detonation speed. Section 4 discussesthe validity of Reynolds’ algorithms and the order of convergence of the iterative scheme. Section 4 also1Zeldovich, von Neumann, Doering detonation structure1gives a comparison between results obtained with the Matlab implementation and results created withShepherd’s CV and ZND programs, the corresponding Fortran implementations [7].2 Governing EquationsFor simplicity, we assume that our shock wave has no volume and is simply a plane passing throughour mixture. We also assume that our mixture does not begin reacting until after it has reached thepost-shock state, allowing us to treat the post-shock calculation as a nonreacting fluids problem. Theshock jump conditions given below are the governing equations for our system. These expressions are inthe fixed-shock frame of reference .• Continuityρ1w1= ρ2w2• Momentum ConservationP1+ ρ1w21= P2+ ρ2w22• Energy Conservationh1+w212= h2+w222We also assume that our gas mixture is sufficiently rarefied such that the ideal gas equation of state,which isP v =RWTapplies. This equation of state introduces temperature into the parameter space, but enthalpy is afunction strictly of temperature for and ideal gas, i.e. h = f (T ).If we eliminate the post-shock velocity from the conservation equations, we can derive the Rayleighline,P = P1− ρ21w21(v − v1) (1)and the reactant Hugoniot,P = 2h − h1v + v1+ P1(2)This simplified representation of the system is depicted in Figure 1. T he initial and post-shock states areindicated by the squares. The post-shock state properties are the quantities of interest in our problem.We will find these properties by solving the above equations simultaneously with the equation of state.3 Numerical MethodsIn this sec tion, we will first outline Shepherd’s algorithm [7] and then outline the successful algorithmsthat Reynolds’ employs in his STANJAN code [5]. STANJAN implements many algorithms but wewill focus on two specific algorithms: the method for finding the post-shock state for a specified shockvelocity and the method for determining the Chapman-Jouguet velocity. These two methods althoughsimilar solve slightly different problems.20 0.25 0.5 0.75 1 1.25 1.5051015202530Pressure (atm)specific volume (m3/kg)Rayleigh LineReactant HugoniotInitial StatePost-shock StateFigure 1: The Rayleigh line and reactant Hugoniot for stoichiometric hydrogen and air traveling at theChapman-Jouguet detonation velocity with initial pressure of 1 atm and initial temperature of 300 K.3.1 Shepherd’s AlgorithmThe original implementation of Shepherd’s algorithm in Fortran was called SJUMP and iterates on thedensity ratio to solve the following se