Franc ¸ois Fraysse Richard Saurel

François Fraysse

Damien Furfaro

Emmanuel Lapebie

Multiscale multiphase modeling of detonations in condensed energetic materials

Keywords: hot spots, hyperbolic, relaxation, subscale, non-conservative, shock to detonation transition, multi-dimensional propagation

Introduction

Shocks and detonations in heterogeneous materials differ widely of similar phenomena in gas mixtures as thermal and mechanical disequilibrium are present among the phases with scales much larger than molecular ones. Chemical decomposition phenomenon is different as well, the ignition being governed by local effects (hot spots) resulting of these disequilibria. Material equations of state are obviously also very different to those of gases. Shock initiation of solid explosives is a long lasting issue. It is well known that solid temperature increase resulting of shock compression is not enough, in most situations, to achieve ignition. Part of the shock energy is focused in specific zones where local heating happens, resulting in hot spots appearance, that are considered responsible for the solid material ignition. There are several possible causes of hot spots:

-Shear bands effects, present for example in composite explosives [START_REF] Bowden | Hot spots on rubbing surfaces and the detonation of explosives by friction[END_REF][START_REF] Frey | The initiation of explosive charges by rapid shear. Army ballistic research lab[END_REF][START_REF] Saurel | A preliminary study of shock initiation of composite solid energetic material by friction mechanism[END_REF]. -Small gas pores (micrometer sizes), considered as the main ignition cause for pressed explosives [START_REF] Bowden | Studies in the transition from deflagration to detonation in granular explosives -III -Proposed mechanisms for transition and comparison with other proposals in the literature[END_REF][START_REF] Khasainov | Two-phase visco-plastic model of shock initiation of detonation in high density pressed explosives[END_REF][START_REF] Kang | A thermomechanical analysis of hot spot formation in condensed-phase, energetic materials[END_REF][START_REF] Saurel | A preliminary study of shock initiation of composite solid energetic material by friction mechanism[END_REF]. -Big gas pores and cavities, considered responsible for the ignition of highly heterogeneous explosives [START_REF] Frey | The initiation of explosive charges by rapid shear. Army ballistic research lab[END_REF]) such as civil ones, having large amount of pores (20-30% volume). -Crack propagation, intergranular friction and heterogeneities are other possible sources of hot spots [START_REF] Gilbert | Computational analysis of hot-spot formation by quasi-steady deformation waves in porous explosive[END_REF], Chakravarthy et al., 2013). The present paper considers pressed explosives only, where small gas pores are present with typical gas volume fractions of the order of 1%. Most existing models addressing shock to detonation transition have empirical basis. The most popular ones are the [START_REF] Lee | Phenomenological model of shock initiation in heterogeneous explosives[END_REF] model, the so-called Forest-Fire (FF) model [START_REF] Forest | Burning and detonation[END_REF] and the JTF model [START_REF] Johnson | Shock-wave initiation of heterogeneous reactive solids[END_REF], Tang et al., 1985). These models are widely used as they correctly predict shock to detonation transition once their parameters are calibrated. However, they have limitations:

-They usually contain many parameters: 4 and 13 parameters respectively for the FF and JTF models respectively. -They sometimes use flow variables difficult to compute, such as the shock pressure and temperature. It means that additional equations have to be added to the flow model to transport these variables. Also, shock wave detection has to be achieved, posing non trivial issues [START_REF] Menikoff | Reactive burn models and ignition & growth concept[END_REF]. -The various kinetic parameters mentioned previously are determined for given thermodynamics and given flow model. Most time the JWL EOS [START_REF] Lee | Adiabatic Expansion of High Explosives Detonation Products[END_REF]) is used for detonation products and the [START_REF] Cochran | Shock initiation and detonation models in one and two dimensions[END_REF] EOS is used for the solid reactant. Consequently the reaction rate parameters are dependent to the EOS couple and mixture model, based in most flow models on pressure and temperature equilibrium among the phases. In addition to these restrictions, these models are predictive only in the state variables range where the coefficients have been determined. Attempt to build physically based models is a long lasting issue [START_REF] Bowden | Studies in the transition from deflagration to detonation in granular explosives -III -Proposed mechanisms for transition and comparison with other proposals in the literature[END_REF][START_REF] Khasainov | Two-phase visco-plastic model of shock initiation of detonation in high density pressed explosives[END_REF][START_REF] Kang | A thermomechanical analysis of hot spot formation in condensed-phase, energetic materials[END_REF][START_REF] Saurel | A preliminary study of shock initiation of composite solid energetic material by friction mechanism[END_REF]. [START_REF] Khasainov | Two-phase visco-plastic model of shock initiation of detonation in high density pressed explosives[END_REF] and [START_REF] Kang | A thermomechanical analysis of hot spot formation in condensed-phase, energetic materials[END_REF] are local models that consider the dynamics of isolated pores to determine the characteristic time scales of hot spots ignition. The first model addressing coupling between pore scale and wave propagation at macro-scale is due to [START_REF] Saurel | A preliminary study of shock initiation of composite solid energetic material by friction mechanism[END_REF]. Limitations appeared however:

-Too many evolution equations were present, making the practical use difficult.

