Combustion Modelling of Solid Propellants
Combustion Modelling of Solid Propellants
Combustion Modelling of Solid Propellants
July, 2019
DT-16
INFORMA:
Que la referida Tesis Doctoral, ha sido realizada por Dª. CARMEN LÓPEZ MUÑOZ,
dentro del Programa de Doctorado ENERGÍAS RENOVABLES Y EFICIENCIA
ENERGÉTICA, dando nuestra conformidad para que sea presentada ante el Comité de
Dirección de la Escuela Internacional de Doctorado para ser autorizado su depósito.
Ciencias
Ciencias Sociales y Jurídicas
X Ingeniería y Arquitectura
INFORMA:
Ciencias
Ciencias Sociales y Jurídicas
X Ingeniería y Arquitectura
Among all energetic materials, the combustion of granulated, composite AP/HTPB and double-base
solid propellants are studied in thesis. The improvement on the design and safety in industrial
processes which use solid propellants are of major importance in industry of energetic materials.
Unexpected detonation of granular solid propellants is a key safety issue. Moreover, ignition and
subsequent combustion of granular solid propellants are also complex mechanisms of high importance
in the manufacturing process of many energetic materials, not only because already mentioned safety
issues, but also in due in the design of propulsion systems. Therefore, and in order to properly
understand this problem, a model for the characterisation of the detonation process of granular solid
propellants under shock tube conditions is developed. To be able to reach this objective, a two-phase
model, which considers the conservation equations of mass, momentum and energy and constitutive
relations of mass generation, gas-solid particle interaction, interface heat transfer and particle-particle
stress is defined. The study performed in this thesis, improve the understanding and modelling of the
deflagration-to-denotation (DDT) phenomenon in granular beds of solid propellants. Two different
models have been used to calculate pressure, temperature and porosity distributions. The first model
considered a modification of the particle momentum governing equation in order to prevent the
porosity from reaching values below the minimum value of compaction defined for packed beds of
spherical particles. The second model studied does not consider the porosity limiter in the momentum
conservation equation but represents the limitation of the porosity directly in the code by preventing
the porosity to reach values below a minimum. In addition, the results are obtained by employing
several numerical schemes and are compared against those available in the bibliographical sources
in order to assess their effectivity to predict the early stages of this transient combustion process.
The results show that developed second model accurately represents the physical behaviour of the
propellant combustion for all variables of interest becoming a predictive tool for the characterisation
of granulated solid propellants.
Modelling the combustion of composite propellants and double-base propellants is a key problem
that has focused the interest of several researches in many industrial fields such as chemical engineer-
ing, aerospace engineering or safety in industrial processes. Regarding safety, not only gun tubes, but
also rocket motors experiment pressure build-up which could increase the chamber pressure, causing
the explosion of the rocket motor being of high importance to develop models which could represent
accurately the combustion process of these materials. Modelling the combustion of these propellants
has been also of high interest for researchers in terms of optimising chemical composition to fulfil
industrial requirements, improve propellant design, development and testing activities. The analysis
of the burning rate of these materials, which means their decomposition and their subsequent com-
bustion, is one of the main objectives of combustion modelling. Therefore, in this work, a transient
multidimensional numerical model to describe the combustion of composite and double-base solid
propellants is presented. The kinetics of the model is described by considering firstly, a change of
phase of the solid propellant from condensed to gas and secondly, a reduced chemistry scheme which
defines simplified chemical reactions to represent the combustion itself. To couple both processes,
mass and species conservation, as well as temperature continuity are imposed in the burning surface
in which the burning rate will represent a key factor. Moreover, an energy balance is also applied at
the burning surface which represents that heat that the gas transfers to the burning surface is invested
3
firstly, in raising the surface temperature to produce the phase change and secondly, in warming the
condensed phase by conduction. The results obtained in the combustion modelling of both, composite
and double-base solid propellants, are compared against experimental test and results present in the
existing literature. An absolute error between 0.4 mm/s and 1 mm/s is obtained being in the order
of magnitude of the experimental error. In addition, the pressure oscillations produced in a rocket
motor and the depressurisation suffered in a combustion chamber, as well as at the exit of a barrel,
are very interesting technological problems which have not been deeply studied in the literature.
Therefore, the results obtained in the numerical modelling of transient combustion of composite and
double-base solid propellants are presented in this thesis too. To perform this analysis, a pressure
ramp is applied during a short interval of time at the open side of the already defined geometries. The
developed tool confirms its robustness for the prediction of composite and double-base propellants
combustion behaviour in multidimensional scenarios with transient environmental conditions such as
rocket engines and/or base bleed units in ballistics applications.
4
Resumen
En esta tesis se han estudiado, de entre todos los materiales energéticos, la combustión de propulsantes
sólidos de tipo granulado, composite y doble-base. Tanto la mejora del diseño, como de los procesos
de seguridad y las instalaciones que usan propulsantes solidos son asuntos de gran importancia en la
industria de materiales energéticos.
La detonación inesperada de propulsantes sólidos de tipo granular es clave en materia de seguridad.
Además, la ignición y la subsecuente combustión de dicho propulsante es un mecanismo complejo
también de gran importancia en los procesos de fabricación de muchos materiales energéticos, no
solo por los problemas de seguridad ya mencionados, sino también por el diseño de los sistemas de
propulsión. Por lo tanto, en este trabajo se ha desarrollado un modelo para caracterizar el proceso
de detonación de propulsantes sólidos granulados bajo condiciones de tubo de choque. Para ser
capaz de alcanzar este objetivo, se ha definido un modelo bifásico que considera las ecuaciones de
conservación de masa, momento y energía, así como los términos fuente de generación de masa,
interacción sólido-gas, calor interfacial y tensión intergranular. El estudio llevado a cabo en esta tesis
mejora la comprensión y modelado del fenómeno de deflagración-detonación en propulsantes sólidos.
Las distribuciones de presión, temperatura y porosidad se han obtenido mediante dos modelos. El
primero considera la modificación de la ecuación de conservación de la cantidad de movimiento de
la fase sólida a fin de evitar que la porosidad alcance un valor por debajo del mínimo permitido
para partículas esféricas. Por otro lado, el segundo modelo limita la porosidad directamente en
el código impidiendo que esta variable alcance valores por debajo del mínimo valor físico en lugar
de considerar esta limitación modificando la ecuación de cantidad de movimiento. Además, se han
empleado distintos esquemas numéricos en la resolución del modelo a fin de poder comparar los
resultados obtenidos entre ellos y con los datos disponibles en la bibliografía. Esto permitirá estudiar
su efectividad para predecir las primeras fases de del proceso de combustión. Los resultados obtenidos
muestran que el segundo modelo desarrollado representa de forma precisa el comportamiento físico
de la combustión del dicho propulsante convirtiéndose por tanto en una herramienta predictiva para
la caracterización de la combustión de propulsantes sólidos de tipo granular.
El modelado de la combustión de propulsantes tipo composite y doble-base es un problema clave
que ha ocupado el interés de muchos investigadores en diversos ámbitos industriales como la ingeniería
química, la ingeniería aeroespacial o la seguridad en procesos industriales. En cuanto a seguridad,
los motores cohete pueden experimentar incrementos de presión en la cámara de combustión dando
lugar a la explosión del motor cohete. Por tanto, es de gran importancia desarrollar modelos que
puedan representar de forma precisa el proceso de combustión de estos materiales. El modelado de la
combustión también ha sido de gran interés para investigadores con la idea de optimizar la composición
química de los propulsantes para satisfacer requisitos industriales, así como para mejorar el diseño, el
desarrollo y las actividades de test del propulsante. El análisis de la velocidad de quemado de estos
materiales, que representa su descomposición y su subsecuente combustión, es uno de los objetivos
de más interés en el modelado de la combustión. Por lo tanto, en este trabajo se ha presentado un
modelo multidimensional transitorio para describir el proceso de combustión de propulsantes sólidos
de tipo composite y doble base. La cinética química del modelo considera primero un cambio de fase
del propulsante de fase condensada a fase gaseosa y, a continuación, un esquema químico reducido
que define reacciones químicas simplificadas para modelar la combustión en sí misma. Con el fin de
5
acoplar ambos procesos, se han impuesto en la superficie de quemado las condiciones de conservación
de masa y especies, así como la continuidad de la temperatura, expresiones en las cuales la velocidad
de quemado es determinante. Además, en la superficie de quemado se ha establecido un balance
de energía que representa la inversión del calor que el gas transfiere a la superficie de quemado en
dos procesos: el primero, en elevar la temperatura de la superficie produciendo el cambio de fase
el propulsante y el segundo, en calentar la fase condensada mediante conducción. Los resultados
obtenidos en la combustión de ambos propulsantes se comparan con resultados experimentales y
obtenidos en la bibliografía. El error absoluto obtenido varía entre 0.4 mm/s y 1 mm/s siendo del
orden de magnitud del error experimental. Además, las oscilaciones de presión que se producen en un
motor cohete y la despresurización que se sufre en la cámara de combustión y a la salida del barril,
son problemas de alto interés tecnológico que no han sido estudiados en profundidad en la literatura.
Por tanto, los resultados obtenidos al modelar la combustión transitoria de ambos propulsantes se
han presentado también en esta tesis. Para realizar este último análisis, se ha aplicado una rampa de
presión durante un corto intervalo de tiempo a la salida de la cámara de combustión de las geometrías
previamente utilizadas. Los resultados confirman la robustez de la herramienta desarrollada para la
predicción del comportamiento de la combustión de propulsantes de tipo composite y doble-base en
escenarios multidimensionales con condiciones ambientales transitorias como es el caso de motores
cohete y unidades base-bleed en aplicaciones balísticas.
6
Acknowledgements
The research of this thesis has been performed due to the financial support of the Spanish Ministry
of Economy and Competitiveness and the European Commission throughout the ALMENA project
(RTC-2016-5194-8) obeying to the interest of the company Expal Systems S.A.
I would like to thank to Dr. José Ramón García Cascales and Dr. Francisco Javier Sánchez
Velasco for giving me the opportunity to be part of this project as well as for their guidance and
dedication. To them and to Dr. Ahmed Bentaib I also owe the opportunity to make my pre-doctoral
stay in the Institut de Radioprotection et de Sûreté Nucléaire (IRSN) of Paris.
I would like to dedicate a special mention to my college Francisco Nicolás Perez, from whom I
had the pleasure to learn a great deal and share many moments during these years.
I am extremely thankful to my family and my beloved ones for supporting and encouraging me
as well as for their love and patience. Specially to my parents, for providing me access to education,
for relying on me and, above all, for their help beyond what is possible to imagine.
Finally, I would like to thank all my friends. Particularly, I would like to express all my love and
gratitude to those whom I can sincerely call family, for listening and understanding me every day
without exception.
7
Contents
1 Introduction 21
1.1 Energetic materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.2 Technological context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.3 Thesis scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.4 Thesis objectives and structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.4.1 Thesis objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.4.2 Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
9
CONTENTS
3 Numerical methods 63
3.1 Rusanov numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2 MacCormack and MacCormack-TVD numerical scheme . . . . . . . . . . . . . . . . . 66
3.3 Lax-Wendroff scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
10
CONTENTS
5 Conclusions 211
5.1 Granulated propellant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
5.2 Composite AP/HTPB propellant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
5.3 Double-base homogeneous propellant . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214
Nomenclature 219
Bibliography 230
11
List of Figures
4.1 Influence of Reynolds number in the Nusselt number for several values of porosity 0.6
(above-left), 0.5 (above-right), 0.4 (below-left) and 0.3 (below-right). . . . . . . . . . 74
4.2 Distributions of pressure, porosity and gas velocity for Rusanov (left) and MacCormack-
TVD (right) with Friction 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.3 Distributions of pressure, porosity and gas velocity for Rusanov (left) and MacCormack-
TVD (right) with Friction 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.4 Distributions of pressure, porosity and gas velocity for Rusanov (left) and MacCormack-
TVD (right) with Friction 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.5 Influence of Reynolds number in the friction coefficient for several values of porosity
0.6 (above-left), 0.5 (above-right), 0.4 (below-left) and 0.3 (below-right). . . . . . . . . 78
4.6 Distributions of pressure obtained with Rusanov (above) and MacCormack-TVD (be-
low) numerical schemes for Nusselt 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.7 Distributions of pressure obtained with Rusanov (above) and MacCormack-TVD (be-
low) numerical schemes for Nusselt 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.8 Distributions of pressure obtained with Rusanov (above) and MacCormack-TVD (be-
low) numerical schemes for Nusselt 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.9 Pressure distribution obtained with Rusanov (left) and MacCormack TVD (right) for
20 µs, 40 µs and 60 µs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.10 Porosity distribution obtained with Rusanov (left) and MacCormack TVD (right) for
20 µs, 40 µs and 60 µs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.11 Gas temperature distribution obtained with Rusanov (left) and MacCormack TVD
(right) for 20 µs, 40 µs and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.12 Particle temperature distribution obtained with Rusanov (left) and MacCormack TVD
(right) for 20 µs, 40 µs and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
13
LIST OF FIGURES
4.13 Gas velocity distribution obtained with Rusanov (left) and MacCormack TVD (right)
for 20 µs, 40 µs and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.14 Porosity distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . 87
4.15 Pressure distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . 88
4.16 Gas temperature distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . 88
4.17 Particle temperature distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . 89
4.18 Flame front propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.19 Porosity distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.20 Pressure distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . 91
4.21 Gas temperature distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . 91
4.22 Particle temperature distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . 92
4.23 Flame front propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.24 Porosity distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . 93
4.25 Pressure distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . . . . . . 94
4.26 Gas temperature distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . . . 94
4.27 Particle temperature distribution for 20, 40 and 60 µs. . . . . . . . . . . . . . . . . . 95
4.28 Flame front propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.29 Geometry of AP/HTPB strand configuration. . . . . . . . . . . . . . . . . . . . . . . . 103
4.30 Comparison between calculated and experimental burning rate at 4 MPa (left) and 7
MPa (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.31 Geometry of AP/HTPB sandwich configuration. . . . . . . . . . . . . . . . . . . . . . 105
4.32 Comparison of calculated and experimental burning rate. . . . . . . . . . . . . . . . . 106
4.33 Velocity distribution of gas phase in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1
MPa, 0.5 MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.34 Velocity distribution of gas phase in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1
MPa, 0.5 MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
4.35 Velocity distribution of gas phase in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1
MPa, 0.5 MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.36 AP distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
4.37 AP distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.38 AP distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
4.39 HTPB distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
4.40 HTPB distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.41 HTPB distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
4.42 Burning rate distribution along the burning surface at t = 0.015 s and t = 0.065 s at
70 MPa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4.43 Mass fraction distribution after burning surface at t = 0.015 s and t = 0.065 s at 7
MPa for AP (left) and HTPB (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
4.44 Temperature distribution for t = 0.04 s along x-direction (left) and temperature in
condensed phase before burning surface at t = 0.03 s for 7 MPa to 0.5 MPa (right). . 118
4.45 Temperature distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.46 Temperature distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
4.47 Temperature distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
14
LIST OF FIGURES
4.48 Heat flux distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
4.49 Heat flux distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
4.50 Heat flux distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
4.51 Pressure distribution in t = 0.005 s, t = 0.005000005 s, t = 0.005000009 s, t = 0.0051
s, t = 0.0052 s, t = 0.0053 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . 126
4.52 Pressure distribution against time for P1 = (515 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . 127
4.53 Zoom of pressure distribution against time for P1 = (515 · 10 6 , 85·10 6 , 85·10 6 ). . . 127
4.54 AP mass fraction distribution in t = 0.005 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
4.55 AP mass fraction distribution in t = 0.01 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
4.56 AP mass fraction distribution in t = 0.02 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.57 AP mass fraction distribution in t = 0.005 s, t = 0.0051 s and t = 0.0052 s (above-
below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.58 HTPB mass fraction distribution in t = 0.005 for 0.5 MPa, the pressure ramp and 7
MPa boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.59 HTPB mass fraction distribution in t = 0.01 for 0.5 MPa, the pressure ramp and 7
MPa boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.60 HTPB mass fraction distribution in t = 0.02 for for 0.5 MPa, the pressure ramp and
7 MPa boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . 132
4.61 HTPB mass fraction distribution in t = 0.005 s, t = 0.0051 s and t = 0.0052 s (above-
below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.62 AP mass fraction (%) distribution against time for P1 = (515 · 10 6 , 85·10 6 , 85·10 6 ). 133
4.63 HTPB mass fraction (%) distribution against time for P1 = (515·10 6 , 85·10 6 , 85·10 6 ).
133
4.64 Temperature distribution in t = 0.005 s, t = 0.0051 s and t = 0.0052 s (above-below). 134
4.65 Temperature distribution in t = 0.01 s for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.66 Temperature distribution in t = 0.02 s for for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.67 Heat flux distribution for ramp pressure boundary condition in t = 0.005 s, t = 0.0051
s and t = 0.0052 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.68 Heat flux distribution in t = 0.01 s for 0.5 MPa, the pressure ramp and 7 MPa boundary
conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.69 Heat flux distribution in t = 0.02 s for for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4.70 Burning rate distribution along the burning surface at t = 0.021 s for fixed 7 MPa,
fixed 0.5 MPa and pressure ramp boundary conditions. . . . . . . . . . . . . . . . . . 138
4.71 Burning rate distribution along the burning surface at t = 0.081 s for fixed 7 MPa,
fixed 0.5 MPa and pressure ramp boundary conditions. . . . . . . . . . . . . . . . . . 138
4.72 Temperature at the burning surface at t = 0.021 s for fixed 7 MPa, fixed 0.5 MPa and
pressure ramp boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.73 Temperature at the burning surface at t = 0.081 s for fixed 7 MPa, fixed 0.5 MPa and
pressure ramp boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.74 Comparison of calculated and experimental burning rate. . . . . . . . . . . . . . . . . 151
4.75 Velocity distribution of gas phase in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2
MPa, 1.1 MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
15
LIST OF FIGURES
4.76 Velocity distribution of gas phase in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2
MPa, 1.1 MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
4.77 Velocity distribution of gas phase in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2
MPa, 1.1 MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
4.78 Fuel distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
4.79 Fuel distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
4.80 Fuel distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below) (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
4.81 Oxidiser distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
4.82 Oxidiser distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
4.83 Oxidiser distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
4.84 DR1 distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
4.85 DR1 distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
4.86 DR1 distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
4.87 DR2 distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
4.88 DR2 distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
4.89 DR2 distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
4.90 Final products distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
4.91 Final products distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
4.92 Final products distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
4.93 Temperature in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-
below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
4.94 Temperature in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-
below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172
4.95 Temperature in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-
below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
4.96 Heat flux in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-below).174
4.97 Heat flux in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-below).175
4.98 Heat flux in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-below).176
4.99 Pressure distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t
= 0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . . . . . . . 178
4.100 Pressure distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . . . . . 179
4.101 Pressure distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . . . . . 179
4.102 Pressure distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . . . . . 180
4.103 Pressure distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . . . . . 180
4.104 Pressure distribution for P4 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . . . . . 181
4.105 Pressure distribution for P3 = (950 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . . . . . 181
16
LIST OF FIGURES
4.106 Pressure distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary conditions
at t= 0.015 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
4.107 Pressure distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary conditions
at t = 0.02 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
4.108 Fuel distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t =
0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . 184
4.109 Oxidiser distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t =
0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . 185
4.110 DR1 distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t =
0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . 186
4.111 DR2 distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t =
0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . 187
4.112 Products distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t
= 0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . . . . . . . 188
4.113 Fuel mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.015 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . 189
4.114 Fuel mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.02 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
4.115 Oxidiser mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa bound-
ary conditions at t = 0.015 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . 190
4.116 Oxidiser mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa bound-
ary conditions at t = 0.02 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . 190
4.117 Fuel and oxidiser distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . 191
4.118 Fuel and oxidiser distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . 191
4.119 Fuel and oxidiser distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . 192
4.120 Fuel and oxidiser distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . 192
4.121 Fuel and oxidiser distribution for P3 = (950 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . 193
4.122 Fuel and oxidiser distribution for P3 = (950 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . 193
4.123 DR1 and DR2 distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . 194
4.124 DR1 and DR2 distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . 194
4.125 DR1 and DR2 distribution for P3 = (950 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . 195
4.126 DR1 and DR2 distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . 195
4.127 DR1 and DR2 distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . 196
4.128 DR1 and DR2 distribution for P3 = (950 · 10 6 , 85·10 6 , 85·10 6 ). . . . . . . . . . . 196
4.129 Products mass fraction distribution for P1 (above-left), P2 (above-right) and P3 (below).197
4.130 Temperature at the burning surface after applying depressurisation. . . . . . . . . . . 197
4.131 Global temperature distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t =
0.01003, t = 0.001004 s and t = 0.01005 s (above-below). . . . . . . . . . . . . . . . . 199
4.132 Global temperature distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.015 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . 200
4.133 Global temperature distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.02 s (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . 200
4.134 Temperature distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ) from t = 0.01 s to t
= 0.02 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
4.135 Temperature distribution for P1 = (510 · 10 6 , 85·10 6 , 85·10 6 ) from t = 0.01 s to t
= 0.01001 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
4.136 Temperature distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ) from t = 0.01 s to t
= 0.02 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
4.137 Temperature distribution for P2 = (750 · 10 6 , 85·10 6 , 85·10 6 ) from t = 0.01 s to t
= 0.01001 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
4.138 Temperature distribution for P3 = (950 · 10 6 , 85·10 6 , 85·10 6 ) from t = 0.01 s to t
= 0.02 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
17
LIST OF FIGURES
A I.1 Heat flux distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa (above-below). . . . . 221
A I.2 Heat flux distribution in t = 0.02 s at 2 MPa, 1 MPa, 0.5 MPa (above-below). . . . 222
A I.3 Heat flux distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
A III.1 Granulated propellant code flowchart . . . . . . . . . . . . . . . . . . . . . . . . . . 227
A IV.1 AP/HTPB and double-base propellants code flowchart . . . . . . . . . . . . . . . . . 229
18
List of Tables
2.1 Pressure, gas and particle temperature and ignition temperature data for the con-
sidered initiation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
19
20
Chapter 1
Introduction
21
CHAPTER 1. INTRODUCTION
warheads, etc. These developments are also used for commercial applications such as the High
Altitude Fuel, Air Re-generating Composition etc. Another very well known application of explosives
are space applications, such as space rockets and indian satellite launch vehicles which are used in
meteorology, weather forecasting, satellite communication, etc.,. Other less known applications are
nuclear applications, such as nuclear weapons, agriculture, medical industry (nitroglycerine tablets
for example), food industry, civil engineering, automobile industry and oil and gas industry [2].
22
1.4. THESIS OBJECTIVES AND STRUCTURE
causing the explosion of the gun tube [3]. Particularly, it is known that unexpected detonation of
granular solid propellants is a key safety issue. Moreover, ignition and subsequent combustion of
granular solid propellants are complex mechanisms key in manufacturing process of many energetic
materials, not only because already mentioned safety issues, but also in due in the design of propulsion
systems. The penetration of the hot gases into the voids promote the convective heating of the
granules as well as the granule compaction and the rapid pressurisation of the system. In case the
ignition is not correctly controlled, the process can generate compressive waves which compact the
un-reacting region of the grain bed. Under these conditions, and as it has been already explained, if
the strong pressure waves rise, they might rapidly collapse the bed leading to a compressive reaction
and a violent explosion. Therefore, and in order to properly understand this problem, a model for the
characterisation of the detonation process of granular solid propellants under shock tube conditions
has been developed. To be able to reach this objective, a two-phase model, which considers the
conservation equations of mass, momentum and energy and constitutive relations for mass generation,
gas-solid particle interaction, interface heat transfer and particle-particle stress, has been defined.
Modelling the combustion of composite propellants and double-base propellants is a key problem
that has focused the interest of several researches in many industrial fields such as chemical engin-
eering, aerospace engineering or safety in industrial processes. Regarding safety, not only gun tubes,
but also rocket motors, experiment pressure build-up which could increase the chamber pressure,
causing the explosion of the rocket motor [3]. As a result, it is of high importance to develop models
which could represent accurately the combustion process of these materials. Modelling the combus-
tion of composite and double-base propellants has been also of high interest for researchers in terms
of optimising chemical composition to fulfil industrial requirements, to improve propellant design
as well as for development and testing activities [8]. Analysing the burning rate of this materials,
which means their decomposition and their subsequent combustion, is one of the main objectives of
combustion modelling [8]. Therefore, a transient multidimensional numerical model to describe the
combustion of composite AP/HTPB and double-base solid propellants is presented in this work. The
model considers the coupling of both condensed and gas phases by using mass, energy and species
conservation, as well as temperature continuity in the burning surface. The kinetics of the model is
described a change of phase of the propellant from condensed to gas and afterwards its combustion
by a simplified chemistry scheme.
23
CHAPTER 1. INTRODUCTION
The combustion of the solid propellant will be performed by considering firstly, a change of phase
of the propellant from condensed to gas and secondly, two simplified chemical reactions to represent
the combustion itself. To couple both processes, several conditions shall be imposed in the burning
surface in which the burning rate will represent a key factor. The results obtained for the combustion
of both, composite and double-base propellants, must compared against experimental test and results
present in the literature. Finally, the developed tool shall confirm its robustness for the prediction of
composite and double propellants combustion behaviour in multidimensional scenarios with transient
environmental conditions such as rocket engines and/or base bleed units in ballistics applications.
24
Chapter 2
Throughout this chapter, a state of the art of the literature related to the combustion of solid
propellant will be presented. Firstly, the classical models to characterise transient combustion of
solid propellants will be summarised briefly. Afterwards, more specific models used in this theses for
the study of granulated high energetic propellants, AP-HTPB composite propellants and double-base
propellants will be described.
@T @T @2T Q̇ Q̇rad
+ rb ↵c = + (2.1)
@t @x @x2 ⇢c C c ⇢c Cc
where the last term corresponds to radiation heat which is usually neglected in the majority of
the studies referred in the revision of Kuo et al. [9]. To solve this equation, initial and boundary
conditions depending on the problem will be needed. To set boundary conditions, it will be necessary
to understand the physics of the burning process at the interface between solid and gas phase which
can be seen in Figures and (Fig. 3 and Fig. 4 from [9]). The boundary conditions at the burning
surface are the mass flux, energy and species continuity.
Gas phase system of equations consists, according to Kuo et al. [9] in a set of equations formed
by continuity, energy and species mass fraction governing equations which can be written as,
@ (↵⇢g ) !
+ r· (↵⇢g !
u g) = 0 (2.2)
@t
✓ ◆ N
!
@h @h @p @ @T @ X @Yi
⇢g + ⇢g !
ug = g + ⇢g hi D i (2.3)
@t @x @t @x @x @x i=1
@x
25
CHAPTER 2. STATE OF THE ART
✓ ◆
@Yi @Yi @ @Yi
⇢g + ⇢g !
ug ⇢g D i = !˙ i (2.4)
@t @x @x @x
where i = 1, 2, 3..., N state for N different species.
The system is formed by equations which are coupled between them. In addition, gas phase
equations are liked to solid phase equation due to the balance at the interface. Therefore, as happened
for solid phase, boundary and initial conditions are needed.
In order to simplify the problem, Kuo et al. [9] described several models which can be classified
into two main groups: quasi-steady flame models and unsteady flame models. quasi-steady flame
models, which are considered the easiest to be applied, are divided into three different approaches
will will be detailed hereafter.
@T p02 @ 2 T
T (p, x) = T0 (p, x)|p=p0 + p0 + + ... (2.6)
@p p=p0 2! @p2 p=p0
✓ ◆ ✓ ◆ ✓ ✓ ◆◆ ✓ ◆
@T0 @p @T0 p0 @ @T @2 0 @T
+ rb + r0 = ↵c 2 T0 + p (2.8)
@p p=p0 @t @x p=p0 2 @x @p p=p0 @x @p p=p0
According to [10], the equation form of the non-steady burning rate can be described as,
✓ ◆
n ↵ dp
r = r0 1 + (2.9)
p r02 dt
where r0 = apn is the steady burning rate expression, n is the exponent of the steady burning
rate expression, ↵ is the solid phase thermal diffusivity and is the dynamic burning rate coefficient,
h n/m
i
1 (p/pm 0)
= ⇣ ⌘n/m (2.10)
p 1
p0 2 + m 2H
where m is the power of propellant pyrolysis, H is the dimensionless magnitude of surface heat
due to propellant decomposition, and p/p0 is the dimensionless pressure.
The main advantage of these models is the obtainment of a closed solution for the burning rate
with the instantaneous pressure, pressurisation rate and propellant properties.
26
2.1. CLASSICAL MODELS TO CHARACTERISE TRANSIENT COMBUSTION OF SOLID
PROPELLANTS
27
CHAPTER 2. STATE OF THE ART
Z1
@T ⇢c r b C g x
g = Qf e ⇢g "dx
˙ (2.11)
@x x=0 g
0
where "˙ is the generation product species and Qf is the heat generation term.
The definition of the generation product species as a value uniformly distributed in the reaction
zone is done by Krier et al. [11] as,
@T g !Qf
⇢c rb Cg
(xf xi )
g = e g (2.13)
@x x=0 ⇢c r b C g
where ! is the average value of the dimensionless frequency and xi is initial coordinate and xf is
flame coordinate.
According to Krier et al. [11], since in the propellant pyrolysis the temperature is the main
parameter of the sublimation process, the pyrolysis law can be written as,
m
r = b (Ts T1 ) (2.14)
where Ts is the surface temperature, T1 is the temperature far from the burning surface and m
and b are the exponent and the burning rate constant respectively in the pyrolysis law.
28
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
29
CHAPTER 2. STATE OF THE ART
@ (↵⇢g ) @
+ (↵⇢g ug ) = c (2.16)
@t @x
✓ ◆
@ (↵⇢g ug ) @ 2 @pg @ 4 @ug 3(1 ↵) Pw ⌧ w
+ ↵⇢g ug = ↵ + µ↵ FD (2.17)
@t @x @x @x 3 @x rp A
@ (↵⇢g Eg ) @
+ (↵⇢g Eg ug ) =
@t @x ✓ ◆
@ @ @T 2qt 3(1 ↵)
(pg ↵ug ) + k↵ + [⇢p ṙhf hc (Tg Tp )] (2.18)
@x @x @x rt rp
where Pw is wetted perimeter and ⌧w is the shear stress on the tube wall, which will be neglected.
The mass transfer c is defined as,
3(1 ↵)
c = ṙ⇢p (2.19)
rp
where ṙis the combustion law,
ṙ = apng (2.20)
The drag force FD for a non-fluidized bed of spheres as,
⇢
⇢g ug |ug | 150 (1 ↵)
FD = + 1.75 (2.21)
6↵ Rep
where Rep is the Reynolds number based on particle diameter. Finally, the heat transfer correla-
tion hc is expressed as,
✓ ◆0.7
kg ↵⇢g |ug | rp
hc = 0.65 P r0.33 (2.22)
rp µg
To solve the problem, Kuo et al. [12] assumed neglected the heat loss to the tube wall, the
friction force between the gas and the tube wall, the work of viscous stress from the velocity gradient,
dissipative heat, molecular heat conduction from gas to gas, imperfection of gas equation of state,
dependence of the temperature from the specific heat, dynamic burning effect and erosive burning
effect.
@ (↵⇢g ) @
+ (↵⇢g ug ) = c (2.23)
@t @x
@ (↵⇢g ug ) @ @pg @
+ ↵⇢g u2g = ↵ + (↵⌧xx ) FD As (2.24)
@t @x @x @x
30
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
@ (↵⇢g Eg ) @
+ (↵⇢g Eg ug ) =
@t @x ✓ ◆
@ @ @ @T
(pg ↵ug ) + (↵⌧xx ug ) + k↵ + As [⇢p ṙhf hc (Tg Tps )] (2.25)
@x @x @x @x
where ⌧xx is the shear stress in the normal direction. The mass transfer c is defined as,
c = As ṙ⇢p (2.26)
where
3(1 ↵)
As = (2.27)
rp
and the combustion law is that of equation (2.20). The remaining source terms are considered
the same as in Kuo et al. [12].
@ (↵⇢g ) @
+ (↵⇢g ug ) = g (2.28)
@t @x
@ug @ug @ (↵pg ) @ ⇣ 2
⌘
↵⇢g + ug = + ↵⇢g (ug um ) g (ug u m ) + FD (2.29)
@t @x @x @x
@Eg @Eg @ ⇣ ⇣ 2
⌘⌘
g
↵⇢g + ug = ug ↵pg + ↵⇢g (ug um ) + g (Echem Eg )
@t @x @x
@ @ h ⇣
2
⌘i
+ (↵⇢g Eg (ug um )) Eg (ug um ) ↵pg + ↵⇢g (ug um )
@x @x
N ukg
+ ↵ 2 (Tg Tp ) (2.30)
dp
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢g up ) = p (2.31)
@t @x
@up @up @ ⇣ 2
⌘
(1 ↵) ⇢p + up = (1 ↵) ⇢p (up um ) p (up um ) FD (2.32)
@t @x @x
@Ep @Ep @ ⇣ ⇣ 2
⌘⌘
p
(1 ↵) ⇢p + up = up (1 ↵) ⇢p (up um ) + p (Echem Ep )
@t @x @x
@ @ h ⇣
2
⌘i
+ ((1 ↵) ⇢p Ep (up um )) (up um ) (1 ↵) ⇢p (ug um )
@x @x
N ukg
+ (1 ↵) (Tg Tp ) (2.33)
d2p
where the density and the average velocity are described as,
⇢p (1 ↵) up + ↵⇢g ug = ⇢m um (2.34)
31
CHAPTER 2. STATE OF THE ART
⇢p (1 ↵) + ↵⇢g = ⇢m (2.35)
and the mass transfer,
p = g (2.36)
The total energy is written as,
ṙ = b + apng (2.39)
In order to close the equations heat transfer and drag force are needed. Van Tassel and Krier
presented in their work a summary of several friction coefficient and Nusselt number correlations.
The gas is assumed to be ideal therefore, pressure may be determined by means of,
p g = ⇢g R g T g (2.40)
This model is afterwards modified by Krier et al. [22] to analyse the flame-spreading and com-
bustion of small solid propellant grains. The problem could be described as the combustion of solid
propellant in a cylindrical chamber with the initiator in the end.
@ (↵⇢g ) 1 @
+ (A↵⇢g ug ) = As ⇢p ṙ (2.41)
@t A @x
@ (↵⇢g ug ) 1 @ g @ (↵pg A) g @
+ A↵⇢g u2g + (A↵⌧xx ) =
@t A @x A @x A @x
gpg @A
gFD As + As ṙ⇢p up + (2.42)
A @x
@ (↵⇢g Eg ) 1 @ 1 @ 1 @ 1 @
+ (A↵⇢g Eg ug ) + (Apg ↵ug ) = + (↵⌧xx ug A) (qA↵)
@t A @x " AJ
# @x A @x A @x
Up2 As u p F D pg @(↵A)
+ As ⇢p ṙ hchem + As ht (T Tps ) Q̇w
2gJ J AJ @t
(2.43)
@ ((1 ↵) ⇢p ) 1 @
+ ((1 ↵) ⇢g up A) = As ⇢p ṙ (2.44)
@t A @x
32
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
@ ((1 ↵) ⇢p up ) 1 @ g @
+ A (1 ↵) ⇢p u2p (A (1 ↵) ⌧xx ) =
@t A @x A @x
⌧wp Pwp
As ṙ⇢p up + gFD As g (2.45)
Q
where As is the specific surface of the granular propellant, A is the cross section of the cylindrical
combustion chamber, J is the heat to work conversion factor, ht is the total heat transfer coefficient,
hchem is the enthalpy of the propellant gas at flame temperature, q is the rate of conduction heat
transfer per unit area Tps is the particle surface temperature, Q̇w is the rate of heat loss to the
chamber wall, ⌧wp is the shear stress between the camber wall and particle phase and Pwp is the
wetted perimeter between the chamber wall and the particle phase.
For the particle energy governing equation, instead of using the same formulation as for the gas
phase, they use the heat transfer equation as,
✓ ◆
DTp ↵p @ 2 (rTp )
= (2.46)
Dt p r @r2
To solve this system of equation the solid is assumed to be incompressible and the gas is modelled
by the equation of Noble-Abel as,
⇢g ( g 1)eg
pg = (2.47)
1 ⇢g ⌘g
A very interesting difference respect the previous models is the usage of the expression of Lenoir
and Robbillard [24] as,
✓ ◆
ṙ⇢p
ṙ = apng + K e hc e ⇢g |ug up |
(2.48)
where hc is the heat transfer coefficient, Ke is the erosive burning constant and is the erosive
burning exponent.
To solve the equation system (2.41)-(2.46) heat transfer and drag force expressions are needed
and defined as,
µg (ug up ) (1 ↵)
FD = fpg (2.49)
12grp ↵
where the friction coefficient is calculated as,
⇢ ⇣ ⌘0.87
Re
fpg = 276.23 + 5.05 1 ↵p
(2.50)
Rep
1< 1 ↵ < 24000
and,
✓ ◆0.7
kg ↵⇢g |ug up | rp
hc = 0.65 P r1/3 (2.51)
rp µg
valid for Reynolds up to 50000 and porosity around 0.37.
33
CHAPTER 2. STATE OF THE ART
@ (↵⇢g ) @
+ (↵⇢g ug ) = g (2.52)
@t @x
@ (↵⇢g ug ) @ @↵
+ ↵ ⇢g u2g + pg = pg + g (ug up ) FD (2.53)
@t @x @x
@ (↵⇢g Eg ) @ g
+ (↵ug (⇢g Eg + pg )) = g (Echem Eg ) Q̇ (2.54)
@t @x
@ ((1↵) ⇢p ) @
+ ((1 ↵) ⇢g up ) = p (2.55)
@t @x
@ ((1 ↵) ⇢p up ) @
+ (1 ↵) ⇢p u2p = +FD (2.56)
@t @x
@ ((1
↵) ⇢p Ep ) @ p
+ ((1 ↵) up ⇢p Ep ) = p (Echem Ep ) + Q̇ (2.57)
@t @x
where the mass transfer of gas and solid phases is,
p = g (2.58)
and it is calculated as equation (2.19). The solid is assumed to be incompressible and the gas is
model by the equation of Noble-Abel (equation (2.47)). Burning rate expression used by the authors
is the already mentioned equation (2.20). heat transfer and drag force are calculated as,
kg
hc = Nu (2.59)
dp
µg (ug up )
FD = ff g (2.60)
d2p
where the Nusselt number and the friction coefficient are defined with Denton [25] and Ergun [26]
correlations respectively as,
N u = 0.58Re0.7
p Pr
0.33
(2.61)
2 ⇢
(1 ↵) Rep
fpg = 150 + 1.75 (2.62)
↵2 1 ↵
Finally, to solve the problem, Beckstead et al. [4] considered that the control volume has a small
portion already initiated, which means that it has high pressure and temperature.
@ (↵⇢g ) @
+ (↵⇢g ug ) = g (2.63)
@t @x
@ug @ug @ (↵pg ) @ ⇣ 2
⌘
↵⇢g + ug = + ↵⇢g (ug um ) g (ug um ) FD (2.64)
@t @x @x @x
34
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
@eg @eg @
↵⇢g + ug = (ug ↵pg ) + g (eg egchem )
@t @x @x
@ 3(1 ↵)
+ (↵⇢g eg (ug um )) h (eg ep ) (2.65)
@x rp c v
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢g up ) = p (2.66)
@t @x
@up @up
(1 ↵) ⇢p + up =
@t @x
@ ((1 ↵) ⌧p ) @ ⇣ 2
⌘
+ (1 ↵) ⇢p (up um ) + (up u m ) + FD (2.67)
@x @x
@ep @ep @
(1 ↵) ⇢p + up = (um (1 ↵) ⌧p ) g (epchem ep )
@t @x @x
@ 3(1 ↵)
+ ((1 ↵) ⇢p ep (up um )) + h (eg ep ) (2.68)
@x rp C v g
where the mass transfer g has the expression of equation (2.19), is the relation of specific heats,
and eg = Cvg Tg .
In this case, the relation used for the burning rate is similar to (2.39) and can be written as,
⌧p = p g + p s (2.70)
where ps is defined by the following empirical equation,
@ps @(1 ↵)
(1 = mB 0 (1 ↵)m
↵) (2.71)
@x @x
where m and B are empirical coefficients.
The expressions used in this case for the heat transfer and drag force are the same ones used by
Beckstead et al. [4].
@ (↵⇢g ) !
+ r· (↵⇢g !
u g) = c + ign (2.72)
@t
@ (↵⇢g !
u g) ! ⇣ !⌘ ! !
+ r· ↵⇢g !
ug ⌦! u g + ↵p I + pg r (1 ↵) = ign
!
u ign + c
!
up FD (2.73)
@t
35
CHAPTER 2. STATE OF THE ART
✓ ✓ ◆◆
@ (↵⇢g Eg ) ! ! p !
+ r· ↵⇢g u g Eg + + pg r [(1 ↵) !
u p] =
@t ⇢g
✓ ◆
pg |!
u p| 2 ! !
c Qex + + u pF D Qign ign (2.74)
⇢p 2
@ ((1 ↵) ⇢p ) !
+ r· ((1 ↵) ⇢p !
u p) = c (2.75)
@t
@ ((1 ↵) ⇢p !
u p) ! ⇣ !⌘ ! !
+ r· (1 ↵) ⇢p !
up⌦!
u p + (1 ↵) pp I pg r (1 ↵) = c
!
u p + F D (2.76)
@t
@Hs ! !
+ u p · r· (Hs ) = q˙t (2.77)
@t
@dq ! !
+ u p · r· (dq ) = ṙ (2.78)
@t
where is the thermal diffusivity of the solid phase, q̇t is the interfacial heat flux, Qex is the
combustion heat, ign , Qign , ! u ign , mass transfer, heat and velocity of ignition respectively, dp is
the average distance of all propellant grains which is burned due to combustion and hs is specific
enthalphy of the propellant grains and it is defined as,
hs = cv,p Tp (2.79)
being cv,p is the specific heat at constant volume for the solid phase.
Energy for gas and solid phase follow equations (2.37) and (2.38). Equations (2.20) and (2.47)
are used for the evaluation of the burning rate and the characterisation of the thermodynamic state
of the gas phase.
@ (↵⇢g ) @
+ (↵⇢g ug ) = g (2.80)
@t @x
@ug @ug @ (↵pg )
↵⇢g + ug = g (ug up ) FD (2.81)
@t @x @x
@Eg @Eg @
↵⇢g + ug = (ug ↵pg ) g Eg Egchem
@t @x @x
!
u2g u2p
+ g + ug up + FD (ug up ) (2.82)
2 2
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢g up ) = g (2.83)
@t @x
36
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
@up @up @ ((1 ↵) ⌧p )
(1 ↵) ⇢p + up = + FD (2.84)
@t @x @x
@Ep @Ep @⌧p
(1 ↵) ⇢p + up = up (1 ↵) + g Ep + Epchem + (2.85)
@t @x @x
where Eg = Cvg Tg and Ep = Cvp Tp . ⌧p is defined as the resistance that the solid does to avoid
compaction and it has the following form,
✓ ◆
K 1 1
⌧p = (2.86)
(1 ↵) (1 ↵c ) (1 ↵)
where K is the bulk modulus for the propellant.
is the interface heat transfer defined as,
6(1 ↵)
= hc (Tg Tp ) (2.87)
dp
k
where hc = dpg N u is the heat transfer coefficient and Nusselt is calculated according to Denton
[25] equation (2.61).
Drag force FD is calculated as,
2
µg (ug up ) (1 ↵)
FD = ff g (2.88)
d2p ↵2
where ff g is the friction coefficient calculated with Ergun [26] expression (2.62) and Reynolds
number is calculated as,
dp ↵⇢g |ug up |
Rep = (2.89)
µg
To finalise with the source terms used by Krier and Kezerle mass transfer expression needs to be
detailed. In this case, equation (2.19) is used and the burning rate is defined as,
✓ ◆m ✓ ◆n
Tp pg
ṙ = b (2.90)
Tp0 p0
where m, n and b are assumed known constants and Tp0 and p0 are initial temperature of the
solid phase and pressure respectively.
@ (↵⇢g ) @
+ (↵⇢g ug ) = (2.91)
@t @x
@ (↵⇢g ug ) @ @pg
+ ↵⇢g u2g = ↵ + up FD (2.92)
@t @x @x
!
@ (↵⇢g Eg ) @ u2p
+ (↵ug (⇢g Eg + pg )) = Egchem + FD u p Q̇ (2.93)
@t @x 2
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢p up ) = (2.94)
@t @x
37
CHAPTER 2. STATE OF THE ART
@ ((1 ↵) ⇢p up ) @ @pp
+ (1 ↵) ⇢p u2p = (1 ↵) u p + FD (2.95)
@t @x @x
!
@ ((1 ↵) ⇢p Ep ) @ u2p
+ ((1 ↵) up (⇢p Ep + pp )) = Epchem + FD up + Q̇ (2.96)
@t @x 2
where gas and solid energies are defined with equations (2.37) and (2.38). Noble-Abel equation
of state for gas is assumed and pp is defined as the particle-particle stress as,
( ⇣ ⌘
K 1 1
↵ ↵c
pp = (1 ↵) (1 ↵c ) (1 ↵)
(2.97)
0 ↵ > ↵c
where ↵c is the critical porosity above which the particles do not touch and K is the particle
stress proportionality constant. This expression is similar to the one from Krier and Kezerle [27] but
considering zero the value of particle stress in case of having porosity values under the critical one.
Remaining source terms are mass transfer, drag and convective heat transfer. Mass transfer has
the same expression detailed for previous authors (2.19) but burning rate is defined as,
✓ ◆m
Tp n
ṙ = (bpg ) (2.98)
Tp0
Drag force is calculated as equation where the friction coefficient has the expression defined by
Kuo and Nydegger [30] written as,
( ✓ ◆0.87 )
Rep
fpg = 276 + 5.05 (2.99)
1 ↵
valid for Reynolds up to 1500. Reynolds is defined with equation (2.89). Gas viscosity is defined
in this case as,
N u = 0.65Re0.65
p P r0.33 (2.101)
One interesting aspect of the work provided by Hoffman and Krier [28], which is not considered
by other authors, is the problem of reaching a porosity under the minimum physical values. In their
work, Hoffman and Krier [28] remarked that the set of equations (2.91) - (2.96) does not prevent
compaction below the minimum porosity that for spheric particles is calculated to ↵min = 0.2595
(see Annex of their work). The approach used to solve this issue consisted in introducing a function
which depends on the porosity in the momentum equation of the solid phase as follows,
@ ((1 ↵) ⇢p !
u p) @ @pg
+ (1 ↵) ⇢p u2p = up + (1 f (↵)) FD (1 ↵) (2.102)
@t @x @x
38
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
@ (↵⇢g ) @
+ (↵⇢g ug ) = (2.104)
@t @x
@ (↵⇢g ug ) @ @ (↵pg ) @ ⇣ 2
⌘
+ ↵⇢g u2g = + u m FD + ↵⇢g (ug um ) (2.105)
@t @x @x @x
@ (↵⇢g Eg ) @ @ @
+ (⇢g Eg ↵ug ) = (um ↵pg ) + (↵⇢g Eg (ug um )) + Egchem FD u m Q̇ (2.106)
@t @x @x @x
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢g up ) = (2.107)
@t @x
@ ((1 ↵) ⇢p up ) @
+ (1 ↵) ⇢p u2p =
@t @x
@ (⌧p (1 ↵)) @ ⇣ 2
⌘
u m + FD + (1 ↵) ⇢p (up um ) (2.108)
@x @x
@ ((1 ↵) ⇢p Ep ) @
+ ((1 ↵) up ⇢p Ep ) =
@t @x
@ @
(um (1 ↵) ⌧p ) + ((1 ↵) ⇢p Ep (up um )) + Epchem + FD um + Q̇ (2.109)
@x @x
where the definition of the average velocity, um , is that introduced in equation (2.34).
On the other hand, the theory of the separated flow is represented by the system of partial
differential equations (2.91)-(2.109) already used by Hoffman and Krier [28].
Regarding the source terms, Denton [25] equation (2.61) is used for the expression of the Nusselt
number and Ergun [26] and Kuo and Nydegger [30] expressions for friction coefficients are used
depending on the Reynolds number. In addition, as in the work done by Hoffman and Krier [28], not
only Noble-Abel equation of state of gases but also same expression for the particle-particle stressed
and mass transfer and are used.
@ (↵⇢g ) @
+ (↵⇢g ug ) = (2.110)
@t @x
@ (↵⇢g ug ) @
+ ↵⇢g u2g + ↵pg = up FD (2.111)
@t @x
39
CHAPTER 2. STATE OF THE ART
!
@ (↵⇢g Eg ) @ chem
u2p
+ (↵ug (⇢g Eg + pg )) = E + FD u p Q̇ (2.112)
@t @x 2
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢p up ) = (2.113)
@t @x
@ ((1 ↵) ⇢p up ) @
+ (1 ↵) ⇢p u2p + (1 ↵) pp = FD up (2.114)
@t @x
!
@ ((1 ↵) ⇢p Ep ) @ u2p
+ ((1 ↵) up (⇢p Ep + pp )) = + FD up + Q̇ (2.115)
@t @x 2
N u = 0.65Re0.65
p P r0.33 (2.117)
Finally pp is called the solid phase stressed and it is expressed as a function of the gas pressure
and equilibrated stress,
pp = pg + pe (2.118)
The equilibrated stress is function of the porosity and it has the following expression depending
of the phase of compaction = 1/(1 ↵),
8 4G0 ( 0 )
>
> 3 (n 1) elastic phase ( 0 > 1)
< h io
e 2 2G0 2G0 ( )
p = 3Y 1 Y ( 0 ) + ln Y(
0
1) elastic plastic phase ( 1 > 2) (2.119)
>
> ⇣ ⌘
: 2 Y ln plastic phase ( > 1.0)
3 ( 1) 2
where Y and G0 are the yield strength and the shear modulus respectively and their values are
found in Table 1 of their work and 1 and 2 have the following form,
2G0 0 + Y
1 = (2.120)
2G0 + Y
2G0 0
2 = (2.121)
2G0 + Y
and 0 = (1 ↵0 ) where ↵0 correspond to initial porosity.
40
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
@ (↵⇢g ) @
+ (↵⇢g ug ) = (2.122)
@t @x
@ (↵⇢g ug ) @
+ ↵⇢g u2g + ↵pg = up FD (2.123)
@t @x
@ (↵⇢g Eg ) @
+ (↵ug (⇢g Eg + pg )) = Ep FD u p Q̇ (2.124)
@t @x
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢p up ) = (2.125)
@t @x
@ ((1 ↵) ⇢p up ) @
+ (1 ↵) ⇢p u2p + (1 ↵) pp = FD up (2.126)
@t @x
@ ((1 ↵) ⇢p Ep ) @
+ ((1 ↵) up (⇢p Ep + pp )) = Ep + FD up + Q̇ (2.127)
@t @x
⇥ ⇤
@ (1 ↵) /r3 @ ⇥ ⇤
+ (1 ↵) up /r3 = 0 (2.128)
@t @x
@ (1 ↵) @ (1 ↵) ↵ pp0 pg0 3(1 ↵)ṙ
+ up (1 ↵) = pp pg (1 ↵) (2.129)
@t @x µc 1 ↵0 r
where µc is the compaction viscosity and r is the solid particle radius. Mass transfer source term is
expressed as equation (2.19) where the burning rate is equation (2.20) and Noble-Abel state equation
is used for gas phase. Particle stress has the following expression,
⇢p0
pp = ( p 1)cvp ⇢p Tp (2.130)
p
where p is the Tait equation parameter, the subscript ‘0’ indicates undisturbed conditions and
is the non-ideal solid parameter.
Mass transfer, drag force and heat transfer expressions are needed to solve the system of equations
being,
3(1 ↵)
= apm
g ⇢p (2.131)
r
(1 ↵)↵
FD = (ug up ) (2.132)
r
(1 ↵)↵
Q̇ = (Tg Tp )h (2.133)
r1/3
where a and m are burning constant and burn index respectively. and h coefficients are defined
as positive constants values. Note that in this model, the particle radius does not have a constant
41
CHAPTER 2. STATE OF THE ART
value as it varies while the propellant is burning. In this model, Powers defines the gas and particle
total energies with the following expressions,
u2g
Eg = c vg T g + (2.134)
2
u2p ⇢p0
Ep = c v p T p + + +q (2.135)
2 p
@ (↵⇢g ) @
+ (↵⇢g ug ) = (2.136)
@t @x
@ (↵⇢g ug ) @ @↵
+ ↵⇢g u2g + ↵pg = up FD + a pg (2.137)
@t @x @x
@ (↵⇢g Eg ) @ @↵
+ (↵ug (⇢g Eg + pg )) = Ep FD u p Q̇ + a u p pg (2.138)
@t @x @x
@ ((1 ↵) ⇢p ) @
+ ((1 ↵) ⇢p up ) = (2.139)
@t @x
@ ((1 ↵) ⇢p up ) @ @↵
+ (1 ↵) ⇢p u2p + (1 ↵) pp = FD up a pg (2.140)
@t @x @x
@ ((1 ↵) ⇢p Ep ) @
+ ((1 ↵) up (⇢p Ep + pp )) =
@t @x
@↵ (1 ↵) ↵
Ep + FD up + Q̇ a u p pg b (pp fp ) [pp fp (pg fg )] (2.141)
@x µc
@ (1 ↵) @ (1 ↵) ↵ 3(1 ↵)ṙ
+ up (1 ↵) = [pp fp (pg fg )] (2.142)
@t @x µc r
@n @
(1 c) + [nup ] = (1 c )F (pg , ⇢g , pp , ⇢p , (1 ↵)) (2.143)
@t @x
where n = 3(1 ↵)
4⇡r 3 , parameter is used to determine if the formulation corresponds to the model
of Baer Nunziato [34] or the one from Powers, Stewart and Krier [37]and fg and fp are the stress of
gas and solid being,
fg = 0 (2.144)
pp0 pg0
fp = (1 ↵) (2.145)
1 ↵0
respectively.
Solid pressure has the same expression as equation (2.130) and gas pressure is modelled by using
Noble-Abel equation.
Drag force and heat transfer are defined as,
42
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
FD = B(ug up ) (2.146)
Q̇ = C(Tg Tp ) (2.147)
where B and C are strictly positives function of the flow variables defined as,
(1 ↵)↵
B= (2.148)
r
(1 ↵)↵
C= h (2.149)
r1/3
Mass transfer was written as,
3(1 ↵)
= apm
g ⇢p H(Tp Tign ) (2.150)
r
where H(Tp Tign ) is a step function included in the combustion model to prevent combustion
until an ignition temperature is reached. Same expressions of gas and solid energy than the ones used
by Powers in his dissertation [33] are used in this model.
6(1 ↵)
Q̇ = hc (Tg Tp ) (2.153)
dp
The convective heat transfer coefficient depends on the Nusselt number in the following manner,
kg
hc = Nu (2.154)
dp
where the Prandtl number can be defined as,
cpg µg
Pr = (2.155)
kg
and gas conductivity as,
✓ ◆
15Ru 4Mg cvg 3
kg = µg + (2.156)
4Mg 15Ru 5
43
CHAPTER 2. STATE OF THE ART
The difference between the convective heat transfer expressions found in the literature comes from
the definition that each author did of the Nusselt number. In this review, three expression will be
considered and analysed calling them Nusselt 1, Nusselt 2 and Nusselt 3 to simplify the reference
afterwards. Nusselt 1 is the expression proposed by Gelperin and Einstein [38] for fluidised beds that
appeared in the work by van Tassel and Krier [21] and Krier et al. [22] as,
⇣ ⌘
N u = 2 + 0.4Re2/3
p Pr
1/3
(2.157)
The second one (called hereafter Nusselt 2) was used by many authors in the existing literature,
among them, Beckstead et al. [4], Krier and Gokhale [7], Krier and Kezerle [27] and Gokhale and
Krier [31],
N u = 0.58Re0.7
p Pr
0.33
(2.158)
To conclude, the last correlation, called from now on Nusselt 3, that used by Hoffman and Krier
[28, 29], Butler et al. [39] and Butler and Krier [32],
N u = 0.65Re0.65
p P r0.33 (2.159)
On the other hand, the expression of the drag force considered in the bibliography is the following,
2
µg (ug up ) (1 ↵)
FD = ff g (2.160)
d2p ↵2
where ff g is the friction coefficient. Depending on the author, the expression of the friction
coefficient changes its form.
Beckstead et al. [4], Krier and Gokhale [7] and Gough [16] used the following form of the friction
coefficient which will be called hereafter “Friction 1”,
⇢
Rep
fpg = 150 + 1.75 (2.161)
1 ↵
According to Gokhale and Krier [31] this correlation is only valid when Reynolds number has a
value below 1500. Therefore, they proposed an alternative expression valid for values of Reynolds
number up to 15000 which has been also used by Hoffman and Krier [28, 29] and will be called from
now on “Friction 2”.
( ✓ ◆0.87 )
Rep
fpg = 276 + 5.05 (2.162)
1 ↵
To conclude, Butler et al. [39] presented in his work an expression, called hereafter “Friction
3”, developed by Wilcox, Jones and Krier [40] which provide them good results for high Reynolds
numbers,
⇢ ⇣ ⌘0.88
Re
fpg = 150 + 3.89 1 ↵p
(2.163)
103 < Rep < 2 · 105
44
2.2. MULTIDIMENSIONAL MODELS TO CHARACTERISE COMBUSTION OF
GRANULATED PROPELLANTS
Beckstead et al. [4] proposed a slope initiation in which the 10% of the tube is considered
initiated region. From this length , the temperature and pressure decrease until reaching 1 bar
and 295K conditions. Krier and Gokhale [7] however, assumed that the ignited region was 15%
of the total length. Krier and Kezerle [27] proposed to have 10% of the domain with high values
of pressure and gas temperature and 3% of the domain with particle temperature 100K above the
ignition temperature. Hoffman and Krier [28] studied two initiation conditions. On the first hand,
one similar to the one of Krier and Kezerle [27] but changing the ignition temperature. On the other
hand, they assume an exponential distribution of the gas pressure and a linear distribution of gas
and particle temperature during 9% of the domain. Finally, the first initiation of Hoffman and Krier
[28] has been reproduced considering the control volume at 1 bar pressure.
The values of pressure, gas and particle temperature for all initiation can be seen in Figure 2.4.
To be able to refer to each initiation type easily, each one of them has been named as Test and a
number.
Table 2.1: Pressure, gas and particle temperature and ignition temperature data for the considered
initiation.
Figure 2.4: Initial conditions initiation of pressure, particle temperature and gas temperature.
45
CHAPTER 2. STATE OF THE ART
46
2.3. MODELS TO CHARACTERISE COMPOSITE PROPELLANTS COMBUSTION
special attention may be drawn to the work of Moriñigo et al. [50] in which the combustion dynamics
of HTPB with oxygen in a wind tunnel is studied by comparing qualitatively the experimental data
obtained with numerical simulations performed with LES (Large Eddy Simulation). In that study,
they propose a reduced chemical model with 2 steps and another model with 6 steps. In addition, Ko-
hga [51] provides experimental results of burning AP/HTPB composite propellant with and without
the use of additives for low and high pressures and AP content from 20% to 80% and his results are
afterwards used by Cao et al. [47]. Ramakrishna et al. [52] modelled the behaviour of sandwich
propellant combustion and validate their results against the experimental ones provided by by Price
[53] and Price et al. [54]. Experimental results in gas temperature for AP/HTPB are obtained by
Fitzgerald and Brewster [55] who made images the combustion of AP/HTPB propellant by using
infrared and ultraviolet emission simultaneously.
Finally, it is worth remarking that, in case of AP/HTPB propellants, several authors as Miccio [56]
and Favale and Miccio [57] have pointed out that it is mandatory to consider the two-dimensional
(2D) and/or the three-dimensional (3D) aspects of the problem due the involvement of diffusion
flames. This remark was also done by other authors such Knott and Brewster [58] which also develop
a two-dimensional model for the combustion of heterogeneous propellants. In addition, Zhou et al.
[59] used the model developed on previous works from Jackson and Buckmaster [60] to model the
combustion of particles of AP embedded in a binder as done by Favale and Miccio [57]. Therefore, it
seems of major importance to develop multidimensional models to study the burning of this type of
propellants, especially in the case of transient combustion.
To end this subsection, a detailed review of the models used from the literature to develop the
one established in this work is performed afterwards.
The model is defined with a set of equations for both phases coupled with the conditions at
47
CHAPTER 2. STATE OF THE ART
the boundary. Assuming ideal gas behaviour and constant thermophysical properties, the partial
differential system of equations for the gas phase can be written as,
@Yi,g
ṁ00 = ⇢g Dg r2 Yi,g + Yi ,g (2.164)
@y
@Tg
ṁ00 cp,g = k g r2 Tg + T,g (2.165)
@y
being equations (2.164) and (2.165) species and energy equations respectively where, i stands for
the i-gas species, ṁ00 is the mass flux, kg is the thermal conductivity of the gas, Dg is gas diffusivity,
Yi ,g is the sum of source terms corresponding to species generation or consumption and T,g is the
source term related to the heat of reaction being both depending of Arrhenius reaction kinetics. Zero
gradient boundary conditions are assumed for temperature and species variables in x direction and
fixed in y direction.
Solid propellant is modelled with an equation for energy and an equation for species such as,
@Tc
ṁ00 cp,c = k c r2 Tc + T,c (2.166)
@y
@Yc
ṁ00 = Yc (2.167)
@y
where Yc is the mass fraction of the condensed phase, kc is the thermal conductivity of the
condensed phase, Yc is the sum of source terms corresponding to species generation or consumption
of the condensed phase and T,c is the source term related to the energy release in the condensed
phase. Zero gradient boundary conditions are assumed for temperature in x direction and fixed to
initial temperature in y direction. For condensed phase concentration the boundaries assumed for
the wall are considered value equal to 0 and 1 in y direction.
To couple both systems, boundary conditions at the burning surface need to be established guar-
antying mass and energy conservation as temperature continuity,
Tg | s = Tc | s (2.169)
48
2.3. MODELS TO CHARACTERISE COMPOSITE PROPELLANTS COMBUSTION
⇣ ⌘
E2
2
R2 = A2 p YHT P B YO2 e RuTg
(2.174)
where A is the reaction rate coefficient and E is the activation energy of reactions 1 and 2.
And the burning rate of the condensed phase is define as,
Ec
Rsolid = ⇢c Be RT (2.175)
where Ec is the activation energy of the solid and B is the burning coefficient.
The model is defined as a system of differential equations for the gas and the condensed phase
being,
@Q @ @
+ (E E v ) + (F F v ) = S (2.176)
@t @x @y
where x e y are the axial and radial coordinates, the vector Q is defined as,
T
Q = y [⇢, ⇢u, ⇢v, ⇢e, ⇢Yi ] (2.177)
being = [1, 0] depending if the problem is considered axisymmetric or two-dimensional, S the
source terms, E and F are the flux vectors and E v andF v the diffusive flux vectors.
According to the authors, when treating the combustion of AP/HTPB propellant, several facts
need to be considered, firstly the mass loading of AP is much higher than HTPB, secondly the AP is
more reactive than HTPB and can provide exothermic reaction without the presence of HTPB and
the last one, the size of AP particles have a decisive role in the burning behaviour. Therefore, in the
combustion of the composite propellant, the authors considered AP degradation as the controlling
factor and HTPB burning is collaborating only in the degradation into gas products in the sublimation
process. As a result, only the AP will be considered in the chemical kinetics,
49
CHAPTER 2. STATE OF THE ART
50
2.3. MODELS TO CHARACTERISE COMPOSITE PROPELLANTS COMBUSTION
✓✓ ◆ ◆
@ (Tg )
+ ⇢g ( !
g
⇢g q ·r) T = r rT + T (2.189)
@t cg
@ (Yi )
⇢g + ⇢g ( !
q ·r) Yi = r (⇢g DrYi ) + Yi i = 1, 2, ..., n 1 (2.190)
@t
@ (Tc ) @ (T )
⇢c c c + ṁ = r ( c rT ) (2.191)
@t @y
where c could state for AP or HTPB, ! q = (q1 , q2 ) is the velocity vector and T and C are the
energy and species source respectively. Gas phase is considered to be ideal.
The reaction mechanism is very similar to that of Knott and Brewster [44] being,
51
CHAPTER 2. STATE OF THE ART
! !
g rT = c rT ṁQs (2.197)
y=0+ y=0
ṁ = ⇢g (ug + rb ) = ⇢c rb (2.198)
@
ṁYs,0 = ṁYi |y=0+ + D Yi (2.199)
@y y=0+
where, D is the diffusivity, Qs is he surface heat release and the burning rate is defined as,
Eb
r b = Ab e R u T s (2.200)
where b can be binder or oxidiser depending on the zone considered.
52
2.4. MODELS TO CHARACTERISE DOUBLE-BASE PROPELLANTS COMBUSTION
@ (Tc )
⇢c c = r ( c r·Tc ) (2.205)
@t
where Qg is the heat of reaction of R1 or R2 depending on the subscript, subscript c stands for
condensed phase which could be AP and HTPB depending on the zone considered, Sm is the mass
source term, S!
u g is the momentum source term which can be defined as,
ṁAcell ⇢c ṙAcell
Sm = = (2.206)
Ccell Vcell
S! ! (2.207)
u g = Sm u g
Energy, species and mass conservation as well as temperature continuity are established as bound-
ary conditions in the burning surface as follows,
Tg = Tc (2.208)
ṁ = ⇢g (|!
u g | + ṙ) = ⇢c ṙ (2.209)
c @Yi
ṁ Yi |0+ = ṁ Yi |0 (2.210)
cp @y 0+
@Tg @Tc
g = c = ⇢c ṙQc (2.211)
@y 0+ @y 0+
where Qc is the heat release of phase change. As for the the far field boundaries considered are
the same as in other works from previous authors.
The chemical kinetics of the combustion process is exactly the same defined by Ye et al. [46] in
equations (2.192) and (2.193) and the reaction rates of R1 and R2 as equations (2.194) and (2.195).
The burning rate, however will be calculated as the product of the burning rates calculated for
each propellant, AP or HTPB being,
1 ↵ ↵
ṙ = ṙAP ṙHT P B (2.212)
where ṙAP and ṙHT P B are the burning rates of AP or HTPB respectively which are calculated
with equation (2.200).
53
CHAPTER 2. STATE OF THE ART
rocket motor but also to study the effect that the turbulence has in the flame structure. The model of
equations is based in that of Roh et al. [70, 66]. Wu et al. [71] developed an aerothermochemical or
erosive model to reproduce erosive combustion. In their work the authors are focused on detailing the
combustion modelling considering the chemical kinetics and the diffusive effects of the gaseous phase.
Tseng and Yang [72] studied the combustion of homogeneous double-base propellant developing a
model and solving it using k ✏ turbulence. In their results they observed that the turbulence seemed
not to have a relevant effect in chemistry combustion but in friction forces. Closing, Arkhipov et al.
[73] presented several methods to describe the results of experimental studies of erosive combustion
of solid propellants under transonic and supersonic flux.
Regarding experimental test, it is relevant to highlight the test performed by several authors such
Robbins and Keys [74] who made a review of experimental results for propellants which contained 12%
or 13% of NC quantifying the effect that this mass fraction, as well as the habitual additives, have in
the burning rate. In addition, Long et al. [75, 76] made an experimental study about the combustion
of several double-base propellants containing nitrogen heterocyclic nitroamines (RDX, TNAD, HMX
and DNP). The experimental apparatus consists in a chimney strand burner with several observation
windows. Aoki and Kubuta [77] also used a chimney-type strand burner with four transparent widows
to analyse the burning rate of twelve different kinds of propellant previously prepared with different
compositions. Finally, Markov et al. [78] studied the combustion of a double-base propellant in a
semi-closed volume. They also presented a short theoretical model in order to interpret and define
several experiments performed in steady and non-steady conditions.
Finally, the work done by the authors from literature which have been used to develop the model
to represent combustion of double-base propellant in this thesis will be described n the following
subsection. The final model itself will be afterwards explained in detail in subsection 4.3.3.1 from
Chapter 4.
54
2.4. MODELS TO CHARACTERISE DOUBLE-BASE PROPELLANTS COMBUSTION
In this model, Wu et al. [71] considered N O2 as oxidiser (O) and CH2 O as fuel (F). In addition,
they consider as a result of the surface sublimation, another group of species that they called delayed
species (DR1). Both, fuel and oxidiser, once they are in gas form, they form a second group of delayed
species called DR2. Both delayed species contribute in the formation of the final products.
To model the rate of production of gas from the solid phase they used a simple Arrhenius expression
of the form,
Ea,s
r b = As e R u Ts (2.213)
Once the fuel, oxidiser and first group of reaction species are in gas phase, the chemical process
is described by the authors as,
!˙ F,ch = Af z ⇢2g e Ru Te
(2.217)
WO
⇣ ⌘
YeDR1
2 Ea,DR1
( " #)
e
@H e
@H 1 @ ⇣µ⌘ e ⇢
@H ⇣ µ ⌘ u2 /2
@e
m
⇢e
u + ⇢e
v = m r + µef f (2.223)
@x @r r @r Sc ef f @r P r ef f @r
✓ ◆ ✓ ◆2
@k @k 1 @ µt @k @e
u
⇢e
u + ⇢e
v = m rm µ + + µt ⇢✏ (2.224)
@x @r r @r C1 @r @r
✓ ◆ ✓ ◆2
@✏ @✏ 1 @ µt @✏ @e
u ✏ ✏2
⇢e
u + ⇢e
v = m rm µ + + C3 µt C4 ⇢ (2.225)
@x @r r @r C2 @r @r k k
55
CHAPTER 2. STATE OF THE ART
where u, u e, and u00 represent the instantaneous, the mass weighted mean and the fluctuating
velocities as,
e ⌘ ⇢u/⇢
u and e + u00
u⌘u (2.226)
and p, pe, p’ represent the instantaneous, the time mean and the fluctuating pressures as,
p ⌘ p + p0 (2.227)
Pressure follows ideal gas equation,
k2
µt = cµ ⇢ (2.229)
✏
As explained for AP/HTPB, to solve the problem, boundary conditions need to be applied in the
boundary solid-gas interface which are written by the authors as,
!
⇣ ⌘ @ Yei
⇢e e
v Yi ⇢s rb Yi,s ⇢D =0 (2.230)
g @r
@T @Tc
= s = ⇢c rb [(cp cc ) (Ts Ts,ref ) + Qc,ref ] (2.231)
@r g @r s
where Qc,ref is the heat release at a certain temperature Ts,ref , cp is the mass averaged specific
heat of gas species, cc is the heat capacity of solid propellant and D is the diffusion term.
The authors use the k-epsilon turbulence model to solve the equation system and discuss the
results accordingly.
dT d2 T Ec P M Eg,N O2
⇢p1 Rb Cp p = Qd ⇢p1 Yaldx Bc e RT + Qr,N O2 BN O2 ,c YN O2 e RT +
dx dx2 RT
✓ ◆2 E
2 PM g,aldx
Qr,aldx Baldx,c Yaldx e RT (2.235)
RT
56
2.4. MODELS TO CHARACTERISE DOUBLE-BASE PROPELLANTS COMBUSTION
where aldx is related to complex aldehyde reaction, B is the pre-exponential factor, E is the
activation energy, Qr is the heat release, Qd is the heat of endothermic degradation, ⇢p1 is the
density of the propellant, Cp is the specific heat capacity of propellant assumed constant, P is the
pressure and M is the average gas molar mass.
On the one hand, the reactions of the condensed phase are defined as,
N O2 ! N O + (2.237)
N O2 ! N O + (2.239)
dT d2 T P M Eg,N O2
ṁCg g 2
= Qr,N O2 Ag,N O2 YN O2 e RT +
dx dx RT
✓ ◆2 E
2 PM g,ald Eg,carb
+ Qr,ald Ag,ald Yald e RT + Qr,carb Ag,carb Scarb ⌘YN O ⇢carb P M e RT +
RT
✓ ◆2
PM Eg,N O
Qr,N O Ag,N O YN2 O e RT (2.243)
RT
dYN O2 P M Eg,N O2
ṁ = Ag,N O2 YN O2 e RT (2.244)
dx RT
✓ ◆2 E
dYald 2 PM g,ald
ṁ = Ag,ald Yald e RT (2.245)
dx RT
dYN O MN O P M Eg,N O2
ṁ = Ag,N O2 YN O2 e RT
dx MN O 2 RT
✓ ◆2
2 PM Eg,N O Eg,carb
Ag,N O YN O e RT Ag,carb Scarb ⌘YN O ⇢carb P M e RT (2.246)
RT
57
CHAPTER 2. STATE OF THE ART
According to the authors, and as it was also explained by Wu et al. [71] the propellant remains
stable in the preheated zone until the surface reaches the temperature of melting, which means that
the temperature is sufficiently high to start the solid degradation.
The problem analysed is a two-dimensional combustion chamber loaded with a double-base pro-
pellant grain and connected to an exhaust nozzle like the one depicted in Figure 2.11.
Figure 2.11: Schematic Diagram of solid-propellant rocket motor. Tseng and Yang [72].
To model it, a set of partial differential equation is defined in vector form as,
@Q @ @
+ (E Ev ) + (F F v) = S (2.248)
@t @x @y
58
2.4. MODELS TO CHARACTERISE DOUBLE-BASE PROPELLANTS COMBUSTION
where x and y are the axial and radial coordinates, the vector Q is defined as,
⇥ ⇤
E = ⇢u, ⇢u2 + p, ⇢uv, u(⇢e + p), ⇢uYi (2.250)
⇥ ⇤
F = ⇢v, ⇢uv, ⇢v 2 + p, v(⇢e + p), ⇢vYi (2.251)
2 3
0
6 ⌧xx 7
6 7
Ev = 6
6 ⌧xy 7
7 (2.252)
4 u⌧xx + v⌧xy qx 5
⇢ûi Yi
2 3
0
6 ⌧xy 7
6 7
Fv = 6
6 ⌧yy 7
7 (2.253)
4 u⌧xy + v⌧yy qy 5
⇢v̂i Yi
2 3
0
6 0 7
6 7
S=66 0 7
7 (2.254)
4 0 5
!˙ i
where the subscript i stands for ith-speciess, ⌧xy , ⌧xx , ⌧yy are the shear and normal stresses and
!˙ stands for the reaction rate.
The model developed by these authors, as it was done by Wu et al. [71], considers that the
degradation of the solid phase will give fuel (CH2 O), oxidiser (N O2 ) and a first group of delayed
reaction species. Once this three components are in gas phase, fuel and oxidiser react to produce
the second group of delayed species an both of them produce the final products. These products are
observed in R1, R2 and R3 defined by Wu et al. [71] which correspond to equations (2.214), (2.215)
and (2.216). Rate of production of the species is given by equations (2.217), (2.218) and (2.219).
In order to be able to couple both phases, and as it was also defined by other authors in the
existing literature, boundary conditions are needed at the burning surface . These conditions are
mass, species and energy balances,
(⇢v)g = ⇢c rb (2.255)
✓ ◆
@Yi
(⇢vYi )g ⇢D = ⇢s rb Yi,s (2.256)
@y g
✓ ◆ " N
# ✓ ◆ " N
#
@T X @T X
g + ⇢ Yi hi (v + v̂i ) = s ⇢c rb Cs (Ts Tref ) + Yi h0f,i (2.257)
@y g i=1
@y s i=1
g
↵ p Ac e
rb2 = T0 Qs
(2.258)
1 Ts 2cs Ts
59
CHAPTER 2. STATE OF THE ART
Ea,c
where ↵p is the thermal diffusivity of the condensed phase, Qs is the heat release, = Ru Ts and
Ts is the burning surface.
⇥ ⇤T
E = ⇢u, ⇢u2 + p, ⇢uv, u(⇢e + p), ⇢uYi (2.261)
⇥ ⇤T
F = ⇢v, ⇢uv, ⇢v 2 + p, v(⇢e + p), ⇢vYi (2.262)
2 3T
0
6 ⌧xx 7
6 7
Ev = 6
6 ⌧xy 7
7 (2.263)
4 u⌧xx + v⌧xy qex 5
q ix
2 3T
0
6 ⌧ 7
6 xy 7
Fv = 6
6 ⌧ yy
7
7 (2.264)
4 u⌧xy + v⌧yy qey 5
q iy
2 3
0
6 0 7
6 7
S=6 6 0 7
7 (2.265)
4 0 5
!˙ i
where subscript i stands for ith-species which varies from 1 to N-1, ⌧ represent the viscous stress
and qex and qey are the thermal diffusion term which stands for the contribution of heat conduction
and mass diffusion processes. The species diffusion terms (qi ) are approximated by Fick’s Law.
As mentioned, the kinetic mechanism is described in [70] so,
60
2.4. MODELS TO CHARACTERISE DOUBLE-BASE PROPELLANTS COMBUSTION
EN O
2
!˙ N O2 = YN O20 !˙ c ⇢g YN O2 p0.39 BN O2 e RT (2.269)
In addition to the gas system of equations of the gas phase, several equations are considered for
the gas phase being, mass, energy and species concentration,
ṁ = ⇢c rb (2.270)
@T @T @2T
⇢c C c + ṁCc = c + q̇c (2.271)
@t @y @y 2
@Yi @Yi
⇢c + ṁ = !˙ i (2.272)
@t @y
where q̇c is the net effect of the endothermic decomposition and the exothermic reaction in the
condensed phase.
The coupling between both phases is done by applying mass, energy and species balance at the
interface so,
(⇢v)g = ⇢c rb (2.273)
✓ ◆ " N
# ✓ ◆ " N
#
@T X @T X
g + ⇢ Yi hi (v + v̂i ) = c ⇢c r b C c T s Ts + Yi h0f,i (2.275)
@y g i=1
@y s i=1
g
The degradation of the solid is defined with an expression similar to that by Tseng and Yang [72].
↵ c Ac e /
rb2 = n o (2.276)
g /ṁ)( @y ) +[Qc +(Cc Cp )(Ts Tref )]·0.5
@T
(
g
Cs T s
61
62
Chapter 3
Numerical methods
The aim of numerical methods is to to solve the PDEs by replacing the continuous problem of Fluid
Dynamics represented by the PDEs (Partial Differential Equations) by a finite set of discrete values.
In order to do that, the domain of the PDEs must be discretised in a finite set of points of volumes
either by a mesh or a grid. In the Finite Difference approach, the discrete values are point values
defined at grid points however, in the Finite Volume approach, the discrete values are average values
over finite volumes.
The multi-dimensional Euler equations used in this work for modelling and describing the physical
problem of combustion of propellants can be written with the following general conservative form,
@U
+ r·H = S(U ) (3.1)
@t
where U is the conserved variables vector, H = (F , G, H) is the tensor of fluxes and S is the
source terms which may or not include non conservative terms and can be defined as,
2 3 2 3 2 3 2 3 2 3
q1 f1 g h1 s1
6 q2 7 6 f2 7 6 g2 7 6 h2 7 6 s2 7
U =6 7 6 7 6 7 6 7 6
4 ... 5 , F = 4 ... 5 , G = 4 ... 5 , H = 4 ... 5 , S = 4 ... 5
7 (3.2)
qm fm gm hm sm
The integration of the homogeneous part of the equation system (3.1) in a control volume ⌦ yields,
Z Z Z Z Z
@
U d⌦ + Hn̂dA = 0 (3.3)
@t ⌦ A
where A is the boundary of the volume ⌦ and n̂ is the normal vector to the surface A depicted in
Figure 3.1.
Considering the first integral as a time rate of change of the averaging of the conserved variables
P
N
U and the boundary A formed by N surfaces so that A = As equation 3.3 can be written as,
s=1
N Z Z
@U 1 X
+ Hn̂dA = 0 (3.4)
@t |⌦| s=1 As
63
CHAPTER 3. NUMERICAL METHODS
U n+1 Un
The time derivative of the conserved variables is defined as j
4t
j
. Therefore, making use of
1
the Rotational property Hn̂ = Ts F Ts U and considering that the surface integral of the flux is
RR 1
approached by As
F n̂dA ⇡ Ts F Ts U As , a finite volume scheme for multiple dimensions in
unstructured grids is obtained,
N
4t X 1
U n+1
j = U nj Ts F Ts U nj As (3.5)
|⌦| s=1
1
where Ts is the rotation matrix, Ts its inverse and As is the area of the sth surface bounding volume
⌦ . Considering the source terms, the system results,
N
n+1 4t X 1
Uj = U nj Ts F Ts U nj As + 4tS(U nj ) (3.6)
|⌦| s=1
where 4t is the time step which is calculated with the following expression,
✓ ◆
4xj
4t = CF L· (3.7)
such that,
0 < CF L 1 (3.8)
CF L stands for the Courant-Friedichs-Lewy coefficient or the Courant number coefficient which
and it is a dimensionless quantity. condition of convergence to solve partial differential equations.
The discretisation of the domain in cells is represented in Figure 3.2.
In order to evaluate the fluxes, several approximations such as the first-order upwind method of
Godunov, Lax-Friedrichs scheme and Lax-Wendroff scheme can be found in the existing literature.
These methods, require the solution of the Riemann problem which can be highly demanding.
Therefore, approximate Riemann solvers are used in this work. Approximate and non-iterative solu-
tions can be used to provide the information for numerical purposes. These information can be
extracted in two ways: one is to find an approximation of the numerical flux employed in the numer-
ical method and the other is to find an approximation to a state and then evaluate the numerical
flux in this state. The approximate Riemann solvers do not need an iteration process since they
approximate the solution for the state required to evaluate the flux.
64
3.1. RUSANOV NUMERICAL SCHEME
Figure 3.2: Discretisation of domain [0, L] into M finite volumes Ii (computing cells) [79].
Among the numerical schemes found in literature, Rusanov and MacCormack numerical schemes
have been chosen in this thesis to solve the problem of propellant combustion. Throughout this
chapter, not only these two, but also Lax-Wendroff scheme, used by the authors already mentioned
in the literature review, will be presented.
The approximation consist in finding, for each interface between two cells, the approximate nu-
merical flow that will be an average value of the numerical flows computed for each cell in the previous
step time, corrected with a factor. This factor is proportional to the difference of conservative vari-
ables of each cell in the previous step time and a positive speed value S + , so the numerical flux at
each interface is computed as,
1
F nj+1/2 = F Ts U nL + F Ts U nR S + (U nR U nL ) (3.9)
2
where R and L stand for right and left respectively, F nj+1/2 is the flux in the interface, F Ts U nL
and F Ts U nR are the convective fluxes evaluated in the previous time step in the left and right cells
65
CHAPTER 3. NUMERICAL METHODS
respectively and S + is a positive value of speed calculated by means of the wave speeds SL and SR
being,
p 4t
U j = U nj F nj+1 F nj + 4tS(U nj ) (3.13)
4x
c 4t
U j = U nj F pj F pj 1 + 4tS(U pj ) (3.14)
4x
1⇣ p c
⌘
U n+1
j = Uj + Uj (3.15)
2
According to Liang et al. [84] the use of a total variation diminishing (TVD) of the MacCormack
scheme will reproduce all flow regimes in an accurate way. Considering the previous equation system
(3.1), the scheme can be written as,
p 4t
U j = U nj F nj F nj 1 + 4tS(U nj ) (3.16)
4x
c 4t
U j = U nj F pj+1 F pj + 4tS(U pj ) (3.17)
4x
1⇣ p c
⌘
U n+1
j = U j + U j + [G(rj+ ) + G(rj+1 )] U nj+1/2 [G(rj+ 1 ) + G(rj )] U nj 1/2 (3.18)
2
where G(x) is a function defined as,
66
3.3. LAX-WENDROFF SCHEME
4U nj 1/2 = U nj U nj 1 (3.23)
The values of rj+ and rj depend on the definition of equations (3.22) and (3.23),
4U nj 1/2 · 4U nj+1/2
rj+ = (3.24)
4U nj+1/2 ·4U nj+1/2
4U nj 1/2 · 4U nj+1/2
rj = (3.25)
4U nj 1/2 ·4U nj 1/2
According to Liang et al. [84] , the use of this scheme will prevent numerical oscillation close to
the sharp gradient regions.
@uni 1 @ 2 un
un+1
i = uni + 4t + 4t2 2i + O 4t3 (3.27)
@t 2 @t
By means of the Cauchy-Kowalewski procedure, the time derivatives can be replaced by space
derivatives,
@q(x,t) @q(x,t)
9
@t = @x
>
=
@q 2 (x,t) 2
2 @q (x,t)
@t2 = @x2
(3.28)
>
@q k (x,t) k @q (x,t) ;
k
@tk
= ( ) @xk
Therefore, equation (3.27) can be written,
@uni 1 @ 2 uni
un+1
i = uni + 4t + 4t2 2 + O 4t3 (3.29)
@x 2 @x2
If the derivatives are approximated by central finite differences (see Figure 3.4)
67
CHAPTER 3. NUMERICAL METHODS
Figure 3.4: Finite differences approximation to the first derivative of a function g(x) [79].
The conventional finite difference method of Lax-Wendroff can also be written in conservative
form.
Considering the one-dimensional homogeneous problem,
@U @F
+ =0 (3.32)
@t @x
A conservative numerical method to solve the previous equation is the scheme in finite volumes,
4t
U n+1
j = U nj F j+1/2 Fj 1/2 (3.33)
4x
where F j+1/2 is the numerical flux which for the two-step version of the Lax-Wendroff method is
written as,
n+1/2
F j+1/2 = F (U j+1/2 ) (3.34)
n+1/2
being U j+1/2 defined with the following expression,
n+1/2 1 1 4t
U j+1/2 = U nj + U nj+1 + F nj F nj+1 (3.35)
2 2 4x
68
Chapter 4
@ (↵⇢g ) !
+ r· (↵⇢g !
u g) = ˙ (4.1)
@t
@ (↵⇢g !
u g) ! ⇣ !⌘ @↵ !
+ r· ↵⇢g !
ug ⌦!u g + ↵pg I = pg + ˙!
up FD (4.2)
@t @x
✓ ✓ ◆◆ ✓ ◆
@ (↵⇢g Eg ) ! pg @↵ |!
u p| 2 !
+ r· ↵⇢g !
u g Eg + = pg +˙ g
Echem + !
u pF D Q̇ (4.3)
@t ⇢g @t 2
@ ((1 ↵) ⇢p ) !
+ r· ((1 ↵) ⇢p !
u p) = ˙ (4.4)
@t
@ ((1 ↵) ⇢p !
u p) ! ⇣ !⌘ @↵ !
+ r· (1 ↵) ⇢p !
up⌦!
u p + (1 ↵) pp I = pp ˙!
up+ FD (4.5)
@t @x
✓ ✓ ◆◆
@ ((1 ↵) ⇢p Ep ) ! ! pp
+ r· (1 ↵) ⇢p u p Ep + =
@t ⇢p
✓ ◆
@↵ |!
u p| 2 !
pg +˙ p
Echem + +!
u p F D + Q̇ (4.6)
@t 2
69
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
pp = pg + ⌧ p (4.7)
where the chemical energy of the gas phase is defined as,
g p
Echem = Echem Echem (4.8)
@ ((1 ↵) ⇢p !
u p) ! ⇣ ! ⌘
+ r· (1 ↵) ⇢p !
up⌦!
u p + (1 ↵) pp I (1 f (↵)) =
@t
@ (↵ (1 f (↵))) ˙! !
pp u p + (1 f (↵)) F D (4.9)
@x
where f (↵) is defined as:
8
>
> 0 ↵ > ↵c
< ⇣ ⌘l
f (↵) = ↵c ↵ ↵min ↵ ↵c (4.10)
>
> ↵c ↵min
: ↵ < ↵min
1
where l is a constant value that considers the properties of the particles framework, ↵min and ↵c
are the minimum and critical porosity respectively.
Therefore, the equation system could be finally written as,
@ (↵⇢g ) !
+ r· (↵⇢g !
u g) = ˙ (4.11)
@t
u g) ! ⇣
@ (↵⇢g ! !⌘ @↵ !
+ r· ↵⇢g !
ug ⌦!u g + ↵pg I = pg + ˙!
up FD (4.12)
@t @x
✓ ✓ ◆◆ ✓ ◆
@ (↵⇢g Eg ) ! pg @↵ |!
u p| 2 !
+ r· ↵⇢g !
u g Eg + = pg +˙ g
Echem + !
u pF D Q̇ (4.13)
@t ⇢g @t 2
@ ((1 ↵) ⇢p ) !
+ r· ((1 ↵) ⇢p !
u p) = ˙ (4.14)
@t
@ ((1 ↵) ⇢p !
u p) ! ⇣ ! ⌘
+ r· (1 ↵) ⇢p !
up⌦!
u p + (1 ↵) pp I (1 f (↵)) =
@t
@ (↵ (1 f (↵))) ˙! !
pp u p + (1 f (↵)) F D (4.15)
@x
✓ ✓ ◆◆
@ ((1 ↵) ⇢p Ep ) ! pp
+ r· (1 ↵) ⇢p !
u p Ep + =
@t ⇢p
✓ ◆
@↵ |!
u p| 2 !
pg +˙ p
Echem + +!
u p F D + Q̇ (4.16)
@t 2
pp = pg + ⌧ p (4.17)
70
4.1. COMBUSTION OF GRANULATED PROPELLANTS
µg (ug up )
FD = fpg (4.19)
d2p
where fpg is the drag coefficient. Butler et al. [39] present in their study an expression developed
by Wilcox, Jones and Krier [40]. The use of that expression provided good result at high Reynolds
numbers.
⇣ ⌘0.88
(1 ↵)2 Rep
fpg = ↵2 150+ 3.89 (1 ↵) (4.20)
10 < Rep < 2 · 105
3
As for the interfacial heat transfer coefficient, the heat expression form is,
Sp
Q̇ = (1 ↵) hc (Tg Tp ) (4.21)
Vp
In the case of spherical particles, Sp /Vp = 6/dp , and therefore the heat is defined as:
Sp
Q̇ = 6(1 ↵) hc (Tg Tp ) (4.22)
dp
The convective heat coefficient depends of the Nusselt number (Nu) with the following expression,
kg
hc = Nu (4.23)
dp
where P r = (cp,g µg )/kg is the Prandtl number and kg is the gas conductivity defined as,
15 Ru 4 cv,g Mg 3
kg = µg ( + ) (4.24)
4 Mg 15 Ru 5
The correlation of the Nusselt number proposed by Hoffman and Krier [28], Butler et al. [39] and
Butler and Krier [32] has the following form:
N u = 0.65Re0.7
p Pr
0.33
(4.25)
The gas mass generation equation assuming that particles are spheres is related to the combustion
velocity and has the following form,
˙ = 6(1 Sp
↵) ṙ⇢p (4.26)
dp
where the combustion law is given by:
71
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
ṙ = apng (4.27)
where a and n are vales that depend on the propellant.
As for the gas equation of state, the equation of Noble-Abel is considered.
⇢g ( g 1)eg
pg = (4.28)
1 ⇢g ⌘ g
In this model, the solid is considered incompressible. However, as presented in equation 4.17
in the model there is a solid pressure defined as the addition of gas pressure and the intergranular
stress (⌧p ), defined as the resistance of the particles to be compactedunder a limit porosity that is the
critical porosity. In this study, the authors consider the expressions proposed by Krier and Kezerle
[27], Hoffman and Krier [28] and Gokhale and Krier [31].
( ⇣ ⌘
K 1 1
↵ ↵c
⌧p = (1 ↵) (1 ↵c ) (1 ↵)
(4.29)
0 ↵ > ↵c
where K is the bulk modulus and ↵c is the critical porosity under which appears the intergranular
stress.
0 1
˙
B ! C
B ✓pg @↵
@x + ˙! up ◆ FD C
B ! ! C
B g | u p| 2
! C
B pg @t + ˙ Echem + 2
@↵
u p F D Q̇ C
B C
S(U ) = B ˙ C (4.33)
B C
B ! C
B (1 f (↵)) pp @↵ ˙! up+ (1 f (↵)) F D C
B ✓ @x ◆ C
@ p |!u p |2 ! ! A
@↵
pg @t + ˙ Echem + 2 + u p F D + Q̇
72
4.1. COMBUSTION OF GRANULATED PROPELLANTS
To solve the system of equation (4.30), Rusanov and MacCormack schemes detailed in chapter 3
are used.
73
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
The behaviour of the variables obtained when using coefficient Friction 2 and 3 (Figures 4.3 and
4.4) is exactly the same that the one with Friction 1.
Figure 4.1: Influence of Reynolds number in the Nusselt number for several values of porosity 0.6
(above-left), 0.5 (above-right), 0.4 (below-left) and 0.3 (below-right).
74
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figure 4.2: Distributions of pressure, porosity and gas velocity for Rusanov (left) and MacCormack-
TVD (right) with Friction 1.
75
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.3: Distributions of pressure, porosity and gas velocity for Rusanov (left) and MacCormack-
TVD (right) with Friction 2.
76
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figure 4.4: Distributions of pressure, porosity and gas velocity for Rusanov (left) and MacCormack-
TVD (right) with Friction 3.
77
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.5: Influence of Reynolds number in the friction coefficient for several values of porosity 0.6
(above-left), 0.5 (above-right), 0.4 (below-left) and 0.3 (below-right).
78
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figure 4.6: Distributions of pressure obtained with Rusanov (above) and MacCormack-TVD (below)
numerical schemes for Nusselt 1.
79
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.7: Distributions of pressure obtained with Rusanov (above) and MacCormack-TVD (below)
numerical schemes for Nusselt 2.
80
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figure 4.8: Distributions of pressure obtained with Rusanov (above) and MacCormack-TVD (below)
numerical schemes for Nusselt 3
81
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
ent x-coordinates depending on the initiation test used. However, Rusanov and MacCormack-TVD
numerical schemes provide very similar results regardless the initiation test used. MacCormack-TVD
presents oscillatory results for all test despite Test 1 which behaves smoothly. It can be concluded
that Test 1 is the one which presents less oscillations in all significant variables and similar results
independently the numerical scheme used. However, high values of pressure and temperature are ob-
tained. As visible in the results, the initiation influences enormously in the results of the simulation.
Due to the difficulty to identify the initiation parameters from bibliography, the results from Hoffman
and Krier [28] at 20 µs are used as initial conditions in order to validate the developed code.
Figure 4.9: Pressure distribution obtained with Rusanov (left) and MacCormack TVD (right) for 20
µs, 40 µs and 60 µs
82
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figure 4.10: Porosity distribution obtained with Rusanov (left) and MacCormack TVD (right) for
20 µs, 40 µs and 60 µs
83
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.11: Gas temperature distribution obtained with Rusanov (left) and MacCormack TVD
(right) for 20 µs, 40 µs and 60 µs.
84
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figure 4.12: Particle temperature distribution obtained with Rusanov (left) and MacCormack TVD
(right) for 20 µs, 40 µs and 60 µs.
85
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.13: Gas velocity distribution obtained with Rusanov (left) and MacCormack TVD (right)
for 20 µs, 40 µs and 60 µs.
86
4.1. COMBUSTION OF GRANULATED PROPELLANTS
4.1.4 Results
In this section, the results obtained once the source terms and the initiation has been defined will be
presented. Firstly, Hoffman and Krier [28] equation system with their modification to the momentum
conservation equation detailed in Section 4.1.1.2 is applied. As will be seen, the system of equations
(4.14)- (4.17) provides a solution with a flame front similar to that obtained by Hoffman and Krier
[28]. However, the values obtained for the temperature of the particle do not match the results
obtained by the authors in the article. Afterwards, the results obtained when limiting directly the
porosity in the equation system (4.4) - (4.7), will provide a displacement of the peak values of pressure
and pressure when time is equal to 60 µs compared those obtained by Hoffman and Krier [28] in their
work.
1
McCormack t = 20 µs
Rusanov t = 20 µs
0.9 McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
0.8 Rusanov t = 60 µs
0.7
0.6
0.5
α
0.4
0.3
0.2
0.1
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
87
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
8
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
7 Rusanov t = 40 µs
McCormack t = 60 µs
Rusanov t = 60 µs
6
5
pg (GPa)
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
4000
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
Rusanov t = 60 µs
3000
Tg (K)
2000
1000
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
88
4.1. COMBUSTION OF GRANULATED PROPELLANTS
12000
McCormack t = 20 µs
11000 Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
10000 McCormack t = 60 µs
Rusanov t = 60 µs
9000
8000
7000
Tp (K)
6000
5000
4000
3000
2000
1000
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
0.06
McCormack
Rusanov
Hoffman and Krier, 1980
0.04
Distance (m)
0.02
0
0 20 40 60 80
Time (µs)
89
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
However, the particle temperature shows high discrepancies between both numerical schemes. In
addition, the magnitude of the variable does not seem to adjust properly to physical conditions since
the gas temperature does not get over 4000 K. In order to reduce particle temperature, the partial
derivative of alpha with time in equation 4.16 is deleted, resulting an equation system as,
@ (↵⇢g ) !
+ r· (↵⇢g !
u g) = ˙ (4.37)
@t
@ (↵⇢g !
u g) ! ⇣ !⌘ @↵ !
+ r· ↵⇢g !
ug ⌦! u g + ↵p I = pg + ˙!
up FD (4.38)
@t @x
✓ ✓ ◆◆ ✓ ◆
@ (↵⇢g Eg ) ! ! p @↵ g |!
u p| 2 ! !
+ r· ↵⇢g u g Eg + = pg +˙ Echem + u pF D Q̇ (4.39)
@t ⇢g @t 2
@ ((1 ↵) ⇢p ) !
+ r· ((1 ↵) ⇢p !
u p) = ˙ (4.40)
@t
@ ((1 u p) ! ⇣
↵) ⇢p ! ! ⌘
+ r· (1 ↵) ⇢p !
up⌦!
u p + (1 ↵) p I (1 f (↵)) =
@t
@ (↵ (1 f (↵))) ˙! !
pp u p + (1 f (↵)) F D (4.41)
@x
✓ ✓ ◆◆
@ ((1 ↵) ⇢p Ep ) ! ! p
+ r· (1 ↵) ⇢p u p Ep + =
@t ⇢p
✓ ◆
|!
u p| 2 !
+˙ p
Echem + +!
u p F D + Q̇ (4.42)
2
pp = pg + ⌧ p (4.43)
The results obtained are detailed in Figures 4.19, 4.20, 4.21, 4.22 and 4.23. The number of cells
and parameter f chosen for calculation are the same as before.
0.96 McCormack t = 20 µs
Rusanov t = 20 µs
0.88 McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
0.8 Rusanov t = 60 µs
0.72
0.64
0.56
α
0.48
0.4
0.32
0.24
0.16
0.08
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
90
4.1. COMBUSTION OF GRANULATED PROPELLANTS
6
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
5 McCormack t = 60 µs
Rusanov t = 60 µs
4
pg (GPa)
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
4000
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
Rusanov t = 60 µs
3000
Tg (K)
2000
1000
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
91
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
5000
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
4000 Rusanov t = 60 µs
3000
Tp (K)
2000
1000
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
0.06
McCormack
Rusanov
Hoffman and Krier, 1980
0.04
Distance (m)
0.02
0
0 20 40 60 80
Time (µs)
92
4.1. COMBUSTION OF GRANULATED PROPELLANTS
Figures 4.19, 4.20 and 4.21 show similar results to those obtained before deleting the temporal
partial derivative in the particle energy equation. However, in Figure 4.22 can be appreciated how
the temperature of the solid phase decrease up to 3000 K. Flame front distribution obtained with
this modification is similar to the previous one. As before, the code has been run with 250 cells and
f = 0.2 for the solution with Rusanov scheme and 100 cell and f = 0.11 for MacCormack-TVD.
0.96 McCormack t = 20 µs
Rusanov t = 20 µs
0.88 McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
0.8 Rusanov t = 60 µs
0.72
0.64
0.56
α
0.48
0.4
0.32
0.24
0.16
0.08
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
93
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
4
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
Rusanov t = 60 µs
3
pg (GPa)
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
6000
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
5000 McCormack t = 60 µs
Rusanov t = 60 µs
4000
Tg (K)
3000
2000
1000
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
94
4.1. COMBUSTION OF GRANULATED PROPELLANTS
3000
McCormack t = 20 µs
Rusanov t = 20 µs
McCormack t = 40 µs
Rusanov t = 40 µs
McCormack t = 60 µs
Rusanov t = 60 µs
2000
Tp (K)
1000
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
x (m)
0.1
McCormack
Rusanov
Hoffman and Krier, 1980
0.08
0.06
Distance (m)
0.04
0.02
0
0 20 40 60 80
Time (µs)
95
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
4.1.5 Conclusions
A code to study the modelling of granular propellants combustion has been developed in C++. To
solve the system of differential equations for both, gas and solid phase, Rusanov and MacCormack-
TVD numerical schemes have been used. In addition, first-order Euler method has been used to solve
the source terms.
Due to the amount of different expression found in literature to evaluate the source terms, there
is a necessity of determining which one is the most suitable for the model developed in this work.
Therefore, a sensitivity analysis is performed to determine the effect that Nusselt number, friction
coefficient and initiation conditions have in the results. Rusanov and MacCormack numerical schemes,
detailed in Chapter 3, have been used to solve the system of equations. The sensitivity analysis
determines on the one hand, which expression of Nusselt number and friction coefficient is the one
which provides more similar results for the numerical schemes used. On the other hand, it is used
to determine how the application of the different initiation conditions found in literature affect to
the variables distributions. As a conclusion, due to the capacity of reproducing very similar results
between both numerical schemes and the good behaviour at a wide range of Reynolds number, the
Nusselt number and friction coefficient presented by Butler et al. [39] in their work will be adopted
afterwards for calculation.
As happened with the closure laws, although all authors consider a small zone of the control
volume already initiated, each one applies different values of pressure and temperature as initial
conditions of the test. High differences in the results can be obtained depending on the initial values
assumed for the variables. Therefore, and in order to validate the model with literature values, the
distributions of pressure, porosity and gas and particle temperatures provided by Hoffman and Krier
[28] at 20 µs, are used as initial conditions for the calculations performed in subsection 4.1.4.
Afterwards, and once initial conditions and source terms are defined, two different models have
been used to calculate pressure, temperature and porosity distributions. The first model studied
is defined by the system of differential equations (4.14) - (4.17) in which the modification of the
particle momentum governing equation done by Hoffman and Krier [28] is considered. According to
the authors, this modification prevents the porosity from reaching values below the minimum value of
compaction that packed beds of spherical particles can reach. However, when this model is applied,
the temperature of the particle phase increases and reaches values up to 8000K. In order to reduce
the particle temperature, a modification of this equation system has been done and the partial time
derivative of the porosity in the particle momentum equation is deleted leading to the system (4.40) -
(4.43). This modification reduces slightly the maximum value of particle temperature but the results
obtained are not matching quantitatively the ones found in the literature. The second model studied
does not consider the porosity limiter from Hoffman and Krier [28] performing the limitation of the
porosity directly in the code by preventing the porosity to reach values below the minimum one.
The magnitude of the values obtained for the main variables of interest when applying this model
(equations (4.4) - (4.7)) are similar to those found in literature. Moreover, the results obtained using
Rusanov scheme agree well with those resulting of applying MacCormack-TVD numerical scheme.
However, the distributions of the variables are displaced in x-direction respect those from literature.
This behaviour can be seen by observing the x-coordinate position of the peak values obtained for
pressure, temperature and porosity variables. These differences could be due on the one hand, to the
initial values of the parameters chosen as initial conditions which highly determine the x-location of
the peak values for all variables and on the other hand, to the lack of all necessary input data in a
single work from bibliography making necessary to collect the values from different works increasing
the difficulty of reproducing the tests available in the bibliography.
Finally, it can be concluded that, despite the difference in the x-location peak values, the last
model considered represents accurately the physical behaviour of the propellant combustion for all
variables of interest becoming a predictive tool for the characterisation of the early stages of the
detonation process of granular solid propellants.
96
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
@ (⇢g Yi ) ! kg
+ r· (⇢g !
u g Yi ) = r2 Y i + i i = 1, 2, ..., n 1 (4.49)
@t g cvg
@ (⇢c Ec )
= k c r2 Tc (4.50)
@t
where the subscripts i stands for the ith-species from 1 to n-1 and j stands for reaction number
from 1 to N. The specific heat at constant volume for a mixture of gases is defined as,
n
X
cvg = Yi cvg,i (4.51)
i=0
97
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Tg = Tc = Ts (4.52)
3. The continuity of mass must be ensured. When the temperature of the surface is higher than
the melting temperature (Tmelt ), the solid changes phase to gas. Otherwise the velocity in the
burning surface must be zero. Therefore, the following relationship is set,
⇢ ⇢c ṙ
ug = ⇢g , Ts > Tmelt (4.54)
0, Ts < Tmelt
where ṙ is the burning rate and i stands for the specie from 1 to n-1. Among all expressions found in
the literature to define the burning rate of solid propellants (see 2.3), that defined by Cai et al. [49],
which considers the burning rate as a function of the pressure, is employed in this case. Therefore,
the burning rates expressions of AP and HTPB which result after applying the formula are,
0 ⇣
E1
⌘1
n1
Qg,1 2 kg (Tg Ts ) ↵AP @ D1 p e RuTg
A Qg,1
ṙ + ṙAP = (4.56)
2cAP Ts AP d⇢AP cAP Ts ⇢AP E1
RuT
Qc,AP
g
Qg,2 2 kg (Tg Ts )
ṙHT PB + ṙHT P B
2cHT P B Ts d⇢HT P B cHT P B Ts
0 ⇣
E2
⌘1
n2
↵HT P B @ D2 p e R u Tg
A Qg,2
= E2
(4.57)
⇢HT P B R T
Qc,HT P B
u g
When the propellant is heterogeneous, each cell will move at a different velocity,ṙAP or ṙHT P B
depending on the content of the cell. However, when the propellant is distributed uniformly, the
burning rate is defined as an average value between both burning rates being function of the mass
fraction of each specie present in the propellant. Therefore, in case of having and homogeneous
AP/HTPB propellant, the burning rate can be defined as,
98
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
p = ⇢g R g T g (4.59)
being Rg the specific constant of gases which for a mixture of gases is as calculated as,
n
X
Rg = Ru Yi /Mg,i (4.60)
i=0
The production rates of AP and HTPB combustion reactions (4.44) and (4.45) are described with
an Arrhenius expression as a function of the mass fraction of the main reactive, the pressure and the
temperature of the gas mixture as,
⇣ ⌘
E1
R1 = D1 pn1 YAP e RuTg
(4.61)
⇣ ⌘
E2
n2
R2 = D2 p YHT P B YO2 e RuTg
(4.62)
where D1 , and D2 are the pre-exponential factors of the reaction rate, and n1 and n2 are the
pressure exponents and E1 and E2 are the activation energies. All this parameters will be considered
input data and will be taken from literature.
99
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
0 1
0
B 0 C
B N C
S 2 (U ) = B
B
P C
C (4.67)
@ j=1 Qg,j Rj A
i
where the numerical flux is computed by solving approximately the Riemann problem at each
interface and the source terms have been evaluated at tn leaving the scheme explicit. The fluxes
are calculated using Rusanov numerical scheme.
2. Next, the ODE
@U
= S 2 (U ) (4.70)
@t
which considers, as mentioned before, the source terms associated to the combustion process,
is solved by using a first-order Euler method,
n+1 n+1
U n+1
j = Uj + 4tchem S 2,j (4.71)
2. Reaction rates R1 and R2 , for AP and HTPB respectively, are calculated with the primitive
values obtained from equation system 4.69.
100
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
3. With these reaction rates, a value of time step associated to the ith-speciess consumed in
jth-reaction is calculated as follows,
⇢g Y i
4tchem,i,j = (4.74)
g,i M
Rj coefi,j Mg,j
where i stands for the ith-speciess and j for reaction and coefi,j was defined as,
being coefi p,j the stoichiometric coefficient of the ith-species in jth-reaction as product, and
coefi r,j the stoichiometric coefficient of the ith-species in jth-reaction as reactive. There will
be as many 4tchem,i,j as consumed species in jth-reactions. From all those 4tchem,i,j the
smallest one is chosen (4tchem ).
4. The chemical time (4tchem ) is compared with that imposed initially (4tchem,0 ) . In case of
being smaller, the calculated 4tchem is used to solve equation (4.71). In opposite case, the
initial chemical time step (4tchem,0 ) will be used to solve the equation.
5. Equation (4.71) is solved and primitive variables are updated.
6. Remaining time is calculated in first iteration as the difference between total numerical time
step step and the calculated chemical time step,
7. The primitive variables calculated after solving Equation (4.71) are used to estimate again the
reaction rates and with them the chemical time step of each species using equation (4.74). If the
smallest chemical time step (4tchem ) is smaller than the remaining numerical time step step
(4tremaining ), 4tchem is used to solve equation (4.71) again. This process is repeated until the
calculated chemical time step is higher than the remaining numerical time step.
8. When this happens, the remaining numerical time step step is the one used to solve equation
(4.71) and the calculation goes to the next convective 4t according to equation 4.72.
U = ⇢c E c (4.80)
and S(U ) as,
S(U ) = kc r2 Tc (4.81)
101
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
102
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
have been considered for temperature, density, mass fractions and pressure along the walls despite
the pressure at the exit of the tube which have been fixed either to 4 MPa or 7 MPa depending on the
reproduced test. On the other hand, velocity have been fixed to zero in the left closed wall, to zero
gradient in the right open wall and slip boundary condition have been considered in upper-bottom
and front-back walls.
Regarding initial conditions, the control volume is assumed to be already initiated. The initial
condition of temperature is described by a linear distribution from the value fixed on the left side,
which can be 220.15 K, 294.15 K or 347.15 K depending on the test considered, to 811 K on the right
side. Initial pressure in the gas domain have been set as the pressure at the right boundary. Finally,
the initial velocity of the gas phase have been set to zero.
The result of the calculation and their comparison with the experimental results from INTA
facilities are plotted in Figure 4.30. Quantitative results, absolute and relative error scan be found
in Tables 4.4 and 4.5.
20 20
OpenFoam2D OpenFoam2D
18 INTA 18 INTA
16 16
14 14
12 12
r (mm/s)
r (mm/s)
10 10
8 8
.
6 6
4 4
2 2
0 0
200 240 280 320 360 400 200 240 280 320 360 400
Temperature (K) Temperature (K)
Figure 4.30: Comparison between calculated and experimental burning rate at 4 MPa (left) and 7
MPa (right).
Table 4.4: Quantitative comparison of calculated and experimental burning rate for 4 MPa.
103
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Table 4.5: Quantitative comparison of calculated and experimental burning rate for 7 MPa.
For both pressures, the difference between numerical and experimental results for 220.15 K is
lower than that obtained either for 347.15 K nor 220.15 K. Therefore, it can be said that the smallest
absolute error is obtained when initial temperature is 220.15 K in the range of pressure studied.
Moreover, the difference between the results obtained with the model and experimental ones when
initial temperature is 220.15 K is higher than for 347.15 K. In addition, differences between code and
experimental results increase with the pressure, obtaining a better agreement for 4 MPa than for 7
MPa being the magnitude of the absolute error less than 1 mm/s for 4 MPa and around 1 mm/s 7
MPa, being in the order of magnitude of the experimental error. However, the relative error, not only
increases considerable for low initial temperatures independently of the pressure of the test, but also
for low pressures up to reach values of 36%. This is due to the relation used to calculate error,
|Vexperimental Vnumerical |
Relative error = (4.82)
Vexperimental
where V is the value of the magnitude whose error is being calculated and |Vexperimental Vnumerical |
is its absolute error. Despite the absolute error is mostly constant independently of the initial tem-
perature, the relative error increases its magnitude drastically. This happens because experimental
value of the burning rate decrease with the initial temperature chosen to run the test and this affects
to the value calculated with expression (4.82).
As a conclusion, it can be said that the results obtained with the numerical model agree with the
experimental ones for the pressure range considered.
104
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
According to Cao et al. [47], the width of the domain is obtained as a function of both AP and HTPB
particle and AP mass fraction as per the following two expressions,
⇢AP ·dAP ↵
= (4.84)
⇢HT P B ·dHT P B (1 ↵)
A strand with 80% of AP content and particle average size of 110 mm was used for calculation,
as it is also considered by Cao et al. [47]. In this case and according to equations (4.84) and (4.83),
the widths obtained are 60 mm for HTPB and 170 mm for the total height. In order to define a 3D
domain, both condensed and gas phase are set 500 mm height as seen in Figure 4.31. This first test
is run in 2D configuration for several pressures in a range from 0.5 MPa to 7 MPa.
Firstly, a comparison between the burning rate obtained using the proposed model and the ex-
perimental results from Kohga [51] is presented in Figure 4.32. The results obtained with the model
previously explained agree well for the considered pressure range and the difference between calcu-
lated and experimental data maintains the same magnitude independently of the pressure considered
which can be seen quantitatively in Table 4.7. The absolute error, remains constant and have a value
around 0.5 mm/s being in the order of magnitude of the experimental error. However, the value of
the relative error increases considerably with the pressure fall to reach up to 22%. This is due to
expression (4.82) used to calculate error as explained in subsection 4.2.3.1.
Besides burning rate, other relevant variables such as gas velocity and AP and HTPB mass fraction
distributions are obtained in Figures from 4.33 to 4.41. The first phenomenon that can be observed
at glance in any of these figures is the decrease of the distance between the left side of the tube and
the burning surface with the increase of the pressure. That means, the higher the pressure is, the
faster the condensed phase changes from solid to gas which in absolute concordance with the results
of Figure 4.32 where the direct relationship between burning rate and increase of pressure is plotted.
Velocity field is depicted in Figures 4.33, 4.34 and 4.35 for time equal to 0.02 s, 0.04 s and 0.06 s
respectively. The behaviour of the gas velocity is opposed to the burning rate, the higher the pressure
is, the lower gas velocity we have. The reason of this behaviour could be the increase of gas density
105
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
with pressure. By comparing the render views for 0.5 MPa of Figures 4.33, 4.34 and 4.35 against
Figures 4.36, 4.37 it can be observed that an increase of AP mass fraction close to the burning surface
leads to a local raise of the gas velocity. Since the increase of AP mass fraction is directly related
with the increase of the burning rate, as it will be explained afterwards analysing Figures 4.43 and
4.42, and the gas velocity at the burning surface is calculated as ug = ⇢⇢cgṙ , a higher result of the local
gas velocity is obtained at the areas where the burning rate is higher too.
10
Kohga, 2011
OpenFoam 2D
9
6
r (mm/s)
5
.
0
0 1 2 3 4 5 6 7 8
p (MPa)
AP mass fraction distributions for time equal to 0.02, 0.04 and 0.06 seconds are represented in
Figures 4.36, 4.37 and 4.38 respectively. The very first visible attribute is the local concentration of
AP in the cells close to burning surface which has not an uniform distribution along the tube unlike
what happens with HTPB (see Figures 4.39, 4.40 and 4.41). This behaviour may be due to the facility
that the first reaction (equation (4.44)) has to take place compared with second one (equation (4.45)).
Since the activation energy which the oxidiser needs to be decomposed is not really high, it could
be deduced that, as soon as condensed AP sublimates into AP in gas phase, the first reaction takes
place and most of the AP is decomposed in P1 (see equation (4.44)). However, binder combustion
reaction (equation (4.45)) not only presents a higher activation energy, but also needs the products
106
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
from the first reaction to take place. This could lead to the conclusion that when HTPB sublimates
to gas, it does not burn immediately and therefore, it has more time to spread along the gas domain
and accumulate higher concentration in the cells than AP.
The figures show how the content of AP and HTPB increases as the pressure decreases as happened
for gas velocity field. The increase of mass fraction of both binder and oxidiser with the fall of pressure
could be assigned to the reaction rates formulation. Equations (4.61) and (4.62) are function of the
pressure therefore, the higher the pressure, the higher reaction rate will be obtained. High values of
R1 and R2 mean that the consumption of both, AP and HTPB. Therefore, it could be said that high
pressures increase the reaction rate leading to a decrease of propellant mass fraction due to mass
fraction conservation equation (4.49).
When the pressure is high, AP in gaseous form is located at the centre of the tube, close to AP
in its condensed form, that means, at the same y-coordinate range values. However, the lower the
pressure, the higher the AP gas is spread increasing the amount of cells containing AP in its gaseous
form close to the burning surface. The increase of the diffusivity as the pressures falls down may
be due to the rise in gas velocity with the reduction of pressure which drags AP throughout the
domain. However, at 0.5 MPa, binder and oxidiser mass fractions present distributions with higher
concentrations of each propellant close to their condensed phases respectively. At low pressures, the
amount of propellant sublimating is very small. In addition, as has been already explained, the
reaction rates of R1 and R2 decrease with the pressure. Therefore, the combustion of AP and HTPB
close to the burning surface is reduced leading to an accumulation of both of them in those cells.
Since in this test, the AP and HTPB are placed in sandwich configuration, it is of high interest
to study the temperature, the burning rate and the mass fraction of both, binder and oxidiser, in
the y-coordinate. Therefore, the regression rate at the burning surface and the AP and HTPB mass
fractions along a cross section plane adjacent to the burning surface for a pressure of 7 MPa are
represented in Figures 4.42 and 4.43 respectively where 0 is the minimum value and 1 the maximum
one. The content of HTPB is mostly constant along y-direction however, the content of AP reaches
its minimum value in the lower or upper sides of the sandwich and its maximum at y = 30 µm where
the interface between AP and HTPB is located, having a symmetric distribution. The burning rate
profile (see Figure 4.42) is similar to that of AP mass fraction unlike the upper and lower layers of
the sandwich where the burning rate does not decrease as abruptly as AP mass fraction does. Both
magnitudes, burning rate and propellant mass fractions decrease with time. The reason why HTPB is
uniformly distributed and AP gas is localised close to its condensed phase has been already explained
before.
107
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.33: Velocity distribution of gas phase in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1
MPa, 0.5 MPa (above-below). 108
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.34: Velocity distribution of gas phase in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1
MPa, 0.5 MPa (above-below).
109
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.35: Velocity distribution of gas phase in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1
MPa, 0.5 MPa (above-below).
110
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.36: AP distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
111
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.37: AP distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
112
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.38: AP distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
113
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.39: HTPB distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
114
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.40: HTPB distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
115
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.41: HTPB distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
116
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
15
time = 0.015 s
14 time = 0.065 s
13
12
11
10
9
r (mm/s)
7
.
0
0.00000 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 0.00014 0.00016
y (m)
Figure 4.42: Burning rate distribution along the burning surface at t = 0.015 s and t = 0.065 s at 70
MPa.
0.2
0.0014 AP time = 0.015 s HTPB time = 0.015 s
AP time = 0.065 s 0.18 HTPB time = 0.065 s
0.0012 0.16
0.14
0.001
Mass fraction
Mass fraction
0.12
0.0008
0.1
0.0006 0.08
0.06
0.0004
0.04
0.0002
0.02
0 0
0.00000 0.00004 0.00008 0.00012 0.00016 0.00000 0.00004 0.00008 0.00012 0.00016
y (m) y (m)
Figure 4.43: Mass fraction distribution after burning surface at t = 0.015 s and t = 0.065 s at 7 MPa
for AP (left) and HTPB (right).
Finally, and before finishing this subsection, the results related to the global temperature and
the heat flux analysis are detailed. A global variable temperature which takes the value of the
condensed phase temperature or the gas temperature depending on the region considered is defined.
This variable is firstly depicted in the x-direction for a pressure range from 0.5 to 7 MPa (see Figure
4.44, left) for time equal to 0.04 s. As observed, there is an abrupt raise of the temperature close to
the burning surface and afterwards it remains stable. The value of constant temperature in gas phase
increases slightly with the fall of pressure fixed at the right boundary. This phenomenon could take
place because, as it was already discussed, the fall of pressure provokes a decrease of the reaction rates
and therefore, there is a decrement in the combustion of AP and HTPB. As a consequence of the
117
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
loss of strength of AP and HTPB combustion reactions, the concentration of the reaction products
(HCl, CO2 , H2 O, N2 ) decreases mildly too. The slight decrement of the reaction products mass
fractions is enough to make the density decrease its value provoking an increment of the temperature.
Moreover, the effect of the temperature decrement with the raise of pressure can be also seen in
Figure 4.44 (right), where the distributions in the y-direction of the condensed phase temperature at
the cells adjacent to the burning surface when time is equal to 0.03 s are plotted for the each pressure
of the considered range.
Its distribution has a shape similar to the one of the regression rate at the burning surface plotted
in Figure 4.42. This behaviour is exactly the one expected since the higher temperature the gas phase
is, the higher will be the heat transferred to the solid along the burning surface.
Global temperature is depicted in colour diagrams as in Figures 4.45, 4.46 and 4.47 for times 0.02,
0.04 and 0.06 seconds. The behaviour observed in Figure 4.44 is repeated, temperature diagrams show
how the decrement of pressure leads to a raise of the temperature far from the burning surface. In
addition, it can be also seen how for low pressures, the temperature in the condensed phase increases
more gradually than for higher ones. The global temperature variable is used afterwards to obtain
the heat flux along the tube. Heat flux has been calculated by considering the heat conductivities of
the different zones and multiplying by the temperature gradient. As expected, the main transference
of heat is located at the burning surface where the solid changes its phase from condensed to gas (see
Figures 4.48, 4.49 and 4.50). The remaining parts of the tube do not present a significant value of
heat flux since in that volume the temperature remains mostly constant as it is already explained
when analysing Figures 4.45, 4.46 and 4.47. Moreover, the greatest amount of energy is obtained in
the central part of the sandwich where the change of phase of AP from condensed into gas phase and
its combustion happen. The asymmetry encountered in the figures is due to the interpolation done
by the software Paraview to obtain the point values in order to provide a smooth result. Cell values
results can be found in Appendix.
4000 2000
7 MPa 7 MPa
3500 5 MPa 1800 5 MPa
3 MPa 3 MPa
2 MPa 1600 2 MPa
3000
1 MPa 1400 1 MPa
Temperature (K)
Temperature (K)
1500 800
600
1000
400
500 200
0 0
0 200 400 600 800 1000 0 40 80 120 160
x (µm) y (µm)
Figure 4.44: Temperature distribution for t = 0.04 s along x-direction (left) and temperature in
condensed phase before burning surface at t = 0.03 s for 7 MPa to 0.5 MPa (right).
118
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.45: Temperature distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below).
119
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.46: Temperature distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below).
120
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.47: Temperature distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5
MPa (above-below).
121
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.48: Heat flux distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
122
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.49: Heat flux distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
123
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.50: Heat flux distribution in t = 0.06 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
124
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
125
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
is applied behaves exactly as 7 MPa fixed boundary one until time equal to 0.005 s. After this
time, although its burning surface is located at the same position as 7 MPa one, AP mass fraction
distribution at this point raises due to the increase of the diffusion of AP with the fall of pressure.
126
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
8
Pressure ramp boundary
0.5 MPa boundary
7 MPa boundary
7
5
Pressure (MPa)
0
0.0045 0.005 0.0055 0.006 0.0065 0.007 0.0075 0.008
time (s)
0.515
Pressure ramp boundary
0.514 0.5 MPa boundary
0.513
0.512
0.511
0.51
0.509
0.508
0.507
Pressure (MPa)
0.506
0.505
0.504
0.503
0.502
0.501
0.5
0.499
0.498
0.497
0.496
0.495
0.0052 0.0054 0.0056 0.0058 0.006 0.0062 0.0064 0.0066 0.0068 0.007
time (s)
127
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.54: AP mass fraction distribution in t = 0.005 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
Figure 4.55: AP mass fraction distribution in t = 0.01 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
128
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.56: AP mass fraction distribution in t = 0.02 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
Figure 4.57: AP mass fraction distribution in t = 0.005 s, t = 0.0051 s and t = 0.0052 s (above-below).
129
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
The behaviour of HTPB in gas phase in case of considering a pressure ramp boundary condition,
as it happens for fixed pressure boundary condition, it is completely opposed to the one of AP. Unlike
AP which is mostly close to the burning surface, HTPB in gaseous form is uniformly distributed
along the gas domain. Since the activation energy which the oxidiser needs to be decomposed is not
really high, as soon as condensed AP sublimates to AP in gas phase, the first reaction takes place and
most of the AP is decomposed in reaction products. However, binder combustion reaction, not only
presents a higher activation energy, but also needs the products from the first reaction to take place.
This could lead to the conclusion that when HTPB sublimates to gas, it does not burn immediately
and therefore, it has enough time to spread along the control volume.
By observing HTPB mass fraction distributions from Figures 4.58, 4.59 and 4.60 the same con-
clusion about the burning rate is obtained. For each time considered, the position of the burning
surface when a fixed boundary condition of 7 MPa is applied is closer to the beginning of the con-
densed phase domain which means that the propellant is burning quicker at high pressures. When
time is 0.005 s, the position of the burning surface is the same for fixed 7 MPa and pressure ramp
boundary conditions since during this time the same value of the pressure is applied at the end of
the tube. However, at 0.01 s the pressure ramp test has already suffered the depression at the end of
the tube and the pressure is established to 0.5 MPa. Therefore, the position of the burning surface
in x-direction for pressure ramp is further from the left side than for 7 MPa fixed boundary test but
closer than for fixed 0.5 MPa boundary condition case, as expected.
As happens with AP mass fraction, HTPB mass fraction increases with fall of pressure. As
explained in subsection 4.2.3.2, the increase of mass fraction of both binder and oxidiser when the
pressure decreases could be assigned to the reaction rates formulation. Since equations (4.61) and
(4.62) are functions of the pressure, the higher the pressure, the higher reaction rate will be obtained.
High values of R1 and R2 mean that the consumption of both, AP and HTPB. Therefore, due to
mass fraction conservation equation (4.49), it could be said that high pressures increase the reaction
rate leading to a decrease of propellant mass fraction in the domain.
It is of high interest to remark the phenomenon observed in Figure 4.63 as result of plotting the
mass fraction of HTPB for the same specific point in the geometry used in Figure 4.62 (P1 ). From
0.0045 s to 0.005 s, the behaviour of the pressure ramp boundary is exactly the same as that of 7
MPa fixed pressure boundary. However, together with the sharp fall of pressure, comes the strong
raise of HTPB mass fraction. Comparing this figure with 4.53, the similarities are easily visible. The
depressurisation at the boundary provokes a fall of pressure close to the burning surface even below
0.5 MPa. This decrement of pressure will lead to a raise of the value of HTPB mass fraction over
the one obtained for 0.5 MPa fixed boundary test. This is the expected behaviour since it has been
already concluded that the fall of pressure brings an increase of propellant in gas phase along the
geometry. The mass fraction values obtained for both tests will come closer as the pressure for the
pressure ramp boundary condition test becomes closer to 0.5 MPa. However, due to the fact that
the mass fraction of HTPB present in the gas phase before reaching equilibrium is higher for the
depressurisation test than for fixed pressure test, the mass fraction of binder obtained along the tube
is higher when applying transient combustion although the final pressure of both tests is the same.
130
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.58: HTPB mass fraction distribution in t = 0.005 for 0.5 MPa, the pressure ramp and 7
MPa boundary conditions (above-below).
Figure 4.59: HTPB mass fraction distribution in t = 0.01 for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
131
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.60: HTPB mass fraction distribution in t = 0.02 for for 0.5 MPa, the pressure ramp and 7
MPa boundary conditions (above-below).
Figure 4.61: HTPB mass fraction distribution in t = 0.005 s, t = 0.0051 s and t = 0.0052 s (above-
below).
132
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
2.5
AP MassFraction ramp boundary
AP MassFraction 0.5 MPa boundary
AP MassFraction 7 MPa boundary
1.5
Mass Fraction (%)
0.5
0
0.0045 0.005 0.0055 0.006 0.0065 0.007 0.0075 0.008
time (s)
Figure 4.62: AP mass fraction (%) distribution against time for P1 = (515 · 10 6
, 85·10 6
, 85·10 6
).
50
HTPB MassFraction ramp boundary
HTPB MassFraction 0.5 MPa boundary
HTPB MassFraction 7 MPa boundary
40
30
Mass Fraction (%)
20
10
0
0.0045 0.005 0.0055 0.006 0.0065 0.007 0.0075 0.008
time (s)
Figure 4.63: HTPB mass fraction (%) distribution against time for P1 = (515·10 6
, 85·10 6
, 85·10 6
).
As in subsection 4.2.3.2, global temperature and heat flux distribution are depicted. Firstly, the
temperature distribution obtained when applying the ramp pressure boundary condition are shown
in Figure 4.64. The immediate effect of the sharp fall of pressure at the end of the tube is the
decrease of the temperature is this location as expected due to ideal gas law as it is visible in time
equal to 0.0051 s. Afterwards, temperature raises again remaining mostly stable in all gas domain
(see Figures 4.64, 4.65 and 4.66). As can be seen, the highest value of temperature reached at the
end of the tube results when pressure applied at the right boundary is that obtained when fixing the
133
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
pressure at the right boundary to 0.5 MPa followed by the depressurisation test and finally, the 7
MPa fixed boundary test. The raise of the temperature field at the right side of of the tube when the
pressure at the boundary falls down has been already explained in section 4.2.3.2. Since the reaction
rates depend on the pressure, the lower the pressure, the lower reaction rate, leading to a decrease
of the mass fractions of the reaction products with the consequent decrease of density and therefore,
an increment of temperature. The temperature obtained when pressure is fixed to 0.5 MPa at the
right boundary is slightly higher than that in the depressurisation test. This happens because before
applying the pressure ramp at the right boundary, the pressure is fixed to 7 MPa and therefore, its
density is higher than in 0.5 MPa test. The difference between the temperature obtained for 0.5
MPa fixed pressure and depressurisation boundary condition tests is narrowed with the time since
the density of both, depressurisation and 0.5 fixed boundary condition tests, comes closer too. In
addition, at low pressures, the evolution of the temperature along the control volume is softer than
at high ones. This effect occurs immediately after applying depressurisation at the boundary (see
Figure 4.64) and it is maintained with time as visible in Figures 4.65 and 4.66.
As expected, the main transference of heat is located at the burning surface where the propellant
changes its phase from condensed into gas (see Figures 4.67, 4.68 and 4.69). The remaining parts of
the tube do not present a significant value of heat flux since in that volume the temperature remains
mostly constant as it is already explained analysing Figures 4.64, 4.65 and 4.66. In addition, the
greatest amount of energy is obtained in the central part of the sandwich where the change of phase
of AP from condensed to gas phase and its combustion happen. These are the same conclusions
than those obtained in subsection 4.2.3.2 therefore, it can be said that the depressurisation of the
ramp it has not a significant impact in the heat flux distribution. However, the magnitude of heat
flux obtained for the pressure ramp test is smaller than that obtained when fixing the pressure to 7
MPa. This may be due to the smoother evolution of temperature along the domain produced at low
pressure against the sharp increase of temperature at high ones.
134
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.65: Temperature distribution in t = 0.01 s for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
Figure 4.66: Temperature distribution in t = 0.02 s for for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
135
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.67: Heat flux distribution for ramp pressure boundary condition in t = 0.005 s, t = 0.0051
s and t = 0.0052 s (above-below).
Figure 4.68: Heat flux distribution in t = 0.01 s for 0.5 MPa, the pressure ramp and 7 MPa boundary
conditions (above-below).
136
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
Figure 4.69: Heat flux distribution in t = 0.02 s for for 0.5 MPa, the pressure ramp and 7 MPa
boundary conditions (above-below).
Finally, burning rate and temperature are plotted along the burning surface for times equal to
0.0021 s (before the pressure ramp is applied) and 0.0081 (after the pressure ramp boundary condition
is applied). As in Figure 4.32 from subsection 4.2.3.1, it is observed how the burning rate decreases
with pressure. Before applying the pressure ramp, the burning rate obtained at the burning surface is
the same as the test cases in which the pressure is fixed to 7 MPa (see Figure 4.70) being the difference
with the one for fixed 0.5 MPa very high. When the pressure ramp boundary condition is applied, the
burning rate decreases reaching a value even lower than that of fixed 0.5 MPa boundary condition (see
Figure 4.71). This phenomenon occurs because applying a sharp depressurisation at the boundary
leads to obtaining lower values of pressure closer to the burning surface than those obtained when
fixing the pressure to 0.5 MPa. Since burning rate is directly related with the pressure, together with
a decrease of pressure below 0.5 MPa a fall of the burning rate below the value obtained for 0.5 MPa
fixed pressure test will come. As the difference between the pressure obtained after depressurisation
and the pressure obtained with fixed 0.5 MPa boundary is not really high, the gap between both
burning rates is mild too.
Regarding the temperature at the burning surface, before applying the pressure ramp at the right
boundary, the maximum value of temperature is that obtained when fixing the pressure of the right
boundary to 0.5 MPa (see Figure 4.72). When the depressurisation is applied, the temperature of
the burning surface for this test increases (see Figure 4.73) until reaching a value very close to that
obtained for 0.5 MPa fixed pressure but lower. As a result, the distribution of the surface temperature
when applying depressurisation is located between those obtained for 0.5 and 7 MPa fixed pressure
boundary conditions tests. The same behaviour is noticed when analysing global temperature variable
distribution in Figures 4.65 and 4.66.
137
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
15
Ramp boundary
0.5 MPa boundary
14 7 MPa boundary
13
12
11
10
9
r (mm/s)
7
.
0
0 2x10-5 4x10-5 6x10-5 8x10-5 0.0001 0.00012 0.00014 0.00016
y (m)
Figure 4.70: Burning rate distribution along the burning surface at t = 0.021 s for fixed 7 MPa, fixed
0.5 MPa and pressure ramp boundary conditions.
15
Ramp boundary
0.5 MPa boundary
14 7 MPa boundary
13
12
11
10
9
r (mm/s)
7
.
0
0 2x10-5 4x10-5 6x10-5 8x10-5 0.0001 0.00012 0.00014 0.00016
y (m)
Figure 4.71: Burning rate distribution along the burning surface at t = 0.081 s for fixed 7 MPa, fixed
0.5 MPa and pressure ramp boundary conditions.
138
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
1400
Ramp boundary
0.5 MPa boundary
7 MPa boundary
1200
1000
800
Ts (K)
600
400
200
0
0 2x10-5 4x10-5 6x10-5 8x10-5 0.0001 0.00012 0.00014 0.00016
y (m)
Figure 4.72: Temperature at the burning surface at t = 0.021 s for fixed 7 MPa, fixed 0.5 MPa and
pressure ramp boundary conditions.
2000
Ramp boundary
0.5 MPa boundary
7 MPa boundary
1800
1600
1400
1200
Ts (K)
1000
800
600
400
200
0
0 2x10-5 4x10-5 6x10-5 8x10-5 0.0001 0.00012 0.00014 0.00016
y (m)
Figure 4.73: Temperature at the burning surface at t = 0.081 s for fixed 7 MPa, fixed 0.5 MPa and
pressure ramp boundary conditions.
139
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
4.2.4 Conclusions
A code to study the multidimensional modelling of unsteady combustion of AP/HTPB propellants has
been developed in OpenFoam. The code employs Rusanov numerical scheme to solve the homogeneous
system of differential equations and afterwards, a first-order Euler methods to solve the ordinary
differential equations including the source terms. The kinetics of the model is described by considering
firstly, a change of phase of the solid propellant from condensed to gas and secondly, a reduced
chemistry scheme which defines simplified chemical reactions to represent the combustion itself. To
couple both processes, mass and species conservation, as well as temperature continuity are imposed
in the burning surface in which the burning rate will represent a key factor. Moreover, an energy
balance is also applied at the burning surface which represents that heat that the gas transfers to
the burning surface is invested firstly, in raising the surface temperature to produce the phase change
and secondly, in warming the condensed phase by conduction.
The first step done to validate the code is to perform a test which simulates the problem of
combustion of a homogeneous AP/HTPB strand in a Crawford bomb strand burner. The simulation
assumes that the tube is insulated and there is no transference of mass or energy in any boundary
despite the right cross-section which open to atmosphere. The simulation has been performed twice
by fixing the pressure in the right cross-section to 4 MPa and 7 MPa. The values obtained for
the burning rate have been compared against experimental data from INTA facilities [86] obtaining
positive results. Once this analysis has been concluded, the configuration of the binder and oxidiser
is changed from a homogeneous strand to a geometry in which a layer of AP is located in between two
layers of HTPB, called sandwich configuration. In this case the condensed control volume presents
two different regions, one containing AP and other of HTPB, which have different thermophysical
properties. As in the Crawford Strand Burner test, the simulation assumes that the tube is insulated
and there is no transference of mass or energy in any boundary despite the right one which is open
to atmosphere. The pressure at the right side of the tube has been fixed for this test and a range of
pressures between 0.5 MPa to 7 MPa have been tested. Finally, the same geometry has been used
to study the combustion of AP/HTPB propellant in sandwich configuration when a depressurisation
ramp is performed. Instead of having a fixed pressure in the open side of the tube, a ramp boundary
condition have been applied in which, in a very short interval of time, the pressure decreases from
7MPa to 0.5 MPa.
Regarding the first test considered in this subsection, it can be said that the results obtained in the
simulation of the burning of an AP/HTPB homogeneous strand in a Crawford Bomb Strand Burner
at 4 MPa and 7 MPa for several initial temperatures agree with the experimental ones. The absolute
error present a range of values between 0.4 mm/s to 1.1 mm/s, being in the order of magnitude of
the experimental error. However, the calculated relative error around 36% for 220.15K and 7 MPa
and decreasing with the decrement of pressure and the increment of temperature to values around
10% at 4 MPa for 294.15 K and 347.15 K. As explained throughout subsection 4.2.3.1, the high
magnitude obtained for the relative error is due to the low value of experimental burning rate for low
temperatures. Therefore, although the absolute error remains around 1 mm/s, or even 0.4 mm/s for
294.15 K and 347.15 K, the decrease of the experimental burning rate raises the relative error.
Once the code have been validated against experimental data, the combustion of an AP/HTPB
propellant in sandwich configuration fixing the pressure at the open side of the tube has been per-
formed for a pressure range from 0.5 to 7 MPa. Distributions of the main parameters of interest such
as, burning rate, gas velocity, mass fraction of AP and HTPB distributions, global temperature and
heat flux have been studied. The main conclusions extracted from the analysis can be summarised
in the following lines:
1. The burning rate results obtained with the developed two-dimensional model of this work
applied to the combustion of AP/HTPB composite propellant in sandwich configuration agree
well with the experimental ones from Kohga [51] for the considered pressure range. The relative
error obtained from the results is less than 10% for a pressure range between 7 MPa and 2 MPa
and below 25% between 0.5 MPa and 2 MPa. However, the absolute error remains constant
140
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
with a value around 0.5 mm/s, being in the order of magnitude of the experimental error. As
it has been already explained, the increase in the relative error with the fall of pressure is due
to the low values of burning rate obtained for the pressure range of 0.5 MPa to 2 MPa.
2. Despite the raise of the burning rate with the pressure, the higher the pressure, the lower the
gas velocity has been obtained. This behaviour may be due to the increase of the gas density
with the pressure. However, it is also observed that an increase of AP mass fraction close to
the burning surface leads to a local raise of the gas velocity. This behaviour is attributed to
the relationship between burning rate and gas velocity at the burning surface (4.54).
3. The behaviour of AP and HTPB mass fraction distribution along the gas domain is completely
opposed. AP mass fraction presents a local concentration in the cells close to burning surface
which has not an uniform distribution along the tube like it happens for HTPB. This behaviour
may be due to the facility that the first reaction has to take place compared with second one.
Since the activation energy that the oxidiser needs to be decomposed is not high, compared with
the one from HTPB, it could be deduced that, as soon as condensed AP sublimates to AP in gas
phase, the first reaction takes place and most of the AP is decomposed in reaction products.
However, binder combustion reaction, not only presents a higher activation energy, but also
needs the products from the first reaction to take place. This could lead to the conclusion that
when HTPB sublimates to gas, it does not burn immediately and therefore, it has enough time
to spread along the gas domain making possible for the cells far from the burning surface to
accumulate higher concentrations of binder than oxidiser.
4. The content of AP and HTPB in gas phase increases with the pressure decrement as happened
for gas velocity field. The increase of mass fraction of both binder and oxidiser with the fall of
pressure could be assigned to the dependency of reaction rates expressions with the pressure.
Therefore, the higher the pressure, the higher reaction rate will be obtained. A large value of
reaction rates means that the consumption of both, binder and oxidiser, will be high. Therefore,
due to mass fraction conservation equation, it could be said that high pressures increase the
reaction rate leading to a consumption of propellant mass fraction.
5. The diffusivity of AP increases at low pressures. The rise in gas velocity at low pressure could
increase the movement of AP throughout the domain. However, it is observed that, at 0.5
MPa, both binder and oxidiser mass fractions present a nonuniform distribution with higher
concentrations in the cells of gas adjacent to their condensed phases respectively. This could
happen because at low pressures, the degradation of the propellant is less intense and, in
addition, R1 and R2 reaction rates of decrease. Therefore, the combustion of AP and HTPB
in the cells close to the burning surface is reduced increasing as a result the mass fractions of
binder and oxidiser in these cells.
6. Analysing the gas species distribution along a cross section plane close to the burning surface
for a pressure of 7 MPa provides an almost constant content of HTPB. However, the content
of AP reaches its minimum value in the lower or upper sides of the sandwich and its maximum
at y = 30 µm where the interface between AP and HTPB is located, having a symmetric
distribution. The burning rate profile is similar to that obtained for AP mass fraction despite
for upper and lower layers of the sandwich where the burning rate does not decrease as abruptly
as AP mass fraction does but softer due to the presence of HTPB which contribute to enhance
the burning rate.
7. The value of constant temperature in gas phase increases slightly with the fall of pressure
fixed at the right boundary. This phenomenon could take place because, as have been already
discussed, the fall of pressure provokes a decrease of reaction rates and therefore, there is a
decrement in the combustion of AP and HTPB. As consequence of the loss of strength of AP and
HTPB combustion reactions, the concentration of the reaction products (HCl, CO2 , H2 O N2 )
141
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
decreases mildly too. The slight decrement of the reaction products mass fractions is enough
to make the density decrease its value provoking an increment of the temperature.
8. The main transference of heat is located at the burning surface where the solid changes its
phase from condensed to gas phase. In the rest of the tube, no significant change have been
appreciated in the heat flux. This behaviour is the one expected since the temperature remains
mostly stable all along the condensed phase and the gas phase presenting the only step in the
burning surface.
Finally, the combustion simulation of AP/HTPB composite propellant in sandwich configuration has
been performed applying a depressurisation ramp at the open side of the tube in order to see how it
influences in the main variables of interest. The conclusions obtained are the following:
1. The pressure ramp induced at the right boundary leads to a decrease of the pressure value
below 0.5 MPa at the cells located in the middle of the tube. Therefore, it can be said, that
when the end of the chamber is put down to a rapid depressurisation the pressure wave evolves
from the end to the beginning of the gas control volume and comes back causing a depression
under 0.5 MPa.
2. As has been remarked in the previous conclusions, the higher the pressure, the higher the
burning rate and therefore, the quicker the burning surface moves to the left side of the tube.
As expected, the x-coordinate of the burning surface has an intermediate value between those
of 7 MPa and 0.5 MPa tests.
3. The fall of pressure provoked by the depressurisation at the boundary leads to an increase of
both AP mass fraction and its diffusivity. As expected, the distribution of oxidiser mass fraction
obtained before applying the ramp boundary condition behaves exactly as that obtained for 7
MPa fixed boundary one. After applying the pressure ramp, AP mass fraction distribution
raises due to the increase of the AP diffusivity with the fall of pressure although the burning
surface is located at the same position as the one of 7 MPa test.
4. HTPB mass fraction also increases with the fall of pressure. Before the depressurisation is
applied, the distribution of HTPB is exactly the same as that obtained in the 7 MPa fixed
pressure boundary test. However, together with the sharp fall of pressure, comes a raise of
HTPB mass fraction. The depressurisation at the boundary provokes a fall of pressure close to
the burning surface even below 0.5 MPa. This decrement of pressure will lead to a raise of the
value of HTPB mass fraction over the one obtained for 0.5 MPa fixed boundary test. The mass
fraction values obtained for both tests will come closer as the pressure for the pressure ramp
boundary condition test becomes closer to 0.5 MPa.
5. The immediate effect of the sharp fall of pressure at the end of the tube is the decrease of
the temperature at this location as could be expected by taking into account ideal gas law.
Afterwards, temperature raises again remaining mostly stable. The highest value of temperature
reached at the end of the tube is that obtained when fixing the pressure at the right boundary
to 0.5 MPa followed by the depressurisation test and finally, the value obtained for the 7 MPa
fixed boundary test. Since the reaction rates depend on pressure, the lower the pressure, the
lower reaction rate. This leads to a decrease of the reaction products mass fractions with the
consequent decrease of density and therefore, an increment of temperature. The temperature
obtained when pressure is fixed to 0.5 MPa at the right boundary is slightly higher than that of
the depressurisation test. This happens because before applying the pressure ramp at the right
boundary, the pressure is fixed to 7 MPa and therefore, its density is higher than in 0.5 MPa
test. The difference between the temperature obtained for the 0.5 MPa fixed pressure test and
the depressurisation test is narrowed with the time since the density of both, depressurisation
and 0.5 fixed boundary condition tests, comes closer too.
142
4.2. COMBUSTION OF COMPOSITE AP/HTPB PROPELLANT
6. The large value of heat transference is located at the burning surface where the propellant
changes its phase from condensed to gas. The greatest amount of energy is obtained in the
central part of the sandwich where degradation of condensed AP takes place. Since, the beha-
viour is the same than that obtained for the fixed pressure boundary condition, it can be said
that the depressurisation of the ramp has not a significant impact in the heat flux. However,
the magnitude of the heat flux obtained for the pressure ramp test is smaller than that obtained
when fixing the pressure to 7 MPa. This may be due to the smoother evolution of temperature
along the domain produced at low pressure against the sharp increase of temperature at high
ones.
7. The burning rate decreases with pressure. Before applying the pressure ramp, the burning
rate obtained at the burning surface is the same for fixed 7 MPa and pressure ramp boundary
conditions. When the pressure ramp boundary condition is applied, the burning rate decreases
reaching a value even lower than that obtained for the 0.5 MPa fixed boundary test. This phe-
nomenon occurs because applying a sharp depressurisation at the boundary leads to obtaining
lower values of pressure closer to the burning surface than those obtained by fixing the pressure
to 0.5 MPa. Since burning rate is directly related with the pressure, together with a decrease of
pressure below 0.5 MPa, a fall of the burning rate below the value obtained for 0.5 MPa fixed
pressure test will come.
8. When the depressurisation is applied, the temperature of the burning surface increases until
reaching a value very close to that obtained for 0.5 MPa but lower. As a result, when applying
depressurisation, the distribution for the surface temperature is located between those obtained
for 0.5 and 7 MPa fixed pressure tests. The same behaviour was encountered when analysing
the distribution of global temperature variable.
143
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Parameter Value
cF 1
cO 1.4
cDR1,DR2 3.4
cP 3.9
144
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
@ (⇢g !
u g) ! ⇣ ! !⌘
+ r· ⇢g u g ⌦ !
ug +pI = 0 (4.89)
@t
✓ ✓ ◆◆ XN
@ (⇢g Eg ) ! p
+ r· ⇢g !
u g Eg + = k g r2 Tg + Qg,j Rj (4.90)
@t ⇢g j=1
@ (⇢g Yi ) ! kg
+ r· (⇢g !
u g Yi ) = r2 Y i + i i = 1, 2, ..., n 1 (4.91)
@t g cvg
@ (⇢c Ec )
= k c r2 Tc (4.92)
@t
where the subscripts i and j stand for the ith-species from 1 to n-1 and for reaction number from
1 to N respectively. The specific heat at constant volume for a mixture of gases is defined with the
following expression,
n
X
cvg = Yi cvg,i (4.93)
i=0
Tg = Tc = Ts (4.94)
3. The continuity of mass must be ensured. When the temperature of the surface is higher than
the melting temperature (Tmelt ), the solid changes phase to gas. Otherwise the velocity in the
burning surface must be zero. Therefore, the following relationship is set,
⇢ ⇢c ṙ
ug = ⇢g , Ts > Tmelt (4.96)
0, Ts < Tmelt
The velocity at which the burning surface of the condensed phase changes into gas is called the burning
rate (ṙ). The formula that Tseng and Yang [72] use to describe the burning velocity is written as,
E a,c
Ea,c
R u T s ↵ c Ac e
R u Ts
2
ṙ = Qc
(4.98)
1 TT0s 2cc Ts
where Ac ,Ea,c , cs and Qc are the pre-exponential factor, activation energy, specific heat and
specific heat release of the condensed phase respectively and their values are provided in Tseng and
Yang [72]. Since the propellant is homogeneous, all points of the burning surface will progress at the
same velocity.
145
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
p = ⇢g R g T g (4.99)
being Rg the specific constant of gases which for a mixture of gases is as calculated as,
n
X
Rg = Ru Yi /Mg,i (4.100)
i=0
The modelling of the kinetic scheme detailed at the beginning of subsection 4.3.1 is done by
calculating firstly, the reaction rates of R1, R2 and R3 which are equations (4.85), (4.86) and (4.87),
and afterwards the rate of production or consumption of each specie due to the combustion process.
The reaction rates have an Arrhenius form however, opposite to the model used for AP/HTPB
composite propellant, they are not function of the gas pressure but the density, the temperature and
the mass fraction of the species involve in the following form,
⇣ ⌘
YF YO E1
R1 = D1 ⇢2g e RuTg
(4.101)
Mg,O
⇣ ⌘
YDR1 E2
R2 = D2 ⇢2g e RuTg
(4.102)
Mg,DR1
⇣ ⌘
YDR2 E3
R3 = D3 ⇢2g e RuTg
(4.103)
Mg,DR2
where D1 , D2 and D3 are the pre-exponential factors, and E1 , E2 and E3 are the activation
energies of the reaction rates.
Once the reaction rates are defined, and in order to solve equation (4.91), the following step is to
calculate the rate of production or consumption of each species involved in the combustion process,
F = R1 (4.104)
cO Mg,O
O = F (4.105)
cF Mg,F
DR1 = R2 (4.106)
DR2 = R3 ( F + O) (4.107)
P = R2 R3 (4.108)
146
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
The first test performed is the combustion of the propellant when the pressure in the open side
147
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
of the tube is fixed in a range from 1 to 9 MPa. The results obtained are compared with the ones
obtained by Tseng and Yang [72] and Aoki and Kubota [77].
In the second test the burning of the same strand when a pressure step boundary condition is
applied at the right side of the control volume is simulated. Firstly the pressure at the right boundary
is fixed to 9 MPa and afterwards, it decreases abruptly until reaching 1 MPa at the open side of the
tube as it is done in subsection 4.2.3.3.
The necessary thermophysical and chemical kinetic properties of double-base propellant considered
as input variables are taken from bibliography and detailed in Tables table 4.10, 4.11 and 4.12 .
148
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Parameter Value
Ys,F 0.1662
Ys,O 0.04601
Ys,DR1 0.4762
Ys,DR2 0
Ys,P 0
The heat of the chemical reactions in the combustion of the double-base propellant is not provided
in the bibliography as is done for AP/HTPB propellant. However, the heat of formation of each
component is known (see Table 4.12). Therefore, the heat of each jth-reaction is obtained as,
X X
Qg,j = h0f,products h0f,reactives (4.111)
To run the test, initial conditions of temperature need to be established. The tube is assumed
initialised therefore, the left side of the domain is set at 294 K (T1 ) and the right side at 1100 K
(T2 ) and the rest of the control volume will have the temperature determined by the linear equation
between these two values. Initial temperatures, as well as reference temperature and melting tem-
perature are summarised in Table 4.15. Melting temperature (Tmelt ) and reference temperature (T0 )
values are obtained from Wu et al. [71]. The value adopted for the temperature at the right side of
the control volume is justified throughout this subsection.
149
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
found in the literature. Aoki and Kubota [77] presented in their work, among other variables, the
burning rate versus the pressure for 12 propellants each one of them with a different content of
nitroglycerine and nitrogen dioxide. As can be seen in their Table 1, EC-1 has a mass fraction of
nitrogen dioxide of 0.466 being the closest that which Tseng and Yang [72] use in their study. Due
to the amount of initial necessary data provided their paper, this propellant has been the one used
to validate the model on double-base propellant.
Moreover, as has explained before, the control volume is assumed already initiated, therefore, in
order to study the effect that the initial temperature distribution has in the burning rate, the results
from Aoki and Kubota [77] and the ones from Tseng and Yang [72] are plotted together with the
burning rates obtained when final temperature is 950K, 1000K, 1100 K and 1300 K in Figure 4.74. As
can be seen, the temperature at which the domain is considered initiated has an effect in the burning
rate. The higher the final temperature is assumed, the better adjust the results of the calculated
burning rate to the results from Aoki and Kubota [77] for high and low pressures. However, when the
temperature is high, the results obtained differ from those of the literature in the range from 2 MPa
to 6 MPa. The quantitative values and their comparison with the results from Aoki and Kubota [77]
together with absolute and relative error depending is described in Tables 4.16, 4.17, 4.18 and 4.19.
The results obtained agree well with those from Aoki and Kubota [77] independently of the initial
temperature chosen for calculation . On the one hand, the absolute error is less than 1 mm/s in each
case except for high pressures where, as seen in Figure 4.74, the curves start to diverge reaching a
value around 1.7 mm/s at 950 K. On the other hand, the relative error present a high value compared
with the absolute one. This is due to the expression used to calculate it, as it as already explained
in subsection 4.2.3.1 when describing the burning rate results obtained for the combustion of an
AP/HTPB homogeneously distributed in a Crawford Strand Burner. When the initial temperature
chosen for the right boundary is high, the calculated burning rate agrees well with the experimental
ones for low and high pressures (see Tables 4.16, 4.17, 4.18 and 4.19).
Analysing in detail the quantitative results obtained in Tables 4.16, 4.17, 4.18 and 4.19 it can be
seen that when the pressure is 9 MPa, the higher the temperature, the lower the absolute error. The
same phenomenon is observed when running the test at 1.1 MPa. In case of having a temperature of
1300 K, the relative and absolute errors obtained at 9.1 MPa and 1.1 MPa, are less than 11% and 0.8
mm/s respectively. However, in this case, the relative error obtained for pressures range from 2 MPa
to 6 MPa increases being even 26 % for 2 MPa. The results obtained for 950 K and 1000 K are very
similar: around 14% of relative error for the results at 9.1 MPa, around 20 % for 1.1 MPa and less
than 9% for intermediate pressures, which means a very good agreement between the calculated and
literature results. When the temperature is increased to 1100 K, not only the relative error for high
pressure but also for low pressure is reduced to 10% and the for intermediate pressures, the relative
error is maintained under 9%. Therefore, the final temperature of 1100 K is the one chosen to present
the results hereafter.
The first variable to be analysed is the velocity field in the gas domain which is represented in
Figures 4.75, 4.76 and 4.77 for 0.01 s, 0.02 s and 0.05 s respectively. The first phenomenon observed is
the big difference between the burning rate for high pressures compared with the one of low pressures
since, for example for time equal to 0.05 (Figure 4.77), the burning surface for 9.1 MPa has arrived
almost to the beginning of the control volume, which means that most of the condensed phase has been
already consumed, while for 1.1 MPa, the burning surface is still at the half of the tube. Therefore,
it could be stated that the higher the pressure is, the quicker the condensed phase changes from solid
to gas form. However, the behaviour of the gas velocity is opposed to the burning rate, the higher the
pressure is, the lower gas velocity we have. This could be due, as explained in AP/HTPB subsection,
to the increment of the gas density with the pressure which leads to a reduction of the gas velocity
field as a result of solving equation (4.89). The velocity fields behaves differently depending on the
pressure fixed at the right boundary. There are 3 mm/s gap between the burning rate of 3.4 MPa
and the one of 9.1 MPa however, the difference in burning rate between 3.4 MPa and 1.1 MPa is
close to 4.5 mm/s (see Table 4.18). This behaviour affects on the first hand, to the position of the
burning surface in the x-coordinate, being almost the same for the first two time steps considered
150
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
and slightly different for 3.4 MPa at 0.05 s. On the other hand, this behaviour affects to how the
velocity of the gas phase increases with time. For the first range of pressures (3.4 MPa to 9.1 MPa)
the velocity magnitude remains almost constant for each of the time steps considered. However, for
2 MPa and 1.1 MPa fixed boundary tests, an increment with time is observed. This phenomenon
could take place because chemical reactions at low pressures need some time to start the reaction and
therefore, for the gas velocity to increase while for high pressures the reaction mechanism is activated
from the very beginning of the combustion process.
15
Calculated Tf = 950 K
Calculated Tf = 1000 K
14 Calculated Tf = 1100 K
Calculated Tf = 1300 K
13 Aoki and Kubota (1982)
Tseng and Yang (1994)
12
11
10
9
r (mm/s)
7
.
0
0 1 2 3 4 5 6 7 8 9 10
Pressure (MPa)
Aoki and
OpenFoam Absolute
Pressure Kubota Relative
Burning rate error
(MPa) Burning rate error (%)
(mm/s) (mm/s)
(mm/s)
9.1 9.447 11.098 1.651 14.879
5.99 8.142 8.917 0.776 8.698
3.4 6.222 6.523 0.301 4.622
2 4.4122 4.483 0.071 1.577
1.1 2.141 2.776 0.635 22.880
Table 4.16: Quantitative comparison of calculated and experimental burning rate for 950 K.
151
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Aoki and
OpenFoam Absolute
Pressure Kubota Relative
Burning rate error
(MPa) Burning rate error (%)
(mm/s) (mm/s)
(mm/s)
9.1 9.611 11.098 1.486 13.40
5.99 8.351 8.917 0.565 6.34
3.4 6.465 6.523 0.058 0.90
2 4.412 4.483 0.070 1.58
1.1 2.243 2.776 0.533 19.21
Table 4.17: Quantitative comparison of calculated and experimental burning rate for 1000 K.
Aoki and
OpenFoam Absolute
Pressure Kubota Relative
Burning rate error
(MPa) Burning rate error (%)
(mm/s) (mm/s)
(mm/s)
9.1 9.905 11.098 1.193 10.75
5.99 8.726 8.917 0.191 2.15
3.4 6.900 6.523 0.376 5.76
2 4.855 4.483 0.372 8.30
1.1 2.478 2.776 0.298 10.73
Table 4.18: Quantitative comparison of calculated and experimental burning rate for 1100 K.
Aoki and
OpenFoam Absolute
Pressure Kubota Relative
Burning rate error
(MPa) Burning rate error (%)
(mm/s) (mm/s)
(mm/s)
9.1 10.347 11.098 0.751 6.772
5.99 9.351 8.917 0.432 4.865
3.4 7.666 6.523 1.146 17.577
2 5.641 4.483 1.161 25.907
1.1 3.061 2.776 0.285 10.264
Table 4.19: Quantitative comparison of calculated and experimental burning rate for 1300 K.
152
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.75: Velocity distribution of gas phase in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa,
1.1 MPa (above-below).
153
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.76: Velocity distribution of gas phase in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa,
1.1 MPa (above-below).
154
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.77: Velocity distribution of gas phase in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa,
1.1 MPa (above-below).
Oxidiser and fuel mass fractions along the control volume are the next variables to be studied.
As can be observed in Figures 4.78, 4.79 and 4.80 independently of the considered pressure, the
highest amount of fuel and oxidiser is located close to to the burning surface decreasing the mass
fraction smoothly with x-coordinate. The maximum value of fuel mass fraction obtained for each
pressure decreases slightly with the time. Figures show how the increase of the time leads to a fall
of the fuel mass fraction at the end of the tube. The same behaviour is observed for the oxidiser
mass fraction, being the expected behaviour since both, fuel and oxidiser, react together in reaction
1 (see equation (4.101)) and need each other for the reaction to take place. Regarding the mass
fraction values obtained for each one of the species analysed, depending on the pressure at which
the boundary is fixed, both fuel and oxidiser present an increase in their maximum values of mass
fraction with the fall of pressure. The same effect was observed when analysing the combustion of
155
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
AP/HTPB composite propellant. Equations (4.101), (4.102) and (4.102), which correspond to the
chemical kinetics reaction rates, depend on the density of the gas. In consequence, the higher the
pressure, the higher the density and therefore, the higher reaction rate will be obtained. A raise of
reaction rates means an increment in the consumption of fuel, oxidiser and intermediate products DR1
and a decrement in opposite case. Therefore, since low pressure means a decrease in the combustion
of fuel and oxidiser, the mass fraction of both will increase.
Figure 4.78: Fuel distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
156
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.79: Fuel distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
157
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.80: Fuel distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below) (above-below).
158
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.81: Oxidiser distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
159
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.82: Oxidiser distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
160
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.83: Oxidiser distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
The mass fraction distributions of the intermediate species DR1 and DR1, as well as the mass
fraction of the final products, are depicted from Figures from 4.84 to 4.92. Figures 4.90, 4.91 and
4.92, show that final products (P) mass fraction increases with the pressure, which means that, the
higher the pressure, the higher reaction rates of R2, R3 or both there will be. On the one hand, the
higher the pressure, the higher reaction rate of reaction 2 and therefore, the higher amount of DR1
will react to become final products leading to a decrease of DR1 mass fraction with the pressure (see
Figures 4.84, 4.85 and 4.86). On the other hand, the mass fraction of the intermediate specie DR2
increases with the reaction of fuel and oxidiser (R1) and decreases when it changes to final products
(R3). Therefore, the amount of mass fraction in the control volume be function of the reaction rates
of both reactions (equations (4.101) and (4.103)). According to Figures 4.87, 4.88 and 4.89, the
161
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
mass fraction of DR2 increases with the pressure despite for the value obtained at 9.1 MPa which
suffers a slight decrement in comparison with that at 5.9 MPa. This could be due to the fact that,
since reaction rates get stronger with the increase of pressure, the reaction rate of R3 provokes a
consumption of DR2 high enough to decrease the content produced by R1.
Figure 4.84: DR1 distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
162
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.85: DR1 distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
163
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.86: DR1 distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
164
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.87: DR2 distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
165
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.88: DR2 distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
166
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.89: DR2 distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa
(above-below).
As general rule, it can be assumed that the final products mass fractions increase with the time.
However, when pressure is fixed between 9.1 MPa to 5.99 MPa, the mass fraction of the final products
decreases slightly from 0.01 s to 0.02 s (hereafter first slot) and presents the raise between 0.02 s and
0.05 s (hereafter second slot). This phenomenon could be explained taking into account that final
products are the result of R2 and R3 and R3 depends strongly on R1. Figures 4.84, 4.85 and 4.86
show how DR1 mass fraction is maintained between the first slot and decreased in the second one
which means that reaction 2 needs some time to increase their effect. Moreover, observing Figures
from 4.78 to 4.83, two behaviours can be observed. On the one hand, at high pressures, the mass
fraction of DR2 increases slightly in the first time slot while decreasing in the second one. On the
one hand, at high pressures, the mass fraction of DR2 increases slightly in the first time slot while
167
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
decreasing in the second one. On the other hand, at low pressures, DR2 mass fraction increases with
the time which means that, during the first slot, reaction 1 produces more intermediate species than
reaction 3 can consume and therefore, there is an increment of DR2 mass fraction. The same is
happening at high pressures for the first time slot. However, between 0.02 s and 0.05 s, the amount of
DR2 mass fraction when pressure is high suffers a decrement. Since during the second time slot and
the pressure fixed is high there is a strong increment of the final products mass, it can be concluded
that, reaction 3 increase its strength with pressure and time.
Figure 4.90: Final products distribution in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1
MPa (above-below).
168
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.91: Final products distribution in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1
MPa (above-below).
169
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.92: Final products distribution in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1
MPa (above-below).
A global temperature variable which takes the condensed phase or the gas value depending on
each domain is defined. As visible in Figures 4.93, 4.94 and 4.95, the location of the burning surface
depends on the pressure fixed at the right boundary. The first conclusion extracted from these figures
is the increment of the temperature with the pressure. Another effect very visible is the temperature
evolution along the control volume. The lower the pressure, the softer is the transition between
minimum and maximum temperature. However, this effect behaves in an opposite manner with the
time evolution. As the time passes, the jump between hot temperature at the gas phase and cold
temperature in the condensed phase becomes more abrupt which means that, the transference of heat
in the solid increases with the fall of pressure and decreases when the time passes. By studying the
already mentioned figures, it can be concluded that there is an increment of the temperature with
170
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
the time. This effect seems to be more visible at low than high pressures.
Figure 4.93: Temperature in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-
below).
171
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.94: Temperature in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-
below).
172
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.95: Temperature in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-
below).
173
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Finally, the heat flux is analysed by depicting its value in Figures 4.96, 4.97 and 4.98. The
behaviour of the heat transfer is in perfect consonance with all previous statements. The heat flux
increases with pressure and with time. The greatest amount of heat is transferred in the burning
surface between the solid and the gas.
Figure 4.96: Heat flux in t = 0.01 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-below).
174
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.97: Heat flux in t = 0.02 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-below).
175
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.98: Heat flux in t = 0.05 s at 9.1 MPa, 5.99 MPa, 3.4 MPa, 2 MPa, 1.1 MPa (above-below).
176
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
177
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
far from point one. Therefore, as happened for the 9.1 fixed boundary condition test, fluctuations
of pressure are less visible. The evolution of pressure with time when ramp boundary conditions are
established and its comparison with the 1.1 MPa and the 9.1 MPa boundary tests are represented in
Figures 4.106 and 4.107.
178
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
11
Pressure ramp boundary
1.1 MPa boundary
9.1 MPa boundary
10
7
Pressure (MPa)
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
3
Pressure ramp boundary
1.1 MPa boundary
2.8
2.6
2.4
2.2
1.8
Pressure (MPa)
1.6
1.4
1.2
0.8
0.6
0.4
0.2
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
179
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
11
Pressure ramp boundary
1.1 MPa boundary
9.1 MPa boundary
10
7
Pressure (MPa)
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
3
Pressure ramp boundary
1.1 MPa boundary
2.8
2.6
2.4
2.2
1.8
Pressure (MPa)
1.6
1.4
1.2
0.8
0.6
0.4
0.2
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
180
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
11
Pressure ramp boundary
1.1 MPa boundary
9.1 MPa boundary
10
7
Pressure (MPa)
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
3
Pressure ramp boundary
1.1 MPa boundary
2.8
2.6
2.4
2.2
1.8
Pressure (MPa)
1.6
1.4
1.2
0.8
0.6
0.4
0.2
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
181
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.106: Pressure distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary conditions
at t= 0.015 s (above-below).
Figure 4.107: Pressure distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary conditions
at t = 0.02 s (above-below).
182
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
The same analysis performed for the pressure, will be done for mass fractions of fuel and oxidiser,
intermediate species and products of chemical reaction in order to study how the depressurisation
affects to the combustion process. The first two distributions of species mass fraction analysed are
those corresponding to fuel and oxidiser (Figures 4.108 and 4.109) observing an increase of their
concentrations during the time after applying depressurisation. By studying the other species, it can
be seen how DR1, remains almost constant for the considered time slot (see Figure 4.84) increasing
slightly its concentration too. However, Figure 4.111 shows the decrement of DR2 mass fraction
from 0.01 s to 0.01005 s. Finally, products mass fraction distribution is plotted in Figure 4.112
where the fall of the specie concentration after applying depressurisation is observed too. How the
depressurisation affect the mass fraction of the gas species can be explained with the same argument
of subsection 4.3.3.1. The magnitude of the mixture density is smaller for low than for high pressures.
Since the reaction rates which represent the chemical kinetics of the combustion process depend on the
density, with the fall of pressure, the reaction rates decrease too. As a result, when depressurisation
is applied, not only less fuel and oxidiser react to form intermediate product DR2, but also less DR1
and DR2 change to become final products. The sublimation of condensed phase and the decrease of
reaction rates are the reason of the increment of fuel, oxidiser and intermediate products DR1. As a
result, there creation of new final products will fall down and part of them will be lost by the exit
of the tube concluding in a reduction of their mass fraction. The same phenomenon is suffered by
intermediate specie DR2 just after depressurisation until time is equal to 0.001004 s. At this time
pressure stabilised to 1.1 MPa (see Figure 4.99) and therefore, the mass fraction of DR2 remains
almost constant too. In order to study in detail the evolution of fuel and oxidiser mass fractions, not
only immediately after depressurisation, but also once pressure stability is reached, fuel and oxidiser
mass fraction distribution for pressure ramp test have been plotted together with the ones of 1.1 MPa
and 9.1 MPa fixed boundary condition test for 0.015 and 0.02 s in Figures 4.113, 4.114, 4.115 and
4.116. The maximum values of fuel and oxidiser mass fraction remain almost constant between the
considered interval of time. As is done for pressure, mass fraction distributions are also studied at
the three points of interest and the results are depicted from Figures 4.117 to 4.122. Figure 4.117 and
its zoom in time (Figure 4.118) show very clearly how the mass fraction of both, fuel and oxidiser for
depressurisation test behaves exactly the same. At the beginning their values are the same that the
ones obtained for 9.1 MPa fixed boundary condition however, when depressurisation is applied, both
mass fractions start to increase until reaching the ones obtained at 1.1 MPa fixed boundary condition.
As observed when analysing the pressure, fuel and oxidiser mass fraction values in the first point of
interest when pressure at the boundary is fixed to 1.1 MPa fluctuate slightly and being less visible
for point two and three. Afterwards, intermediate species DR1 and DR2, and final products are also
represented in these points obtaining the same evidences, both variables remain almost stable during
the interval between 0.01 s and 0.02 s for the fixed boundary condition test. In the case of applying a
ramp boundary condition, the mass fraction of DR1 starts with the value of 9.1 MPa fixed pressure
and it increases when depressurisation is applied until reaching the value of 1.1 MPa fixed pressure
test. However, DR2 mass fraction starts with the same percentage than the test of 9.1 MPa, and
with the depressurisation, instead of increasing its value, it falls down until reaching the amount of
DR2 obtained in 1.1 MPa fixed boundary condition test. Finally, to end this analysis, final products
mass fraction at the three considered points are plotted in three graphs collected in Figure 4.129. As
happened with the other species, no visible differences are observed during the selected interval of
time in case of fixing the pressure at the boundaries but a slight decrement for 9.1 MPa test. The
behaviour for the ramp test is exactly the one expected, the mass fraction of products start with a
value close to 9% which corresponds to 9.1 MPa fixed boundary condition. When depressurisation is
applied, chemical reaction lose their strength and as a result, less amount of products are obtained.
The value of 9% will fall gradually until reaching that obtained for 1.1 MPa fixed boundary condition.
Before analysing temperature distribution and heat flux, it is of high interest to study how the
surface temperature varies when strong pressure step is applied at the right boundary. As expec-
ted, when the depressurisation is applied, the temperature decreases sharply from the value of 9.1
fixed boundary condition test and starts to fluctuate slightly to finally reach stability around 650 K
183
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
184
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
185
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
186
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
187
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
188
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.113: Fuel mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.015 s (above-below).
Figure 4.114: Fuel mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.02 s (above-below).
189
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.115: Oxidiser mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t = 0.015 s (above-below).
Figure 4.116: Oxidiser mass fraction distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t = 0.02 s (above-below).
190
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
40
Fuel Pressure ramp boundary
Fuel 1.1 MPa boundary
Fuel 9.1 MPa boundary
Oxidizer Pressure ramp boundary
35 Oxidizer 1.1 MPa boundary
Oxidizer 9.1 MPa boundary
30
25
Mass fraction (%)
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
40
Fuel Pressure ramp boundary
Fuel 1.1 MPa boundary
Oxidizer Pressure ramp boundary
Oxidizer 1.1 MPa boundary
35
30
25
Mass fraction (%)
20
15
10
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
191
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
40
Fuel Pressure ramp boundary
Fuel 1.1 MPa boundary
Fuel 9.1 MPa boundary
Oxidizer Pressure ramp boundary
35 Oxidizer 1.1 MPa boundary
Oxidizer 9.1 MPa boundary
30
25
Mass fraction (%)
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
40
Fuel Pressure ramp boundary
Fuel 1.1 MPa boundary
Oxidizer Pressure ramp boundary
Oxidizer 1.1 MPa boundary
35
30
25
Mass fraction (%)
20
15
10
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
192
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
40
Fuel Pressure ramp boundary
Fuel 1.1 MPa boundary
Fuel 9.1 MPa boundary
Oxidizer Pressure ramp boundary
35 Oxidizer 1.1 MPa boundary
Oxidizer 9.1 MPa boundary
30
25
Mass fraction (%)
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
40
Fuel Pressure ramp boundary
Fuel 1.1 MPa boundary
Oxidizer Pressure ramp boundary
Oxidizer 1.1 MPa boundary
35
30
25
Mass fraction (%)
20
15
10
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
193
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
80
DR2 Pressure ramp boundary
75 DR2 1.1 MPa boundary
DR2 9.1 MPa boundary
DR1 Pressure ramp boundary
70 DR1 1.1 MPa boundary
DR1 9.1 MPa boundary
65
60
55
50
Mass fraction (%)
45
40
35
30
25
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
80
DR2 Pressure ramp boundary
75 DR2 1.1 MPa boundary
DR2 9.1 MPa boundary
DR1 Pressure ramp boundary
70 DR1 1.1 MPa boundary
DR1 9.1 MPa boundary
65
60
55
50
Mass fraction (%)
45
40
35
30
25
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
194
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
80
DR2 Pressure ramp boundary
75 DR2 1.1 MPa boundary
DR2 9.1 MPa boundary
DR1 Pressure ramp boundary
70 DR1 1.1 MPa boundary
DR1 9.1 MPa boundary
65
60
55
50
Mass fraction (%)
45
40
35
30
25
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
80
DR2 Pressure ramp boundary
75 DR2 1.1 MPa boundary
DR2 9.1 MPa boundary
DR1 Pressure ramp boundary
70 DR1 1.1 MPa boundary
DR1 9.1 MPa boundary
65
60
55
50
Mass fraction (%)
45
40
35
30
25
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
195
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
80
DR2 Pressure ramp boundary
75 DR2 1.1 MPa boundary
DR2 9.1 MPa boundary
DR1 Pressure ramp boundary
70 DR1 1.1 MPa boundary
DR1 9.1 MPa boundary
65
60
55
50
Mass fraction (%)
45
40
35
30
25
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
80
DR2 Pressure ramp boundary
75 DR2 1.1 MPa boundary
DR2 9.1 MPa boundary
DR1 Pressure ramp boundary
70 DR1 1.1 MPa boundary
DR1 9.1 MPa boundary
65
60
55
50
Mass fraction (%)
45
40
35
30
25
20
15
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
196
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
20 20
Products Pressure ramp boundary Products Pressure ramp boundary
Products 1.1 MPa boundary Products 1.1 MPa boundary
Products 9.1 MPa boundary Products 9.1 MPa boundary
18 18
16 16
14 14
12 12
Mass fraction (%)
8 8
6 6
4 4
2 2
0 0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02 0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s) Time (s)
20
Products Pressure ramp boundary
Products 1.1 MPa boundary
Products 9.1 MPa boundary
18
16
14
12
Mass fraction (%)
10
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
Figure 4.129: Products mass fraction distribution for P1 (above-left), P2 (above-right) and P3
(below).
1000
Pressure ramp boundary
9.1 MPa boundary
950
900
850
800
Ts (K)
750
700
650
600
550
500
0.01000 0.01001 0.01002 0.01003 0.01004 0.01005 0.01006 0.01007 0.01008 0.01009 0.01010
time (s)
197
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Finally, global temperature and heat flux variables are analysed. The variable global temperature,
represents the temperature from condensed to gas domain and is depicted by using colour diagrams
in Figures 4.131, 4.132 and 4.133. The first one of the three figures represents the evolution of the
temperature just after the pressure ramp is applied to the right boundary. The other two figures
compare the distribution of temperature obtained in each one of the three test considered in this
subsection: 1.1 MPa and 9.1 MPa fixed boundary condition and ramp boundary condition. Figure
4.131 shows how the depressurisation provokes immediately a diffusion of the temperature along the
condensed phase, phenomenon already seen in the previous subsection from Figure 4.93 to Figure
4.95. This behaviour can be easily seen thanks to Figures 4.132 and 4.133. The higher the pressure,
the sharper will be the transition between cold condensed and hot gas temperatures.
Another effect that can be appreciated is the sharp fall of the maximum temperature (around 400
K) immediately after the time in which the pressure ramp is started; afterwards the value decrease
smoothly. The differences obtained in the maximum values and the position evolution of the burning
surfaces for each test is stated very clear in Figures 4.132 and 4.133. The burning surface for the
depressurisation test is located between those of 1.1 MPa and 9.1 MPa fixed boundary condition
and the maximum temperature is in between the ones obtained for these two tests too. In order
to observe in detail the behaviour of the temperature, as done for pressure and for mass fraction
distributions, the evolution of the temperature with the time at the three chosen three points has
been represented from Figures 4.134 to Figure 4.139. As happened for the previous variables studied,
the temperature for 1.1 MPa fixed boundary condition test measured at the point close to the burning
surface present some fluctuations which decrease with the x-coordinate. However, and as was observed
when analysing the pressure, this phenomenon does not take place when ramp boundary condition is
applied. As has been already explained, the boundary surface for the pressure ramp test when time
is equal to 0.01 s is located further from the one of 1.1 MPa fixed boundary condition since before
the depressurisation is applied the pressure at the boundary is fixed 9.1 MPa. Figures 4.134, 4.136
and 4.138 show the same conclusion already extracted from Figures 4.132 and 4.133, the temperature
suffers a sharp decrease when depressurisation is applied until reaching a value close to the one of 1.1
MPa boundary condition test.
By plotting the temperature from t = 0.01 s to t = 0.01001, the behaviour of the temperature
just when depressurisation is applied can be studied in detail (see Figures 4.135, 4.137 and 4.139).
As happened for the rest of the variables analysed, the temperature falls abruptly due to the pressure
ramp at the end of the tube. After the fall, the temperature will fluctuate until t = 0.01001 from
which it will experiment a soft increase until reaching the same value as the one obtained for the 1.1
MPa fixed boundary condition test. The behaviour is exactly the same for the three points but the
fluctuations once the stability temperature is reached are more visible in the first point.
The last variable to be studied is the heat flux along the control volume which presents the
expected behaviour after studying the rest of the variables. Figures 4.141 and 4.142 show how the
maximum value obtained for heat flux when applying the pressure ramp test is in between the ones
obtained for 1.1 MPa and 9.1 MPa fixed boundary condition. In addition, the light blue of the
condensed phase, which indicates a transfer of heat in the condensed phase due to the progressive
warming from the left side of the tube to the burning surface, is more visible for the distribution
obtained in the ramp test than the one of 9.1 MPa but less than for 1.1 MPa. The fluctuation of the
heat flux just after depressurisation is very well observed in Figure 4.140. The strong fall of pressure
at the exit provokes a decrease of the heat flux between 0.01 s and 0.0101 s. After that interval of
time, its maximum fluctuates until reaching a magnitude closer to the one obtained for 1.1 MPa fixed
boundary condition test.
198
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
199
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.132: Global temperature distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.015 s (above-below).
Figure 4.133: Global temperature distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary
conditions at t= 0.02 s (above-below).
200
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
1800
Pressure ramp boundary
1700 1.1 MPa boundary
9.1 MPa boundary
1600
1500
1400
1300
1200
1100
Temperature (K)
1000
900
800
700
600
500
400
300
200
100
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
1800
Pressure ramp boundary
1700 1.1 MPa boundary
1600
1500
1400
1300
1200
1100
Temperature (K)
1000
900
800
700
600
500
400
300
200
100
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
201
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
1800
Pressure ramp boundary
1700 1.1 MPa boundary
9.1 MPa boundary
1600
1500
1400
1300
1200
1100
Temperature (K)
1000
900
800
700
600
500
400
300
200
100
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
1800
Pressure ramp boundary
1700 1.1 MPa boundary
1600
1500
1400
1300
1200
1100
Temperature (K)
1000
900
800
700
600
500
400
300
200
100
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
202
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
1800
Pressure ramp boundary
1700 1.1 MPa boundary
9.1 MPa boundary
1600
1500
1400
1300
1200
1100
Temperature (K)
1000
900
800
700
600
500
400
300
200
100
0
0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.02
Time (s)
1800
Pressure ramp boundary
1700 1.1 MPa boundary
1600
1500
1400
1300
1200
1100
Temperature (K)
1000
900
800
700
600
500
400
300
200
100
0
0.01 0.0101 0.0102 0.0103 0.0104 0.0105 0.0106 0.0107 0.0108 0.0109 0.011
Time (s)
203
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
Figure 4.140: Heat flux distribution for t = 0.01000001, t = 0.01001 s, t = 0.01002, t = 0.01003, t
= 0.001004 s and t = 0.01005 s (above-below).
204
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
Figure 4.141: Heat flux distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary conditions
at t= 0.015 s (above-below).
Figure 4.142: Heat flux distribution for fixed 1.1 MPa, ramp and fixed 9.1 MPa boundary conditions
at t= 0.02 s (above-below).
205
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
4.3.4 Conclusions
The code to study the multidimensional modelling of unsteady combustion of AP/HTPB propellant,
which uses Rusanov numerical scheme and first-order Euler method to solve the homogeneous differ-
ential equation system and source terms respectively, has been adapted to simulate the combustion
of double-base propellants. The kinetics of the model is described by considering firstly, a change of
phase of the solid propellant from condensed to gas and secondly, a reduced chemistry scheme which
defines simplified chemical reactions to represent the combustion itself. To couple both processes,
mass and species conservation, as well as temperature continuity are imposed in the burning surface
in which the burning rate will represent a key factor. Moreover, an energy balance is also applied at
the burning surface which represents that heat that the gas transfers to the burning surface is invested
firstly, in raising the surface temperature to produce the phase change and secondly, in warming the
condensed phase by conduction.
Firstly, and in order to validate the code, the problem of a homogeneous double-base strand
combustion when pressure at the right boundary is fixed is performed obtaining positive results. The
same geometry as that used to validate the code for AP/HTPB composite propellants in sandwich
configuration is taken. In this case, the fuel and oxidiser will not be placed in layers, but they will
distribute homogeneously along the condensed phase. Once this analysis is performed, the same
geometry have been used to check how a depressurisation ramp affect to the combustion of double-
base propellants.
For the first part of the analysis, the test has been run for a range of pressures from 1.1 to 9.1
MPa and the results of burning rate has been compared with those found in literature. Afterwards,
the main variables o interest have been analysed observing the following conclusions:
1. The burning rate increases with the pressure. In addition, the adjustment between the results
obtained with the model and the ones from literature depends on the initial distribution assumed
for the temperature. The higher the temperature in the domain, the better adjustment between
the calculated and the literature burning rate will be obtained at low and high values of fixed
pressure but the worse for intermediate values of pressure range. Finally, temperature of 1100
K is the one chosen to present the results due to its better adjustment with bibliography results.
2. The higher the pressure, the quicker the condensed phase changes sublimates. However, the gas
velocity presents its highest value for the lowest fixed pressure boundary condition. This could
be due to he increment of the gas density with the pressure which leads to a reduction of the
magnitude of the gas velocity field.
3. The evolution of the velocity with the time behaves different depending on pressure considered.
For pressures between 3.4 MPa and 9.1 MPa, the maximum value of the velocity field remains
almost constant for the time steps considered. However, in case of considering pressures between
3.4 and 1.1 MPa the maximum value of the velocity experiments an increment for the considered
time steps. This phenomenon could happen because chemical reactions at low pressures need
some time to be able to start and therefore, for the gas velocity to increase its value while at
high pressures the reaction mechanism is activated from the very beginning of the combustion
process.
4. Since the propellant is homogeneously distributed in the condensed phase, not only fuel but also
oxidiser in gas phase are uniformly distributed too, not presenting local concentration of any of
them, at happened when analysing AP/HTPB composite propellant combustion. Independently
of the considered pressure, the highest amount of fuel and oxidiser is located close to the burning
surface decreasing the mass fraction smoothly with x-coordinate. From this behaviour can be
concluded that, when the propellant in condensed phase sublimates at the burning surface, it
progresses along the control volume until reaching the end of the tube.
5. The maximum value of fuel mass fraction decreases slightly with the time. The increase of the
time leads to a fall of the fuel mass fraction at the end of the tube. The same behaviour is
206
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
observed for the oxidiser mass fraction. Both, fuel and oxidiser, present an increment of their
mass fraction maximum values with the fall of pressure, effect also observed when analysing
the combustion of AP/HTPB composite propellant. Reaction rates of chemical kinetics depend
on the density of the gas therefore, the higher the pressure, the higher the density and the
higher reaction rate will be obtained. A raise of the reaction rate means an increment in the
consumption of fuel, oxidiser and intermediate products. Therefore, since at low pressures the
combustion of fuel and oxidiser is lower, the mass fraction of both will increase.
6. Final products mass fraction increases with the pressure, which means that, the reaction rates
of reaction 2 (R2), reaction 3 (R3) or both increase at high pressures. On the one hand, as
expected and explained in the previous conclusion, the mass fraction of DR1 increases with the
fall of pressure. That means that the higher the pressure, the higher reaction rate of reaction
2 and therefore, the higher amount of DR1 will react to become final products. On the other
hand, DR2 is produced in the reaction of fuel and oxidiser (R1) and consumed when it changes
to final products (R3). The overall result is the increase of DR2 mass fraction with the pressure.
However, the maximum value of mass fraction at 9.1 MPa experiments a slight decrement. This
could be due to the fact that, since reaction rates get stronger with the increase of pressure,
at 9.1 MPa, the reaction rate of R3 provokes a consumption of DR2 enough to decrease the
content produced in R1.
7. Final products mass fraction increase with the time. However, when pressure is fixed between
9.1 MPa to 5.99 MPa, the mass fraction of final products decreases slightly from 0.01 s to 0.02
s (first slot) and raises between 0.02 s and 0.05 s (second slot). This could be explained by
taking into account that final products are the result of R2 and R3 and R3 depends strongly
on R1. DR1 mass fraction value is constant in the first slot and decreases in the second one.
This means that reaction 2 needs some time to increase their effect. Moreover, two phenomena
are observed. On the one hand, at high pressures, the mass fraction of DR2 increases slightly
in the first time slot while decreasing in the second one. On the other hand, at low pressures,
DR2 mass fraction increases with the time which means that, during the first slot, reaction
1 produces more intermediate species than reaction 3 can consume and therefore, there is an
increment of DR2 mass fraction. The same is happening at high pressures for the first time
slot. However, between 0.02 s and 0.05 s, the amount of DR2 mass fraction when pressure is
high suffers a decrement. Since during the second time slot and the pressure fixed is high there
is a strong increment of the final products mass, it can be concluded that, reaction 3 increase
its strength with pressure and time.
8. The temperature increases with pressure. Moreover, the lower the pressure, the softer is the
transition between minimum and maximum temperature. However, as the time passes, the
jump between hot temperature in the gas phase and cold temperature in the condensed phase
becomes more abrupt which means that, the transference of heat in the solid increases with the
fall of pressure and decreases when the time passes.
9. Temperature increases with time. This effect seems to be more visible at low than at high
pressures.
10. The heat flux increases with the pressure and with the time. The greatest amount of heat is
transferred in the burning surface between the solid and the gas. Being the expected result to
observe since the variable is defined as the thermal conductivity multiplied by the temperature
gradient. The value along the rest of the control volume remains very low or even zero.
To conclude, as it was done for AP/HTPB composite propellant, a depressurisation ramp is applied
to the right side of the control volume using the same geometry than the one chosen for the previous
test. The ramp is designed to decrease the pressure from 9.1 MPa to 1.1 MPa in 1 · 10 8 s starting
at 0.01 s. Besides the render views provided by ParaView software, in order to observe properly how
the pressure ramp affect to the behaviour of the variables of interest, three points in the geometry
207
CHAPTER 4. COMBUSTION MODELLING AND RESULTS
have been chosen to plot the variables. The first point chosen is located close to the burning surface
of the 1.1 MPa fixed boundary case when time is equal to 0.01 seconds, the second one is located at
the middle of the initial gas domain and the third one almost in the end of the gas zone. The main
conclusions extracted from the simulation are the following:
1. The pressure wave due to depressurisation evolves from the right side of the tube until the
burning surface and afterwards, it comes back to the right side. This movement is repeated
until reaching the equilibrium pressure of 1.1 MPa.
2. The pressure at P1 (x = 510·10 6 m) when pressure is fixed to 1.1 MPa at the right boundary
fluctuates around this value. However, the pressure at the same point for the 9.1 MPa fixed
boundary condition test remains stable. This is be explained due to the distance between the
first point and the burning surface. When the pressure in the boundary is fixed to 1.1 MPa,
the boundary is located close to P1 while when the pressure fixed at the boundary is 9.1 MPa,
the burning surface is located further, closer to the left boundary. Therefore, the further from
the burning surface the pressure is measured, the more stable will be the pressure. The cells
close to the burning surface are the first to receive the sublimated gas species from the burning
surface and, as a result, the fuel, oxidiser and intermediate species combust firstly in these cells.
Therefore, the density of cells located close to the burning surface will fluctuate and that will
influence to the temperature and pressure distributions.
3. When a pressure ramp is applied to the right boundary the domains suffers a sharp decrease
of pressure. Afterwards, the pressure value at P1 (x = 510·10 6 m) oscillates around 1.1 MPa
reaching even lower pressure values. Therefore, a strong pressure jump applied at the exit of
the chamber leads to an immediate depressurisation in the cells close to the burning surface
even lower than the minimum value of the pressure step. Once the pressure reach 1.1 MPa of
stable pressure, the fluctuations observed when applying 1.1 MPa fixed boundary condition are
less visible. The reason of this phenomenon is attributed again to the position of the burning
surfaces for the different tests. Before suffering depressurisation, the pressure of the ramp test
was fixed to 9.1 MPa, as a consequence, when time is equal to 0.01 seconds its burning surface
is located in the same position far from point one.
4. Depressurisation provokes immediately a diffusion of the temperature along the condensed
phase. The higher the pressure, the sharper will be the transition between cold condensed
and hot gas temperatures.
5. Just after depressurisation, fuel and oxidiser mass fraction increase their values. DR1, remains
almost constant but with a slight increment of concentration while DR2 and final products
mass fractions decrease. As has been already explained, the magnitude of the mixture density
is smaller for low than for high pressures. Since the reaction rates which represent the chemical
kinetics of the combustion process depend on the density, with the fall of pressure, the reaction
rates decrease too. As a result, when depressurisation is applied, not only less fuel and oxidiser
react to form intermediate product DR2, but also less DR1 and DR2 change to become final
products. The sublimation of condensed phase and the decrease of reaction rates are the reason
of the increment of fuel, oxidiser and intermediate products DR1. As a result, there creation of
new final products will fall down and part of them will be lost by the exit of the tube concluding
in a reduction of their mass fraction. The same phenomenon is suffered by intermediate specie
DR2 just after depressurisation until time is equal to 0.001004 s. At this time pressure stabilised
to 1.1 MPa and therefore, the mass fraction of DR2 remains almost constant too.
6. Fuel and oxidiser mass fraction distributions behave exactly the same. At the beginning the val-
ues obtained are the same than those obtained for 9.1 MPa fixed boundary condition. However,
when depressurisation is applied, the mass fractions of fuel and oxidiser start to increase until
reaching the values obtained when 1.1 MPa fixed boundary condition is applied. This could
happen because an increase of pressure leads to a gain in strength of reaction 1 and therefore,
208
4.3. COMBUSTION OF DOUBLE-BASE HOMOGENEOUS PROPELLANT
to the fall of fuel and oxidiser mass fraction which react to provide intermediate specie DR2.
When depressurisation is applied, reaction 1 loses strength and as a result, the mass fraction of
fuel and oxidiser increase.
7. Before depressurisation is applied, the values of fuel and oxidiser mass fractions obtained for
the ramp test are the same than those obtained for the 9.1 MPa fixed boundary test. However,
when depressurisation is applied, fuel and oxidiser mass fractions start to increase until reaching
the values obtained for the 1.1 MPa fixed boundary test. In addition, DR2 and final products
increment their concentration while DR1 falls down in order to match the values of the 1.1 MPa
fixed boundary test.
8. At the beginning, the temperature at the surface has the same value than that obtained at 9.1
MPa fixed boundary condition. When the depressurisation is applied, the temperature at the
surface suffers a sharp decrement fluctuating slightly and afterwards, it remains around 650 K.
9. The maximum temperature of the control volume falls sharply immediately after the pressure
ramp at the boundary is applied. Afterwards, it fluctuates until t = 0.01001 s, time from which
it will experiment a soft increase until reaching the same value as that obtained for the 1.1 MPa
fixed boundary condition test.
10. The temperature for 1.1 MPa fixed boundary condition test measured at the point close to the
burning surface presents some fluctuations, phenomenon which does not take place when ramp
boundary condition is applied since its boundary surface when time is equal to 0.01 s is located
further from the one of 1.1 MPa fixed boundary condition test.
11. The maximum value obtained for heat flux when applying the pressure ramp test is in between
those obtained for the 1.1 MPa and the 9.1 MPa fixed boundary tests. The heat flux, as
the other variables, presents fluctuations immediately after the applying depressurisation. The
strong fall of pressure at the exit provokes a decrease of the heat flux between 0.01 s and 0.0101
s. After that interval of time, its maximum fluctuates until reaching a magnitude closer to the
one obtained for 1.1 MPa fixed boundary condition test.
209
210
Chapter 5
Conclusions
The improvement on the design and safety in industrial processes which use solid propellants are
of major importance in industry of energetic materials. That is the reason why, modelling the
combustion of solid propellants is a key problem that has focused the interest of several researchers
over the years.
The main objective of this thesis have been the study of the combustion of granulated, composite
and double-base solid propellant. In order to reach this target a literature review have been performed
for each propellant in the first place and in the second place, numerical codes have been developed
to characterise their combustion. Finally, the models developed have been validated, either with the
results obtained by other researchers of the existing literature or with experimental data.
Specifically, the objectives of the thesis could be classified in two, depending on the solid propel-
lant. The first one is to improve the understanding and modelling of the deflagration-to-denotation
(DDT) phenomenon in granular beds of solid propellants. In order to reach this, a tool which provides
the characterisation of the combustion process has been developed. In addition, the results obtained
have been compared with those available in the existing literature to assess the effectivity of the model
to predict the early stages of transient combustion processes. The second objective is the study of
composite and double-base propellants and the characterisation of their unsteady combustion. In this
case, a multi-dimensional code which solves the set of differential equations for condensed and gas
phases have been programmed in OpenFoam. Finally, the results obtained for the combustion simu-
lation of both propellants have been compared against experimental tests and results from existing
literature. As a result, a tool which predicts composite and double propellant combustion behaviour
in multidimensional scenarios with transient environmental conditions has been developed.
The main conclusions obtained after the modelling of each propellant are summarised throughout
this chapter.
211
CHAPTER 5. CONCLUSIONS
Firstly, a sensitivity analysis is performed to determine the effect that Nusselt number, friction
coefficient and initiation conditions have in the results. The analysis have been performed using both,
Rusanov and MacCormack numerical schemes. Due to the capacity of reproducing very similar results
between both numerical schemes and the good behaviour at a wide range of Reynolds number, the
Nusselt number and friction coefficient presented by Butler et al. [39] in their work will be adopted
afterwards for calculation.
Although the authors consider a small zone of the control volume already initiated, each one
applies different values of pressure and temperature as initial conditions of the test. As a result, high
differences in the results can be obtained depending on the initial values assumed for the variables.
Therefore, and in order to find to validate the model with literature values, the distributions of
pressure, porosity and gas and particle temperatures provided by Hoffman and Krier [28] at 20 µs,
have been used as initial conditions for the calculations performed afterwards.
Once initial conditions and source terms are defined, two different models have been used to
calculate pressure, temperature and porosity distributions. The first model which have been studied
considers the modification of the particle momentum governing equation done by Hoffman and Krier
[28]. According to the author, this modification prevents the porosity from reaching values below
the minimum value of compaction that packed beds of spherical particles can reach. A non-realistic
value of particle temperature has been obtained as a result of applying this model. The second
model which has been studied does not consider the porosity limiter from Hoffman and Krier [28]
performing the limitation of the porosity directly in the code. The magnitude of the values which have
been obtained for the main variables of interest are similar to those found in literature. Moreover,
the results obtained using Rusanov scheme agree well with those resulting of applying MacCormack-
TVD numerical scheme. However, the distribution of the variables are displaced in x-direction respect
those from literature. These differences could be due on the one hand, to the initial values of the
parameters chosen as initial conditions which highly determine the x-location of the peak values for
all variables and on the other hand, due to the lack of all necessary input data in a single work from
bibliography making necessary the values from different works increasing the difficulty of reproducing
the tests available in the bibliography.
Finally, it can be concluded that, despite the difference in the x-location peak values, the last
model considered represents accurately the physical behaviour of the propellant combustion for all
variables of interest becoming a predictive tool for the characterisation of the early stages of the
detonation process of granular solid propellants.
212
5.2. COMPOSITE AP/HTPB PROPELLANT
different regions, one containing AP and other of HTPB, which have different thermophysical prop-
erties. As in the Crawford Strand Burner test, the simulation assumes that the tube is insulated and
there is no transference of mass or energy in any boundary despite the right one which is open to
atmosphere. The pressure at the right side of the tube has been fixed for this test and a range of
pressures between 0.5 MPa to 7 MPa have been tested. The main conclusions extracted from the
analysis can be summarised in the following lines:
1. The burning rate results obtained with the developed two-dimensional model of this work
applied to the combustion of AP/HTPB composite propellant in sandwich configuration agree
well with the experimental ones from Kohga [51] for the considered pressure range.
2. The higher the pressure, the lower the gas velocity is obtained. This behaviour may be due
to the increase of the gas density with the pressure. However, it has also been observed that
an increase of AP mass fraction close to the burning surface leads to a local raise of the gas
velocity. This behaviour is attributed to the relationship between burning rate and gas velocity
at the burning surface.
3. AP mass fraction presents a local concentration in the cells close to burning surface which
has not an uniform distribution along the tube like it happens for HTPB. Since the activation
energy that the oxidiser needs to be decomposed is not high, compared with the one from
HTPB, it could be deduced that, as soon as condensed AP sublimates to AP in gas phase, the
first reaction takes place and most of the AP is decomposed in reaction products. However,
binder combustion reaction not only presents a higher activation energy, but also needs the
products from the first reaction to take place. This could lead to the conclusion that when
HTPB sublimates to gas, it does not burn immediately and therefore, it has enough time to
spread along the gas domain making possible for the cells far from the burning surface to
accumulate higher concentrations of binder than oxidiser.
4. The content of AP and HTPB in gas phase multiplies with the pressure decrement. This could
be assigned to the dependency of reaction rates expressions with the pressure. A large value of
reaction rates means that the consumption of both, binder and oxidiser, will be high. Therefore,
it has been concluded that high pressures increase the reaction rate leading to a consumption
of propellant mass fraction.
5. The diffusivity of AP increases at low pressures. However, it has been observed that at 0.5
MPa, both binder and oxidiser mass fractions present a nonuniform distribution with higher
concentrations in the cells of gas adjacent to their condensed phases respectively. This could
happen because at low pressures, the degradation of the propellant is less intense and, in
addition, R1 and R2 reaction rates of decrease. Therefore, the combustion of AP and HTPB
in the cells close to the burning surface is reduced increasing as a result the mass fractions of
binder and oxidiser in these cells.
6. The value of constant temperature in gas phase increases slightly with the fall of pressure
fixed at the right boundary. The fall of pressure provokes a decrease of reaction rates and
therefore, there is a decrement in the combustion of AP and HTPB. As consequence of the loss
of strength of AP and HTPB combustion reactions, the concentration of the reaction products
(HCl, CO2 , H2 O N2 ) decreases mildly too. The slight decrement of the reaction products
mass fractions is enough to make the density decrease its value provoking an increment the
temperature.
Finally, the combustion simulation of AP/HTPB composite propellant in sandwich configuration has
been performed applying a depressurisation ramp from 7 MPa to 0.5 MPa during a short interval of
time at the open side of the tube in order to see how it influences in the main variables of interest.
The conclusions obtained are the following:
213
CHAPTER 5. CONCLUSIONS
1. The pressure ramp induced at the right boundary leads to a decrease of the pressure value below
0.5 MPa at the cells located in the middle of the tube. Therefore, it has been concluded that
when the end of the chamber is put down to a rapid depressurisation the pressure wave evolves
from the end to the beginning of the gas control volume and comes back causing a depression
under 0.5 MPa.
2. The fall of pressure provoked by the depressurisation at the boundary leads to an increase of
both AP mass fraction and its diffusivity.
3. Together with the sharp fall of pressure, comes a raise of HTPB mass fraction. Since the
depressurisation at the boundary provokes a fall of pressure close to the burning surface even
below 0.5 MPa, the value of HTPB mass fraction increases over the one obtained for the 0.5
MPa fixed boundary test. The mass fraction values obtained for both tests come closer as the
pressure for the pressure ramp boundary condition test becomes closer to 0.5 MPa.
4. The immediate effect of the sharp fall of pressure at the end of the tube is the decrease of
the temperature at this location as could be expected by taking into account ideal gas law.
Afterwards, temperature raises again remaining mostly stable. The highest value of temperature
reached at the end of the tube is that obtained when fixing the pressure at the right boundary
to 0.5 MPa followed by the depressurisation test and finally, the value obtained for the 7 MPa
fixed boundary test. Since the reaction rates depend on the pressure, the lower the pressure, the
lower reaction rate. This leads to a decrease of the reaction products mass fractions with the
consequent decrease of density and therefore, an increment of temperature. The temperature
obtained when pressure is fixed to 0.5 MPa at the right boundary is slightly higher than that
of the depressurisation test. This happens because before applying the pressure ramp at the
right boundary, the pressure is fixed to 7 MPa and therefore, its density is higher than in 0.5
MPa test.
5. When depressurisation is applied, the burning rate decreases reaching a value even lower than
that obtained for the 0.5 MPa fixed boundary test. Sharp depressurisation at the boundary
leads to obtaining lower values of pressure closer to the burning surface than the ones obtained
by fixing the pressure to 0.5 MPa. Since burning rate is directly related with the pressure,
together with a decrease of pressure below 0.5 MPa, a fall of the burning rate below the value
obtained for 0.5 MPa fixed pressure test will come.
214
5.3. DOUBLE-BASE HOMOGENEOUS PROPELLANT
3. The evolution of the velocity with the time behaves different depending on pressure considered.
For pressures between 3.4 MPa and 9.1 MPa, the maximum value of the velocity field remains
almost constant for the time steps considered. However, in case of considering pressures between
3.4 and 1.1 MPa the maximum value of the velocity experiments an increment for the considered
time steps. This phenomenon could happen because chemical reactions at low pressures need
some time to be able to start and therefore, for the gas velocity to increase its value while at
high pressures the reaction mechanism is activated from the very beginning of the combustion
process.
4. The maximum value of fuel mass fraction decreases slightly with the time. The increase of the
time leads to a fall of the fuel mass fraction at the end of the tube. The same behaviour is
observed for the oxidiser mass fraction. Both, fuel and oxidiser, present an increment of their
mass fraction maximum values with the fall of pressure, effect also observed when analysing
the combustion of AP/HTPB composite propellant. Reaction rates of chemical kinetics depend
on the density of the gas therefore, the higher the pressure, the higher the density and the
higher reaction rate will be obtained. A raise of the reaction rate means an increment in the
consumption of fuel, oxidiser and intermediate products. Therefore, since at low pressures the
combustion of fuel and oxidiser is lower, the mass fraction of both will increase.
5. Final products mass fraction increases with the pressure, which means that, the reaction rates of
reaction 2 (R2), reaction 3 (R3) or both increase at high pressures. On the one hand, the mass
fraction of DR1 increases with the fall of pressure. That means that the higher the pressure,
the higher reaction rate of reaction 2 and therefore, the higher amount of DR1 will react to
become final products. On the other hand, DR2 is produced in the reaction of fuel and oxidiser
(R1) and consumed when it changes to final products (R3). The overall result is the increase
of DR2 mass fraction with the pressure. However, the maximum value of mass fraction at at
9.1 MPa experiments a slight decrement. This could be due to the fact that, since reaction
rates get stronger with the increase of pressure, at 9.1 MPa, the reaction rate of R3 provokes a
consumption of DR2 enough to decrease the content produced in R1.
6. Final products mass fraction increases with the time. However, when pressure is fixed between
9.1 MPa to 5.99 MPa, the mass fraction of final products decreases slightly from 0.01 s to 0.02
s ( first slot) and raises between 0.02 s and 0.05 s ( second slot). This could be explained by
taking into account that final products are the result of R2 and R3 and R3 depends strongly on
R1. DR1 mass fraction value is constant in the first time slot and decreases in the second one.
This means that reaction 2 needs some time to increase their effect. Moreover, two phenomena
are observed. On the one hand, at high pressures, the mass fraction of DR2 increases slightly
in the first time slot while decreasing in the second one. On the other hand, at low pressures,
DR2 mass fraction increases with the time which means that, during the first slot, reaction
1 produces more intermediate species than reaction 3 can consume and therefore, there is an
increment of DR2 mass fraction. The same is happening at high pressures for the first time
slot. However, between 0.02 s and 0.05 s, the amount of DR2 mass fraction when pressure is
high suffers a decrement. Since during the second time slot and the pressure fixed is high there
is a strong increment of the final products mass, it can be concluded that, reaction 3 increase
its strength with pressure and time.
7. The temperature increases with pressure. Moreover, the lower the pressure, the softer is the
transition between minimum and maximum temperature. However, as the time passes, the
jump between hot temperature in the gas phase and cold temperature in the condensed phase
becomes more abrupt which means that, the transference of heat in the solid increases with the
fall of pressure and decreases when the time passes.
8. The heat flux increases with the pressure and with the time. The greatest amount of heat is
transferred in the burning surface between the solid and the gas. Being the expected result to
215
CHAPTER 5. CONCLUSIONS
observe since the variable is defined as the thermal conductivity multiplied by the temperature
gradient. The value along the rest of the control volume remains very low or even zero.
To conclude, and as it was done for AP/HTPB composite propellant, a depressurisation ramp is
applied to the right side of the control volume using the same geometry than the one chosen for the
previous test. The ramp is designed to decrease the pressure from 9.1 MPa to 1.1 MPa in 1 · 10 8
s starting at 0.01 s. Three points in the geometry have been chosen to plot the variables. The first
point chosen is located close to the burning surface of the 1.1 MPa fixed boundary case when time is
equal to 0.01 seconds, the second one is located at the middle of the initial gas domain and the third
one almost in the end. of the gas zone. The main conclusions extracted from the simulation are the
following:
1. The pressure wave due to depressurisation evolves from the right side of the tube until the
burning surface and afterwards, it comes back to the right side. This movement is repeated
until reaching the equilibrium pressure of 1.1 MPa.
2. The pressure for the 1.1 MPa fixed boundary test measured at P1 (x = 510·10 6 m) fluctuates
around 1.1 MPa. The cells close to the burning surface are the first to receive the sublimated
gas species from the burning surface and, as a result, the fuel, oxidiser and intermediate species
combust firstly in these cells. the density of cells located close to the burning surface will
fluctuate and that will influence to the temperature and pressure distributions.
3. Depressurisation provokes immediately a diffusion of the temperature along the condensed
phase. The higher the pressure, the sharper will be the transition between cold condensed
and hot gas temperatures.
4. The reaction rates which represent the chemical kinetics of the combustion process depend
on the density, with the fall of pressure, the reaction rates decrease too. As a result, when
depressurisation is applied, not only less fuel and oxidiser react to form intermediate product
DR2, but also less DR1 and DR2 change to become final products. The sublimation of condensed
phase and the decrease of reaction rates are the reason of the increment of fuel, oxidiser and
intermediate products DR1. As a result, there creation of new final products will fall down and
part of them will be lost by the exit of the tube concluding in a reduction of their mass fraction.
The same phenomenon is suffered by intermediate specie DR2 just after depressurisation until
time is equal to 0.001004 s. At this time pressure stabilised to 1.1 MPa and therefore, the mass
fraction of DR2 remains almost constant too.
5. The maximum temperature of the control volume falls sharply immediately after the pressure
ramp at the boundary is applied. Afterwards, it fluctuates until t = 0.01001 s, time from which
it will experiment a soft increase until reaching the same value as that obtained for the 1.1 MPa
fixed boundary condition test.
6. The heat flux present fluctuations immediately after the applying depressurisation. The strong
fall of pressure at the exit provokes a decrease of the heat flux between 0.01 s and 0.0101 s.
After that interval of time, its maximum fluctuates until reaching a magnitude closer to the one
obtained for 1.1 MPa fixed boundary condition test.
216
Future work
Several aspects could be considered as future work of this thesis, not only to complete the understand-
ing of combustion of solid propellants but also to improve the results and decrease the differences
obtained between experimental and numerical ones.
The first idea to be considered is the influence of AP size in the burning rate. This effect was taken
into account by Kohga [51] who studied the influence of AP size in the burning rate by comparing
the results obtained when using particles diameters of 4 mm and 110 mm. In addition, Jeppson and
Beckstead [42], studied this effect by choosing AP grains with diameters from 20 mm to 200 mm and
plotting the burning rate results for each one of them. Kubota [3] also mentioned the influence that
the AP particle size has on the burning rate. The three authors remarked that, considering the
same propellant composition, the smaller the particle size, the higher burning rate. Therefore, it
is considered of high importance to study the effect of the particle size together with the pressure,
not only in the burning rate, but also in the temperature and AP mass fraction distributions. This
analysis will be performed by varying the parameters of equations (4.83) and (4.84). By considering
a fixed AP composition (a), for each chosen AP particle diameter (dAP ), equation (4.84) will provide
the height of HTPB (dHT P B ). By introducing both values in equation (4.83), the total height of the
control volume, (L), can be obtained. Moreover, the burning rate results could be compared with the
ones available in the literature, such as the ones from Jeppson and Beckstead [42] and from Kohga [51].
Another aspect considered as future work is the study of the influence of the propellant composition in
the burning rate, this could be done, not only for the AP/HTPB in sandwich configuration, but also
for the homogeneous strand and the double-base solid propellant. The results could be compared
with the ones presented by Kohga [51] who also studied the influence of AP mass fraction in the
burning rate.
Finally, and in order to improve the results obtained, it could be interesting to include the release
of heat due to radiation and the emissivity of the surfaces which have been assumed negligible when
developing the mathematical model. In the present work, only the heat due to conduction has been
considered and since the experiments are done in a closed combustion chamber, reflecting walls could
be approximated as black body or grey body of high emissivity. This consideration could reduce the
error between experimental and numerical results.
217
218
Nomenclature
c Specific heat
cv Specific heat at constant volume
cp Specific heat at constant pressure
cg Speed of sound
dAP Diameter of AP particles
dHT P B Diameter of HTPB particles
d Normal distance between cells centre and considered surface
dp Granulated particles diameter
Dj Reaction rate constants of jth-reaction
Da Damköhler number
Ec Activation energy of condensed phase
Ej Activation energy of jth-reaction
E Energy
FD Drag force
k Thermal conductivity
Mg Molecular mass
nj Exponential pressure factor of jth-reaction
Nu Nusselt number
p Pressure
Pe Peclet number
Pr Prandtl number
Q̇ Heat transfer
Qg Heat of reaction
Qc Heat of decomposition of condensed phase
Rj Reaction rates of jth-reaction
Re Reynolds number
Rg Specific constant of gases
Ru Universal constant of gases
ṙ Burning rate
Sp Granulated particles surface
T Temperature
t Time
u Velocity
Vp Granulated particles volume
Y Mass fraction
219
Greek Symbols
↵ Porosity
˙ Mass transfer
Heat capacity ratio
⌘g Co-volume
µg Gas viscosity
Subscripts
c Condensed
g Gas
melt Melting
p Propellant
s Surface
220
Appendix I Heat Flux detailed figures
Cell values of the heat flux can bee seen in Figures A I.1and A I.2.
Figure A I.1: Heat flux distribution in t = 0.02 s at 7 MPa, 5 MPa, 3 MPa (above-below).
221
Figure A I.2: Heat flux distribution in t = 0.02 s at 2 MPa, 1 MPa, 0.5 MPa (above-below).
222
Figure A I.3: Heat flux distribution in t = 0.04 s at 7 MPa, 5 MPa, 3 MPa, 2 MPa, 1 MPa, 0.5 MPa
(above-below).
223
224
Appendix II Dimensionless numbers
In order to understand and characterise the combustion problem, a nondimensionalisation has been
introduced. According to Margolis and Williams [87] the first step should be define a characteristic
velocity defined as the speed of propagation Uc = dxs /dt , where xs is the coordinate of the burning
surface. Therefore, the characteristic velocity chosen in this work will be the burning speed which
depends on the pressure of the selected test. The remaining expressions [87] for dimensionless variables
are,
where the subindex a stands for dimensionless. T0 , pg,0 are the initial unburned temperature and
pressure respectively and therefore, ⇢g,0 is the density at this pressure and temperature conditions.
Lc is the characteristic length of the problem which in this case is 170 µm. The characteristic time
of the problem (tf ) has been calculated with the following expression,
g
tf = (5.2)
⇢g cvg Uc2
Besides the nondimensionalisation of the problem variables, it is also of high interest to calculate
the typical non dimensional numbers in order to characterise the properties of the fluid. Characteristic
diffusion and combustion time, as well as Reynolds, Peclet and Damköhler number has been calculated
and their results are presented afterwards. The expressions of combustion and diffusion times have
been taken from Cai et al. [45] as,
˙ td = L2 /D
tc = ⇢/!, (5.3)
Peclet, Reynolds and Damköhler numbers have the following expressions, .
u g L c ⇢g ug Lc
Re = , Pe = , Da = td /tc (5.4)
µ D
For AP/HTPB sandwich propellant, the results obtained have been summarised in Table A II.1
and Table A II.2. Dimensionless values of gas and condensed phase temperatures are almost the same
independently of the pressure set therefore,
Tc,g 2700 Tc,g 370
Ta g = = = 8, Tac = = = 1.233 (5.5)
T0 300 T0 300
Since the combustion chamber is open, the pressure of the domain is the same as the initial
pressure. As a result, the dimensionless pressure is pga = 1.
225
OpenFoam
Pressure
Burning rate tf (µs) ta xa
(MPa)
(mm/s)
7 7.52 1.7880 0.0006 7.3800
5 6.61 2.0282 0.0005 6.4870
3 5.53 2.1823 0.0004 5.4271
2 4.82 1.7975 0.0002 4.7303
1 3.89 1.4319 0.0001 3.8176
0.5 3.19 1.0070 0.0001 3.1306
Pressure
ug (m/s) ug,a ⇢g (kg/m3 ) ⇢g,0 (kg/m3 ) ⇢g
(MPa)
7 1.8 239.3617 9.2063 80.7614 0.1140
5 2.0 302.5719 6.5740 57.6867 0.1140
3 2.5 452.0796 3.9103 34.6120 0.1130
2 3.4 705.3942 2.5666 23.0747 0.1112
1 5.5 1413.8817 1.2313 11.5373 0.1067
0.5 9.5 2978.0564 0.5868 5.7687 0.1017
Pressure 8 8 4
tc1 · 10 tc2 · 10 td ·10
(MPa)
7 1.3170 4.5967 70.0160
5 1.6911 5.9025 49.9966
3 2.45168 8.5570 29.7389
2 3.2637 11.3914 19.5198
1 5.2446 18.3056 9.3645
0.5 8.3723 29.2221 4.4629
Pressure
Re Da1 Da2 Pe
(MPa)
7 41.2 531644.5 152318.7 74.1
5 32.7 295647.4 84704.5 58.8
3 24.3 121302.5 34753.7 43.7
2 21.7 59809.0 17135.6 39.0
1 16.8 17855.4 5115.6 30.3
0.5 13.9 5330.6 1527.2 24.9
226
Appendix III Granulated propellant
code structure
227
228
Appendix IV AP/HTPB and
double-base propellants code
structure
229
230
Bibliography
[2] J.P. Agrawal. High Energy Materials. Propellants, Explosives and Pyrotechnics. Wiley-VCH,
2010.
[3] N. Kubota. Propellants and Explosives. Thermochemical Apects of Combustion. Wiley-VCH,
2002.
[4] M.W. Beckstead, N.L. Peterson, D.T. Pilcher, B.D. Hopkins, and H. Krier. Convective combus-
tion modeling applied to deflagration-to-detonation transition of HMX. Combustion and Flame,
30:231–241, 1977.
[5] A. Harrland and I.A. Johnston. Review of Solid Propellant Ignition Models Relative to the
Interior Ballistic Modelling of Gun Systems. Weapons Systems Division Defence Science and
Technology Organisation, 2012.
[6] K.L. Kosanke, B.J. Kosanke, B.T. Sturman, and R.M. Winokur. Encyclopedic dictionary of
pyrotechnics and related subjects. Pyrotechnic reference series, no. 5. Whitewater, Colo. : Journal
of Pyrotechnics, Inc., 2012.
[7] H. Krier and S. Gokhale. Modeling of Convective Mode Combustion through Granulated Pro-
pellant to Predict Detonation Transition. AIAA Journal, 16(2):177–183, 1978.
[8] M.W. Beckstead, K. Puduppakkam, P. Thakre, and V. Yang. Modeling of combustion and
ignition of solid-propellant ingredients. Progress in Energy and Combustion Science, 33(6):497–
551, 2007.
[9] K.K. Kuo and M. Summerfield. Fundamentals of Solid-Propellant Combustion, volume 90 Pro-
gress in Astronautics and Aeronautics. American Institute of Aeronautics and Astronautics, Inc.,
1984.
[10] H. Krier. Solid Propellant Burning Rate During a Pressure Transient. Combustion Science and
Technology, 5(1):69–73, 1972.
[11] H. Krier, W.A. Sirignano, M. Summerfield, and J.S. Tien. Nonsteady burning phenomena of
solid propellants - Theory and experiments. AIAA Journal, 6(2):278–285, 1968.
[12] K.K. Kuo, R. Vichnevetsky, and M. Summerfield. Theory of Flame Front Propagation in Porous
Propellant Charges under Confinement. AIAA Journal, 11(4):444–451, 1973.
[13] J.H. Koo and K.K. Kuo. Transient Combustion in Granular Propellant Beds. Part I. Theoret-
ical Modeling and Numerical Solution of Transient Combustion Processes in Mobile Granular
Propellant Beds. Technical report, Pennsylvania State University Park Dept. of Mechanical
Engineering, 1977.
231
BIBLIOGRAPHY
[14] N.C. Markatos. Modelling of two-phase transient flow and combustion of granular propellants.
International Journal of Multiphase Flow, 12(6):913–933, 1986.
[15] D.B. Spalding. Numerical computation of multi-phase fluid flow and heat transfer. In Recent
Advances in Numerical Methods in Fluids, pages 139–167. 1980.
[16] P.S. Gough and F.J. Zwarts. Modeling Heterogeneous Two-Phase Reacting Flow. AIAA Journal,
17(1):17–25, 1979.
[17] R.A. Oton-Martinez, G. Monreal-Gonzalez, J.R. Garcia-Cascales, F. Vera-Garcia, F.J.S. Velasco,
and F.J. Ramirez-Fernandez. An approach formulated in terms of conserved variables for the
characterisation of propellant combustion in internal ballistics. International Journal for Nu-
merical Methods in Fluids, 79(8):394–415, 2015.
[18] H.W. Sandusky and W.L. Elban. Quasi-static compaction of porous propellant beds. II. Exper-
iments and application of lattice compaction model to cannon propellants. Powder Technology,
89(3):219–229, 1996.
[19] D.E. Kooker, S.L. Howard, and L.M. Chang. Convective ignition of a granular solid propel-
lant bed: Influence of propellant composition. Symposium (International) on Combustion,
26(2):2033–2040, January 1996.
[20] K.K. Kuo and M. Summerfield. Theory of Steady-State Burning of Gas-Permeable Propellants.
AIAA Journal, 12(1):49–56, 1974.
[21] W.F. van Tassell and H. Krier. Combustion and flame spreading phenomena in gas-permeable
explosive materials. International Journal of Heat and Mass Transfer, 18(12):1377–1386, 1975.
[22] H. Krier, S. Rajan, and W.F. Van Tassell. Flame-Spreading and Combustion in Packed Beds of
Propellant Grains. AIAA Journal, 14(3):301–309, 1976.
[23] K.K. Kuo, J.H. Koo, T.R. Davis, and G.R. Coates. Transient combustion in mobile gas-permeable
propellants. Acta Astronautica, 3(7):573–591, 1976.
[24] J.M. Lenoir and G. Robillard. A mathematical method to predict the effects of erosive burning
in solid-propellant rockets. Symposium (International) on Combustion, 6(1):663–667, 1957.
[25] W. Denton. The Heat Transer and Flow Resistance for Fluid Flow Through Randomly Packed
Spheres. American Society of Mechanical Engineers, 1951.
[26] S. Ergun. Fluid flow through packed columns. Journal of Chemical Engineering Progress, 48,
No. 2:89–94, 1952.
[27] H. Krier and J.A. Kezerle. A separated two-phase flow analysis to study deflagration-to-
detonation transition (DDT) in granulated propellant. Symposium (International) on Com-
bustion, 17(1):23–34, 1979.
[28] S.J. Hoffman and H. Krier. Fluid Mechanical Processes of Deflagration to Detonation Transition
in Beds of Porous Reactive Solids., volume Technical report AAE -80_2. Aeronautical and
Astronautical Engineering Department, University of Illinois at Urbana-Champaign, 1980.
[29] S.J. Hoffman and H. Krier. Fluid Mechanics of Deflagration-to-Detonation Transition in Porous
Explosives and Propellants. AIAA Journal, 19(12):1571–1579, 1981.
[30] K.K. Kuo and C.C. Nydegger. Cold flow resistance measurement and correlation in a packed
bed of WC-870 spherical propellants. Journal of Ballistics, 2:1–26, 1978.
[31] S.S. Gokhale and H. Krier. Modeling of unsteady two-phase reactive flow in porous beds of
propellant. Progress in Energy and Combustion Science, 8(1):1–39, 1982.
232
BIBLIOGRAPHY
[32] P.B. Butler and H. Krier. Analysis of deflagration to detonation transition in high-energy solid
propellants. Combustion and Flame, 63(1):31–48, 1986.
[33] J.M. Powers. Theory of detonation structure for two-phase materials. Thesis. University of
Illinois at Urbana-Champaign, 1988.
[34] M.R. Baer and J.W. Nunziato. A two-phase mixture theory for the deflagration-to-detonation
transition (ddt) in reactive granular materials. International Journal of Multiphase Flow,
12(6):861–889, 1986.
[35] J.M. Powers, D.S. Stewart, and H. Krier. Theory of two-phase detonation Part I Modeling.
Combustion and Flame, 80(3):264–279, 1990.
[36] J.M. Powers, D.S. Stewart, and H. Krier. Theory of two-phase detonation Part II: Structure.
Combustion and Flame, 80(3):280–303, 1990.
[37] J. Powers, D.S. Stewart, and H. Krier. Two-phase steady detonation analysis. NASA STI/Recon
Technical Report N, 1988.
[38] N.I. Gelperin and V.G. Einstein. Heat transfer in fluidized beds. Fluidization, Academic Press.
New York, 1971.
[39] P.B. Butler, M.F. Lembeck, and H. Krier. Modeling of shock development and transition to
detonation initiated by burning in porous propellant beds. Combustion and Flame, 46:75–93,
1982.
[40] S.F. Wilcox, D.P. Jones, and H. Krier. Technical Report UILU-ENG 80-0501, 1980.
[41] M.W. Beckstead, R.L. Derr, and C.F. Price. A model of composite solid-propellant combustion
based on multiple flames. AIAA Journal, 8(12):2200–2207, 1970.
[42] M. Jeppson, M. Beckstead, and Q. Jing. A kinetic model for the premixed combustion of a
fine AP/HTPB composite propellant. In 36th AIAA Aerospace Sciences Meeting and Exhibit,
Aerospace Sciences Meetings. American Institute of Aeronautics and Astronautics, 1998.
[43] C. Guirao and F.A. Williams. A model of ammonium perchlorate deflagration between 20 and
100 atm. AIAA Journal, 9(7):1345–1356, 1971.
[44] G.M. Knott and M.Q. Brewster. Modeling the combustion of propellant sandwiches. Combustion
Science and Technology, 174(4):61–90, 2002.
[45] W. Cai and V. Yang. A model of AP/HTPB composite propellant combustion. In 38th Aerospace
Sciences Meeting and Exhibit, Aerospace Sciences Meetings. American Institute of Aeronautics
and Astronautics, 2000.
[46] R. Ye, Y.G. Yu, and Y.J. Cao. Analysis of Micro-scale Flame Structure of AP/HTPB Base
Bleed Propellant Combustion. Defence Technology, 9(4):217–223, 2013.
[47] Y.J. Cao, Y.G. Yu, and R. Ye. Numerical analysis of AP/HTPB composite propellant combustion
under rapid depressurization. Applied Thermal Engineering, 75:145–153, 2015.
[48] N.D. Vo, M.Y. Jung, D.H. Oh, J.S. Park, I. Moon, and M. Oh. Moving boundary modeling for
solid propellant combustion. Combustion and Flame, 189:12–23, 2018.
[49] W. Cai, P. Thakre, and V. Yang. A Model of AP/HTPB Composite Propellant Combustion in
Rocket-Motor Environments. Combustion Science and Technology, 180(12):2143–2169, 2008.
[50] J.A. Moríñigo and J. Hermida-Quesada. Evaluation of reduced-order kinetic models for HTPB-
oxygen combustion using LES. Aerospace Science and Technology, 58:358–368, 2016.
233
BIBLIOGRAPHY
[55] R.P. Fitzgerald and M.Q. Brewster. Infrared imaging of AP/HTPB laminate propellant flames.
Combustion and Flame, 154(4):660–670, 2008.
[56] F. Miccio. Numerical modeling of composite propellant combustion. Symposium (International)
on Combustion, 27(2):2387–2395, 1998.
[57] G. Favale and F. Miccio. Modeling unsteady and perturbed combustion of heterogeneous com-
posite propellants. Aerospace Science and Technology, 12(4):285–294, 2008.
[58] G.M. Knott and M.Q. Brewster. Two-dimensional combustion modeling of heterogeneous solid
propellants with finite Peclet number. Combustion and Flame, 121(1):91–106, 2000.
[59] X. Zhou, T.L. Jackson, and J. Buckmaster. A numerical study of periodic sandwich propellants
with oxygenated binders. Combustion Theory and Modelling, 7(2):435–448, 2003.
[60] T.L. Jackson and J. Buckmaster. Heterogeneous Propellant Combustion. AIAA Journal,
40(6):1122–1130, 2002.
[61] W. Cai. Two-phase flow interactions and combustion of AP/HTPB composite propellant in
rocket motors with acoustic oscillations. PhD thesis, 2001. OCLC: 55721452.
[62] G. Lengelle. Thermal degradation kinetics and surface pyrolysis of vinyl polymers. AIAA Journal,
8(11):1989–1996, 1970.
[63] M.M. Ibiricu and F.A. Williams. Influence of externally applied thermal radiation on the burning
rates of homogeneous solid propellants. Combustion and Flame, 24:185–198, 1975.
[64] G. Lengelle, J. Duterque, C. Verdier, A. Bizot, and J.F. Trubert. Combustion mechanisms
of double base solid propellants. Symposium (International) on Combustion, 17(1):1443–1451,
1979.
[65] A. Bizot and M.W. Beckstead. A model for double base propellant combustion. Symposium
(International) on Combustion, 22(1):1827–1834, 1989.
[66] T.S. Roh, S. Apte, and V. Yang. Transient combustion response of homogeneous solid propellant
to acoustic oscillations in a rocket motor. Symposium (International) on Combustion, 27(2):2335–
2341, 1998.
[67] S.Y. Hsieh and V. Yang. A Preconditioned Flux-Differencing Scheme for Chemically Reacting
Flows at all Mach Numbers. International Journal of Computational Fluid Dynamics, 8(1):31–
49, 1997.
[68] M.Q. Brewster, M.J. Ward, and S.F. Son. Simplified Combustion Modeling of Double Base
Propellant: Gas Phase Chain Reaction Vs. Thermal Decomposition. Combustion Science and
Technology, 154(1):2–30, 2000.
234
BIBLIOGRAPHY
[69] S. Apte and V. Yang. Unsteady flow evolution and combustion dynamics of homogeneous solid
propellant in a rocket motor. Combustion and Flame, 131(1):110–131, 2002.
[70] T.S. Roh, I.S. Tseng, and V. Yang. Effects of acoustic oscillations on flame dynamics of homo-
geneous propellants in rocket motors. Journal of Propulsion and Power, 11(4):640–650, 1995.
[71] X. Wu, M. Kumar, and K.K. Kuo. A comprehensive erosive-burning model for double-base
propellants in strong turbulent shear flow. Combustion and Flame, 53(1):49–63, 1983.
[72] I. Shih Tseng and Vigor Yang. Combustion of a double-base homogeneous propellant in a rocket
motor. Combustion and Flame, 96(4):325–342, 1994.
[73] V.A. Arkhipov, V.E. Zarko, I.K. Zharova, A.S. Zhukov, E.A. Kozlov, D.D. Aksenenko, and
A.V. Kurbatov. Solid propellant combustion in a high-velocity cross-flow of gases (review).
Combustion, Explosion, and Shock Waves, 52(5):497–513, 2016.
[74] F.W. Robbins and T. Keys. The Burning Rate Behavior of Pure Nitrocellulose Propellant
Samples:. Technical report, Defense Technical Information Center, Fort Belvoir, VA, 1993.
[75] Q.L. Yan, X.J. Li, Y. Wang, W.H. Zhang, and F.Q. Zhao. Combustion mechanism of double-base
propellant containing nitrogen heterocyclic nitroamines (I): The effect of heat and mass transfer
to the burning characteristics. Combustion and Flame, 156(3):633–641, 2009.
[76] Q.L. Yan, Z.W. Song, X.B. Shi, Z.Y. Yang, and X.H. Zhang. Combustion mechanism of double-
base propellant containing nitrogen heterocyclic nitroamines (II): The temperature distribution
of the flame and its chemical structure. Acta Astronautica, 64(5):602–614, 2009.
[77] I. Aoki and N. Kubota. Combustion Wave Structures of High- and Low-Energy Double-Base
Propellants. AIAA Journal, 20(1):100–105, 1982.
[78] V.N. Marshakov and B.V. Novozhilov. Transient modes of double-base propellant combustion
in a semiclosed volume. Russian Journal of Physical Chemistry B, 5(1):45–56, 2011.
[79] E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics - A Practical Intro-
duction. Springer, 2009.
[80] V.V. Rusanov. Calculation of interaction of non-steady shockwaves with obstacles,. National
Research Council of Canada, Ottawa, 1962. OCLC: 12207123.
[81] R.W. MacCormack. The effect of viscosity in hypervelocity impact cratering. In 4th Aerodynamic
Testing Conference, Aerodynamic Testing Conference, Cincinnati,OH,U.S.A., 1969. American
Institute of Aeronautics and Astronautics.
[82] J.K. Lee and T.K. Kim. Application of TVD-McCormack Scheme to Analysis of Dam-Break
Problems. Journal of Korea Water Resources Association, 36(3):365–374, 2003.
[83] P. Lax and B. Wendroff. Systems of conservation laws. Communications on Pure and Applied
Mathematics, 13(2):217–237, 1960.
[84] D. Liang, R.A. Falconer, and B. Lin. Comparison between TVD-MacCormack and ADI-type
solvers of the shallow water equations. Advances in Water Resources, 29(12):1833–1845, 2006.
[85] D. Gidaspow. Multiphase Flow and Fluidization. Academic Press, 1994.
[86] Ensayos de la velocidad de combustión. Technical report, Ministerio de Defensa, Instituto
Nacional de Ténica Aeroespacial (INTA).
[87] S.B. Margolis and F.A. Williams. Structure and Stability of Deflagrations in Porous Energetic
Materials. Technical Report SAND99-8458, Sandia National Labs., Albuquerque, NM (US);
Sandia National Labs., Livermore, CA (US). US Department of Energy (US), 1999.
235