Rafik Absi
email: r.absi@ebi-edu.com

Charlène Le Gall
email: clegall@elcimai.com

Stéphane Marchandon
email: s.marchandon@ebi-edu.com

Pierre-Antoine Hammad
email: pa.hammad@ebi-edu.com

Interprétation des profils de concentration pour des particules de sable en suspension de différentes tailles

Keywords:

La forme des profils verticaux de concentration des sédiments en suspension sur un fond de rides, qui peut être concave ou convexe, dépend de la taille des sédiments. Pour améliorer la prédiction de ces profils, nous considérons d'une part, le modèle de longueur de mélange finie avec différentes descriptions du mélange turbulent et d'autre part, l'équation d'advection-diffusion classique basée sur l'approximation de Fick avec une description appropriée du paramètre  (rapport entre diffusivité et viscosité turbulente). L'analyse détaillée des données de mesures a permis d'identifier deux couches distinctes, une première proche du fond où le profil de concentration est concave et une seconde plus éloignée où le profil est convexe. La limite de passage d'une couche à l'autre dépend de la taille des sédiments.

1

Introduction Le transport sédimentaire côtier, résultant de l'hydrodynamique (houles, marées, courants), est un phénomène critique qui provoque l'évolution de la ligne de côte et des fonds marins. La compréhension de ce phénomène a un intérêt environnemental et économique fort.

Les modèles actuellement utilisés pour le calcul de la concentration des sédiments en suspension sur un fond de rides présentent quelques lacunes. En effet, ils ne permettent pas de prédire les profils de concentration concave et convexe qu'on rencontre pour des particules de sables de différentes tailles. Ces profils ont été observés in situ et en laboratoire [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF] sur un fond constitué de rides. Cette faiblesse au niveau de ces modèles est étroitement liée au problème fondamental de l'interaction entre les particules de sable et l'écoulement turbulent.

L'objectif de notre communication est de contribuer à améliorer nos capacités de calcul des profils de concentration des sédiments en suspension. Pour cela, nous utiliserons une modélisation analytique du mélange turbulent avec d'une part l'équation d'advection-diffusion classique basée sur l'approximation de Fick avec une nouvelle fonction pour le paramètre  (rapport entre diffusivité et viscosité turbulente) et d'autre part le modèle de longueur de mélange finie [START_REF] Nielsen | Turbulent diffusion of momentum and suspended particles: A finite-mixing-length-theory[END_REF]).

2

Modélisation des profils de concentration En écoulement oscillant (cas d'une houle), les sédiments en suspension sont soumis à un flux de sédimentation, descendant, et à un flux net ascendant dû au mélange turbulent. Le flux s q  peut s'écrire en fonction de la concentration :

s ) ( w z c q s    
; où s w est la vitesse de chute des sédiments et z est la distance par rapport au fond. Pour établir le profil moyen de concentration c(z) des sédiments à l'équilibre, il suffit d'écrire que le flux de sédimentation s q  compense exactement le flux de mélange m q  . L'équilibre des flux impose la relation 0

     m s q q , d'où :   0   z c w q s m (1)

2.1

Modèle classique basé sur l'approximation de Fick L'écriture du flux de mélange m q  peut être simplement obtenue en écrivant une proportionnalité entre m q  et le gradient de concentration, ou loi de Fick, de même Actes de Xèmes Journées Nationales Génie côtier Génie civil, Sophia-Antipolis, 14-16 octobre 2008, ISBN 2-9505787-9-9, Tome 1, pp. 9-18.

forme que celle obtenue dans le cas de la diffusion moléculaire :

dz dc D q t m   (2)
Où D t est la diffusivité turbulente. Soit, à partir des écritures (1) et [START_REF] Absi | On the effect of sand grain size on turbulent mixing[END_REF], nous obtenons l'équation d'advection-diffusion classique : (3)

2.2

Modèle de longueur de mélange finie Dans le modèle de longueur de mélange finie proposé par NIELSEN & TEAKLE (2004), le flux de mélange ascendant s'écrit [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF] Nous pouvons remarquer que nous obtenons l'équation (3) à partir de l'équation [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF], en négligeant le premier terme de [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF]

      2 / 2 / m m m m l z c l z c w q      (
. 2 exp( 1 ) . 2 exp( ² 1 2 2 / 1 1 2     C C C D U k f        (9) où, h z   ,
  ) exp( ) exp( 1 w l m f m L z z w C U w          (10) où :   f l m U z w   ; l z un niveau à partir duquel la vitesse de mélange (et donc k) commence à décroitre (ABSI, 2008a), et 1 / C h L w 
. NIELSEN & TEAKLE (2004) ont pris la vitesse de mélange donnée par l'éq. (10) et une longueur de mélange qui croît linéairement

z z l m    ) ( 1 . Leur viscosité turbulente s'écrit alors   ) exp( w l m m m t L z z w z w l         (
             w L z m e z z l   ) ( 2 (12)
Les équations (10) et (12) permettent d'écrire la viscosité turbulente comme suit :

                          h z C f L z f t e z U e z U w 1 2 2 2 2      (13) 4
Calcul des profils de concentration pour les sédiments en suspension NIELSEN & TEAKLE (2004) ont comparé les solutions obtenues à partir des équations ( 7) et (3) pour le cas d'un écoulement oscillant (houle) sur un fond de rides pour deux types de particules de sable, des petites particules s w =0.65cm/s et des plus grosses s w =6.1cm/s. La figure (1) présente les profils de concentration moyenne, les données de mesures sont celles de [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF]. Les traits pleins représentent les solutions de l'équation [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF] avec l m = z, =1, z l =0.005m, w m (z l )=0.025m/s, L w =0.022m, tandis que les traits en tiretspointillés représentent les solutions de l'équation [START_REF] Absi | Analytical solutions for the modeled k-equation[END_REF]. La figure 3 représente le profil de concentration obtenu après résolution de l'équation [START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF] avec le profil concave 4 m l ; cette figure montre que le profil de concentration c(z) est convexe pour les particules fines. Il apparaît donc essentiel de distinguer les particules fines des grosses, nous avons considéré donc deux descriptions différentes du mélange turbulent.

En effet, le profil de longueur de mélange convexe est capable de prédire un profil de concentration concave et inversement le profil de longueur de mélange concave est capable de prédire un profil de concentration convexe. Dans cette partie, nous avons supposé que la taille des sédiments affecte le mélange turbulent. Le mélange des particules qui passe par le produit m m l w . semble différent (le produit m m l w . pour les petites est différent du produit pour les grosses).

4.2

Profils obtenus avec le modèle classique Il est possible aussi d'avoir une bonne prédiction des profils de concentration aussi bien pour les petites particules que pour les grosses à partir de l'équation (3), avec la même description pour

  h z C C f h z C f C t t e z U e z U e D                    1 1 2 2 0
f Dt U     2 0  et    C C h h Dt   1 2 .
Les résultats de l'équation (3) avec

Interprétation et Discussions

La concavité/convexité du profil de concentration est liée à Comme la limite entre les deux couches dépend de la taille des particules de sable, nous avons suggéré dans un premier temps que le mélange turbulent dépend aussi de la taille des particules de sable. Cependant, ceci ne semble pas être la seule explication. En effet, proche du fond, les profils sont liés aussi au dépôt, le phénomène de sédimentation prédomine. Le profil concave serait donc lié au mélange turbulent, alors que le profil convexe serait plus lié au phénomène de sédimentation qui est plus important près du fond. C'est pourquoi la taille des particules impacte fortement les paramètres des profils et notamment la hauteur qui définit le seuil entre les deux types de profils : concave et convexe.

Conclusion

Pour améliorer la prédiction des profils de concentration des sédiments en suspension sur un fond de rides pour des sédiments de différentes tailles, nous avons considéré d'une part le modèle de longueur de mélange finie et de l'autre l'équation d'advection-diffusion classique basée sur l'approximation de Fick.

Le modèle de longueur de mélange finie a nécessité deux descriptions différentes pour le mélange turbulent, tandis que l'équation d'advection-diffusion classique a nécessité une fonction adéquate pour le paramètre  . La forme des profils de concentrations a été interprétée par une relation mathématique entre une première proche du fond où le profil de concentration est concave et une seconde plus éloignée où le profil est convexe. La hauteur par rapport au fond de la limite de passage d'une couche à l'autre dépend de la taille des particules de sable.

6

Références bibliographiques 1 ABSI, R. (2005). Comment on Turbulent diffusion of momentum and

  de mélange et m l la longueur de mélange. Le développement en série de Taylor donne : (2004) ont écrit l'équation (6), en se limitant aux deux premiers termes du développement, sous la forme d'une équation différentielle ordinaire du troisième ordre :

  et en décrivant le mélange turbulent non par le produit La résolution des équations différentielles (3) et[START_REF] Mcfetridge | Sediment suspension by non-breaking waves[END_REF] nécessite la connaissance du mélange turbulent, qui passe soit par la diffusivité turbulente t D (pour l'éq. 3'éq. 7). Ce produit représente la viscosité turbulente t  , qui est reliée à la diffusivité turbulente par le paramètre  : modèle k - ). Le profil de vitesse de mélange la solution analytique générale de l'équation de l'énergie cinétique de la turbulence obtenue Actes de Xèmes Journées Nationales Génie côtier Génie civil, Sophia-Antipolis, 14-16 octobre 2008, ISBN 2-9505787-9-9, Tome 1, pp. 9-18.dans ABSI (2008a). La solution particulière, qui correspond à l'équilibre entre production et dissipation de l'énergie cinétique turbulente, est similaire d'un coté à la fonction semi-empirique de NEZU & NAKAGAWA (1993) et de l'autre à la relation empirique utilisée par NIELSEN & TEAKLE (2004) qui s'appuie sur les mesures de DU TOIT & SLEATH (1981) sur un fond de rides.

)

  11) Le choix de NIELSEN & TEAKLE pour m w et m l donne un profil de vitesse non logarithmique. Si on suppose un profil de vitesse logarithmique, la longueur de mélange s'écrit selon ABSI (2008b)

  Actes de Xèmes Journées Nationales Génie côtier Génie civil, Sophia-Antipolis, 14-16 octobre 2008, ISBN 2-9505787-9-9, Tome 1, pp. 9-18.