-The macro-scale model was based on solid-gas mixture Euler-type flow model where microscale motion was absent. Micro-macro scale motion coupling has been addressed by van Wijngaarden (1968) in the context of bubbly liquids with a model involving a single velocity. Extension to two-phase mixtures in velocity A C C E P T E D M A N U S C R I P T 3 disequilibrium was addressed in [START_REF] Drumheller | Transient wave propagation in bubbly liquids[END_REF], [START_REF] Bedford | Theories of immiscible and structured mixtures[END_REF] and [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF]. In the present work micro inertia is neglected, as assumed in [START_REF] Baer | A two-phase mixture theory for the deflagration-to detonation transition (DDT) in reactive granular materials[END_REF] type models. Here the symmetric variant of [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF] is considered and the pressure relaxation process is used to model pore collapse. The visco-plastic heat deposited at pore surfaces is determined through a simple approach based on 'configuration' pressure (Passman et al. 1984, Baer and[START_REF] Baer | A two-phase mixture theory for the deflagration-to detonation transition (DDT) in reactive granular materials[END_REF]) and radial pore collapse velocity, determined through a pressure relaxation process. This method considerably simplifies the approach of [START_REF] Khasainov | Two-phase visco-plastic model of shock initiation of detonation in high density pressed explosives[END_REF], [START_REF] Kang | A thermomechanical analysis of hot spot formation in condensed-phase, energetic materials[END_REF] and [START_REF] Saurel | A preliminary study of shock initiation of composite solid energetic material by friction mechanism[END_REF] as no radial integration normal to the pore interface is needed. Relaxation terms combined with the 'granular equation of state' summarize 3D subscale motion and associated dissipation. Efforts in the same direction of simplification were done formerly by [START_REF] Gonthier | Predictions for weak mechanical ignition of strain hardened granular explosive[END_REF], [START_REF] Hamate | Ignition and detonation of solid explosives: a micromechanical burn model[END_REF] and possibly other researchers. From pore collapse, heat is deposited at pore surface and hot spots are created. These hot spots may be subcritical or supercritical depending on heat losses. Dominant heat loss occurs in the solid phase, between the (hot) pore interface and the (cold) solid core. Its determination usually needs resolution of the heat equation and this is problematic as the solid temperature varies as a consequence of shock propagation and compressibility effects as well as subscale heat diffusion. A method is developed herein to compute heat losses without subscale resolution. It uses on one hand the time varying volume average solid temperature, determined from the solid phase balance energy equation, and the solid core temperature that varies as a consequence of compressibility and shocks. Its precise evolution is determined through a non-conventional energy equation, or more precisely through two-phase shock jump relations. Having in hand the mean solid temperature and the solid core temperature, the interfacial heat flux (heat loss) is determined from an approximate profile expressing correctly the physics of heat transfer. At this point, both subscale heat production and heat loss are computed without subscale resolution. These heat sources enter in the definition of the interface temperature that now involves all relevant contributions reported in shock ignition sensitivity studies: pore size, shock pressure, thermal conductivity, subscale plastic stress, pressure relaxation rate (related to pore collapse velocity). Transient evolution of these variables directly enters in the interface temperature formulation. This temperature is used as criterion to initiate solid deflagration, modeled by Vieille's law. Collective effects of subscale deflagrations are summarized at global scale, in the two-phase model, through the specific interfacial burning surface. At short times, total pore surface is used in the global mass transfer rate to compute the burning surface. After some time, the fluidization limit of the granular bed is reached and burning continues through the outer grains surface. Last the model is extended to deal with quasi-steady detonation propagation in explosive cylinders of variable diameters. To this end, modification of the decomposition kinetics with burning area increase is addressed. With this modification good agreement with experimental measurements of detonation velocity as a function of charge diameter is obtained. The full model is thus able to deal with both shock to detonation transition, occurring in pressure ranges of 50-100 kbar, and detonation propagation regime, occurring in pressure ranges of 300-1000 kbar. The paper is organized as follows. The macro-scale flow model is presented in Section 2, as an extension of the [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF] model, this one being itself a modification of the [START_REF] Baer | A two-phase mixture theory for the deflagration-to detonation transition (DDT) in reactive granular materials[END_REF] model. The various closure relations are provided in the same section. Section 3 analyses pore collapse dynamics and hot spot creation through relaxation effects of the former model. Section 4 analyses the extra evolution equations used to determine the solid core temperature. This temperature is used to compute the solid-gas interface temperature as well as the interfacial heat flux. Section 5 addresses validations of shock to detonation transition in 1D for PBX9501 explosive. A sensitivity analysis of the flow model to the various parameters is done in the same section. Section 6 deals with extra extensions and validations to render the model able to deal with detonation propagation in multi-dimensions. Validations against experimental data of velocity of detonation as a function of charge diameter are addressed in the same section. Conclusions are given in Section 7.

A C C E P T E D M

A N U S C R I P T

Multiscale Multiphase Model

The flow model of [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF] is considered as symmetric variant of the [START_REF] Baer | A two-phase mixture theory for the deflagration-to detonation transition (DDT) in reactive granular materials[END_REF] two-phase flow model. It is complemented by extra equations used to compute pore and grain radii, important for interfacial area determination. Extra equations are also added to determine solid core temperature, important for the computation of conductive heat loss at pore scale. For the sake of simplicity the equations are given in 1D.

Solid phase

ss I s g I m u ( ) tx             (2.1) ss ( ) ( u) m tx       (2.2) 2 ' s s s I D I ( u) ( u p) F mu t x x                (2.3) ' s s s I I D I I s g I ( E) ( ( E p)u) u F u ( ) mE t x x                     
(2.4)

ss ss ( n) ( nu) ( ) n tx         (2.8) Gas phase gg ( ) ( u) m tx         (2.9) 2 g g g ' I D I ( u) ( u p) F mu t x x                (2.10) g g g ' I I D I I s g I ( E) ( ( E p)u) u -F u ( ) mE t x x                     (2.11) gg gg ( n) ( nu) ( ) n tx        
(2.12)

In this system

    k k k k k k k k , , u , p ,e , E , B , n k s,g   
represent respectively the volume fractions, the densities, the velocities, the pressures, the internal energies, the total energies, the 'granular' or 'configurational' energies, the specific number of grains and pores per unit mass, The subscale physics in this model is expressed through equations (2.1), (2.5), (2.6), (2.7), (2.8), (2.12) and will be examined in detail later. The sources terms 0 s p and 0 s v represent the Hugoniot pole reset, important when multiple shocks are under consideration. These source terms are examined in Section 4. The fragmentation rates s n and g n are important as well and are directly connected to the burning surface. The interfacial area model is detailed later in the present section. Mass transfer is represented by m (negative for the solid phase). The total energies of the phases are defined as,

2 k k k k k 1 E e B ( ) u 2     . (2.13)
The granular pressure are defined as, .14) and exists only in the solid phase. In the gas phase g 0  .

k k k k k dB d      , ( 2
The 'granular energy' s B , represents the specific energy stored in the deformed layers around grains.

Following [START_REF] Saurel | Modeling dynamic and irreversible powder compaction[END_REF] 1) ln( 1)

(1 ln(1 ))( ) (1 ) ln(1 ) if 1 B ( ) 0 otherwise                          
(2.15) where 0 s

 represents the threshold value of the solid phase from which the granular pressure and energy depart from zero. Parameters gr a , gr n and 0 s  are material dependent. The equation of state (2.14) is determined from pressed granular beds experiments [START_REF] Kuo | Intergranular stress particle-wall friction and speed of sound in granular propellant beds[END_REF]. The present work considers condensed materials with low porosity, but conceptually the same equation of state is considered, with a smaller range of variation of the volume fraction. It is worth to mention that conductive heat exchanges are absent in equations (2.4) -(2.11). As will be shown later, heat deposition at the interface is considered through the interstitial power ' I s g ()     that induces heat diffusion from the interface to the solid phase. After pore collapse and solid ignition, consideration of solid -gas heat exchange becomes of minor importance compared to mass transfer effects and is therefore neglected.

The following interfacial variables render the model symmetric with respect to the phases indexes:

s s g g g s s I s g s g Z u Z u u sgn Z Z x Z Z              , s s g g ' I s g Z u Z u u ZZ    , (2.16) g s s g s g s I g s s g s g Z Z Z Z sgn (u u ) Z Z x Z Z               , g s s g ' I s g Z Z ZZ      ,
where

k k k Z c 
represents the acoustic impedance of the phase k. k c represents the conventional sound speed of phase k. The symmetric formulation of the BN model has some advantages. As these formulas correspond to local 'granular Riemann problem' solutions [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF], fluid-fluid interfaces, fluid-granular media interfaces as well as permeable interfaces are handled properly by the non-conservative terms

kk II u, xx     (2.17)
Although not directly linked to hot spots ignition, it is worth to mention that computation of material interfaces can be achieved with System (2.1-2.12). In [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flows[END_REF] stiff mechanical relaxation was used to enforce interface conditions, while in Le [START_REF] Metayer | Modelling evaporation fronts with reactive Riemann solvers[END_REF] nonconservative terms only (2.17) were used in the same aim. In this context there is no need to use interface reconstruction, Level Set, or any Front Tracking method. At interfaces, when a volume

A C C E P T E D M A N U S C
( ) s ( ) s ( ) u s ( ) u s t x Z Z Z 1 ( ) sgn Z (u u ) (u u ) ( ) T x x Z Z Z Z ZZ Z Z ( ) sgn Z (u u ) x x ZZ 1 T                                                                                g 2 2 g s g s s g s g s g s g s s g s I Z (u u ) ( ) Z Z Z Z 1 m g g (T T )s B ( )                                       (2.18)
where k k k k g h T s  denotes the Gibbs free energy of phase k. Except regarding the entropy production due to mass transfer that is a solid-gas irreversible reaction, all terms are symmetric. Last, the model is hyperbolic with 7 wave speeds (instead of 6 in the original BN model), this feature having importance regarding numerical resolution [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible nonequilibrium two-phase flows[END_REF].

Interfacial area

The specific interfacial area through which mass, momentum and energy exchanges occur at the pore external surface is determined from equation (2.12). g n represents the number density of pores per unit mass of gas. Obviously,

gg ( n) N  
, where g N represents the specific number of pores per unit volume. From the solid volume fraction determined by equation (2.1) and number density g n there is no difficulty to determine the sizes of inclusions that vary with both time and space, as a result of compressibility, dilatation and mass consumption. Considering spherical pores, their radius g R is determined by:

3 g g g 4 RN 3    .
(2.19)

The specific interfacial area g I A is deduced as:

g 2 I g g A 4 R N  . (2.20)
During the ignition and growth stages, combustion occurs through the pore surface and continues through the outer grains surface. The switch is based on the fluidization limit of the granular bed, estimated to occur (quite arbitrarily) at 90% ( SW s 0.9 

). This parameter is one of the set of adjustable parameters of the flow model. Continuity of the combustion kinetics, between inner pore combustion to outer grain surface, is enforced at the volume fraction switch SW s  by imposing equality of the interfacial areas. Indeed, as the local mass transfer rate, modelled through Vieille's law (Eq. 2.22) is considered constant during both combustion modes, the burning surface must be continuous to preserve continuity of the various flow variables. This is done through the following identity,

sg 2 s I s s I s 3 A 4 R N A R      ,
that defines the grains radius when the combustion regime changes,

A C C E P T E D M A N U S C R I P T 7 SW sg s SW s R R 1    .
Once the grains radius is determined from the relation above it is transported with the help of (2.8).

Interfacial area creation

To predict both detonation ignition and growth as well as quasi-steady 2D propagation in explosive charges of varying diameters, dynamic fragmentation and interfacial area creation have been addressed. As the shock pressure varies between 50 kbar during the ignition stage to 500 kbar during the propagation regime, an interfacial area creation factor (IACF) has been introduced as a function of the maximum pressure recorded by the explosive reactant. This maximum pressure corresponds to the shock pressure during the quasi-steady propagation stage and to the local pressure during the ignition stage. The IACF function is shown in Figure 1, and obeys the following formula: f represents the maximum interfacial area amplification factor. The interfacial area used in the mass transfer rate, during both pore and grains combustion is thus modified as:

max 2s cp max s1 f (p ) c e  (2
k k IACF max 2 max I s k k s I A f (p )4 R N f (p )A    .
During the first stage corresponding to pore combustion The IACF is aimed to model ignition of pores of decreasing sizes when the shock pressure increases. Indeed, various sizes of pores are present in an explosive sample, as reported for example in [START_REF] Mang | An optical microscopy and small-angle scattering study of porosity in thermally treated PBX9501. Shock Compression of Condensed Matter Massoni[END_REF].

It also models grains fragmentation as a function of the maximum pressure seen by a given volume of solid reactant.

A C C E P T E D M A N U S C R I P T

The use of max s p as the maximum pressure between the local one and the shock pressure enables evolution of the IACF function from SDT to the quasi-steady detonation propagation regime. Indeed, at early stages of ignition, the solid material subjected to projectile impact sees the shock pressure computed with (2.7). At later times max s p takes the local pressure value, as during SDT, the pressure in the post-shock reaction zone reaches higher values than the shock pressure. Finally, when quasisteady conditions are reached max s p corresponds to the Neumann spike pressure, transported by equation (2.7). Note that equation (2.7) needs a shock detection algorithm that will be described in Section 4. The interfacial area creation model involves both activation of small pores and fragmentation of grains. Their sizes reduction follows approximately the inverse rate of the interfacial area creation.

Pressure relaxation

Following [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF] the pressure relaxation parameter is a consequence of local Riemann problem solutions around pores and grains. When the mixture is made of compressible ideal fluids it reads,

Ig s g A Z Z   .
In the present context, an accommodation factor acc  is needed as an elastic-plastic solid is present.

The effective pressure relaxation parameter is thus corrected by an accommodation coefficient,

I g acc s g A Z Z     , acc
 being one of the adjustable parameters of the flow model.

Drag exchanges

The specific drag force D F appearing in equations (2.3), (2.4), (2.10) and (2.11) is given by,

s D I D g s 1 F A f m (u u ) 2      
, where D f is estimated as the acoustic (or pressure) drag coefficient [START_REF] Chinnayya | Modelling detonation waves in heterogeneous energetic materials[END_REF],

sg D sg ZZ f ZZ   .
This estimate seems appropriate in the present context rather than correlations of porous media, determined in quasi-steady low speed flow conditions. Moreover, the present modeling doesn't introduce extra parameters and seemed accurate enough as will be observed with computed results.

Mass transfer

Mass transfer is modeled through Vieille's law with same parameters for both pore and grain burning modes:

k IACF I I s m A R  , v n s v mix R a p  , (2.22)
where a v and n v are material dependent constants and mix p is the mixture pressure defined as During pore collapse, the interface temperature is determined by a subscale model, as will be examined in the next section.

The interfacial density appearing in (2.1) follows [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF] and is given by,

2 2 gg ss sg I 2 2 gg ss s g s g C C c dB c dB dd             .
This estimate was derived for flows with mass transfer in the stiff mechanical relaxation limit, a situation quite close to the present context. Last, the 'interfacial' energy is given by [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF]:

gs ss I s I s sI uu p E E (u u ) 2          .

Equations of state

The solid phase thermodynamics is modeled by the Cochran-Chan (CC) equation of state, whereas the Jones-Wilkins-Lee (JWL) equation of state is used for detonation products. Details are given in Appendices A, B and C. Equipped with these relations, System (2.1-2.12) is able to model flows of reactive granular materials. But hot spots creation through pores collapse need subscale modeling, as addressed in the next section.

Pore collapse and hot-spot ignition

The aim is to model shock ignition of low porosity energetic materials (typically of the order of 1%). In the present section focus is done on ignition resulting of subscale pore collapse. The sequence of events considered in the ignition and growth process is summarized in Figure 2 as follows:

-The energetic material is initially at rest and modeled as a multiphase mixture initially in mechanical and thermal equilibrium. -A shock wave is emitted by the impact of a projectile on the explosive. It propagates in the multiphase mixture. -Immediately after shock passage the mixture is in pressure, velocity and temperature disequilibrium. -Pore collapse occurs as a consequence of pressure disequilibrium.

-Pore interface temperature increases as a consequence of local dissipation. It results in solid material ignition if the interface temperature exceeds a given threshold limit. Pore collapse being a consequence of pressure disequilibrium effects, System (2.1-2.12) is consequently well suited to model it, as examined hereafter.

Configuration energy deposition

Let us examine System (2.1-2.12) at a given point of space after shock propagation. During pore collapse and pressure relaxation, the system to examine reduces to:

s sg () t        , ss I s g ( ) (e B) () t           , gg I s g ( ) (e B) () t           , (3.1) s ( ) cst1   , g ( ) cst2   , s ( u) cst3   , g ( u) cst4   , s ( n) cst5   , g ( u) cst6   .
Heat exchange between the gas phase and the interface has been neglected compared to the interstitial power

I s g
()     as a consequence of low thermal conductivity of the gas phase [START_REF] Jackson | Pore collapse in an energetic material from the micro-scale to the macro-scale[END_REF]. Heat exchange in the solid phase doesn't appear explicitly in the energy equation of phase s but is implicitly present, as examined in the next paragraph. The interface condition (3.9) shows balance of the interstitial power deposited on the interface during pore collapse with the heat exchange between the interface and the bulk solid. The sum of these two effects being zero, they are absent in the solid balance energy equation. Thanks to the interfacial pressure estimate, When pressure relaxation is stiff, the production term becomes negligible and each phase evolves isentropically [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF]. But in the present context, finite rate pressure relaxation is considered.

During time increment t  , the configuration energy change per volume unit reads:

      c s s s s s Q ( ) B (t t) B (t)        
Dividing this quantity successively by the number of pores per unit volume g N , the pore surface and the time increment, it becomes,

      s s s s s c c 2 p g int g ( ) B (t t) B (t) Q q S N t 4 R N t            (3.3)
This term corresponds to the local specific power released during pore collapse. It can also be understood as,

A C C E P T E D M A N U S C R I P T 11 3 int g c s s s s s s cs 2 2 2 p g int g int g s int g 4 RN Q ( ) B ( ) dB 1 3 q ( ) S N t 4 R N t 4 R N d t 4 R N t                     
, where int R denotes the pore radius, as shown in Figure 3.

After simplifications it becomes,

int cs R q t     (3.4)
This last formula provides the power of 'configuration forces' released at the interface (deposited heat) as the product of the interfacial velocity

() t        , transforms to, int s g int g R ( ) R t3        .
(3.5)

The deposited heat alternatively reads,

int s c s g g R q ( ) 3        , ( 3.6)
c q represents the heat source in the pore collapse mechanism. This source is in competition with heat loss, whose determination is addressed in the following paragraph.

Heat loss

In order to avoid subscale numerical resolution of the heat equation in the solid phase (a non-trivial issue when pressure wave propagates in compressible materials), let us consider an approximate method to determine heat loss. Heat loss between the interface and the gas phase are omitted as a result of low thermal conductivity of the gas phase and heat diffusion in the solid phase only is considered. Let us assume a temperature profile in the solid phase of the form [START_REF] Baer | A model for interface temperatures induced by convection heat transfer in porous materials (SAND-88-1073)[END_REF]:

n c I c int r T T (T T ) R         if int ext R r R     (3.7) c T T  otherwise   int R 
represents the thermal boundary layer thickness (to be determined), and n a parameter which determines the shape of the temperature profile (Figure 3). Dependence of the solution to this parameter will be addressed later by comparing computed profiles to one-dimensional heat equation solutions.

The 'core temperature' c T corresponds to the temperature of the solid phase outside the thermal boundary layer, far from the interface. Its determination will be addressed later with the help of equations (2.5)-(2.6).

A C C E P T E D

n1 Ic n int (T T ) T n ( r) r ( R )          
At the pore surface it becomes,

int Ic rR int (T T ) T n r ( R )           .
It can be checked easily that this derivative is negative. The heat flux at the solid surface thus reads,

int Ic s s r s r rR int (T T ) T q u n u r ( R )              , ( 3.8)
where r u is the polar vector.

Before solid ignition, mass ejection from the solid-gas interface is absent. Therefore the energy interface condition is, s r c q .u q 0   

(3.9) where c q is the heat deposited at the interface. The conductive heat flux between the interface and the gas phase has been neglected as a consequence of low thermal conductivity compared to the solid one. The interface condition (3.9) is local and expresses balance of plastic dissipation with heat conduction in the solid phase. Inserting (3.7) and (3.8) in (3.9) the interface temperature is obtained as:

int s I c int s g gs R T T ( R ) ( ) 3n           (3.10)
Formula (3.10) is an interesting interface temperature model:

• It increases when the solid pressure ( s  ) increases, in agreement with pressure sensitivity of explosives. As the solid core temperature c T increases identically with pressure, the same tendency is expressed by the two terms appearing in the sum in (3.10).

• It increases when the pore collapse is stiff (large  ), in agreement with pore collapse viscosity models [START_REF] Saurel | A preliminary study of shock initiation of composite solid energetic material by friction mechanism[END_REF].

• It increases when the thermal boundary layer int ( R )  increases too. So it increases versus time as the thermal boundary layer is time dependent.

• It decreases when the solid conductivity is high, as heat deposited at the interface is diffused in the solid phase.

• It decreases for pores of small radius int R , in agreement with experimental observations (Borne andBeaucamp, 2002, Borne et al., 2008).

Thus, formula (3.10) contains most observations of solid explosives sensitivity reported in the literature.

However, for practical use, (3.10) needs solid code temperature c T and thermal boundary layer thickness int ( R ) 

. Their determination is addressed in the following section.

Subscale heat transfer

In this section we address computation of the thermal boundary layer thickness that appears in the temperature profile (3.7) and interface temperature model (3.10). The method is then validated against unsteady solutions of the heat equation.

The solid core temperature c T appearing in the same formulas is addressed later.

Thermal boundary layer thickness

The thermal boundary layer thickness is determined on the basis of:

-The solid mean temperature T knowledge, determined from the balance equation (2.4), as done usually in multiphase flow computations.

-The solid core temperature c T , assumed known at this level.

-The approximate temperature profile (3.7). The approximate temperature profile (3.7) is inserted in the volume average solid temperature definition: 3 is linked to the internal pore radius int R and the gas volume fraction g  through the relation:

R R ext ext 22 33 33 ext int V R R int int ext int 1 1 3 T TdV T4 r dr Tr dr 4 V (R R ) (R R ) 3          (4.1) In this expression, ext R shown in Figure
3 int g ext R R

   

Inserting (3.7) the following result is obtained: . In that situation the approximate temperature profile is inappropriate and long term solution has to be considered. Its determination is addressed in Appendix D. As the temperature profile (3.7) contains a parameter (the exponent n) a sensitivity analysis of the solution to this parameter is needed as well as comparisons with solutions of the 1D heat conduction equation.

          2 23 int ext ext 2 ext c int int int int int 33 2 ext int 0 2 int int s int int s g gs R RR R TR R R 3 3R 3(R ) T (R R ) R R 2 ( R ) ( R ) ( ) 3 n n 2 n 1 n 3                                                          

Determination of the exponent n of the temperature profile

The temperatures obtained with the short and long term profiles presented formerly are compared to the numerical solution of the heat equation in the solid hollow sphere. To compare the solutions, the 'configuration heat flux' c q is taken constant.

Let us consider the following configuration corresponding to slightly porous PBX9501 explosive with 5.14 GPa incident shock pressure:

 . The practical use of the reduced heat conduction model needs the solid mean temperature as input. It is usually provided by the multiphase flow model through the solid phase balance energy equation (2.4). Here an isolated hollow sphere of solid is considered instead. In this simplified situation, the hollow sphere of solid mean temperature is given by the following differential equation,

p pore c Vsolid d C TdV S q dt   , resulting in, 2 int c 33 ext int p Rq d3 T dt (R R ) C  
, where  and p C represent respectively the solid density and specific heat at constant pressure.

As in the present test problem all parameters in the right hand side of this equation are constant, its integration is immediate:

2 int c 0 33 ext int p Rq 3 T(t) T t (R R ) C  
The interface temperature and interfacial heat flux determined from the 1D heat equation, ).

ps T 1 T C r² t r² r r            
Good agreement between the heat equation solution and the reduced model is observed for large exponents n (=4, 5). Even if some differences are reported in the dimensionless interfacial heat flux at early times, the deviation does not exceed a few percent. Only weak dependence to the parameter n is observed, provided it is large enough. Consequently the reduced heat model is considered as an approximate method free of parameter. The exponent n5  is used in all computations shown in the forthcoming section.

The last issue with subscale heat transfer is related to the solid core temperature determination.

Solid core temperature determination

During shock and pressure wave propagation the solid temperature varies. Close to grains interfaces the various energy exchanges induce temperature variations, summarized in the balance solid energy equation (2.4). But far from interfaces, a core zone is present and evolves adiabatically during some short duration. This duration is function of the thermal boundary layer evolution, determined previously. Knowledge of the solid core temperature is needed only during this stage. The core temperature is determined from the phase Hugoniot relation [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF]:

0 0 0 0 0 ss s s s s s s s s pp e (p , v ) e (p , v ) (v v ) 0 2      , ( 1 v   ) (4.3)
Given pre-(denoted by superscript 0) and post-shock values, the solid phase core internal energy is determined. The use of (4.3) needs knowledge of the Hugoniot pole (4.5) Otherwise, the method that follows has to be used and needs shock detection. Shock detection is a non-trivial issue [START_REF] Menikoff | Shock Detector for SURF model[END_REF]Shaw, 2010, Menikoff, 2016). In the following we use the simplified algorithm that follows. A shock is assumed to be present in the solid phase if, for a given cell i,

  n 1 n s,i s,i p 1 p     (4.6)
where  is a parameter. Numerical experiments show that 0.01  is a fair choice. When a shock is present, or more precisely when the current computational point belongs to the shock layer, pole 00 ss (p , v ) 4.4) contain relaxation terms, expressing Hugoniot pole reset. These relaxation terms are not expressed herein for the sake of simplicity, as they involve shock detection that is more a computational issue than a modeling one.