Figure 1 .

 1 Figure 1. Profil de concentration moyenne sur un fond de rides pour des petites particules (o) w s =0.65cm/s et des grosses particules (×) w s =6.1cm/s (symboles mesures). Traits pleins: solution de (7) avec l m = z (=1) et z l =0.005m, w m (z l )=0.025m/s, L w =0.022m. tirets-pointillés, solution de l'éq. (3) (Figure d'après NIELSEN & TEAKLE, 2004).

  le modèle de longueur de mélange finie Comme précisé dans (ABSI, 2005), certaines équations de longueur de mélange permettent d'améliorer les prédictions du modèle de longueur de mélange finie. La solution (9) pour l'énergie cinétique de la turbulence permet de développer des équations pour la longueur de mélange. Nous avons proposé d'autres équations pour la longueur de mélange. Une première, noté 3 m l , qui s'écrit selon ABSI (et une autre obtenue à partir d'une extension de l'hypothèse de similitude de Von Karman qui s'écrit selon ABSI (2006) La figure 2 représente le profil de concentration obtenu après résolution de l'équation (7) avec le profil convexe 3 m l . Cette figure montre que le profil de concentration c(z) obtenu est concave pour les Actes de Xèmes Journées Nationales Génie côtier Génie civil, Sophia-Antipolis, 14-16 octobre 2008, ISBN 2-9505787-9-9, Tome 1, pp. 9-18. grosses particules. Le profil vertical de 4 m l est concave.

Figure 2 .

 2 Figure 2. Profil de concentration moyenne sur un fond de rides pour des particules grosses (w s =6.1cm/s). Courbe : solution de l'éq. (7) avec l'éq. (14) ; (=1) ; L w =0.05m , z l =0.005m; w m (z l )=0.025m/s. ×: mesures.

Figure 3 .

 3 Figure 3. Profil de concentration moyenne sur un fond de rides pour des particules fines (w s =0.65cm/s). courbe : solution de l'éq. (7) avec l'éq. (15) ; (=1) ; L w =0.022m , z l =0.005m ; w m (z l )=0.025m/s. o: mesures.

  petites et aux grosses particules. La viscosité turbulente est la même quelle que soit la taille des particules, cependant la diffusivité turbulente t D , qui est donnée par L'effet de la taille des particules sera pris en compte par le paramètre  .Nous avons proposé (ABSI, 2008b) pour le paramètre  une fonction empirique, à partir des expériences de GRAF et[START_REF] Graf | Suspension flows in open channels: experimental study[END_REF], dépendant de la distance z équation (13) pour t  , la diffusivité turbulente s'écrit donc :

tD

  donné par l'équation (18) sont présentés dans la figure4.

Figure 4 .

 4 Figure 4. Profil de concentration moyenne sur un fond de rides pour des petites particules (o) w s =0.65cm/s et des grosses (×) w s =6.1cm/s (mesures (McFETRIDGE & NIELSEN, 1985)). Courbes : solutions de (3) avec Dt  =0.025m/s. Trait plein : Dt h =0.24m. Tirets : Dt h =0.022m.

w

  la croissance/décroissance de la diffusivité turbulente est liée à dz dD t . Pour comprendre le lien entre concavité/convexité du profil de concentration et croissance/décroissance de la diffusivité turbulente, nous avons besoin dest toujours positif, l'équation (19) montre que Journées Nationales Génie côtier Génie civil, Sophia-Antipolis, 14-16 octobre 2008, ISBN 2-9505787-9-9, Tome 1, pp. 9-18. dz dD t ont le même signe. Par conséquent une diffusivité turbulente croissante entraîne un profil de concentration concave tandis qu'une diffusivité turbulente décroissante entraine un profil de concentration convexe. Notre analyse fine des mesures de McFETRIDGE & NIELSEN 1985 (voir figure 5) montre que pour les grosses particules il y a deux zones, une près du fond où le profil c(z) est convexe et une autre plus éloignée où c(z) est concave. Pour les particules fines, nous n'observons qu'un seul profil convexe. L' analyse d'autres données expérimentales sur un fond plat confirme qu'il y a deux couches, une première proche du fond où le profil c(z) est convexe et une seconde où c(z) est concave plus éloignée. Nous avons remarqué également que le seuil entre les deux types de profils, la limite de passage (à partir du fond) d'une couche à un autre, est plus importante pour les fines particules que pour les grosses (ABSI, 2006).

Figure 5 .

 5 Figure 5. Profil de concentration moyenne sur un fond de rides. , sable fin w s =0.65 cm/s ; , sable plus gros w s =6.1 cm/s.

  analyse fine a permis d'identifier deux couches distinctes,

  avec h pouvant représenter soit la hauteur d'eau soit l'épaisseur de la

	couche limite.	D 2 , 1 C et  sont des coefficients, et	U	f	est la vitesse de
	frottement. La vitesse de mélange s'écrit donc :


Rafik Absi
email: r.absi@ebi-edu.com

ANALYTICAL METHODS FOR VELOCITY DISTRIBUTION AND DIP-PHENOMENON IN NARROW OPEN-CHANNEL FLOWS

Keywords: Open-channel flows, free surface, velocity distribution, dip phenomenon, aspect ratio, channel slope

In this study, we present analytical tools for velocity distribution in steady uniform open-channel flows. Simple dip-modified laws are presented. An ordinary differential equation based on an analysis of the Reynolds-averaged Navier-Stokes (RANS) equations is proposed. Comparisons of predicted velocity profiles with experimental data show good agreement. It is well known that dip-phenomenon depends on aspect ratio Ar (ratio of the channel width b to the water depth h). Results

show that it should depend also on channel slope and friction velocity.

INTRODUCTION

In open-channel flows, the vertical velocity profile is well described by the classical log law in the inner region (ξ < 0.2, where ξ is the ratio of the distance from the bed y to the water depth h). However, in the outer region (ξ > 0.2) this law deviates from experimental data. This deviation is accounted for by adding the Coles' wake function. In 2D openchannel flows, the log-wake law appears to be the most reasonable extension of the log law [START_REF] Nezu | Turbulence in open-channel flows[END_REF]. However, in narrow open-channels (i.e. aspect ratio Ar smaller than 5, where Ar is the ratio of the channel width b to the water depth h) the maximum velocity appears below the free surface (i.e. the velocity dip phenomenon) involving deviation from the logwake law. This phenomenon is related to secondary currents which appear in 3D open-channel flows. Coles' wake function is unable to represent this behavior since it predicts a velocity which increases with distance from the bed.

In order to predict the velocity dip phenomenon, empirical relations were proposed. However, these empirical laws need a certain number of parameters such as dip position. Based on an analysis of the Reynolds-averaged Navier-Stokes (RANS) equations, Yang et al. [START_REF] Yang | Velocity distribution and dip-phenomenon in smooth uniform open channel flows[END_REF] proposed a dip-modified log law which consists of two logarithmic distances, one from the bed (i.e. the log law), and the other from the free surface. This law presents the advantage that it contains only one parameter of dip-correction α. The dip-modified log law reverts the classical log law for α = 0.

In our study, we will point out that even if the dip-modified log law [START_REF] Yang | Velocity distribution and dip-phenomenon in smooth uniform open channel flows[END_REF] predicts dip-phenomenon, it presents important difference with experimental data. In most cases, it is not possible to improve velocity profiles by adjusting the parameter α. The aim of this study is to improve the prediction of velocity distribution with dip phenomenon in open channel flows.

The shortcoming of the dip-modified log law will be analyzed and a simple dip-modified log-wake law, which allows improving predicted velocity profiles, will be presented. An ordinary differential equation, based on an analysis of the Reynolds-averaged Navier-Stokes (RANS) equations, will be proposed. Predicted profiles will be compared with experimental data of measured velocities in open-channels.

MODEL EQUATIONS

For steady uniform open-channel flows, using the continuity equation, the Reynolds-averaged Navier-Stokes (RANS) momentum equation becomes in the streamwise direction x (1) where x, y and z are respectively the streamwise, vertical and lateral directions and U, V and W the three corresponding mean velocities (u, v and w are the turbulent fluctuations), ν the fluid kinematic viscosity, g the gravitational acceleration, and θ is the angle of the channel bed to the horizontal axis [START_REF] Absi | Comments on `Turbulent velocity profile in fully-developed open channel flows[END_REF] [2].

DIP-MODIFIED LAWS

Dip-modified log law

From an analysis of the Reynolds-averaged Navier-Stokes (RANS) equations and by assuming a parabolic profile for the eddy viscosity, Yang et al. [START_REF] Yang | Velocity distribution and dip-phenomenon in smooth uniform open channel flows[END_REF] proposed the following law [START_REF] Absi | An ordinary differential equation for velocity distribution and dip phenomenon in open channel flows[END_REF] where u * is the friction velocity, h the water depth, κ the von Karman constant (≈ 0.41) and y 0 the distance at which the velocity is hypothetically equal to zero. This law is referred as the dip-modified log law. It presents the advantage that it contains only one parameter α and it reverts the classical log law for α = 0.

Figure [START_REF] Absi | Comments on `Turbulent velocity profile in fully-developed open channel flows[END_REF] presents comparisons of predicted velocity profiles obtained by dip-modified log law (Eq. 2) for different values of coefficient α with experimental data [START_REF] Lyn | Turbulence and Turbulent Transport in Sediment-Laden Open-Channel Flows[END_REF]. This figure shows that Eq. ( 2) with α given by empirical formula α = 1.3 exp (-Ar/2) proposed by Yang et al. (solid lines) is able to predict a small deviation from log law profile (dip phenomenon). We notice that, it is not possible to improve predicted velocity profiles by adjusting the parameter α; since by increasing α, the deviation from experimental data increases (dash-dotted lines, with α = 0:5) and by decreasing the value of α predicted profiles approach log law profiles (dashed lines).

A simple dip-modified log-wake law

Figure [START_REF] Absi | Comments on `Turbulent velocity profile in fully-developed open channel flows[END_REF] shows that log law profiles (dashed lines) are valid only in the inner region (ξ = y/h < 0.2). In the outer region (ξ > 0.2) log law deviates from experimental data. In 2D open-channel flows, this deviation is accounted for by adding the Coles wake function as [START_REF] Guo | Application of the modified log-wake law in open-channels[END_REF] where П is the Coles parameter expressing the strength of the wake function. It was found experimentally that П increases with the Reynolds number Re in zero-pressure-gradient boundary layers, attaining an asymptotic value of П = 0.55 at high Re. Examination of the log-wake law in the outer layer of 2D open-channel flows [START_REF] Nezu | Turbulence in open-channel flows[END_REF], shows that П increases from zero with the friction Reynolds number Re * and becomes nearly constant П ≈ 0.2, for Re * > 2000. In 2D open-channel flows, the log-wake law appears to be the most reasonable extension of the log law.

However, in 3D open-channel flows with secondary currents, the log-wake law is unable to predict velocity dipphenomenon (figure 2). A suitable simple law is possible by adding to log law both Coles' wake function and the term linearly proportional to the logarithmic distance from the free surface as Figure ( 2) presents comparisons of predicted velocity profiles obtained by our simple dip-modified log-wake law (Eq. 4) with experimental data of Lyn [START_REF] Lyn | Turbulence and Turbulent Transport in Sediment-Laden Open-Channel Flows[END_REF] for an aspect ration Ar equal to 4.68 but for different channel slope S = 0.00296 for C3 (figure 2.a) and S = 0.00401 for C4 (figure 2.b). Our profiles (solid lines) are also compared to log law profiles (dashed lines, α = 0 and П = 0) and log-wake law profiles (dash-dotted lines, α = 0 and П = 0.2). Comparison with experimental data shows reasonable agreement. Velocity profiles (solid lines) are obtained by dip-modified log-wake law with П = 0.2 but with different coefficients α. This indicates that coefficient α is not dependent only on aspect ration Ar, as proposed by Yang et al. [START_REF] Yang | Velocity distribution and dip-phenomenon in smooth uniform open channel flows[END_REF], but should depend also on channel slope S.

Modified log-wake law

Guo and Julien (2008) proposed a modified log-wake law as (5) where δ is the dip distance from the bed. Eq. ( 5) fits better experimental data. However, it needs a certain number of parameters such as δ the dip position. This seems to exclude the possibility of its adoption in predictive and practical applications.

AN ORDINARY DIFFERENTIAL EQUATION

From the analysis of the Reynolds-averaged Navier-Stokes (RANS) equations and by using a more appropriate approximation for the eddy viscosity in accordance to log-wake law, we obtain the following ordinary differential equation (ODE) [START_REF] Absi | An ordinary differential equation for velocity distribution and dip phenomenon in open channel flows[END_REF] (6) where , . Eq. ( 6) shows the interest that for α = 0 and П = 0 it gives the log law and for α = 0 it gives the log-wake law.

CONCLUSION

In this study, we proposed an ordinary differential equation based on an analysis of the Reynolds-averaged Navier-Stokes (RANS) equations. Comparisons of predicted velocity profiles with experimental data show good agreement. Our results show that coefficient α is not dependent only on aspect ration Ar, but should depend also on channel slope S and friction velocity.

Figure 1 .

 1 Figure 1. Comparison of dip-modified log law profiles (Eq. 2) [6] for different values of coefficient α with experimental data; dashed lines: α = 0 log law profile, solid lines: α = 1.3 exp (-Ar/2) proposed by Yang et al. [6], dash-dotted lines: α = 0.5; Symbols: measurements [4] for (a) C4, Ar = 4.68 and S = 0.00401, (b) C3, Ar = 4.68 and S = 0.00296.

Figure 2 .

 2 Figure 2. Comparison of our simple dip-modified log-wake law profiles (Eq. 4, solid lines) with experimental data; Symbols: measurements [4], dashed lines: log law profile (α = 0 and П = 0), dash-dotted lines: log-wake law profiles (α = 0 and П = 0.2), solid lines: simple dip-modified log-wake law (Eq. 4) with (a) α = 0.22 and П = 0.2 (C4 for Ar = 4.68 and S = 0.00401), (b) α = 0.27 and П = 0.2 (C3 for Ar = 4.68 and S = 0.00296).

Figure ( 3 )Figure 3 .

 33 Figure (3) presents comparisons of predicted velocity profiles obtained from the proposed ODE (Eq. 6, solid lines) with experimental data of Lyn (1986) [4] (fig. 3.a, 3.b) for an aspect ration Ar = 4.68 and Sarma et al. [7] (fig. 3.c, 3.d) for an aspect ration Ar = 4.98. In these figures, α and П are chosen to give the best fit. Predicted profiles show good agreement. This figure confirms that the parameter α is not only a function of aspect ratio but should include effect of channel slope and friction velocity.

ACKNOWLEDGEMENTS

The author would like to thank J. Guo for providing the experimental data files.


SNAPPING BICEPS FEMORIS TENDON: A DYNAMIC SONOGRAPHIC EVALUATION

Keywords: Knee, sonography, dynamic, snapping, biceps femoris Knee, biceps femoris, snapping, dynamic, sonography Introduction

Snapping of the distal arms of the biceps femoris tendon may explain pain and discomfort on the lateral aspect of the knee. Dynamic sonography has been recently developed in order to assess transient pathologic phenomenons. In our two cases, this modality was able to confirm the diagnosis and document which of the main arms was involved in the process. This modality represents a cheap, accessible and useful tool in the evaluation of the disease.

Case 1

A 44 year-old male, regular amateur cyclist, was reffered to our hospital to assess the association of a right knee pain and repeated snapping on the postero-lateral aspect of the knee, especially during sport activities. Symptoms had exacerbated progressively to the point that his activities were limited. There was no history of trauma. On physical examination, knee motion and stability remained in normal ranges. An audible snap and pain could be generated in the vicinity of the fibular head by exerting 70-90° flexion and extension manoeuvres. At static US examination [START_REF] Tubbs | Descriptive anatomy of the insertion of the biceps femoris muscle[END_REF](8)(9)(10)(11)(12)(13)(14)(15) MHz linear array transducer, ATL HDI 5000, Philips Medical Systems) during knee extension, axial examination showed, above the level of the fibular head, the bifurcation of the LHBFT. Echotexture of the two distal arms was normal and no bursitis was encountered. The patient was then placed in left lateral decubitus while the probe was firmly positioned in the axial plane along the fibular head (figure 1). As the patient was asked to flex the knee above 90 degrees, the anterior arm was shown to travel backward over the fibular head and distal lateral collateral ligament insertion. At subsequent extension, the tendon was demonstrated to spring over the same structure through a sudden move, thus producing postero-lateral pain and snap that were easily recognized by the patient (video and figures 1 and 2). This manoeuvre could be repeated with good reproductibility. Presence of two distal arms of the LHBFT was confirmed by MRI.

Case 2:

A 25 years old male, with irregular sports activities, was referred for a recurrent lateral knee pain along the bulge of the fibular head. For the last year, symptoms had become significant while the patient was following an intensive long distance running program, and tended to completely fade during periods of inactivity. If snapping of the knee was only intermittent during such efforts, it could be reproduced when the knee was flexed around 90 degrees at physical examination. At US (12 MHz linear array transducer, iU 22, Philips Medical System), both direct and tibial arms of the LHBFT were seen, without evidence of tendinosis or bursitis. During dynamic study, the body of the tibial arm was gliding along the lateral aspects of the distal lateral collateral ligament and underlying fibular head, thus producing a sudden snap. Presence of this anterior arm that was attached to the proximal tibia was also confirmed by MRI. As the patient accepted to lower the intensity of his activities, symptoms improved and surgery was not required.

Discussion

Lateral snapping of the knee may be regarded as a rare condition if one considers no wide serie is available in the literature. Case reports have involved numerous structures in lateral snaps, such as foreign bodies, fabellas (9), lateral meniscus (10), the proximal Other studies have emphasized the role of impingements between tendons and neighbouring bony structures, such as the popliteal tendon and the edge of the popliteal groove (14, 15), the iliotibial band and the lateral femoral condyle (16) and the distal arms of the LHBFT and fibular head (1,2,3,4,[START_REF] Bagchi | Partial fibular head resection for bilateral snapping biceps femoris tendon[END_REF][START_REF] Kristensen | Snapping knee from biceps femoris tendon[END_REF]17).

In all studies in the literature, involvement of the LHBFT in a lateral snapping process of the knee was only clinically suspected but was not confirmed before the time of open surgery. Moreover, despite the fact that both direct (3,[START_REF] Bagchi | Partial fibular head resection for bilateral snapping biceps femoris tendon[END_REF] and tibial arms of the LHBFT (1,2,4,[START_REF] Kristensen | Snapping knee from biceps femoris tendon[END_REF]17) may be separately responsible for the snap, the one that was brought to glide against the fibular head was never documented preoperatively.

In our cases, involvement of the tibial arm of the LHBFT and subsequent mechanism were clearly identified by US. Use of dynamic US to explain snapping phenomenons has been widely explored in the recent literature (18,19,20,21,22). At the level of the knee, this modality has only been cited in orthopedic studies in the field of the "snapping pes syndrome", in which the semi tendinosus and gracilis tendons were involved (23, 24), or in a case of snapping fabella in a prosthetic knee (9).

In previously cited surgical works, various anatomical abnormalities were proposed to explain how distal arms of the LHBFT were brought to snap against the fibula. Such hypothesis included an excessive prominence of the fibular head [START_REF] Bagchi | Partial fibular head resection for bilateral snapping biceps femoris tendon[END_REF]4), a two distal bifurcation of the distal heads of the LHBFT (2), a tear of an aponevrotic reflected arm of the LHBFT (17), a too anterior insertion on the fibula of the direct head of the LHBFT (3) and an exclusive attachment of the biceps femoris onto the proximal tibia (1,[START_REF] Kristensen | Snapping knee from biceps femoris tendon[END_REF]. Apart from the latter, none of those causes were properly evaluated at the time of examination for our two patients.

In conclusion, we report, for the first time, the use of sonography in the field of the snapping biceps femoris phenomenon at the knee. As shown in our two cases, US may, thanks to its dynamic capabilities, represent a unique, reliable, cheap and accessible tool to confirm the diagnosis and help in surgical planning.

Introduction:

Numerous structures may be responsible of a lateral snapping at the level of the knee.

Among them, distal arms of the long head of the biceps femoris (LHBFT) have been reported to impinge against the lateral aspect of the fibular head in the orthopedic literature (1,2,3,4,[START_REF] Bagchi | Partial fibular head resection for bilateral snapping biceps femoris tendon[END_REF][START_REF] Kristensen | Snapping knee from biceps femoris tendon[END_REF]. In addition, anatomical studies have emphasized a normal bifurcation of the distal LHBFT into two main arms, the direct one reaching the tip of the fibular head while an anterior one, also named "tibial" arm, ends up onto the anterolateral aspect of the proximal tibia through a lateromedial and slightly downward trajectory ( 1b). Here we report two cases in which dynamic ultrasonography (US) was useful to confirm the involvement of the tibial arm of the LHBFT in lateral snapping of the knee.

Case 1

A 44 year-old male, regular amateur cyclist, was reffered to our hospital to assess the association of a right knee pain and repeated snapping on the postero-lateral aspect of the knee, especially during sport activities. Symptoms had exacerbated progressively to the point that his activities were limited. There was no history of trauma. On physical examination, knee motion and stability remained in normal ranges. An audible snap and pain could be generated in the vicinity of the fibular head by exerting 70-90° flexion and extension manoeuvres. At static US examination (7-15 MHz linear array transducer, ATL HDI 5000, Philips Medical Systems) during knee extension, axial examination showed, above the level of the fibular head, the bifurcation of the LHBFT. Echotexture of the two distal arms was normal and no bursitis was encountered. The patient was then laid supine on the controlateral side placed in left lateral decubitus while the probe was firmly positioned in the axial plane along the fibular head (figure 1). As the patient was asked to flex the knee above 90 degrees, the anterior arm was shown to travel backward over the fibular head and distal lateral collateral ligament insertion. At subsequent extension, the tendon was demonstrated to spring over the same structure through a sudden move, thus producing postero-lateral pain and snap that were easily recognized by the patient (video Presence of two distal arms of the LHBFT was confirmed by MRI.

Case 2:

A 25 years old male, with irregular sports activities, was referred for a recurrent lateral knee pain along the bulge of the fibular head. For the last year, symptoms had become significant while the patient was following an intensive long distance running program, and tended to completely fade during periods of inactivity. If snapping of the knee was only intermittent during such efforts, it could be reproduced when the knee was flexed around 90 degrees at physical examination. At US (12 MHz linear array transducer, iU 22, Philips Medical System), both direct and tibial arms of the LHBFT were seen, without evidence of tendinosis or bursitis. During dynamic study, the body of the tibial arm was gliding along the lateral aspects of the distal lateral collateral ligament and underlying fibular head, thus producing a sudden snap. Presence of this anterior arm that was attached to the proximal tibia was also confirmed by MRI. As the patient accepted to lower the intensity of his activities, symptoms improved and surgery was not required.

Discussion

Lateral snapping of the knee may be regarded as a rare condition if one considers no wide serie is available in the literature. Case reports have involved numerous structures in lateral snaps, such as foreign bodies, fabellas (9), lateral meniscus (10), the proximal Other studies have emphasized the role of tendinoosseous impingements impingements between tendons and neighbouring bony structures, between such as the popliteal tendon and the edge of the popliteal groove (14, 15), between the iliotibial band and the lateral femoral condyle (16) and between the distal arms of the LHBFT and fibular head (1, 2, 3,

4, 5, 6, 17).

In all studies in the literature, involvement of the LHBFT in a lateral snapping process of the knee was only clinically suspected but was not confirmed before the time of open surgery. Moreover, despite the fact that both direct (3,[START_REF] Bagchi | Partial fibular head resection for bilateral snapping biceps femoris tendon[END_REF] and tibial arms of the LHBFT (1, 2, 4, 6, 17) may be separately responsible for the snap, the one that was brought to glide against the fibular head was never documented preoperatively.

In our cases, involvement of the tibial arm of the LHBFT and subsequent mechanism were clearly identified by US. Use of dynamic US to explain snapping phenomenons has been widely explored in the recent literature (18,19,20,21,22). At the level of the knee, this modality has only been cited in orthopedic studies in the field of the "snapping pes syndrome", in which the semi tendinosus and gracilis tendons were involved (23, 24), or in a case of snapping fabella in a prosthetic knee (9).

In previously cited surgical works, various anatomical abnormalities were proposed to explain how distal arms of the LHBFT were brought to snap against the fibula. Such hypothesis included an excessive prominence of the fibular head [START_REF] Bagchi | Partial fibular head resection for bilateral snapping biceps femoris tendon[END_REF]4), a two distal (3) and an exclusive attachment of the biceps femoris onto the proximal tibia (1,[START_REF] Kristensen | Snapping knee from biceps femoris tendon[END_REF]. Apart from the latter, none of those causes were properly evaluated at the time of examination for our two patients.

In conclusion, we report, for the first time, the use of sonography in the field of the snapping biceps femoris phenomenon at the knee. As shown in our two cases, US may, thanks to its dynamic capabilities, represent a unique, reliable, cheap and accessible tool to confirm the diagnosis and help in surgical planning.

  joint (11), ganglion cysts (12) or other type of intraarticular tumours (13).

5 -

 5 Bagchi K, Grelsamer RP. Partial fibular head resection for bilateral snapping biceps femoris tendon. Orthopedics 2003; 26(11):1147. 6-Kristensen G, Nielsen K, Blyme PJ. Snapping knee from biceps femoris tendon. Acta Orthop Scand 1989;60(5): 621. 7-Tubbs RS, Caycedo FJ, Oakes WJ, et al. Descriptive anatomy of the insertion of the biceps femoris muscle. Clin Anat 2006;19(6):517. 8-Terry GC, LaPrade RF. The biceps femoris muscle complex at the knee. Its anatomy and injury patterns associated with acute anterolateral-anteromedial rotatory instability. Am J Sports Med 1996;24(1):2. 9-Segal A, Miller TT, Krauss ES. Fabellar snapping as a cause of lateral knee pain after total knee replacement: assessment using dynamic sonography. Am J Roentgenol 2004; 183:352. 10-Youm T, Chen AL. Discoid lateral meniscus: evaluation and treatment. Am J Orthop 2004;33(5):234. 11-Cazeneuve JF, Bracq H. Chronic instability of the proximal tibio-fibular articulation: hemi-long biceps ligamentoplasty by the Weinert and Giachino propos of 3 cases. Rev Chir Orthop Reparatrice Appar Mot 1998;84(1): 84. 12_ Liu PC, Chen CH, Huang HT, et al. Snapping knee symptoms caused by an intraarticular ganglion cyst. Knee 2007;14(2):167. 13-Mine T, Ihara K, Taguchi T, et al. Snapping knee caused by intra-articular tumors. Arthroscopy 2003;19(3):E21. 14-Krause DA, Stuart MJ. Snapping popliteus tendon in a 21-year-old female. J Orthop Sports Phys Ther 2008;38(4):191. 16-Orchard JW, Fricker PA, Abud AT, et al. Biomechanics of iliotibial band friction syndrome in runners. Am J Sports Med 1996;24:375. 17-Bansal R, Taylor C, Pimpalnerkar AL. Snapping knee: an unusual biceps femoris tendon injury. The knee 2005;12:458. 18-Neustadter J, Raikin SM, Nazarian LN. Dynamic sonographic evaluation of peroneal tendon subluxation. Am J Roentgenol 2004;183(4):985. 20-Bureau NJ, Beauchamp M, Cardinal E, et al. Dynamic sonography evaluation of shoulder impingement syndrome. Am J Roentgenol 2006;187:216.

  V, Cardinal E, Bureau NJ. Musculoskeletal sonography: A dynamic tool for usual and unusual disorders. Am J Roentgenol 2007;188:W63. 21-Deslandes M, Guillin R, Cardinal; et al. The snapping iliopsoas tendon: new mechanism using dynamic sonography. Am J Roentgenol 2008;190:576. 22-Guillin R, Cardinal E, Bureau NJ. Sonographic anatomy and dynamic study of the normal iliopsoas musculotendinous junction. Eur Radiol 2009;19(4):995. 23-Bollen SR, Arvinte D. Snapping pes syndrome: a report of four cases. J Bone Joint Surg Br 2008;90(3):334. 24-Bae DK, Kwon OS. Snapping knee caused by the gracilis and semitendinosus tendon. A case report. Bull Hosp Jt Dis 1997;56(3):177.

  words or sentences in red Added words or sentences in blue Abstract: Snapping of the distal arms of the biceps femoris tendon may explain pain and discomfort on the lateral aspect of the knee. Dynamic sonography has been recently developed in order to assess transient pathologic phenomenons. In our two cases, this modality was able to confirm the diagnosis and document which of the main arms was involved in the process. This modality represents a cheap, accessible and useful tool in the evaluation of the disease.

  This manoeuvre could be repeated with good reproductibility.

  joint (11), ganglion cysts (12) or other type of intraarticular tumours (13).

  distal heads of the LHBFT (2), a tear of an aponevrotic reflected arm of the LHBFT (17), a too anterior insertion on the fibula of the direct head of the LHBFT

5 -

 5 Bagchi K, Grelsamer RP. Partial fibular head resection for bilateral snapping biceps femoris tendon. Orthopedics 2003; 26(11):1147. 6-Kristensen G, Nielsen K, Blyme PJ. Snapping knee from biceps femoris tendon. Acta Orthop Scand 1989;60(5): 621. 7-Tubbs RS, Caycedo FJ, Oakes WJ, et al. Descriptive anatomy of the insertion of the biceps femoris muscle. Clin Anat 2006;19(6):517. 8-Terry GC, LaPrade RF. The biceps femoris muscle complex at the knee. Its anatomy and injury patterns associated with acute anterolateral-anteromedial rotatory instability. Am J Sports Med 1996;24(1):2. 9-Segal A, Miller TT, Krauss ES. Fabellar snapping as a cause of lateral knee pain after total knee replacement: assessment using dynamic sonography. Am J Roentgenol 2004; 183:352. 10-Youm T, Chen AL. Discoid lateral meniscus: evaluation and treatment. Am J Orthop 2004;33(5):234. 11-Cazeneuve JF, Bracq H. Chronic instability of the proximal tibio-fibular articulation: hemi-long biceps ligamentoplasty by the Weinert and Giachino propos of 3 cases. Rev Chir Orthop Reparatrice Appar Mot 1998;84(1): 84. 12_ Liu PC, Chen CH, Huang HT, et al. Snapping knee symptoms caused by an intraarticular ganglion cyst. Knee 2007;14(2):167. 13-Mine T, Ihara K, Taguchi T, et al. Snapping knee caused by intra-articular tumors. Arthroscopy 2003;19(3):E21. 14-Krause DA, Stuart MJ. Snapping popliteus tendon in a 21-year-old female. J Orthop Sports Phys Ther 2008;38(4):191. 16-Orchard JW, Fricker PA, Abud AT, et al. Biomechanics of iliotibial band friction syndrome in runners. Am J Sports Med 1996;24:375. 17-Bansal R, Taylor C, Pimpalnerkar AL. Snapping knee: an unusual biceps femoris tendon injury. The knee 2005;12:458. 18-Neustadter J, Raikin SM, Nazarian LN. Dynamic sonographic evaluation of peroneal tendon subluxation. Am J Roentgenol 2004;183(4):985. 20-Bureau NJ, Beauchamp M, Cardinal E, et al. Dynamic sonography evaluation of shoulder impingement syndrome. Am J Roentgenol 2006;187:216.

  V, Cardinal E, Bureau NJ. Musculoskeletal sonography: A dynamic tool for usual and unusual disorders. Am J Roentgenol 2007;188:W63. 21-Deslandes M, Guillin R, Cardinal; et al. The snapping iliopsoas tendon: new mechanism using dynamic sonography. Am J Roentgenol 2008;190:576. 22-Guillin R, Cardinal E, Bureau NJ. Sonographic anatomy and dynamic study of the normal iliopsoas musculotendinous junction. Eur Radiol 2009;19(4):995. 23-Bollen SR, Arvinte D. Snapping pes syndrome: a report of four cases. J Bone Joint Surg Br 2008;90(3):334. 24-Bae DK, Kwon OS. Snapping knee caused by the gracilis and semitendinosus tendon. A case report. Bull Hosp Jt Dis 1997;56(3):177.

  Figure 1a:Axial sonographic view in the vicinity of the tip of the fibular head. Knee is in full extension. The tibial arm (TA, white arrows) of the LHBFT is encountered in front of the attachment of the lateral collateral (LCL) ligament on the fibular head (FH).

Figure 1b :

 1b Figure 1b : Lateral view of the knee in extended position. Lateral collateral ligament (LCL), iliotibial tract (ITT) and tibial (TA) and direct (DA) arms of the LHBFT are shown.

Figure 2a :

 2a Figure 2a:Axial sonographic view during flexion/extension of the knee between 70 and 90 degrees. The tibial arm (TA, white arrows) of the LHBFT is shown to travel over the fibular head, to glide in a sudden move along the distal attachment of the lateral collateral ligament (LCL) on the tip of the fibular head (FH) during extension.

Figure 2b :

 2b Figure 2b:During extension from flexion of the knee, the tibial arm of the LHBFT travels backward over the fibular head, thus producing an audible snap (orange star). Lateral collateral ligament (LCL), iliotibial tract (ITT) and tibial (TA) and direct (DA) arms of the LHBFT are shown.

  Figure 1aAxial sonographic view in the vicinity of the tip of the fibular head. Knee is in full extension. The tibial arm of the LHBFT (TA, white arrows) is encountered in front of the attachment of the lateral collateral ligament (LCL) on the fibular head (FH). Anterior (A) and posterior (P) sides are mentioned. 245x190mm (150 x 150 DPI)

John Wiley & SonsJournal of Clinical Ultrasound


C Gruau

D Picart

R Belmas

E Bouton

F Delmaire-Sizes

J Sabatier

H Trumel

Ignition of a confined high explosive under low velocity impact

Keywords: PBX safety, Steven test, ignition criterion ---

In the framework of low velocity impact studies, dedicated to safety analyses of plastic bonded explosives (PBX), we propose a new numerical tool, designed for restituting the ignition of a HMX (high melting point explosive) based composition. Major results are the use of a concrete-like constitutive law for the PBX and an efficient implementation of an ignition criterion. We also put forward two variants of classical Steven tests, which enable us to visualize either a dot ignition or an unusual ring-shaped ignition. It is shown that the calculation tool is able to restitute accurately both results.

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

Introduction

High explosive parts confined in pyrotechnic devices can unintentionally ignite, when submitted to impacts, even at low velocities. Ignition can then evolve to more violent and undesired events, such as deflagration or detonation. Prediction of such transitions being hardly tractable, we are focusing our attention on the ignition stage, which is the very first event that must be correctly restituted, in the framework of safety analyses.

The considered high explosive is composed of a poly-dispersed distribution of HMX grains, with a mean diameter of 200 µm, which are mixed with a few percentage of a polymeric binder. An isostatic compaction process is used to reduce the composition porosity at few percentages. Its mechanical response is similar to the behavior exhibited by PBX-9501, another explosive composition.

The low velocity impacts of interest may occur during transport, storage, assembling or disassembling of devices, leading to a great number of accident scenarios. Thus, safety analyses of pyrotechnic structures involve many parameters, like target configuration (shape, materials and boundary conditions) and mechanical loadings (characteristics of the projectile and its movement law), which prevents us from dealing only with experimental techniques. The aim of this paper is to describe a numerical tool designed for the ignition, applied to this kind of accidental configurations.

In order to validate such a numerical tool, an experimental database is needed. The so-called "Steven test", introduced by Chidester and co-workers [START_REF] Chidester | A frictional work predictive method for the initiation of solid high explosives from low-pressure impacts[END_REF] was modified to simplify the boundary conditions and to limit strain localization [START_REF] Delmaire-Sizes | Low-velocity impact tests on an HMX-based explosive[END_REF]. This test is devoted to the determination of a safety threshold (in terms of impact velocity), under which no violent reaction is observed. Figure 1 shows three pictures taken during one of these tests. Deflagrations and detonations induced by low velocity impacts are not the results of shock-to-detonation transitions, since pressures generated by these impacts are too low. Another mechanism should be invoked; however such a model is beyond the scope of this paper. We actually consider that the material ignites when a sufficient mechanical energy is locally dissipated in the vicinity of defects within the PBX microstructure. Investigations are being made to identify and model the hot spot mechanisms involved in high explosives during ignition. Among all local mechanisms proposed by Field et al. [START_REF] Field | Hot-spot ignition mechanisms for explosives and propellants[END_REF], an ignition by low velocity impact is probably caused by frictional heating of preexisting or induced closed micro cracks [START_REF] Dienes | Frictional hot-spots and propellant sensitivity[END_REF]. Such a frictional heating at the microstructure scale is assumed to result from a macroscopic plastic shear strain under a high level of macroscopic pressure. In other words, frictional heating leads to ignition when the levels and the durations of pressure and plastic shear rate are sufficient. In the literature, several ignition criteria have been proposed, relying on the macroscopic pressure p and on the macroscopic plastic shear strain rate pl γ & . Partom [START_REF] Partom | A threshold criterion for impact ignition[END_REF] used the instantaneous quantity pl p γ & , compared with a threshold value of 2.35 10 13 Pa.s -1 . Browning [6] and co-workers [7] employed the time-to-ignition ig t , in a first formula , derived from a constant thermal flux assumption. Here, we adopt a criterion based on the work of Browning and Scammon [START_REF] Browning | Microstructural model of ignition for time varying loading conditions[END_REF], which integrates the history of pressure and plastic shear rate, for varying flux conditions:

1 d ) ( ) ( 1 0 3 2 * * * = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ∫ - τ τ γ τ τ pl t n n ig ig p p t t c & (1)
where * c is the threshold value, * t is a characteristic time, * p is a characteristic pressure, p is the positive part of the pressure and the exponent n comes from the HMX decomposition kinetics. When (1) is met, ignition is supposed to occur.

Criterion (1) relies on the computation of the macroscopic pressure and the macroscopic plastic shear rate, inside the high explosive, seen as a continuum. In the history of the Steven test, several material models and several simulation softwares have been used to compute stresses and strains, inside the target. At the beginning, 2D computations were run with a double Jones-Wilkins-Lee equation of state, in conjunction with a reaction model [START_REF] Chidester | Low amplitude impact testing and analysis of pristine and aged solid high explosives[END_REF].Other 2D computations [START_REF] Scammon | Low amplitude insult project: structural analysis and prediction of low order reaction[END_REF] were also used with a statistical crack mechanics model [START_REF] Dienes | A unified theory of flow, hot spots, and fragmentation with application to explosive sensitivity[END_REF]. Then, 3D computations were performed, whilst taking a perfect plasticity model for the material [START_REF] Vandersall | Experimental and modeling studies of crush, puncture, and perforation scenarios in the Steven impact test[END_REF] and [START_REF] Vandersall | Investigation of Steven impact test using a transportation hook projectile with gauged experiments and 3D modeling[END_REF].

In fact, a low velocity impact requires a more complex modeling of the PBX behavior, than more violent events (like shocks that induce a pressure dominant state of stress). In the framework of low velocity impacts, the hydrodynamics assumption cannot be made and a tensorial approach of the relationship between stress and strain is adopted, like in impact studies on concrete [START_REF] Bazant | Large-strain generalization of microplane model for concrete and application[END_REF], rock [START_REF] Grange | Modelling dynamic fragmentation of rocks[END_REF] or ceramics [START_REF] Espinosa | A numerical investigation of penetration in multilayered material/structure systems[END_REF] structures. Furthermore, the material is very complex itself, involving many difficult issues like: pressure dependence,

T P I R C S U N A M D E T P E C C A
coupled with rate dependence, plastic compressibility and anisotropy. Unfortunately, all these ingredients influence the mechanical response of the material under a low velocity impact. Here, we use a concrete-like constitutive law, which involves a hypoelasticity model, a pressure dependent and non-associated plasticity with isotropic hardening. Dilation, dissymmetry between tensile and compression states of stress, altogether with isotropic damage, are taken into account. Some features of the material behavior have not been accounted for yet, like anisotropy due to oriented microcracking or plastic deformation of grains, rate dependent plasticity, a cap yield surface describing compaction at high pressure and a transition between brittle and ductile behavior (strain softening has been taken into account, but in an artificial manner, by a negative hardening evolution).

Abaqus/Explicit was chosen to run computations. This choice is motivated by several reasons. Firstly, we have to deal with continua (thanks to the confinement of the PBX, there is neither pulverization nor fragmentation), whose behaviors exhibit strong nonlinearities under finite strains (but no excessive distortion). So, a Lagrangian finite element method is convenient for our purpose. Secondly, even for low velocity impacts, the strain rate can overcome 10 5 s -1 (fast transient dynamics) and the contact between projectile and target can be strongly non linear. So, we focus our attention on an explicit time integration. Thirdly, preliminary tests on low velocity impacts show that Abaqus Explicit exhibits a better energy conservation, a stronger stability and broader functionalities than hydrocodes [START_REF] Gruau | Simulation of PBX when submitted to low velocity mechanical impacts[END_REF].

When associated to the concrete-like constitutive model and the ignition criterion (1), an Abaqus computation of the target 3 test case (see Figure 15 for a description) indicates that ignition happens at the PBX sample center (Figure 2). The remainder of this paper is devoted to the description of this numerical tool design and is organized in four sections. In section 2, material models are described and the parameter determination is validated against experimental measurements. Building, scaling and implementing of the ignition criterion are detailed in section 3 and appendix 1. Section 4 is devoted to numerical and experimental results obtained for different impact test variants. Finally, in section 5 a summary of the results is drawn and further investigations in progress are emphasized, while a second appendix is dedicated to the HMX decomposition model.

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
International Journal of Impact Engineering 4

Constitutive modeling

In order to take into account many (but not all) influential ingredients on the response of low velocity impacts, the behavior of our PBX is modeled by an elastic-plastic-damage law. This law based on the work of Lubliner and coworkers [START_REF] Lubliner | A plastic-damage model for concrete[END_REF], further improved by Lee and Fenves [START_REF] Lee | Plastic-damage model for cyclic loading of concrete structures[END_REF], about concrete materials. This model, called "concrete damaged plasticity" is natively implemented in Abaqus. Here we propose a brief review of the simplified version in use for the PBX. A description of other material models involved in the target and the confinement is also given. Then, confrontations are made between numerical results and experimental measurements.

PBX modeling

In this model, an isotropic damage is taken into account via a scalar parameter d :

( ) eff d σ σ - = 1 ( 2
)
where σ is the Cauchy stress tensor and eff σ is the effective (undamaged) stress tensor. In our case, the evolution of d depends on the equivalent plastic deformation, but only during tensile loadings:

∫ = = t pl t pl t pl t r d r d 0 max , , , with ) ( ε ε ε & (3)
where max , pl ε& is the maximal eigenvalue of the plastic strain rate tensor pl D . The curve ( )

t pl t d , ε
is given by an experimental characterization.

In expressions (3), the scalar variable r takes a continuous value between 0 and 1, quantifying the proportion of strictly positive eigenvalues of σ , via Then, the elastic strain is assumed to be small enough for the additive decomposition to be adopted. Elastoplasticity is written in terms of the hypoelastic Jaumann stress rate: ) ( :

pl eff K D D σ - = ∇ (4)
where K is the Hookean (isotropic) elastic stiffness tensor and D is the total strain rate tensor.

The yield function f is an extension of the linear Drucker-Prager criterion, using a friction angle ϕ , as follows:

c Y t Y c Y eff eff eff eff p q f , , , max , max , 3 tan 1 3 tan 1 3 tan 1 tan σ ϕ σ σ σ ϕ σ ϕ ϕ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - - = (5)
where

) ( 3
1 eff eff Trace p σ - =
is the effective pressure and ) ( :

) ( 2 3 I σ I σ eff eff eff eff eff p p q + + =
is the effective von Mises stress ( I being the identity matrix). The friction angle is determined from several triaxial compression tests (quasi-static), performed at different confining pressures. A value of 20° agrees well with pressures between 40 MPa and 230 MPa (Figure 3). At zero pressure, an initial yield stress of 5 MPa is considered. Strength loss is reached at a tensile pressure of about -13 MPa.

A major difference between (5) and the linear Drucker-Prager criterion is the influence of the third invariant, via the maximal eigenvalue denoted by σ , which enables us to deal with the great dissymmetry between compression behavior (rather ductile) and tensile behavior (brittle). Strain softening in tension is given by ( )

T P I R C S U N A M D E T P E C C A
t pl t Y ,
, ε σ , while strain hardening in compression is given by:

∫ - = t pl c pl c pl c Y r 0 min , , , , ) 1 ( with )
( ε ε ε σ & (6)
where min , pl ε& is the minimal eigenvalue of pl D . Both strain-stress curves are plotted on Figure 4.

Unfortunately, our model for compression cannot take the (experimentally observed) brittle-to-ductile transition into account. So, two assumptions are made. Firstly, strain hardening is identified at a 200 MPa confining pressure (which is the order of magnitude in an impact test). Secondly, a strain softening effect is introduced after a compressive plastic strain of 20 %, with a slope of -53 MPa. .

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

International Journal of Impact Engineering 6

A non-associated plastic flow is considered via the hyperbolic Drucker-Prager plastic potential, using a dilation angle ψ :

( ) ψ ψ σ tan tan ) 0 ( 2 , 2 eff t Y eff p e q g - + = ( 7
)
where e is the eccentricity and ) 0 ( ,t Y σ is the initial value of the yield stress in tension. The dilation angle is determined so as to reproduce the transversal response, observed on the triaxial compression tests. Figure 5 shows both longitudinal and transversal curves for different confining pressures (from 25 MPa to 800 MPa). A low value of 1° leads to a good agreement between the computed transversal curves and the experimental ones, at finite strain.

Materials parameters

The material parameters used for the PBX are summarized in Table 1. The dilation angle is lower than the friction angle, ensuring that plastic dissipation is positive. Other materials are involved in impact tests. The projectile is made of stainless steel and the confinement can be composed of: stainless steel, Plexiglas (PMMA), Teflon (PTFE) and transparent glass (Figures 15 and16). Glass and Plexiglas are not supposed to fail during the tests; an isotropic elastic model is adopted. Steel and Teflon are ductile, but they are also not supposed to fracture during the tests; an elastoplastic model is used, with a perfect plasticity for steel (a Johnson-Cook failed to give more accurate results) and with linear strain hardening up to 20 MPa for Teflon (Table 2). Friction between PBX and steel is modeled by a constant friction coefficient (without threshold). Other contact interactions (especially between PBX and glass or Plexiglas) are assumed to be frictionless.

Comparisons between experiments and numerical results

Confrontation between experimental and numerical results, obtained with the aforementioned models and parameters, is performed on impact tests, on instrumented targets. Like in [START_REF] Scammon | Low amplitude insult project: structural analysis and prediction of low order reaction[END_REF], [START_REF] Chidester | Low amplitude impact testing and analysis of pristine and aged solid high explosives[END_REF] and [START_REF] Vandersall | Investigation of Steven impact test using a transportation hook projectile with gauged experiments and 3D modeling[END_REF], instrumentation consists in strain and pressure gauges. We also use post-mortem observations made on soaked targets, especially to quantify the permanent strain and fracture.

In what concerns the test case of Figures 1 and2, the target contains a 100 mm diameter and 13 mm thick PBX sample, within a 150 mm diameter and 35 mm thick cylindrical steel confinement. A hemispherical steel projectile of 1.2 kg is launched at 77 m/s. The impacted target was equipped with a strain gauge placed in the center of the confinement rear face.

During the first 100 μs (Figure 6), numerical and experimental strains increase in the same way (except a slight delay for the numerical one). Then, the numerical curve shows a plateau around 0.25 %, which is nearly the value of the first experimental peak. Unlike the experimental result, the numerical one does not oscillate, except after 300 μs, but with an opposite phase.

T P I R C S U N A M D E T P E C C A
A central pressure gauge had also been placed between the front confinement and the explosive. Numerical results show a first plateau above 100 MPa, while the experimental gauge records two peaks under 100 MPa (Figure 7). Then, pressure rises in both cases, but the third experimental peak at 60 μs is not reproduced numerically. The largest calculated value is under 500 MPa and is attained before 270 μs, although the largest experimental value is about 550 MPa and is attained at 180 μs. The experimental pressure decrease after 200 μs is quite well restituted (though delayed), up to 300 μs where a violent reaction makes the experimental gauge signal increase suddenly. Another target, corresponding to a non reactive test case, was soaked, sliced and photographed (target 2, hemispherical projectile, according to Figure 15). Here, we use this information in order to measure the permanent strain and the fracture underwent by the PBX sample (Figure 8). Consequent profile and thickness sizes along the radius agree well between computations and measures.

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

International Journal of Impact Engineering

This confrontation between numerical and experimental results does not perfectly validate the current state of the material modeling and deserves further investigations on rheological characterization. However, agreement between the curves of Figures 6 and7 is not so bad until ignition (before 100 μs and under a pressure of 300 MPa), which does not prevent us from performing relevant safety studies.

Ignition criterion

Even in a nearly adiabatic framework (since the PBX conductivity is low and since the mechanical stimulus is shorter than 1 ms), the macroscopic temperature rise, due to plastic dissipation, does not exceed 40 K (except at contact between PBX and confinement, where a plastic strain localization occurs). Thus, a PBX thermal ignition due to a low velocity impact cannot be explained at the macroscopic scale, but rather at the microscopic one. As stated in the introduction, the heating mechanism privileged in our modeling is a local friction between HMX grains. The Browning & Scammon ignition criterion [START_REF] Browning | Microstructural model of ignition for time varying loading conditions[END_REF], described in this section (without innovation, except in what concerns its computation), is based on the dissipation generated by friction between HMX grains. Macroscopic pressure and plastic shear rate are transferred to a simple microstructure representation. Then, the criterion is used to compute whether the local heating, due to these quantities, is sufficient for the material to ignite. This criterion can take into account ignition through localization at adiabatic shear bands, like in [START_REF] Dey | Shear band formation in plastic bounded explosives (PBX)[END_REF], since frictional heating between grains is compatible with shear banding. In fact, the Browning & Scammon criterion can take into account the high level of plastic shear strain inside these bands, whose size is comparable with the coarse grain size.

Criterion building

Although a wide range of grain sizes and the presence of a polymeric binder should be taken into account for our material, the microstructure is assumed to be a 3D regular sphere packing of identical HMX grains (Figure 9). A sample is subjected to a compressive macroscopic pressure ) (t p and a macroscopic plastic shear rate

) (t pl γ &
. Only the plastic part of shear rate is considered, since micro-surfaces are not supposed to slide during elastic loadings and unloadings. Positive pressure ) (t p is transferred at the grain scale into a contact force ) (t F c , applying between two adjacent grains, which is proportional to ) (t p . Thus, there exists a constant area 1 a , such that:

) ( ) ( 1 t p a t F c = (8)
T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS
International Journal of Impact Engineering 10 Considering a Hertzian contact between two adjacent grains (Figure 10), the contact surface radius

) (t R c is proportional to ) ( 3 / 1 t F c
. So, there exists a constant 2 a (m.Pa -1/3 ), such that:

3 / 1 2 ) ( ) ( t p a t R c = (9)
Fig. 10. Hertzian contact and frictional dissipation.

Furthermore, the mean contact pressure

) (t p c
on the contact surface is proportional to ) (

) ( 2 t R t F c c . A constant 2 2 1 3 - = a
a a (Pa 2/3 ) can be found, such that:

3 / 1 3 ) ( ) ( t p a t p c = (10)
Then, the average heat flux ) (t ϕ , induced by frictional dissipation, is assumed to be proportional to

) ( ) ( t t p pl c γ & . So, a constant 4
a (m.Pa 2/3 ) can be defined such that:

) ( ) ( ) ( 3 / 1 4 t t p a t pl γ ϕ & = (11)
Now, an ignition criterion can be expressed as follows:

1 ) ( d ) ( ) ( 0 2 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ∫ ig t c c ig R R t K τ τ τ ϕ τ τ ( 12
)
where K is a kernel that depends on the HMX chemical decomposition kinetics.

To determine the kernel K , we consider a stationary case, namely

ϕ ϕ = ) (t and c c R t R = ) (
. In this case, ignition criterion (12) can be expressed as follows:

c c t c ig R R R t K ig ϕ τ τ 1 d 2 0 2 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ∫ (13)
In this context, Browning and Scammon [START_REF] Browning | Microstructural model of ignition for time varying loading conditions[END_REF] proved that a power law relationship can be established between [START_REF] Tarver | Chemical kinetic modeling of HMX and TATB laser ignition tests[END_REF] and [START_REF] Tarver | Thermal decomposition models for HMX-based plastic bonded explosives[END_REF] (see appendix 2). This power law introduces a new constant 5 a (J.m -3 .(m 2 .s -1 ) n ):

1 2 5 - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = n c ig c R t a R ϕ (14)
Here, we will take a value of 0.447 for n (appendix 2). Engineering 11 Combining ( 13) and ( 14), we obtain the following equation for K :

T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS International Journal of Impact
n c ig c t c ig R t a R R t K ig - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ∫ 1 2 5 2 0 2 1 dτ τ (15) which solution is n x a x K - = 6 ) (
, where Thus, ( 12) can be rewritten as:

6 0 1 2 1 d ) ( ) ( ) ( a R t ig t n c n ig = - ∫ - - τ τ ϕ τ τ (16)
and combining with ( 9) and ( 11), yields to:

7 0 3 / 2 d ) ( ) ( ) ( a p t ig t pl n n ig = - ∫ - τ τ γ τ τ & ( 17
)
with a constant threshold 7 a (in s -n .Pa 2n/3 ). In order to remove the units of quantities elevated to a non-integer power, we introduce a characteristic time * t and a characteristic pressure * p . Then, with a threshold value * c (dimensionless), ( 17) can be recast into the following dimensionless form:

1 d ) ( ) ( 1 0 3 2 * * * = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ∫ - τ τ γ τ τ pl t n n ig ig p p t t c &
An efficient and accurate computation relies on a Prony series approximation (appendix 1). The relative error of the numerical scheme ( 23)-( 24) is quantified in the following section, after a needful signal filtering. . With such a raw signal, integral [START_REF] Dey | Shear band formation in plastic bounded explosives (PBX)[END_REF] would inherit much numerical noise, which would prevent us from determining a discriminating threshold value for the ignition criterion.

Threshold value

Hence, we apply a low pass filter to ) ( ~τ f , so as to obtain a filtered signal ) (τ f , used in the criterion computation (previous section). A third order Bessel filter was chosen, because of a sufficiently regular gain and a predictable phase shift (which is nearly the cutting frequency inverse). This filter relies on the following scheme:

⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ - + + - = - + + - = - - + - - - - = )) ( ) ( ( 2 ) ( ) ( )) ( ) ( ( 2 ) ( ) ( )) ( ) ( ( 15 ) ( 15 ) ( ) 6 1 ( ) ( 3 2 t t f t f t t t f t f t t f t f t t t f t f t t f t f t t t f t t t f t t f c c c δ δ δ δ δ δ δ δ ω δ δ ω δ δ ω (18)
where t δ is the time increment and c ω is the cutting frequency (at Engineering 12 Fig. 11. Filtering of the signal at an integration point: the noise is highly reduced. Fig. 12. Consequence on the criterion evolution, at the same integration point: the noise is highly reduced, but an offset of nearly 5 μs is introduced. This filter works well in smoothing the signal field (Figure 13). These time and spatial noise reductions facilitate the determination of the ignition time and location. Once applied to the filtered signal, the approximate criterion agrees well with the exact one (Figure 14). Except at early times, where a high accuracy is not required, the relative error decreases rapidly beneath 5 %. Thus, the scheme ( 23)-( 24) shows a good behavior. Fig. 14. Evolution of relative error between exact and approximated criterion at the same integration point. Now, the quantity ) (t c can be computed at each integration point and if its value is greater than * c , ignition occurs.

0 = t we consider 0 = f , 0 = f & and 0 = f & T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS International Journal of Impact
T P I R C S U N A M D E T P E C C A
However, the threshold value * c remains to be calibrated. For the sake of simplicity, we consider a characteristic time * t of For a non reactive case, the maximum value of ) (t c (in space and time) is a lower bound for the threshold. Meanwhile, for a reactive case, the maximum value of ) (t c (in space and time) is an upper bound for the threshold (Table 3). Finally, a value of 4.868 10 4 for * c is compatible, according to the upper and lower bounds obtained with our four configurations. , which is greater than 1 when ignition occurs.

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

Confrontations between simulations and experiments

Previous impact test configurations were dedicated to a velocity threshold determination. However, their design prevents us from determining the ignition location and time of reactive cases. So, the computed ignition time of 98 μs for the test case 2, with a small ring-shaped (or a cap-shaped) ignition zone, cannot be validated experimentally, with the target 3 configuration. That is the reason why two new configurations with a transparent visualization window were designed (Figure 16), allowing fast video recording. The first configuration, with the round nosed projectile, leads to the observation of a punctual ignition (next section). The second configuration, with a nearly flat projectile and a puncher, was designed so as to obtain a ring-shaped ignition area (section 4.2). Experiments and simulations were performed with both configurations, so as to evaluate the capabilities of the numerical tool.

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

Impact test with rear visualization

The former configuration involves a hemispherical nosed projectile, launched at 103 m/s. The first lightening event is recorded between 50 μs and 60 μs (Figure 17). This experimental lightening event is due to local ignition, which is roughly located at PBX center and near the PBX rear face. The timestamp of 60 μs is considered as an upper bound of the computed ignition time, since the transition to lightening is not modeled. 60 μs is a classical time to ignition for low velocity impacts (contrary to shock ignition in 1 microsecond), even shorter than already observed time to ignition in references [START_REF] Chidester | A frictional work predictive method for the initiation of solid high explosives from low-pressure impacts[END_REF], [START_REF] Scammon | Low amplitude insult project: structural analysis and prediction of low order reaction[END_REF], [START_REF] Chidester | Low amplitude impact testing and analysis of pristine and aged solid high explosives[END_REF] and [START_REF] Vandersall | Experimental and modeling studies of crush, puncture, and perforation scenarios in the Steven impact test[END_REF].

The computation of this test case indicates an ignition at 53 μs (Figure 18), which is coherent with the 60 μs experimental upper bound (the filter delay of 5 μs being taken into account). Furthermore, the computed ignition zone is punctual and located at the center of the PBX rear face. In fact, no integration point lies exactly on the axis and, even in the vicinity of the axis, plastic shear strain is not negligible. Ignition occurs near the axis because the pressure level is very high at this place, with a strictly positive plastic shear strain. At the projectile radius, shear strain is higher, but the pressure level is not high enough to lead to ignition.

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

Impact test with a penetrating part

As for the classical tests, the previous configuration is a non penetrating impact, leading to a single dot ignition. Actually, safety analyses of pyrotechnic structures may involve scenarios with penetrating projectiles. We thus derived another configuration, with a hole pierced in the front confinement, in order to insert a cylindrical puncher (Figure 16 on the right). This puncher penetrates the PBX sample, when impacted by a flat nosed projectile, launched at 110 m/s. In this case, shear strain is enhanced and the pressure level remains high.

The numerical results show that ignition takes place before 55 μs (Figure 19), not only around the puncher corner (where the stimulus is quite severe, excluding an artifact), but also along a ring-shaped zone on the PBX rear face. The ring diameter is 12.6 mm. The corresponding experimental test was then performed, with a fast video recording (Figure 20). The very first lightening event appears between 53 μs and 60 μs. Again, the ignition time is correctly restituted. This experimental result validates the computed ring-shaped ignition on the PBX rear face. However, the experimental diameter is evaluated to be of nearly 25 mm, which is twice as the computed diameter. Furthermore, the fact that numerical simulation raises two separate ignition zones can be explained by the lack of a characterization about compressive strain softening in the material model. Without any experimental measurement, strain localization may not be well reproduced, so that a shear band connects both ignition zones.

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

Conclusion

In this paper, we propose two extensions of the classical Steven impact test, in order to validate a numerical tool for ignition under a low velocity impact. Major improvements are the ability to visualize the ignition location, the ignition time and the ability to generate shear induced ignition via a puncher. This numerical tool relies on an existing ignition model and its threshold value. This criterion integrates the macroscopic pressure and the macroscopic plastic shear strain rate, using an accurate and efficient approximation by Prony series. Computation of these macroscopic quantities involves a concrete-like constitutive model for the PBX material, with pressuredependent plasticity and dilation. Triaxial compressive tests, at different confining pressures, are used to determine the friction and the dilation angles.

In this framework, instrumented targets are used to compare the calculated and the measured strain and pressure evolutions. Meanwhile, post-mortem analyses make it possible to correlate the final target dents. Impact tests are also used to calibrate the ignition model threshold value and to compare computed ignition time and location to measurements. This confrontation is quite fair, even if remaining discrepancies are still observed and should deserve further work on the constitutive models.

The numerical tool is not a definitive one, great progresses remains to achieve. Work in progress is focusing on material characterization and modeling, as an extension of [START_REF] Gratton | Mechanical behaviour of a viscoelastic plastic granular material: experimental procedure and modelling[END_REF] and [START_REF] Frachon | Mechanical modelling and integration of a viscoelastic plastic material, sensitive to hydrostatic pressure -bend and brezilian tests calculations[END_REF], which involves difficulties to attain both high pressures and high strain rates, in finite strain. Furthermore, our PBX is also modeled at the microscopic scale, so as to improve the estimation of the energy dissipation generated by friction of closed microcracks, as initiated in [START_REF] Vivier | Studying and modelling a pressed HMX-based energetic material[END_REF] and [START_REF] Vivier | On the stored and dissipated energies in heterogeneous rate-independent systems[END_REF]. Developments about shear banding (localization limiter) and macroscopic fracture (discontinuity tracking) are also under investigation.

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

Appendix 1: Ignition criterion computing

In order to determine whether ignition occurs or not in the calculation, the following quantity has to be computed:

τ τ τ d ) ( ) ( ) ( 0 ∫ - - = t n f t t c (19)
where ) (τ f is a signal function to be integrated with the kernel

n t - -) ( τ .
Unfortunately, since n is not an integer, neither a simplifying Laplace transform, nor a Taylor development can lead to a simple incremental scheme. Furthermore, computing ) (t c at each time step of the simulation, would involve the storage of all the signal history and an expensive numerical integration, from the initial time. A discretization of kernel

n t - -) ( τ
could reduce the computational cost, but not the storage.

That is the reason why, following the work of [START_REF] Sabatier | Past and some recent CRONE group applications of fractional differentiation -Part 2[END_REF] and [START_REF] Sabatier | Fractional integration for energetic material ignition prediction[END_REF], kernel

n t - -) ( τ
is replaced by a Prony series, so as to reduce both computational cost and storage. This Prony series approximation consists in considering:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - + ≅ - ∑ ∫ ∫ = - N k t k k k t n f t A t f A C f t 1 0 0 0 d ) ( )) ( exp( ) ( d ) ( ) ( τ τ τ ω ω τ τ τ (20)
where C is a constant, k A and k ω are suitable amplitudes are frequencies. In practice, a few terms are required ( N = 10 is a usual value). Furthermore, in the right hand side of ( 20), the new kernels

( ) ( ) τ ω - - t k exp
can now be expanded in Taylor series, which removes the need of storing all the signal history. This Taylor development relies on the following time derivative identity:

τ τ τ ω ω τ τ τ ω d ) ( )) ( exp( ) ( d ) ( )) ( exp( 0 0 ∫ ∫ - - - = - - ∂ t k k t k t f t t f f t (21)
According to [START_REF] Tarver | Chemical kinetic modeling of HMX and TATB laser ignition tests[END_REF], with the following notation:

τ τ τ ω d ) ( )) ( exp( ) ( 0 ∫ - - = t k k f t C t c ( 22
)
it comes:

∑ = + ≅ N k k k k t c A t f A C t c 1 0 ) ( ) ( ) ( ω (23)
Now, using [START_REF] Tarver | Thermal decomposition models for HMX-based plastic bonded explosives[END_REF] with a time increment t δ , we can consider:

( ) Engineering 20 and denoting the Laplace transform by L , we use the following identity:

) ( 1 ) ( ) ( t t c t t f t C t c k k k δ δ ω δ - - + ≅ (24
τ τ τ τ d ) ( ) ( ) ( ) ( 0 ∫ +∞ - - - = f t H t t c n (25) T P I R C S U N A M D E T P E C C A ARTICLE IN PRESS International Journal of Impact
( ) ( ) ( ) ( ) ) , ) 1 ( ( , ), ( 1 1 1 τ ω τ ω τ τ τ τ - - Γ = - = - - - - - - - t n L t H L L t H t n n n (26)
where Γ is the classical Gamma function.

The Prony series is introduced by considering the following approximation:

∏ + + ≅ = - N k k k n C 1 ' 0 1 1 1 ω ω ω ω ω (27)
where 0

C is a constant to be determined and the N principal frequencies are computed recursively by:

l N n l h k N l h k ω ω ω ω ω ω ω ω 2 1 1 1 with ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = - ( 28
)
and the auxiliary frequencies are given by

k N n l h k ω ω ω ω - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = 1 ' .
After a partial fraction expansion, it comes:

∑ ∏ = = + + = + + N k k k N k k k A A 1 0 1 ' 1 1 1 ω ω ω ω ω ω (29)
where the N corresponding amplitudes are: So, we obtain the following approximation:

∑ ∑ = = - - - - + - ≅ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ - + + ≅ - - N k k k k N k k k n t A C t A C t A C A C L t H t 1 0 1 0 1 )) ( exp( ) ( , 1 ) ( ) ( τ ω ω τ δ τ ω ω τ τ (31)
where δ is the Dirac function.

( ) ∑ = - - = N k k k k n t A t C 1 exp ω ω (32) T P I R C S U N A M D E T P E C C A ARTICLE IN

Appendix 2: HMX ignition and exponent computation

The exponent n , appearing in the ignition criterion building (section 3.1), has to be determined by fitting a power law relationship between the ignition time ig t and a constant flux ϕ , imposed on a micro-surface of constant radius c R .

Following the work of Browning and Scammon [START_REF] Browning | Microstructural model of ignition for time varying loading conditions[END_REF], we establish this relationship by considering the global multistep chemical kinetic model for the thermal decomposition of HMX, proposed by Tarver [START_REF] Tarver | Chemical kinetic modeling of HMX and TATB laser ignition tests[END_REF] and [START_REF] Tarver | Thermal decomposition models for HMX-based plastic bonded explosives[END_REF]. The experimental data for the Tarver model in [START_REF] Lubliner | A plastic-damage model for concrete[END_REF] involves laser induced ignition of HMX, with ignition time from 80 milliseconds to 2 seconds. Obviously, there is a 10 3 factor between this lowest ignition time and our ignition timing. However, in the work of [START_REF] Henson | Ignition chemistry in HMX from thermal explosion to detonation[END_REF], measurements have been made about ignition times from 2 nanoseconds to 30 000 seconds. It is shown in [START_REF] Henson | Ignition chemistry in HMX from thermal explosion to detonation[END_REF] that the Tarver decomposition model interpolates well the ignition times between 2 nanoseconds and 300 micro-seconds, which is in agreement with our observed ignition timing.

This model relies on four reactions and five chemical species. The reaction sequence is:

β-HMX → δ-HMX δ-HMX → intermediate solid intermediate solid → intermediate gas 2 intermediate gas → final gas (A) (B) (C) (D) Each reaction source term A q & , B q & , C q & and D q & is based on an Arrhenius kinetic law: ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = T R E Z Q q I I I m I I exp α & (
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - = T R E Z T R E Z T R E Z T R E Z T R E Z T R E Z T R E Z T R E Z D D D E D D D C C C D C C C B B B C B B B A A A B A A A A exp exp exp exp exp exp exp exp 2 2 α α α α α α α α α α α α α (34)
Furthermore, the mass conservation equation for the chemical species is satisfied:

1 = + + + + E D C B A α α α α α (35)
Temperature T is computed by solving the heat flux equation with Abaqus/Standard (convection and diffusion):

D C B A p q q q q T T C & & & & & + + + = Δ -λ ρ (36)
where ρ is the HMX density, p C the HMX specific heat and λ the HMX conductivity. In fact, this conductivity is low (0.4 W/m/K) and it takes several milliseconds for the heat to really diffuse out.

In preliminary studies, such a thermal decomposition model has been used to compute, with a good accuracy, times-toexplosion for PBX confined in an ODTX apparatus (One Dimension Time to eXplosion [START_REF] Tarver | Thermal decomposition models for HMX-based plastic bonded explosives[END_REF]). For sufficiently high initial temperatures at the PBX boundary, these explosion times correspond to a very high rising of the temperature, which leads to a value for the source term

D C B A q q q q & & & & + + +
greater than 180 MW.m -3 . Therefore, in our computations, the ignition time ig t is the timestamp at which equality is reached between the source term and the threshold value of 180 MW.m -3 .

Contrary to [START_REF] Browning | Microstructural model of ignition for time varying loading conditions[END_REF], this thermal decomposition model is now applied at the grain scale, as illustrated in Figure 21. Axisymmetrical computations are performed with a HMX grain of 200 μm (diameter) and four different contact surface radii c R : 12.5, 25, 50 and 100 μm. Thanks to symmetry reasons, only half of the upper grain section is represented (Figure 21). Various constant heat fluxes ϕ are imposed on the contact surface: 10 7 , 10 8 and 10 9 W.m -2 . According to equation (11), with a constant 4 a of 77.7 m.Pa 2/3 for HMX, a flux of 10 7 W.m -2 corresponds, for instance, to a 0.1 MPa pressure and a 2 780 s -1 plastic shear rate, while 10 9 W.m -2 corresponds to a 500 MPa pressure and a 16 250 s -1 plastic shear rate.

According to [START_REF] Browning | Microstructural model of ignition for time varying loading conditions[END_REF], we study the evolution of fluence where 8 a is a dimensionless constant and exponent n is 0.447. This power law can be expressed as follows:

Fig. 1 .

 1 Fig. 1. High speed recording of an impact test (the PBX sample is confined inside the steel target): picture before impact, during impact and at the beginning of reaction (the configuration of the target is given on Figure 15, target 3, projectile velocity = 77 m/s).

Fig. 2 .

 2 Fig. 2. Computational result of the target 3 test case: the left hand side of equation (1) is denoted by "Browning" in the legend.

  are adopted: stands for positive part and for absolute value). Its use in the computation of d enables stiffness recovery in compression (for such stress states, r decreases and ( )

Fig. 3 .

 3 Fig. 3. Friction angle determination: ultimate yield stresses obtained with triaxial compression tests (in blue) and the considered initial yield function (in red) are plotted against the effective pressure.

Fig. 4 .

 4 Fig. 4. Isotropic strain hardening curve in tension and compression (unconfined).

Fig. 5 .

 5 Fig. 5. Longitudinal and transversal behaviors for different confining pressures. Comparisons between experimental results and numerical approximation.

Fig. 6 .

 6 Fig. 6. Comparison between strain evolutions.

Fig. 7 .

 7 Fig. 7. Comparison between pressure evolutions.

Fig. 8 .

 8 Fig. 8. Comparison between observed and calculated post-mortem sample (target 2, hemispherical projectile, following the Figure 15 numbering)

Fig. 9 .

 9 Fig. 9. Schematics of the microstructure and notations.

  thermal decomposition of HMX, based on the global multistep chemical kinetic model proposed by Tarver

Fig. 13 .

 13 Fig. 13. Smoothing of the signal field at an intermediate time step: spatial noise is reduced.

Fig. 15 .

 15 Fig. 15. Four axi-symmetrical impact test configurations (deformed geometries at 100 μs) are used for determining an ignition threshold.

Fig. 16 .

 16 Fig. 16. Two new axi-symmetrical impact test configurations (deformed geometries at 100 μs) are used to determine the ignition time and location.

Fig. 17 .

 17 Fig. 17. High speed recording of the rear face (target 4 test case).

Fig. 18 .

 18 Fig. 18. Computational result at 53 μs of the target 4 test case.

Fig. 19 .

 19 Fig. 19. Computational result at 55 μs of the new penetrating reactive test (target 5).

Fig. 20 .

 20 Fig. 20. High speed recording of the target 5 test case.

Finally

  33) where = I A, B, C or D, α is the mass fraction of the left hand side specie, I m the reaction order, I Q the heat of reaction per unit of volume, I Z the frequency factor, I E the activation energy, R the gas constant and T the temperature. The chemical kinetic parameters for the HMX (coarse grains) decomposition model are listed in Table4, according to[START_REF] Tarver | Thermal decomposition models for HMX-based plastic bonded explosives[END_REF].

Fig. 21 .

 21 Fig. 21. Different HMX grains with an imposed heat flux on a surface with a radius of, respectively, 12.5, 25, 50 and 100 μm (a mesh size of 250 nm is required on those surfaces, other surfaces being considered as adiabatic).

  ig t , with different radii c R . In practice, this evolution is investigated after an appropriate scaling:

  time, Q being a global reaction energy per unit of volume. So as to compute a power law interpolation, a log-log graph is plotted between the scaled fluence 2 1 π π and the scaled time 2 π (Figure22).

Fig. 22 .

 22 Fig. 22. The log of scaled fluence is quite well linear with the log of scaled time: a power law interpolation is relevant.

From

Table 1 .

 1 Material parameters for the PBX.

	parameters	value
	density	1800 kg.m -3
	Young modulus	4 GPa
	Poisson's ratio	0.4
	friction angle	20°
	dilation angle	1°
	eccentricity	0.1
	initial yield stress	5 MPa
	compressive strain softening slope	-53 MPa

Table 2 .

 2 Parameters for other materials.

	parameters steel	Plexiglas glass	Teflon

Table 3 .

 3 Upper and lower bounds for the threshold value.

	configurations	highest non reactive velocities	lowest reactive velocities	threshold lower bounds	threshold upper bounds
	target 1 projectile 1		77 m/s	84 m/s 4.281 10 4 4.925 10 4
	target 2 projectile 1		62 m/s	76 m/s 4.867 10 4 8.693 10 4
	target 2 projectile 2		76 m/s	81 m/s 4.302 10 4 4.871 10 4
	target 3 projectile 1		61 m/s	77 m/s 4.066 10 4 7.523 10 4
	Figures 2, 18 and 19 display the scalar field	1 * c	c	( t	)

Table 4 .

 4 parameters of HMX decomposition.

	reaction	reaction order I m	heat of reaction I Q (MJ.m -3 )	1 ln -s I Z 1	activation energy I E (kJ.mol -1 )
	A	1	-76	48.1	203
	B	1	-460	48.7	220
	C	1	1 000	37.8	185
	D	2	10 000	28.1

International Journal of Impact Engineering

& & & & & & & & & & & &

& ).Empirically, a cutting frequency of 200 kHz leads to a correct filtering of the signal (Figure11), with a reasonable phase shift of nearly 5 μs on the criterion evolution (Figure12).

s and a characteristic pressure * p of 1 Pa. Then, an experimental database, involving several reactive and non reactive tests, is used to determine a lower bound and a upper bound for * c . This database relies on four configurations, with different confinements and different projectiles (Figure15).

& & & &

Acknowledgments

Authors wish to thank gratefully M. Biessy, E. Espinet, C. Calleau and other contributors from the reactive experiment facility team, J. Belair, M. Leblanc and S. Veyre for their work on Abaqus, J.-L. Brigolle, P. Lambert, M. Labrunie, Y. Sadou and P. Rey for their work on PBX. Reviewer's helpful suggestions are also deeply acknowledged.