The overall computational algorithm that couples micro and macro scale models during pore collapse is summarized in Figure 5.

Computed results

In this section one-dimensional numerical experiments are considered to compute ignition and growth of a detonation wave under shock interaction, as illustrated in Figure 2. Quasi-steady propagation regime in 2D, with a study of charge diameter effects is addressed subsequently. All computations are done with the Furfaro and Saurel (2015) HLLC solver, embedded in a Godunov type scheme, related details being given in the same reference.

Shock to detonation transition (SDT) in 1D

Careful experimental analysis of shock to detonation transition in PBX9501 explosive has been done in [START_REF] Gustavsen | Shock initiation of new and aged PBX 9501 measured with embedded electromagnetic particle velocity gauges[END_REF] with the configuration depicted in Figure 6. A 2.5cm long explosive cylinder made of PBX9501 with 0.017 porosity filled with air is impacted by a projectile in vistal (high density aluminium) of width 1.1cm at 817m/s resulting in a 51.4 kbar shock wave. Lagrangian gauges are settled in the explosive to record local velocity. These gauges are located at abscissas 0.51 mm, 1.01 mm, 1.5 mm, 2.02 mm, 2.49 mm, 3 mm, 3.48 mm, 3.99 mm, 4.47 mm and 4.98 mm of the impact surface. The EOS parameters used for PBX9501 solid and gas phases are given in Appendixes A, B and C. The vistal projectile is modeled as a fluid through the stiffened gas EOS with parameters given in Table 1. In particular, the vistal EOS parameters have been determined from the Hugoniot given in [START_REF] Gustavsen | Shock initiation of new and aged PBX 9501 measured with embedded electromagnetic particle velocity gauges[END_REF]. Vistal The physical properties of PBX9501 used in the present study are given in Table 2. Thermal conductivity is taken from [START_REF] Gibbs | LASL Explosive property data[END_REF], pore radius from [START_REF] Mang | An optical microscopy and small-angle scattering study of porosity in thermally treated PBX9501. Shock Compression of Condensed Matter Massoni[END_REF] and ignition temperature is taken from [START_REF] Henson | The phase transition in the energetic nitramine octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine: Thermodynamics[END_REF]. Vieille's law parameters are adjustable parameters of the model as their determination is accessible in the kilobar range but hardly in detonation conditions. Granular EOS parameters are obtained through quasi-steady granular beds compaction experiments [START_REF] Kuo | Intergranular stress particle-wall friction and speed of sound in granular propellant beds[END_REF][START_REF] Saurel | Modeling dynamic and irreversible powder compaction[END_REF]. Extrapolated data from this last reference have been used in the present work. The pressure relaxation rate accommodation factor and fluidization limit are arbitrary parameters in the present study. However, their determination can be addressed through mesoscale computations [START_REF] Chakravarthy | Analysis of microstructure-dependent shock dissipation and hot-spot formation in granular metalized explosive[END_REF]. SDT of PBX9501 is considered first. The mixture flow velocity is recorded at the various Lagrangian gauges and compared to the experimental result shot #1133 of [START_REF] Gustavsen | Shock initiation of new and aged PBX 9501 measured with embedded electromagnetic particle velocity gauges[END_REF] in Figure 7, without (top) and with (bottom) interfacial area creation factor (IACF). Three distinct stages are visible:

PBX9501

 First a velocity plateau is visible, especially on the first two gauges. Only slight velocity increase happens, characteristic of the pore collapse stage and related ignition where the reaction is weakly exothermic.  The second part corresponds to the first growth of the detonation stage. The combustion occurs at the internal pore surface and increases stiffly as visible on gauges 3, 4, 5 and 6. The combustion surface increases during this stage resulting in rapid growth of the reaction. Then the fluidization limit is reached and combustion continues by the external grain's surface.  Typical velocity profiles of detonations gradually appear at signals recorded by gauges 7, 8, 9

and 10, where both pore collapse and inner pore combustion happen immediately after shock propagation and then followed by grain's combustion. SDT computations are in good agreement with the experiments when grain fragmentation is taken into account through IACF, as well as small pores activation when the shock strength increases. It is clear that in the absence of the IACF, large discrepancies are visible. Indeed, the first gauges match well the experiments indicating that Vieille's law parameters are well chosen. This stage is mainly governed by inner pore combustion, and corresponds to the ignition of the biggest ones. When the shock strength increases, pores of smaller size become activated as well (through the IACF function), increasing the total burning surface. Same surface increase is considered for grains that become smaller as a consequence of fragmentation effects when the shock strength increases. Computed (x-t) diagram for the leading wave of the PBX9501 test case is compared to the experimental records of [START_REF] Gustavsen | Shock initiation of new and aged PBX 9501 measured with embedded electromagnetic particle velocity gauges[END_REF]. The corresponding comparison is given in Figure 8 showing very good agreement. Convergence of the computations is checked by considering mesh refinement, as shown in Figure 9.

The computed steady reaction zone has a width of approximately 120µm handled with 20 computational cells with the coarsest mesh.

Interface temperature sensitivity analysis

Figure 10 shows the various temperatures involved in the subscale model: interface, outer pore boundary and solid average temperature when the same explosive sample is subjected to the same impact and incident shock wave. Mass transfer has been removed to simplify the illustrations. While the outer pore boundary and average temperatures of the solid phase have closed evolutions during pore collapse, the interface temperature is strongly different, as expected. This is the result of heat deposition at the interface during pore collapse. Short term / long term switch can be observed on the temperature profile of the outer pore boundary showing continuity of the solution. 2.

This set of parameters yields an interfacial temperature slightly greater than the ignition threshold of 550K, in agreement with the minimum shock pressure where ignition happens (23 kbar). From this configuration, the various parameters are varied one by one in order to demonstrate their individual influence on the interfacial temperature whose expression is recalled hereafter:

0 int s I c int s g gs R T T ( R ) ( ) 3n          
In Figure 12, same computations are done with smaller initial pores radius ( int R 0.3 m  ). In Figure 13, the thermal conductivity is increased to In Figure 14, the accommodation factor in the pressure relaxation parameter is lowered to acc 0.005  . As expected, variations in these parameters have important consequences to the results. The former three parameters are consequently first order ones and particularly important for the explosive ignition threshold, as already known in the literature.

The pressure relaxation accommodation parameter acc 0.01  is chosen such that the material becomes insensitive for lower shock pressures. A test has been performed using parameters from Table 2 and a lower velocity impact of 320m/s resulting in a 19 kbar shock. Results are displayed in Figure 15. 2 and lower impact velocity (320m/s). Change of this parameter renders the explosive insensitive in the present impact conditions.

Steady detonation

We now check model's ability to compute quasi-steady detonation regime in 1D. PBX9501 is considered and associated pressure profiles are shown in Figure 16 at various times. We now address multi-D computations with the same flow model to check its ability to compute front curvature effects and velocity of detonation as a function of charge diameter.

2D explosive cylinder computations

Two-dimensional axisymmetric configuration is now considered to examine the ability of the flow model to deal with multi-D propagation regime. The pressure level during propagation is about ten times higher than during the ignition stage. There is thus no guarantee that the approach developed for ignition be also able to deal with propagation. In the author's knowledge, it is rare to achieve both ignition and steady propagation with a single kinetics and flow model. The test under consideration consists in several explosive cylinders of decreasing diameter connected to a primary cylinder as shown in Figure 18. The explosive is initially at atmospheric state and surrounded by air in the same condition. The first cylinder is impacted by a vistal projectile at 3km/s. A detonation wave forms and propagates in the next cylinder. As the explosive and detonation products are in contact with the ambient air, strong rarefaction waves appear. They result in local pressure decrease which in turn decreases the shock strength and pore collapse dynamics. In the limit, the pores interface temperature becomes lower than the ignition temperature, resulting in detonation failure.

In order to capture accurately the reaction zone, a highly anisotropic mesh composed of quadrangles is used. The cylinder under consideration is meshed with 6800 points in the x-direction, corresponding to a minimum of 10 points in the reaction zone whereas the y-direction is resolved with 20 points only. The total number of elements for the whole computational mesh is about 1 500 000 quadrangles.

All parameters of the flow model, including those of the IACF are the same as before. Detonation propagation in small two cylinders of diameters 1.667mm and 1.538mm, connected to a larger one, is shown in Figure 19 at times 2.5µs and 4.3µs. It can be seen on the 1.538mm cylinder snapshot at 4.3µs that the reaction zone thickens suddenly resulting in detonation failure. 9.12 10 u ( ) 8802 1 1.9 10 , u mm.µs 1 0.48

             
where  is the inverse charge diameter. Good agreement is obtained with the reference IACF parameters used also for SDT computations of Section 5.1. Computed critical diameter is 1.667mm whereas the experimental one is 1.52mm. Without IACF, computed critical diameter is 10mm and the velocity of detonation is about 4000m/s, in very bad agreement with the experiments, showing the importance of dynamic burning surface modelling. The present flow model is thus able to deal with both ignition stage and propagation regime using the same set of parameters, all having physical meaning or at least physical representation.

Conclusion

A coupled multiscale multiphase flow model has been built on the basis of [START_REF] Baer | A two-phase mixture theory for the deflagration-to detonation transition (DDT) in reactive granular materials[END_REF] model. The subscale model uses pore collapse and pressure relaxation to determine local heating source responsible for hot spot creation. These hot spots may be sub or super critical depending on heat loss between the hot interface and cold solid core. Associated heat transfer is determined with a subgrid model based on grain's core temperature, determined itself from two-phase Hugoniot relations and volume average temperature. The corresponding subgrid model is free of parameter. The multiscale flow model is validated against shock to detonation transition experiments of pressed PBX9501, with low initial porosity (about 1%). It is then illustrated on multi-D detonation computations showing its capabilities for the computation of velocity of detonation as a function of charge diameter. Such capability is rendered possible with an extension dealing with burning surface creation versus shock strengthening. When the explosive is correctly characterized by, -condensed phase and detonation products EOS; -initial size of pores (the largest ones being responsible for ignition), the overall flow model involves 6 adjustable parameters all having physical meaning.

These parameters are:

 The solid thermal conductivity, not necessarily known in detonation conditions.  Two Vieille's law parameters appearing in Equation (2.22). These parameters can be determined in the kilobar range, but difficultly in detonation conditions.  The fluidization limit of the granular bed.  Two parameters in the IACF function. Numerical predictions with the flow model are quite sensitive to these parameters, especially when both detonation ignition and propagation regimes have to be reproduced. As a perspective, the same flow model will be extended to hot spots creation by intergranular friction. Cochran-Chan parameters for PBX9501 material are summarized in the Table A.1: Appendix B -Detonation products: JWL EOS The Jones-Wilkins-Lee (JWL) EOS [START_REF] Lee | Adiabatic Expansion of High Explosives Detonation Products[END_REF] and have been determined with the Cheetah 2.0 thermochemical code (Fried et al., 1998) To render compatible the Cochran-Chand and JWL EOS during mass transfer special care is needed in the computation of various energy constants. This issue is addressed in the appendix that follows.

A C C E P T E D

M A N U S C R I P T         ks ks e e v P v,e P v v   where,   1 E 1 E 12 0 1 Re f 2 Re f ks vs Re f s 1 Re f 2 Re f A v A v vv e v 1 1 C T e 1 E v 1 E v                                      ,    

Appendix C: Heat of reaction -Determination of the solid reference energy

In the frame of [START_REF] Kamlet | The chemistry of detonations. 1 -A simple method for calculating detonation properties of CHNO explosives[END_REF] In mass units, it becomes, . Therefore, two relations are needed. The solid core temperature is no longer provided by the Hugoniot as thermal boundary layer has penetrated the entire solid phase.

(R R ) R nT R R ( ) 3 Y Z R R 3 T nY YZ (R R ) T Z T 3 R R T YZ                     (D.4)
The long terms interface temperature and fluxes are thus determined for given n and T .

The

(R R ) Y R R 3 R 2R (R R ) (R R ) Z n 1 n 2 n 3                   

  . At mechanical equilibrium sg    and sg uu  . The rate at which mechanical equilibrium is reached is controlled by the relaxation parameter  and drag force D F .

  being based on the volume fraction SW s  . This fragmentation law thus requires two parameters, a threshold pressure max p (approximately the Neumann spike pressure) and a maximum value max f .

Figure 1 :

 1 Figure 1: Interfacial area creation factor as a function of pressure.

  is computed from the pores or grains radii depending on the combustion mode, based on SW s  , and scaled by the IACF. Mass transfer is activated when the interface temperature exceeds a given threshold limit, A C C E P T E D M A N U S C R I P T I ign TT  .

Figure 2 :

 2 Figure 2: Schematic representation of pore collapse, ignition and growth of reactions in a weakly porous energetic material.

  result being obtained for the gas phase). Equation (3.2) indicates that entropy is mainly produced in the solid phase during pore collapse, as a consequence of large acoustic impedance s Z (present in the ratio

Figure 3 .

 3 Figure 3. Pore geometry and associated definitions. Schematic representation of the thermal boundary layer thickness   int R 

  (4.2) Equation (4.2) is a non-linear function of  , which is solved numerically (the Ridders method is used A C C E P T E D M A N U S C R I P T in practical computations). Solving (4.2) for given exponent n provides  and the interface temperature I T through (3.10) as a function of the solid core temperature c T and the mean solid temperature T .It is worth to mention that two instances may occur with the determination of  :

  Estimated post shock solid core temperature:

,

  and the reduced model solutions, based on (4.2) and (D.4) expressed in dimensionless variables, are compared in Figure 4. The following dimensionless variables are: represents the dimensionless temperature,  the dimensionless radius,  the dimensionless time,  the dimensionless heat flux and ign T 550 K  an estimate of the explosive (PBX9501) ignition temperature. D represents the thermal diffusivity (

Figure 4 :

 4 Figure 4: Comparison of the heat equation and reduced model solutions for various values of the exponent n of the dimensionless interface temperature (3.7 and D.1) and associated heat flux. The interfacial heat flux deposited at the interface during pore collapse is taken constant to a quite realistic value in the present context ( 8 2 c q 6.9 10 W / m 

  If the configuration studied involves propagation of a single shock equations (4.4) are not necessary and are replaced by, volume of the solid phase).

Figure 5 :

 5 Figure 5: Schematic representation of the coupling procedure between the macro-scale model at left and subscale one at right.

Figure 6 :

 6 Figure 6: Experimental configuration of SDT studies for PBX9501 explosive.

  PBX9501 physical parameters used in the present study.

Figure 7 :

 7 Figure 7: Velocity signals during SDT of PBX9501. Computed results are shown with lines while Gustavsen et al.(1999) records are in symbols. Top: without IACF and bottom: with IACF.

  The use of the IACF function is a rough way to consider complicated effects of interfacial area creation. The parameters used with the IACF function (2.22) in the present computations are:

Figure 8 :

 8 Figure 8: Leading wave trajectory in the (x-t) diagram for PBX9501 SDT. Comparison between computed and experimental results of Gustavsen et al (1999).

Figure 9 :

 9 Figure 9: Mesh convergence of SDT computations for the PBX9501 test case.

Figure 10 :

 10 Figure 10: Computed temperatures used in the present model for the PBX9501 example, in the absence of mass transfer. Short term interface temperature uses the core temperature computed from the Hugoniot relation (4.3). Long term temperature corresponds to the situation the thermal boundary layer has reached the outer pore boundary.

Figure 11 :

 11 Figure 11: Reference interfacial temperature computed at various lagrangian gauges. Computations done with reference parameters of Table2.

Figure 12 :

 12 Figure 12: Computed interfacial temperature at various lagrangian gauges with smaller initial pore radius int R 0.3 m  . Change of this parameter renders the explosive insensitive in the present impact conditions.

Figure 13 :

 13 Figure 13: Computed interfacial temperature at various lagrangian gauges with normal initial pore radius int R 0.5 m  and increased solid heat conductivity 11 s 0.6 W.m K 

Figure 14 :

 14 Figure 14: Computed interfacial temperature at various lagrangian gauges with normal initial pore radius, normal initial heat conductivity 11 s 0.4538 W.m K 

Figure 15 :

 15 Figure 15: Computed interfacial temperature at various lagrangian gauges with reference parameters from Table2and lower impact velocity (320m/s). Change of this parameter renders the explosive insensitive in the present impact conditions.

Figure 16 :Figure 17 :

 1617 Figure 16: Computation of quasi-steady state detonation regime of PBX9501. Gas pressure profiles are plotted at various times showing the steady detonation regime.

Figure 18 :

 18 Figure 18: Initial configuration of 2D cylinder test computations. The first cylinder of inverse diameter 1 50m   is used to initiate and propagate a detonation wave in the cylinder under interest, of smaller size.

Figure 19 :

 19 Figure 19: Transmission of a detonation from a large cylinder to two smaller ones (top, diameter=1.667mm and bottom, diameter=1.538mm). Computed solid volume fraction profiles at times 2.5µs (left) and 4.3µs (right). At time 2.5 µs the detonation propagates in the large cylinder and is not yet transmitted to the small one. At time 4.3µs it is transmitted to small explosive cylinders and continues propagation on the cylinder of 1.667 mm diameter while it fails in the second one. The computed detonation velocity is recorded in each explosive cylinder and compared to the experimental relationship of Campbell and Engelke (1976):   3 21 det det

  Note that in this relation, the CJ speed is 8802 m/s, slightly lower than the value computed with the Cheetah thermochemical code (8857 m/s) and used in the JWL-EOS. This gives error bar estimate for the present computations. It also induces a deviation, as the JWL EOS parameters are used in the present computations. Velocity of detonation comparison is shown in Figure20, where both experimental data and computations (with and without IACF) are shown. Numerical simulations have been done with cylinders of inverse diameters 300, 350, 400, 450, 500, 550 and 600 1 m  . Cylinders of inverse diameters 60, 70, 80, 90 and 100 1 m  have been used for numerical computations without IACF.

Figure 20 :

 20 Figure 20: Velocity of detonation as a function of charge diameter. Experimental data of Campbell and Engelke are shown with lines. Computed results are shown with cross symbols and use the same IACF parameters as before: 10 max Neumann P P 6.10 Pa  and max f 450  in Equation (2.22). The lower curve shows computed results with the same model but without IACF.

  Hugoniot curves are shown in Figure A.1.

Figure A. 1 :

 1 Figure A.1: Comparison of theoretical Cochran-Chan and experimental (Marsh, 1980) Hugoniot curves for condensed PBX9501. The temperature dependence as a function of pressure for this explosive has been taken equal to the one of HMX (right graph of Figure A.1) as PBX9501 is made of 95% HMX.

  be recovered by taking the differences of the enthalpies given by the JWL and CC EOS at standard state: term solution of subscale heat transfer Long term solution has to be considered when the boundary layer exceeds the external pore radius ( solution, the unknowns are the solid/gas interface temperature I T and the solid temperature c T at the symmetry surface ext rR 

  is provided by the interface energy condition (3.9) that reads with (D.1), a direct extension of (3.10), where  has been replaced by ext R .The second relation is deduced from the mean solid temperature definition where the long term temperature profile (D.1) is inserted:

  , transported by (4.4), has to be used in (4.3). More precisely, it has to be used in (4.3) if the heat conduction boundary layer computed from has not exceeded the outer pore boundary (

			R	ext	).
	To summarize, if the computed point doesn't belong to the shock layer (4.6 false) and if 
	from (4.2) is greater than the outer pore boundary (		R	ext	) then,
	0n s,i s,i pp 	,
	0n s,i s,i vv 	,				(4.7)

meaning that the Hugoniot pole is reset to the current state. The new Hugoniot pole is transported by (4.4) and used for another shock computation and subscale heat transfer determination. It means that equations (2.5)-(2.6) and (

Table 1 :

 1 Vistal stiffened gas EOS parameters.

  is widely used to model the thermodynamics of detonation products. In Mie-Gruneisen type form it reads,

															PBX9501
															10 A 1.687 10 1 
															10 A 1.742 10 2 
															1 E 4.1 
															2 E 3.1 
															s 	0.7
															vs C 1444.2 J / kg / K 
															v	re f	43 5.476 10 m / kg  
															Re f T		298 K
	  P v,e	  e e v kg v   g     kg P v	,
	where,
	  e v kg		Re f 1 R Av	e	 R v v 1   0 Re f 2 Bv e R 	 R v v 2	Re f			Re f g   Re f kv v v  		g	and   e v,T e v C T   kg vg 	,
	  P v Ae kg 	 R v v 1     Ref Pv k1 Be 	 R v v 2	Ref		Ref v  v   k 	1  g	and   P v,T P v   g vg kg CT v  	.
	The CJ state appears in the computation of some constants,
	v	CJ	Re f CJ D v     Re f CJ v P	2	,
	and ,
	  CJ P v k1  g vg CJ CJ CJ CT k P vv v Re f CJ             	1  g	,

  . Parameters for PBX9501 are given in the Table B.1: PBX9501 Table B.1: JWL-EOS parameters of PBX9501 detonation products.

			11 A 10.3 10 
			9 B 10.97 10 
				1 R 4.75 
			2 R 1.074 
			g 	0.374
		vg C 2508 J / kg / K 
	v	re f	43 5.476 10 m / kg  
			CJ P		34.7 GPa
			CJ D		1 8835 m.s 
			CJ T		4043 K

  approximation the PBX9501 decomposition reaction reads,

	4 0.9394C H N O 0.0376C H N O 8 8 8 10 14.6 0.37 	3.42		8 14 0.0112C H N O 4 10		7 12 0.0117C H N O 4 10
		2 1.846056CO 4.18068H O 3.810356N 2.459044C 2 2   
	The associated heat of reaction	h 	0 react / mol	reads,
	h 	0 react /mol		1.846056h	0 CO 22 0 H O 4.18068h 		h	0 PBX9501	1764KJ / mol  	.

  long term flux reads:

									2
								ext	int
						ext	int
					22
					ext	ext	ext	int	ext	int
	q	s	 	s	int   Ic r s rR ext (T T ) int T u n r (R R )       	u	r	.

Acknowledgements. Part of this work has been funded by CEA Gramat and DGA that are gratefully acknowledged. The authors are also very grateful to the anonymous referees who helped to improve significantly the quality of the manuscript.