P Vannucci
email: paolo.vannucci@uvsq.fr

F Masi

I Stefanou

A comparative study on the effects of blast actions on a monumental structure

Keywords: blast actions, fast dynamics, architectural monuments

The problem of evaluating the effects of an explosion inside a monumental building is addressed; namely, we focus on the determination of the overpressure due to the shock wave. Three models are compared: a thermodynamic approach, based on the JWL model and two empirical approaches, based on the CONWEP method or the recommendations TM5-1300, both of them approaches proposed by USACE.

In particular, a detailed procedure for the use of TM5-1300 is proposed, making use of new analytical expressions of the physical parameters of a hemispherical blast fitting the experimental data of Kingery and Bulmash. The novelty of these expressions is that they fit accurately the experimental data over the entire variation range of the scaled distance.

In the paper, we show the existence of reflected waves focus, a phenomenon that alters sensibly the pressure field in case of non-typical geometries, like vaults or domes. We show that for this kind of structures the assessment of the blast effects is adequately estimated by the JWL method, while in the case of more typical, planar geometries, namely of flat walls, the assessment of the overpressure is sufficiently estimated also by TM5-1300.

Introduction

The most widely used methods for the evaluation of the actions and effects of a blast on a civil structure are based upon the use of empirical laws, i.e. on laws describing the main physical quantities related to an explosion, that have been obtained as the result of experimental measurements. Such an approach has the advantage of being relatively simple to be used in calculations by hand or in numerical procedures, at least when the geometry of the blast scenario (i.e. the geometry of the constructions affected by the blast) is not too complicated.

The use of empirical laws has been widely studied and coded in different recommendations/rules, mostly proposed by military authorities. In particular, the two most commonly used empirical models are based upon different but related studies of the U.S. Army Corps of Engineers (USACE): the document [USACE, 1986], containing the model CONWEP, and the Technical Manual TM5-1300, [USACE, 1990], completed by successive documents, [USACE, 2008]. The Joint Research Center of the European Union has produced in 2013 a Technical Report, [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF], substantially referring to these two last documents and to another Technical Report of the U.S. Army, [START_REF] Kingery | Air blast parameters from TNT spherical air burst and hemispherical burst[END_REF]. In [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF] all the empirical laws of [USACE, 1990] are reproduced using S. I. units.

In particular, TM5-1300 recommendations contain several indications for the calculation of the physical parameters of a blast in different scenarios, for both the cases of an explosion outside or inside a building. Unfortunately, such scenarios are related to relatively simple and recurrent geometries, so such indications/practical computation schemes are not suited in more uncommon situations.

The case of monumental structures is one of those possible unusual scenarios. In fact, the geometry of monuments such as old palaces, cathedrals, temples and so on is more complicated and articulated than simple box-shaped building, like those considered in technical recommendations. In such cases, the scenario cannot be reduced to one of the cases provided by technical rules like TM5-1300.

Another aspect can be of great importance in the case of a blast inside a monument. The peculiar geometry of certain structural elements like vaults, domes, absides can produce a physical phenomenon that is not included in the empirical rules: the reflection and focusing of shock waves. This effect is described in this paper on a simple example representing a possible, common monumental structure, i.e. a barrel-vaulted hall. We show that it can have a great importance in the analysis of the effects of a blast on a monumental structure. Reflected shock waves are not taken into account by empirical methods. However, the effects of reflected waves, especially in the case of internal blasts, can be very important and generally speaking, not negligible.

The use of empirical approaches to the study of the effects of a blast on a monumental structure, especially for the case of internal explosions, may therefore be inappropriate. Approaches based upon a complete, more physical, description of the blast are preferable. Such methods, the most well known being the so-called JWL model, are based upon the laws of thermodynamics to recover the physics of a chemical blast; they allow for a complete description of the phenomenon of the blast, including the propagation of reflected shock waves. The drawback of such approaches is that usually they need heavy numerical simulations.

The study of the effects of a blast on a monumental structure has acquired today great importance: architectural monuments of the past, as symbols of a cultural heritage proper to a nation, to a civilization or to a religion, have been unfortunately too often the objects The present international situation has brought to the attention of people and governments the threat on highly representative monuments, that can be today the target of destructive attacks. It is hence relevant to consider the mechanical problem of the effects of a blast on a monumental structure; rather surprisingly, literature seems poor about this matter. The most part of papers concerning the effects of an explosion on a civil structure regard reinforced concrete structures and usually the geometries considered are really simple, normally a squared building, [START_REF] Remennikov | A review of methods for predicting bomb effects on buildings[END_REF], [START_REF] Ngo | Blast loading and blast effects on structures -An overview[END_REF], [START_REF] Koccaz | Architectural and structural design for blast resistant structures[END_REF], [START_REF] Draganic | Blast loading on structures[END_REF], in some cases bridges, e.g. [START_REF] Birhane | Blast analysis of railway masonry bridges[END_REF]. As mentioned above, this is seldom the case of architectures of the past, often much more articulated than simple squared buildings.

In this paper, we ponder on the best way to study the effects of a blast inside a monumental structure. We compare, to this purpose, different approaches: one based upon the use of JWL so giving a detailed description of the phenomenon, including the reflected waves. Then, three other approaches are used: one using the CONWEP method, [USACE, 1986], another using the approach suggested by the Technical Manual TM5-1300, [USACE, 1990] and finally a last one combining these two last approaches.

The approaches based upon the JWL and CONWEP methods have been used as implemented in finite element codes. For what concerns the method based upon the TM5-1300 recommandations, we propose in this paper an approach, namely a sequence of calculation, which is new and that is based upon a new analytical representation of the empirical data describing the different parameters of a hemispherical blast. This analytical representation is accurate and allows to represent by a unique function each physical parameter (namely, incident and reflected pressure and impulse, arrival time and duration of the positive phase, see Sect. 4). This can lead to a more elegant and straightforward implementation of the aforementioned empirical models in software for structural analyses.

The paper is subdivided into six Sections: the first one, Sect. 2, is devoted to a general description of the mechanical effects of a blast (we do not consider in this paper thermal effects nor the projection of fragments). The JWL model, [START_REF] Jones | The detonation of solid explosives[END_REF], [START_REF] Wilkins | The equation of state of PBX 9404 and LX 04-01[END_REF], [START_REF] Lee | Adiabatic expansion of high explosives detonation products[END_REF], is shortly recalled in Sect. 3, while Sect. 4 is devoted to the empirical models for the calculation of the blast overpressure and to the new analytical representations for the different physical parameters. In Sect. 5 we propose a procedure for the evaluation of the blast overpressure, using precise interpolating functions of the experimental data contained in the report [USACE, 1990], while in Sect. 6 we give a brief account of the model CONWEP, [USACE, 1986]. Finally, a comparison of the three models is given in Sect. 7, in a case study representing schematically the interior of a possible monument. A last Sect. 8 giving the conclusions ends the paper.

The mechanical effects of a blast

An explosion is an extremely rapid and exothermic chemical reaction that lasts just few milliseconds. During detonation, hot gases, produced by this chemical reaction, expand quickly and, for the hot temperatures produced instantaneously, the air around the blast expands too. The result is a blast shock wave, characterized by a thin zone of air propagating spherically much faster than the sound speed, through which pressure is discontinuous.

The shock-wave, travelling along a solid surface, produces an almost instantaneous increase of the air pressure on the surface, that decreases very quickly to the ambient pressure; this is the so-called positive phase of the blast. Then, the pressure decreases further, below the ambient pressure and then increases again to the ambient pressure, but in a longer time; this is the negative phase of the blast, see Fig. 1. The shock wave is the main mechanical effect of a blast on a structure, but not the only one: the hot gases, expanding, produce the so-called dynamic pressure, least in value with respect to the shock wave and propagating at a lower speed, while the impinging shock wave can be reflected by solid surfaces and act again on other surfaces as reflected shock wave.

To better understand all the mechanics of a blast, let us first introduce some quantities, used in the following:

• W : explosive mass;

• R = ||q -o||: distance of a point q from the detonation point o;

• P o : ambient pressure;

• P s : overpressure due to the blast; it is the pressure in the air relative to P o ;

• P so : side-on overpressure peak: the shock-wave peak pressure, relative to P o , measured in the air at q;

• P r : reflected overpressure: the pressure, relative to P o , acting at a point q of a solid surface when hit orthogonally by a shock-wave;

• P ro : peak value of P r ;

• P rα : peak value of the reflected overpressure at a point q of a solid surface atilt of the angle α on the direction of q -o;

• t A : arrival time, i.e. the instant at which the shock-wave peak arrives at q, taking as t = 0 the instant of detonation; of course, t A increases with R, but experimental evidence has shown that it decreases with W , i.e. the velocity of the shock-wave increases with the quantity of explosive;

• t o : positive phase duration; this is the duration of the time interval, starting from t A , during which P s ≥ 0; t o increases with R;

• t o-: negative phase duration; this is the duration of the time interval, starting from t A + t o , during which P s < 0;

• Z: Hopkinson-Cranz scaled distance, [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF], defined as,

Z = R W 1 3 ; (1)
Figure 1: Scheme of the time variation of the pressure due to a blast.

• i s : impulse of the shock-wave positive overpressure, defined as

i s = t A +to t A P s (q, t) dt; (2)
• i r : impulse of the shock-wave positive reflected overpressure, defined as

i r = t A +to t A P r (q, t) dt; (3)
• t rf : fictitious duration of the positive phase of the blast, defined as

t rf = 2 i r P ro ; (4)
• t Aw , t ow , t ow-: scaled durations, obtained dividing the respective durations by W 1 3 ;

• i sw , i rw : scaled impulses, obtained dividing the respective impulses by W 1 3 .

The overpressure P s at a point q decreases with both the time t > t A and R. Generally, the time rate decrease is much greater than the space rate decrease: the blast overpressure is really like a very localized pressure wave that propagates at high speed and whose intensity decreases, like for any other wave, with the travelled distance.

Fig. 1 represents an ideal profile of the overpressure P s (q, t). When the shock wave arrives at q, after t A from detonation, the pressure instantaneously increases, from the ambient pressure P o to a peak P so , which is a strong discontinuity.

For t > t A the overpressure decreases extremely fast, with an exponential rate, until time t A + t o , the end of the positive phase. After t A + t o we have the negative phase: the pressure decreases with respect to P o and then it returns to P o after a time t o-> t o .

During the negative phase, the decrease of the pressure is much lower in absolute value, than the peak pressure of the positive phase. This is why the negative phase is usually neglected in structural analyses, even though it can be important in some specific cases, due to its duration, which is always much longer than the positive phase.

Generally, the arrival time t A decreases, for the same distance R, when the amount of explosive W increases, i.e. the velocity of the shock-wave increases. The peak value P so , on its side, increases with W and decreases with R, while the time duration t o does exactly the opposite, see Fig. 2.

The decrease of the pressure wave, i.e. the function P s (q, t), is an extremely rapid phenomenon; it can be modeled by the Friedlander's equation, [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF]:

P s (q, t) = P so (q) 1 - t -t A t o e -b t-t A to . (5)
With this model, the impulse can be calculated analytically:

i s (q, t) = P so (q) t o b 2 (b -1 + e -b ). (6)
This relation is useful for determining the decay coefficient b if i s is known, e.g. from experimental data or with the method proposed below; a discussion about the value of b is given in Sect. 7.5. The same rate decrease (5) is used also for the reflected overpressure, P r .

P r is the pressure that acts on a surface impinged by the incident overpressure P s . The peak P ro of P r is normally much greater than P so measured at the same point in the absence of any surface. In fact, if we consider air as an ideal linear-elastic fluid, the air particles should bounce back freely from the surface, this giving a P r equal to the double of the incident pressure. But normally, P ro /P so > 2 because a blast is actually a nonlinear shock phenomenon, where the reflection of the particles is hindered by subsequent air particles, with consequently a far higher reflected pressure.

A formula relating the values of P ro and P so for normal shocks is

P ro = 2P so 4P so + 7P o P so + 7P o , (7)
The above equation indicates, on one side, that the ratio P ro /P so is not constant, but depends upon P so and, on the other side, that this ratio can vary between 2 and 8 or more. Of course, it is the value of P ro to be used for structural design. The Friedlander's law is used also for describing the decrease of P r .

Figure 3: The Friedlander law.

Sometimes, the Friedlander's law is replaced by a simpler linear approximation, see Fig. 3. The fictitious positive duration time t rf is then calculated imposing to preserve the same P ro and i r , which gives eq. ( 4).

The reflected pressure becomes, on its turn, an incident pressure for other surfaces, which of course complicates the situation: different incident waves can hit a surface besides the first one originated directly by detonation, all reflected by other surfaces, so giving a time history of the overpressure at a point that can have several successive peaks. In some cases, the effects of reflected waves can produce a peak of the overpressure greater than the one corresponding to the incident wave, especially in the case of a blast inside a building.

Another phenomenon produced by a blast is the dynamic pressure; the air behind the front of the blast wave moves like a wind, but with a smaller velocity. This wind causes a dynamic or drag pressure Q, loading a surface for the whole duration of the wind produced by the blast. Its peak value Q o is less than P so and it is delayed with respect to t A , but Q has a much longer duration (up to 2 ÷ 3 s) than t o (some ms).

The blasts can be of different types: in free-air bursts the detonation occurs in the air and the blast waves propagate spherically outwards and impinge first and directly onto the structure, without prior interaction with other obstacles or the ground. In air bursts the detonation occurs still in the air but the overpressure wave arrives at the structure after having first interacted with the ground; a Mach wave front is created. Surface bursts are the explosions where detonation occurs almost at ground surface: the blast waves immediately interact locally with the ground and they propagate next hemispherically outwards and impinge onto the structure. Finally, in internal blasts the detonation occurs inside a structure: blast waves propagate and interact with the walls, reflected waves are generated and the effects of dynamic pressure due to gas expansion are amplified by the surrounding space.

In the case of explosions in monuments, we are concerned with surface blasts, i.e. with hemispherical overpressure waves. Surface blasts result in much greater shock overpressure than air-blasts, because of the ground effect that reflects and amplifies the overpressure wave.

The case of internal blasts is also important in the study of monuments; nevertheless, in several cases an internal blast can be considered as an external one. This happens for buildings like churches or great halls, where the internal volume is so large to limit the effects of dynamic pressure and reflected waves.

The simulation of a blast can be conducted using different approaches, the most widely used being three: the JWL model, the CONWEP model and the TM5-1300 model. They are detailed in the following sections.

3 The JWL model JWL stands for Jones, Wilkins and Lee, the authors of this model, [START_REF] Jones | The detonation of solid explosives[END_REF], [START_REF] Wilkins | The equation of state of PBX 9404 and LX 04-01[END_REF], [START_REF] Lee | Adiabatic expansion of high explosives detonation products[END_REF]. Basically, JWL is a physically based model using the laws of thermodynamics to recover the physics of a chemical blast.

This model allows, in principle, to obtain a complete description of a blast phenomenon, i.e. including not only the propagation of the shock-wave in a medium, e.g. air, but also its reflection on solid surfaces and the expansion of the hot gases, i.e. the dynamic pressure. The JWL model is implemented in different commercial codes and its use needs the meshing not only of the structure, but also of the air volume involved in the blast.

The JWL model gives the overpressure P s as function of different parameters:

P s = A 1 - ω ρ R 1 ρ 0 exp -R 1 ρ 0 ρ + B 1 - ω ρ R 2 ρ 0 exp -R 2 ρ 0 ρ + ω ρ E m . (8)
In the above equation, A, B, R 1 , R 2 and ω are parameters depending upon the explosive, along with ρ 0 , its density, while ρ is the density of the detonation products and E m is the internal energy per unit mass. In addition, detonation velocity v D and the Chapman-Jouguet pressure p cj need to be specified. All the parameters are derived by fitting experimental results.

The use of JWL model allows a rather precise and complete simulation of the blast phenomenon, but its drawback is the need of discretising, finely, the charge and the fluid domain, that can be very large, besides the structure for the coupled structural analysis. Such multi-physics transient problems, with a strong coupling between fluid and structure dynamics, lead to numerical simulations that can be, in the case of a monument, very heavy, computationally speaking.

So this model, though in principle able to describe precisely the blast event and its mechanics, can be problematic to use in the case of monumental structures, where the fluid volume to discretize is very large.

Empirical models

Because of the drawbacks of the model JWL detailed hereon, empirical methods are more often used in calculations. They offer a good balance between computing cost and precision. These models are based upon the results of experimental tests and model uniquely the effects of the blast, namely the pressure field.

With such models, the characteristic parameters of the explosive serve to calculate the overpressure shock-wave and its propagation speed. The wave propagates spherically from the detonation point o to the elements of the structure. The distance of o from any impact point q on a surface of the structure and the inclination of the perpendicular to the surface with respect to the vector q -o are the only geometric parameters needed by the models.

In fact, these models consider just the incident wave, not the reflected ones, nor the dynamic pressure. So, what is mainly lost with empirical models is the possibility of taking into account for the set of reflected waves that impinge again the wall surfaces.

Nevertheless, the effect of the reflection of the shock-wave by the ground, in the case of a ground-explosion, is taken into account by specific laws, different from those modeling a free-air burst: the models for hemispherical blasts differ from those for the spherical ones, the first ones giving higher values of the overpressure to take into account for the ground reflection and the formation of the so-called Mach stem.

This phenomenon is due to the reflection of the wave pressure by a surface. In general, the overpressure shock-wave due to an air detonation is produced by an incident wave, emanating from the explosive charge, and by a wave reflected, at least, by the ground.

For small incident angles, up to about 40 • , the incident wave is ahead of the reflected wave produced by the surface and typical reflection occurs. However, for larger angles, coalescence between the incident and the reflected wave takes place, creating a Mach stem.

In ground explosions, the interaction between the ground and the blast wave takes place since the beginning, due to the closeness of the detonation point to the ground surface. Instead of the creation of a Mach front at a certain distance from the detonation point, the incident wave is reflected immediately by the ground. This wave coalescence can give much greater pressure values than the normal reflection.

If the ground was a rigid surface, the reflected pressure P r would be twice that of a free-air burst. In practice, a part of the energy is absorbed by the creation of a crater, so P r is less than the double of a free-air burst pressure, say 1.7 ÷ 1.8 times.

As said in the introduction, the two most commonly used empirical models are CONWEP, [USACE, 1986], and the Technical Manual TM5-1300, [USACE, 1990], completed by successive documents, [USACE, 2008], [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF],

All these documents are based upon the experimental results contained in another Technical Report of the U.S. Army, [START_REF] Kingery | Air blast parameters from TNT spherical air burst and hemispherical burst[END_REF], which contains a large collection of numerical data measured during blast tests and universally accepted as reference experiments in blast analysis. Such data concern both the cases of spherical and hemispherical explosions and most important they provide the data not only for the incident pressure peak P so , but also for the reflected one, P ro , which is by far most interesting for design purposes.

In Fig. 4 we show the results of these experiments for what concerns the value of P so and P ro for hemispherical blasts. It is immediate to observe that P ro is always greater than P so for any Z.

The experimental data of Kingery and Bulmash for all the parameters of the positive phase of a hemispherical blast are represented in Fig. 5. The diagrams are referred to the explosion of 1 kg of TNT and concern a distance range from 0.5 to 40 m. As they are given as functions of Z, it is possible to adapt such data also to other cases, by multiplying by W 1 3 the value of the scaled parameters, except for pressure and velocity.

New analytical expressions for the Kingery-Bulmash data

The experimental data of Kingery and Bulmash are the basis of the most common empirical methods. To this purpose, it is important to interpolate the experimental data. Several authors, starting from Kingery and Bulmash, have proposed analytical expressions fitting the data (see [START_REF] Ullah | Review of analytical and empirical estimations for incident blast pressure[END_REF] for a recent review). We consider here only the case of hemispherical blasts.

The original work of Kingery and Bulmash, [START_REF] Kingery | Air blast parameters from TNT spherical air burst and hemispherical burst[END_REF], gives the value of the logarithm base 10 of a physical parameter Y as a linear combination of powers of a function of the logarithm base 10 of Z:

log(Y ) = n j=0 c j U j , U = k 0 + k 1 log(Z). (9)
Just for example, the values of the numerical coefficients n, c j , k 0 and k 1 for some physical parameters are given in Tab. 1.

Swisdak, [START_REF] Swisdak | Simplified Kingery airblast calculations[END_REF] has proposed alternative analytical representations for the experimental data of Kingery and Bulmash, namely of the form

ln(Y ) = n j=0 k j Z j ; (10)
unfortunately, such expressions do not cover the entire range of the scaled distance Z for the experiments, so the coefficients k j are given for different sub-intervals of validity, see Tab. 2.

Several other analytical expressions have been proposed in the literature, for representing the Kingery and Bulmash experimental data or also other data; a rather exhaustive Pr, kPa Pso, kPa ir/W 1/3 , kPa-ms/kg 1/3 is/W 1/3 , kPa-ms/kg 1/3 tA/W 1/3 , ms/kg 1/3 to/W 1/3 , ms/kg 1/3 U, m/ms Lw/W 1/3 , m/kg 1/3

Figure 5: Blast parameters for the positive phase of a hemispherical blast, according to Kingery and Bulmash (from [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF]).

list of them, for the different physical parameters, can be found in [START_REF] Ullah | Review of analytical and empirical estimations for incident blast pressure[END_REF]. However, the most precise expressions, when compared to the experimental results, are those proposed by Swisdak, eq. ( 10).

In view of an empirical approach, our purpose was to dispose of reliable, i.e. as much precise as possible, analytical representations of the Kingery and Bulmash data for P so , P ro , i ow , i rw , t Aw and t ow valid for all the range of the data, i.e. for Z ∈ [0.05, 40] m/kg 1/3 . Therefore, we have limited the comparison just to the original data and to the corresponding diagrams obtained with eqs. (9), and (10). Such comparisons, along with the analytical expressions that we propose for different physical parameters of hemispherical blasts, are presented below. In all the diagrams presented below, the experimental data are represented by a series of blue dots, almost contiguous, while eq. ( 9) is represented by an orange curve and eq. ( 10) by a green one. The proposed analytical expressions are represented by red line; very often, the analytical expressions are so close to the empirical data that they are hard to be distinguished. It is shown that the proposed interpolations fit accurately the experimental data over the entire range of the scaled distance, contrary to the Swisdak [START_REF] Swisdak | Simplified Kingery airblast calculations[END_REF] interpolation.

Incident pressure P so :

P so = exp 0.26473 -1.5168 ln Z -0.079822 ln 2 Z -0.57802 sin(ln Z) -0.228409 sin 2 (ln Z) (1 + 1 2e 10Z ); (11)
In Fig. 6 we show the comparison of the results given by eq. ( 11) with those given by eq. ( 9), eq. ( 10) and the experimental results. Reflected pressure P ro :

P ro = exp 2.03043 -1.80367 ln Z -0.0929389 ln 2 Z -0.877952 sin(ln Z) -0.360291 sin 2 (ln Z) (1 + 1 2e 10Z ). (12)
The comparison of eq. ( 12) with eqs. ( 9) and ( 10) and the experimental data is shown in Fig. 7.

Scaled incident impulse i sw :

i ow = exp -2.11036 -1.06541 ln Z + 0.0733434 ln 2 Z + 0.85168 sin(ln Z) -0.0367758 sin 3 (ln Z) + 0.512623 cos(ln Z) .

(13)

The comparison of eq. ( 12) with eqs. ( 9) and ( 10) and the experimental data is shown in Fig. 8.

Scaled reflected impulse i rw :

i rw = exp -0.110157 -1.40609 ln Z + 0.0847358 ln 2 Z , (14)
The comparison of eq. ( 12) with eqs. ( 9) and ( 10) and the experimental data is shown in Fig. 9.

Scaled arrival time t Aw :

t Aw = exp -0.684763 + 1.42884 ln Z + 0.0290161 ln 2 Z + 0.410838 sin(ln Z) , (15)
The comparison of eq. ( 12) with eq. ( 10) and the experimental data is shown in Fig. 10. Table 2: Coefficients for the fitting functions of Swisdak (10), from [START_REF] Swisdak | Simplified Kingery airblast calculations[END_REF]. ∆Z: range of validity in terms of the scaled distance (1) [m/kg 1/3 ]. The comparison of eq. ( 12) with eq. ( 10) and the experimental data is shown in Fig. 11.

Scaled
∆Z k 0 k 1 k 2 k 3 k 4 k 5 k 6 Pso [kPa] 0.2 ÷ 2
From the diagrams presented in Figs. 6 to 11, it appears clearly that, unlike eqs. ( 9) and (10), eq. ( 11) fits rather well, by a unique set of coefficients, the experimental data over the whole interval of Z, i.e. from the near field to the far field. The fitting of the experimental data obtained by eq. ( 10) is always very precise, but the analytical representations that we propose, represent with a high degree of accuracy the experimental data with the advantage of giving, in all the cases presented here, a unique analytical expression of each physical parameter over the whole range of variation of Z. We hence use the above analytical expression for the method for the computation of P r that we propose in Sect. 5.2.

Influence of the type of explosive

The experimental data and formulae for the blasts are always referred to TNT, used as reference explosive. To assess the effects of a blast produced by another explosive, an equivalent weight W of TNT is computed according to the following formula, that links the weight W e of a chosen explosive to W using the ratio of the heat produced during detonation, [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF]: where H e is the heat of detonation of the explosive and H T N T is that of TNT. The values of the heat of detonation of some explosives are given in Tab. 3. 5 A procedure for computing P r using TM5-1300

W = W e H e H T N T , (17)

Influence of the direction

The values of the incident, P so , and reflected, P ro , pressures and impulses, i sw and i rw respectively, are intended for a normal shock, i.e. when the vector q -o is orthogonal to the impinged surface in q. In such a case, the reflected pressure takes its maximum local value, that decreases when the shock is not orthogonal, i.e. when the vector q -o forms an angle α > 0 with the inward normal to the impinged surface.

The effect of the lack of orthogonality in the shock is taken into account introducing the reflection coefficient c rα , defined as

c rα = P rα P so , ( 18
)
where P rα is the peak of the reflected pressure in q for a surface inclined of the angle α onto the direction of q -o. It is worth noticing that the reflection coefficient is defined as the ratio of the reflected inclined pressure P rα with the incident pressure P so , not P ro .

The value of c rα has been evaluated experimentally, and the results are shown in Fig. 12; this figure has been obtained from the parametrized curves given in [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF]. To remark that c rα depends not only upon α, but also upon P so .

In Fig. 13 we show also the interpolation of the empirical values of i rαw , still obtained from [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF]].

Looking at Fig. 12, one can see that c rα is not a decreasing function of α, as it could be expected. Actually, for α >∼ 40 • , c rα increases reaching a maximum and then it decreases again. This fact is attributed to the creation of the Mach stem cited above. For small values of P so the behavior is more complicated, but this presumably could be the effect of uncertainties and imperfections in the experimental data.

Figure 12: Surface representing c rα (α, P so ).

P r calculation

We detail in the following the sequence of the calculations to obtain, for a given blast, the time history of the reflected pressure P r (q, t) at the time t on a point q of a surface placed at a distance R from the detonation point o and whose inward normal forms an angle α with the vector q -o.

We make the following assumptions: i) the blast occurs at point o and is produced by a mass W e of a given explosive; ii) the blast is hemispherical; iii) the time rate decrease is of the same type for P s and P r and it is ruled by the Friedlander's equation; iv) the calculation of the time history of P s and P r , as well as the impulse, is made pointwise; v) only the positive duration phase is considered (this assumption, usually done, is justified because of the much larger values of the positive pressures in comparison with the negative ones); finally, vi) reflected waves and dynamic pressure are ignored.

Then, the calculation sequence is the following one:

• a time duration t max is fixed;

• t max is subdivided into time intervals dt;

• the impinged surface is discretized into regular patches, whose centroids are points q;

• the equivalent mass of TNT W is calculated using eq. ( 17);

• then, for each point q we calculate:

the distance from the blast point: R = ||q -o||;

the inward normal n to the wall;

the angle α between q -o and n: the scaled distance Z, eq. ( 1);

cos α = q -o R • n; (19)
the normal incident pressure peak P so using the fitted curve of Kingery-Bulmash for hemispherical blasts, eq. ( 11), Fig. 6;

the normal reflected pressure peak P ro using the fitted curve of Kingery-Bulmash for hemispherical blasts, eq. ( 12), Fig. 7;

-the ratio c = P ro P so ; (20)
the angular coefficient c rα , eq. ( 18), by a linear interpolation of the experimental data represented by the surface in Fig. 12;

the ratio c r0 = P r0 P so ,

i.e. P rα for α = 0 • , still linearly interpolating the data of the surface in Fig. 12;

the corrected value of P rα as

P rα = c c rα c r0 P so = c rα c r0 P ro ; ( 22
)
this correction is done to adapt the data of the interpolating surface in Fig. 12 to those of the Kingery-Bulmash fitted curves, Fig. 5, more conservative;

the reduction coefficient

c red = P rα P ro = c rα c r0 ; (23)
the positive normal scaled reflected impulse i rw using the fitted curve of Kingery-Bulmash for hemispherical blasts, eq. ( 14), Fig. 9;

the effective angular reflected impulse i rα as

i rα = c red i rw W 1 3 ; (24)
the value of i rα is voluntarily not calculated interpolating the surface of experimental data, Fig. 13, like for c rα , because it has been observed, numerically, that the procedure described here is conservative;

the scaled arrival time t Aw using the fitted curve of Kingery-Bulmash for hemispherical blasts, eq. ( 15), Fig. 10;

-the arrival time t A as t A = t Aw W 1 3 ; (25)
the scaled positive duration time t ow using the fitted curve of Kingery-Bulmash for hemispherical blasts, eq. ( 16), Fig. 11;

the positive duration time t o as

t o = t ow W 1 3 ; (26)
the fictitious positive duration time t rf , eq. ( 4), as

t rf = 2 i rα P rα ; (27)
check on t o : if t o < t rf then put t o = 1.1 t rf ; this is done to avoid pathological situations, due to the fact that the Kingery-Bulmash curve for t o does not cover low ranges of Z, see Figs. • the time history of the pressure wave can now be calculated:

t = n dt;

-∀q: * if t < t A or t > t A + t o then P r (q, t) = 0; * else, use the Friedlander equation ( 5) with P rα in place of P so to evaluate P r (q, t);

iterate on t until t > t max .

This sequence has been implemented in a program for the formal code Mathematica and applied to the case study presented in Sec. 7.

The CONWEP model

CONWEP is the acronym of CONventional Weapons Effects Programme, a study made by USACE, [USACE, 1986], for the simulation of the effects of a blast produced by conventional, i.e. not nuclear, explosives. The report [USACE, 1986] uses in many parts the same experimental results shown above, but not completely; in particular, a noticeable difference is the way in which the reflected pressure is calculated; in place of using the data represented in Fig. 13, and integrating the effect of the Mach stem, CONWEP makes use of the following law, using circular functions:

P rα = P so (1 + cos α -2 cos 2 α) + P ro cos 2 α if cos α ≥ 0, P so if cos α < 0. ( 29
)
The algorithm CONWEP is today implemented in different finite elements commercial codes, e.g. in LS-DYNA, AUTODYN or ABAQUS. In particular, in ABAQUS, the code we used for numerical simulations with CONWEP, the following assumptions are made: i) for each point q of the impinged walls, P rα is calculated according to eq. ( 29); ii) reflected waves are ignored; iii) the dynamic pressure is neglected; iv) the negative phase is taken into account; v) hemispherical blasts can be modeled.

The assumptions made by CONWEP are similar to those introduced in Sec. 5.2, except for the calculation of P rα and the fact that in CONWEP the negative phase is considered too. This is why a comparison of the two methods is interesting, see Sec. 7.

A final, important remark: neglecting the reflected waves and the dynamic pressure means that there is no difference between internal and external blasts: the only geometric parameters that matter are R and α, regardless the surrounding geometry. For numerical simulations, this is interesting and important because it allows modeling only a part of the building, the closest one to the detonation point, and reasonably neglecting the effects, at least in terms of applied pressures, of more distant parts. Therefore, only a part of the building can be modeled, reducing in this way the calculation cost, but of course reducing in the same time the degree of accuracy compared to more physical approaches, like the JWL, which are able to cover all the aspects of the blast wave.

Comparison of the models on a case study

Even though CONWEP and JWL models are implemented in some commercial finite element codes, this is apparently not the case for the model TM5-1300. This is probably due to the fact that it is much easier to implement the dependence on the direction described by eq. ( 29) than the more accurate of TM5-1300, described in Fig. 12.

As said above, in the case of large buildings, like monuments, to make a complete nonlinear fluid-solid simulation using the model JWL is computationally expensive and unnecessary from an engineering point of view. This is why CONWEP is generally preferred, but its accuracy has to be assessed. The following comparison is hence made with the objective of evaluating the response of the model CONWEP, in comparison with both JWL and TM5-1300. To evaluate the differences in the prediction of the blast parameters, and namely of the reflected pressure between the three models, we consider an explosion in the interior of a building covered by a barrel vault, see Fig. 14. This is intended to simulate the typical structure of a monumental building, like a church or a palace gallery and so on. The dimensions used in this case are: width: 12 m, height of the walls: 12 m, radius of the vault: 6 m, total height: 18 m.

The explosive charge is composed by W = 20 kg of TNT, and detonates at the instant t = 0 at point o = (0, 0). The problem is treated as a planar one, exception made for the JWL model, where a 3D approach is used.

The evaluation of P rα is done in four different ways: in the first one, we use the JWL model; the results are indicated as JWL. Then, we have used the data of TM5-1300 and following modifications, as integrated into [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF], and according to the calculation sequence presented in Sec. 5.2; the results are indicated as TM5-1300. The third way, is the use of the data of TM5-1300 and following modifications, as integrated into [START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF], but evaluating the pressure P rα according to eq. ( 29); the results are indicated as CONWEP*. Finally, we have applied the model CONWEP as implemented in ABAQUS to the above structure; the results are indicated as CONWEP.

The results TM5-1300 and CONWEP* have been obtained using the program we have done in Mathematica.

The results from JWL

The analyses with the JWL model are run with AUTODYN Hydrocode, a code allowing good facilities in modeling the set of explosive, air domain and structure, where the equations of mass, momentum and energy conservation for inviscid flows are coupled with the dynamic equations of solid continua. A Coupled Eulerian-Lagrangian (CEL) method is used: explosive and surrounding air, as fluids, are modeled with an Eulerian frame, while solid walls are identified in a Lagrangian reference.

The 3D geometrical model, taking advantage of symmetries, consists of a quarter of the entire domain. A volume having the transversal section like in Fig. 14 and a depth of 2 m is hence discretized. Ground and planes of symmetry are modeled as reflecting planes to prevent flow of material through them. The ending transversal surfaces are modeled as transmitting planes, i.e. boundary surfaces whereupon the gradients of velocity and stress are put to zero. This approach is used to simulate a far field solution at the boundary, it is only exact for outflow velocities higher than the speed of sound and is an approximation for lower velocities.

A preliminary study of the solution sensitivity to the cells size has been made studying the peak of the reflected pressure at point A of Fig. 14 for five different meshes. For each one of them, the dimensions of the volume elements discretizing the explosive charge and of the surface elements modeling the explosive-air and air-solid interfaces are decreased each time by a factor 1.5 starting from the coarsest mesh.

In Tab. 4 and in Fig. 15 we show the relative percent error, computed with respect to the finest discretization; it is apparent that it remains very small in all the cases. Attention needs to be paid when dealing with a CEL approach: more precisely, an almost equal size of the Eulerian and Lagrangian elements and an average grid size smaller than the thickness of the solid walls must be guaranteed. The simulations have been made with the mesh M2. The results of the numerical simulation are presented in Fig. 16. It is apparent the formation not only of the principal, hemispheric shock-wave, but also of the reflected ones, that produce the Mach stem, well visible in the photogram at t = 14 ms, at the bottom of the waves reflected by the vertical walls. At t ∼ 48 ms, the incident and reflected waves focus just at the top of the vault, giving rise to a very high localized pressure. Then, the reflected waves propagate downwards and upwards and decrease in intensity. Though not well visible in the figures, also the dynamic pressure is taken into account in this analysis.

The focusing of waves in the vault is not a surprising result and it is very similar to what observed by Rayleigh for acoustic waves in the study of the whispering galleries, [J. W. [START_REF] Strutton -Lord Rayleigh | The problem of the whispering gallery[END_REF][START_REF] Strutton -Lord Rayleigh | Further applications of Bessel's functions of high order to the whispering gallery and allied problems[END_REF]: due to its geometry, the vault behaves like a concave (i.e. converging) mirror for the shock waves, which has the tendency to collect the blast energy. As far as it concerns blasts, we can hence say that barrel vaults and domes have a dissatisfactory behavior.

What is clear, is the true complexity of the pressure dynamics, which cannot be predicted in advance and, most important, that is strongly depending upon the geometry and dimensions of the structure. Namely, the focusing of the shock waves described hereon takes place just because of the geometry of the vault and of its dimensions. To this purpose, we have also performed a simulation where the lowest part of the vertical walls, for a height of 7 m, are replaced by transmitting planes, to simulate the presence of openings, like in the case of an aisle. In such simulation, the focus, at the key of the vault, of the reflected waves, gives a sensibly smaller peak of the reflected pressure, reduced by a factor 2.72 with respect to that calculated in the above simulation, see Fig. 17, passing from 0.299 to 0.109 MPa.

It is not possible to recover such a complete description of the pressure history everywhere in the fluid domain, and on the walls, using the other empirical models. Therefore, it is interesting, for the comparison, to have the time-history of the pressure at a given point of the solid boundary. To this purpose, we give in Fig. 18 the diagram of the time variation of the reflected pressure P r at four points of the boundary, indicated in Fig. 14 as points A, at the base of the vertical vault (Z = 2.21 m/kg 1/3 ), B, at 5.63 m from the base of the vertical wall (Z = 3.03 m/kg 1/3 ), C, at the springing of the vault (Z = 5.40 m/kg 1/3 ), and D, the key of the vault (Z = 6.63 m/kg 1/3 ).

It is apparent that the decay phase is not exactly as predicted by the Friedlander's law, and this because of the reflected waves, clearly visible on each diagram as secondary peaks of the curve. It is interesting to notice that the successive peaks are not necessarily decreasing, which confirms the complexity of the interactions and dynamics of the reflected waves. Table 5: Simulation results for models TM5-1300 (P rα 1) and CONWEP* (P rα 2).

# x y R Z α ϕ Prα1 Prα2 t A to t A + to b [m] [m] [m] m/kg 1/3 [ • ] [ • ] [MPa] [ms] [ms] [

7.2

The results from TM5-1300 and CONWEP*

The spatial distributions of the maximum value of P rα as given by TM5-1300 and CON-WEP*, and calculated as specified above and in Sect. 5.2, are shown in Fig. 14. Such distributions are similar, though some differences exist: first of all, CONWEP* gives maximum values of P rα that are almost always less or equal of those given by TM5-1300.

Then, the greatest differences appear for 0 • < α < 90 • ; this can be explained by the fact that the angular variation taken by CONWEP does not take into account for the formation of the Mach stem.

Such occurrence happens for α 40 • and it is clearly indicated by the local increase of P rα in the diagram of TM5-1300, that shows two humps: at midway of the vertical wall and at the springing of the vault where α 40 • in both the cases; this fact can be of a great importance for vaulted structures, because an increase of P rα in the zone between 0 • and 30 • can be very dangerous for the stability of the vault, that normally has on its back a filling with a material like rubble or gravel in order to improve the structural stability of the structure.

The numerical data of the simulations TM5-1300 and CONWEP* are shown in Tab. 5; ϕ is the angle formed by the normal n with the axis y, while b is the coefficient appearing in the Friedlander's law, eq. ( 5).

Observing the results concerning P rα and t o , we see clearly that the peak of the shock wave decreases with the distance R, passing from a maximum of 0.736 MPa for R = 6 m, to a minimum of 0.063 MPa for R = 18 m, while its time duration increases, passing from 7.2 ms to 11 ms.

In Fig. 19 we show the same curves of Fig. 18 but now obtained with the models TM5-1300, red curves, and CONWEP*, green curves. The red and green curves are distinct only for α = 90 • , for the way the values of CONWEP* are calculated.

Table 6: Simulation results for the CONWEP model. The problem is still considered planar and the mesh consists of linear quadrilateral elements with an average size of 0.2 m and a thickness of 0.1 m.

# x y R Z α ϕ Prα t A to t A + to [m] [m] [m] 1/3 [ • ] [ • ] [MPa] [MPa] [ms
The results obtained with this simulation are reported in Tab. 6, while the diagrams of the time variation of the reflected pressure P r at the four points A, B, C and D of the boundary are still shown in Fig. 19.

Comparison

In order to compare the results of the above simulations, something must be mentioned about the results from JWL. On one hand, the multiple reflections of the shock-wave, visible in Fig. 16, completely alter the time variation of P r after the positive phase, to such a point that in some cases the duration of the pure decaying phase can be really difficult to be determined. On the other hand, the multiple reflections add together and can considerably increase the value of P r . Nevertheless, it is still possible, at least for a part of the boundary surface, to focus on just the characteristic elements of the positive phase, namely P rα , t A and t o . Such a comparison, of course, is possible only if the peak P rα has not yet been affected by other reflected shock-waves and if a positive phase is still clearly distinguishable, which depends again on the reflected waves.

To assess the results, we have chosen to consider the relevant physical parameters of the blast: P rα , t A and t o . We have compared them as evaluated at the four points A, B, C and D of Fig. 14; these values are summarized in Tab. 7.

In Fig. 20 we show a comparison of the results for the values of P rα obtained with JWL, TM5-1300, CONWEP* and CONWEP. The latter ones, as already pointed out, use the empirical data for computing P r for α = 0, but a different angular dependence from α, eq. ( 29).

It can be observed that the values given by CONWEP and the code ABAQUS are prac- tically coincident with those calculated using the model CONWEP*, except some oscillations, to be imputed to the finite element approximation. Hence, on one side the computation of P ro is practically the same for the three models, but what changes is the way the effect of the inclination α is taken into account. In particular, except the value α = 0, where the three models give the same P ro , for all the other values of α, CONWEP and CONWEP* underestimate the value of P rα with respect to TM5-1300.

TM5-1300 and JWL give values that are comparable until point C; here, JWL diverges and gives values of P rα that can be considerably greater. This is actually the effect of converging the reflected waves, that increases significantly the value of the overpressure. It is a local phenomenon, essentially depending upon the geometry and dimensions of the structure, as shown above. This same phenomenon is noticed also on the time variation of P rα for point D in Fig. 18, where the peak due to the reflected waves is clearly visible.

It is worth noticing that, besides the convergence of the reflected waves starting from point C, JWL and TM5-1300 give not only comparable values of P rα , but also of its time variation. In Fig. 21 we show the time diagrams of JWL, reduced to the positive phases, presented together with those of TM5-1300. The curves are in a rather good agreement, apart from the peak on the curve of point D, due to reflected waves.

As far as it concerns the time durations, represented in Fig. 22, TM5-1300 and CONWEP* give, of course, the same values, while CONWEP underestimates slightly the durations for small values of Z. Nevertheless, these discrepancies are not significant. The diagrams in Figs. 20 and 22 explain the aforementioned slight differences, namely in the value of t A , appearing in Fig. 19 between the curves obtained with CONWEP and those relative to TM5-1300 and CONWEP*.

For what concerns JWL, the curve of t A is in a very good agreement with those of TM5-1300 and CONWEP until point C, where it diverges, once more due to wave reflection. The curve of t A + t o has been obtained by interpolating the values estimated for points A, B, C and D for the duration t o because, as mentioned above, it is not easy to determine the exact duration of the positive phase for all the points, due to the interaction of the reflected waves. Globally, the four models give values that are comparable, apart the zone of the wave reflections, where the values given by JWL diverge from those of the other models.

The decay coefficient

A point that can be raised, is the assessment of the decay coefficient b in the Friedlander's equation ( 5). We have seen that in the procedure detailed in Sect. 5.2, we propose to calculate b as the solution of eq. ( 28) once P rα , i rα and t o determined using their respective analytical expressions given in Sect. 4.1.

Of course, such an estimation of b should be validated. To this purpose, let us consider again Fig. 21: we observe that the peaks A, B and C as evaluated by JWL and TM5-1300 are in a rather good agreement (this is not the case for peak D, due to the already remarked effect of focused reflected waves, that of course completely invalidate a Friedlander's type decay). The red curves, relative to TM5-1300, decrease slightly more rapidly than those of JWL, in black. This means that b is slightly overestimated by the procedure sketched in Sect. 5.2.

We have then compared the value of b with the data of the literature. A recent paper, [START_REF] Karlos | Analysis of the blast wave decay coefficient using the Kingery-Bulmash data[END_REF], considers the assessment of the decay coefficient b. There, the authors propose a deep analysis of the problem, presenting different laws for b existing in the and in a law of the type (9), for spherical or hemispherical blasts, in both the cases for the incident and reflected overpressure. Also in [START_REF] Ullah | Review of analytical and empirical estimations for incident blast pressure[END_REF] different expressions for b are given. We have then considered such analytical expressions, but we have found a great dispersion of the results, which renders the comparison difficult to be done in a reliable way.

So, we have just compared our results with those given by the law proposed by [START_REF] Karlos | Analysis of the blast wave decay coefficient using the Kingery-Bulmash data[END_REF]. The comparison is shown in Fig. 23 for the reflected pressure, hemispherical blasts. It shows a good agreement between the two estimations of b and, considering the incertitude of the tests on b due to the difficulty of experimental measures, we can conclude on the validity of the way proposed in Sect. 5.2 for evaluating b.

Conclusion

The results of the simple example treated above clearly show that the distribution of the pressures given by TM5-1300 is sensibly similar to that given by the more exact JWL model, apart from those zones where the interaction of the reflected waves alters the overpressure distribution. This phenomenon, as already explained, strongly depends upon the geometry and the dimensions of the structure and cannot be predicted a priori.

Simulations done with the JWL model are hence the best way to predict the effects of an explosion, especially when the analysis concerns monumental structures, which involve non-typical geometries that promote multiple reflections and the concentration of the blast energy to specific places and structural elements.

However, the use of JWL can be impractical in the case of very large structures: the need of discretizing finely not only the structure but, even more, the explosive charge and the volume of air can be really problematic for some studies. That is why the use of empirical models, like TM5-1300 and CONWEP, can be very useful in such cases. If abstraction is made concerning the effect of reflected waves, the results of TM5-1300, though obtained neglecting the dynamic pressure, are close to those of JWL, which confirms that the dynamic pressure can be neglected in simulations.

The comparison made of CONWEP, TM5-1300 and JWL confirms what usually stated in the literature: TM5-1300 is more accurate than CONWEP. Nevertheless, as mentioned above, it is this last model that normally is implemented in commercial finite element codes. Its use must hence be accurately considered, adopting a multiplying factor p f > 1. Such a coefficient, affecting only the reflected pressure P rα and not the other blast parameters, namely t A and t o , is needed to obtain values that are similar to those predicted by TM 5-1300 and JWL; it represents hence, on the average, the ratio between the value of P rα given by JWL and that given by CONWEP. This coefficient must be chosen in function of the problem at hand, namely considering the characteristic dimensions of the building: as apparent from Fig. 20, p f depends upon Z, hence upon R, and upon the angle α. For monuments of large dimensions, we can suggest a value in the interval 1.5 ÷ 2.5, in order to take into account, though indirectly, of the effect of reflected waves.

In this paper, we have also proposed new analytical expressions for the representation of the Kingery and Bulmash experimental data on hemispherical blasts. These expressions are really accurate and have the advantage of representing the different physical parameters of a blast over the entire range of variation of the scaled distance Z.

  of violences and iconoclastic destructions. The examples of the Cathedral of Reims in 1914, the Buddha statues of Bamyan in 2001 and the more recent destructions at Palmyra in 2015 and 2016 or at Mosul in 2017 are just some few recent examples of that.

Figure 2 :

 2 Figure 2: Influence of the distance R on the positive phase of a blast.

Figure 4 :

 4 Figure 4: Results of the Kingery and Bulmash tests for P so and P ro (hemispherical blasts).

22Figure 14 :

 14 Figure 14: Parameters of positive phase of shock hemispherical wave of TNT charges from surface bursts (modified from [9]). Scaled Distance Z = R/W 1/3 0.05 0.05 0.07 0.1 0.2 0.3 0.5 0.7 1 2 3 4 5 6 7 8 10 20 30 40 50 0.01 0.02 0.03 0.05 0.1 0.2 0.3 0.5 1 2 3 5 10 20 30 50 100 200 300 500 1000 2000 3000 5000 10000 20000 30000 50000 100000 200000 300000 500000 1000000

Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 6789 Figure 6: Comparison of experimental and fitted values of P so . 0.1 0.5 1.0 5.0 10.0

Figure 10 :

 10 Figure 10: Comparison of experimental and fitted values of t Aw .

Figure 11 :

 11 Figure 11: Comparison of experimental and fitted values of t ow .

Figure 13 :

 13 Figure 13: Surface representing i rαw (α, P so ).

  5 and 11; -solve numerically the equation b -1 + e -b b 2 P rα t o = i rα (28) to determine the coefficient b of the Friedlander's equation.

Figure 14 :

 14 Figure 14: Case study: layout and P rα distribution for the model TM5-1300 and CON-WEP*.

Figure 15 :

 15 Figure 15: Convergence study for different meshes.

Figure 16 :

 16 Figure 16: Evolution of the shock-wave as predicted by JWL; the color scale represents the overpressure and runs from blue for the minimum value, 0 kPa, to red for the maximum one, 300 kPa

Figure 17 :

 17 Figure 17: Comparison of P rα as obtained by JWL for the cases of vertical openings until the height of 7 m ( curve 1) and of complete vertical walls (curve 2).

Figure 19 :

 19 Figure 19: Time variation of P rα at four points of the solid boundary, as predicted by TM5-1300 (red curves), CONWEP (blue curves) and CONWEP* (green curves).

Figure 20 :

 20 Figure 20: Comparison of P rα as obtained by the four models.

Figure 21 :

 21 Figure 21: Comparison of the positive phases for points A, B, C and D as given by JWL (black) and TM5-1300 (red).

Figure 22 :

 22 Figure 22: Comparison of the arrival time t A and of the positive phase duration t o as obtained by the four models.

Figure 23 :

 23 Figure23: Comparison of the decay coefficient for the reflected pressure, hemispherical blasts, as evaluated by the procedure in Sect. 5.2, red curve, and after[START_REF] Karlos | Analysis of the blast wave decay coefficient using the Kingery-Bulmash data[END_REF], blue curve.

Table 1 :

 1 Coefficients for the fitting functions ofKingery and Bulmash (9), from[START_REF] Unoda | International Ammunition Technical Guideline -IATG 01.80: Formulae for ammunition management (2nd Edition)[END_REF]. ∆Z: range of validity in terms of the scaled distance (1) [m/kg 1/3 ].

		Pso	Pro	isw		irw
		[kPa]	[kPa]	[kPa ms/kg 1/3 ]	[kPa ms/kg 1/3 ]
	∆Z	0.05 ÷ 40	0.05 ÷ 40	0.0674 ÷ 0.955	0.955 ÷ 40	0.05 ÷ 40
	n	8	9	4	7	3
	k 0	-0.214362789151	-0.240657322658	2.06761908721	-1.94708846747	-0.246208804814
	k 1	1.35034249993	1.36637719229	3.0760329666	2.40697745406	1.33422049854
	c 0	2.78076916577	3.40283217581	2.52455620925	1.67281645863	2.70588058103
	c 1	-1.6958988741	-2.21030870597	-0.502992763686	-0.384519026965	-0.949516092853
	c 2	0.154159376846	0.218536586295	0.171335645235	-0.0260816706301	0.112136118689
	c 3	0.514060730593	0.895319589372	0.0450176963051	0.00595798753822 -0.0250659183287
	c 4	-0.0988534365274	0.24989009775 -0.0118964626402	-0.014544526107
	c 5	-0.293912623038 -0.569249436807	-0.00663289334734
	c 6	-0.0268112345019	0.11791682383	-0.00284189327204
	c 7	0.109097496421	0.224131161411		0.0013644816277
	c 8	0.00162846756311	0.0245620259375
	c 9		-0.455116002694

Table 3 :

 3 Heat of detonation of different explosives (from[START_REF] Karlos | Calculation of Blast Loads for Application to Structural Components[END_REF]).

	Type of explosive Heat of detonation
		[MJ/kg]
	TNT	4.10 ÷ 4.55
	C4	5.86
	RDX	5.13 ÷ 6.19
	PETN	6.69
	Pentolite 50/50	5.86
	Nitroglycerin	6.30
	Nitromethane	6.40
	Nitrocellulose	10.60

Table 4 :

 4 Convergence study for different meshes. Relative error on P ro at point A [%] 0. 0.018 0.059 0.106 0.274

	Meshes	M 1 M 2	M 3	M 4	M 5
	Number of elements [10 6 ]	2.23 1.31 0.86 0.59	0.4

  Figure 18: Time variation of P rα at four points of the solid boundary, by JWL.

				800 0.80
					A
				588 0.58
			P (kPa) Pr (MPa)	375 0.37	B		Peak produced by the superposition of the reflected waves
				163 0.16	C	D
				0,00 0 -50 -0.05	25,00 25	50,00 50	75,00 75	100,00 100
							t (ms) time (ms)
												ms]
	1.	6.	0.	6.	2.21042	0.	90.	0.765485	0.765485	5.80175	5.70184	11.5036	3.88204
	2.	6.	0.188558	6.00296	2.21151	1.8	90.	0.760041	0.76395	5.8068	5.70346	11.5103	3.87904
	3.	6.	0.377488	6.01186	2.21479	3.6	90.	0.752468	0.759364	5.82198	5.70837	11.5304	3.87009
	4.	6.	0.567167	6.02675	2.22027	5.4	90.	0.742913	0.751779	5.84741	5.71668	11.5641	3.85531
	5.	6.	0.757976	6.04769	2.22799	7.2	90.	0.731543	0.741284	5.88325	5.72862	11.6119	3.83489
	6.	6.	0.950307	6.07479	2.23797	9.	90.	0.71823	0.728	5.92975	5.74448	11.6742	3.80909
	7.	6.	1.14456	6.10819	2.25028	10.8	90.	0.703064	0.712079	5.98724	5.76463	11.7519	3.77826
	8.	6.	1.34116	6.14807	2.26497	12.6	90.	0.686672	0.693701	6.05612	5.78955	11.8457	3.74281
	9.	6.	1.54054	6.19462	2.28212	14.4	90.	0.668886	0.673071	6.13689	5.81981	11.9567	3.70319
	10.	6.	1.74316	6.24809	2.30182	16.2	90.	0.64965	0.650416	6.23013	5.85607	12.0862	3.65994
	11.	6.	1.94952	6.30877	2.32417	18.	90.	0.629087	0.62598	6.33654	5.89908	12.2356	3.61361
	12.	6.	2.16013	6.377	2.34931	19.8	90.	0.607325	0.60002	6.45693	5.94972	12.4067	3.5648
	13.	6.	2.37557	6.45316	2.37737	21.6	90.	0.584501	0.572801	6.59223	6.00893	12.6012	3.51413
	14.	6.	2.59643	6.5377	2.40851	23.4	90.	0.560758	0.544596	6.74353	6.07776	12.8213	3.46222
	15.	6.	2.82339	6.6311	2.44292	25.2	90.	0.536604	0.515674	6.91207	6.15734	13.0694	3.40968
	16.	6.	3.05715	6.73396	2.48081	27.	90.	0.514591	0.486303	7.09928	6.24889	13.3482	3.35705
	17.	6.	3.29853	6.84692	2.52243	28.8	90.	0.491878	0.456741	7.3068	6.35366	13.6605	3.30483
	18.	6.	3.54839	6.97073	2.56804	30.6	90.	0.469527	0.427236	7.53652	6.47291	14.0094	3.2534
	19.	6.	3.80772	7.10624	2.61796	32.4	90.	0.448879	0.398019	7.79061	6.60789	14.3985	3.20301
	20.	6.	4.0776	7.25443	2.67256	34.2	90.	0.427992	0.369304	8.07157	6.75974	14.8313	3.1537
	21.	6.	4.35926	7.41641	2.73223	36.	90.	0.403807	0.341285	8.38228	6.92942	15.3117	3.10531
	22.	6.	4.65408	7.59345	2.79745	37.8	90.	0.384418	0.314131	8.72609	7.11756	15.8437	3.05738
	23.	6.	4.96363	7.78702	2.86876	39.6	90.	0.367394	0.287992	9.10688	7.32437	16.4312	3.00913
	24.	6.	5.28971	7.99882	2.94679	41.4	90.	0.363427	0.262988	9.52915	7.54941	17.0786	2.95941
	25.	6.	5.63438	8.23081	3.03226	43.2	90.	0.343876	0.239219	9.99818	7.79141	17.7896	2.90667
	26.	6.	6.	8.48528	3.126	45.	90.	0.317226	0.216758	10.5202	8.04807	18.5683	2.84899
	27.	6.	6.38935	8.76492	3.22902	46.8	90.	0.304107	0.195655	11.1025	8.31585	19.4183	2.78408
	28.	6.	6.80566	9.07288	3.34248	48.6	90.	0.270092	0.175939	11.7538	8.58988	20.3436	2.70941
	29.	6.	7.25275	9.41289	3.46774	50.4	90.	0.2356	0.157617	12.4844	8.86397	21.3483	2.62239
	30.	6.	7.73515	9.78941	3.60645	52.2	90.	0.203598	0.140679	13.3068	9.13083	22.4376	2.5206
	31.	6.	8.25829	10.2078	3.76059	54.	90.	0.177162	0.125096	14.2361	9.38266	23.6187	2.40217
	32.	6.	8.82873	10.6746	3.93255	55.8	90.	0.158073	0.110828	15.2906	9.61206	24.9026	2.26615
	33.	6.	9.45449	11.1976	4.12525	57.6	90.	0.141759	0.0978226	16.4929	9.81332	26.3062	2.11303
	34.	6.	10.1454	11.7869	4.34232	59.4	90.	0.127267	0.0860183	17.8711	9.98409	27.8552	1.94507
	35.	6.	10.914	12.4545	4.58828	61.2	90.	0.112877	0.0753469	19.4605	10.127	29.5875	1.76667
	36.	6.	11.7757	13.2161	4.86887	63.	90.	0.0983451	0.0657355	21.3057	10.2512	31.557	1.58446
	37.	5.96239	12.6707	14.0035	5.15893	58.3816	83.5816	0.106085	0.0626654	23.2452	10.3613	33.6065	1.42342
	38.	5.82149	13.4527	14.6582	5.40014	52.5888	75.9888	0.11852	0.0616353	24.88	10.4518	35.3318	1.31115
	39.	5.60197	14.1489	15.2176	5.6062	47.4129	69.0129	0.0852345	0.0607484	26.2908	10.5337	36.8245	1.23022
	40.	5.31949	14.7754	15.7038	5.78534	42.6467	62.4467	0.0750139	0.0599822	27.5269	10.6104	38.1374	1.17038
	41.	4.98432	15.3401	16.1296	5.94219	38.1727	56.1727	0.0684349	0.0593205	28.6161	10.6825	39.2987	1.12546
	42.	4.60409	15.8474	16.5026	6.07963	33.9164	50.1164	0.0656936	0.0587507	29.5755	10.7497	40.3252	1.09142
	43.	4.185	16.2995	16.8282	6.19956	29.8268	44.2268	0.0634732	0.0582626	30.4163	10.8113	41.2276	1.0655
	44.	3.7324	16.6978	17.1098	6.30332	25.8672	38.4672	0.0616696	0.0578481	31.1463	10.8669	42.0131	1.04573
	45.	3.2511	17.0429	17.3502	6.39186	22.0096	32.8096	0.0602101	0.0575007	31.7709	10.9159	42.6868	1.0307
	46.	2.74559	17.335	17.551	6.46586	18.2322	27.2322	0.0590426	0.057215	32.2942	10.9579	43.2522	1.01937
	47.	2.22013	17.5741	17.7138	6.52583	14.517	21.717	0.0581296	0.0569866	32.7191	10.9927	43.7119	1.01097
	48.	1.67885	17.7603	17.8395	6.57213	10.8487	16.2487	0.057444	0.0568122	33.0477	11.02	44.0677	1.00494
	49.	1.12576	17.8934	17.9288	6.60503	7.21431	10.8143	0.0569668	0.0566894	33.2814	11.0396	44.321	1.00089
	50.	0.564835	17.9734	17.9822	6.62471	3.60178	5.40178	0.0566853	0.0566164	33.4213	11.0514	44.4727	0.998567
	51.	0.	18.	18.	6.63126	0.	0.	0.0565922	0.0565922	33.4678	11.0553	44.5232	0.997807

Table 7 :

 7 Comparison of the characteristic parameters evaluated at points A, B, C and D of Fig.14as calculated by JWL, TM5-1300, CONWEP* and CONWEP; R 1 =JWL/TM5-1300, R 2 =JWL/CONWEP, R 3 =TM5-1300/CONWEP.

		JWL TM5-1300 CONWEP* CONWEP R 1 R 2 R 3
				P rα [MPa]
	A 0.746	0.765	0.765	0.771	0.97 0.97 0.99
	B 0.339	0.344	0.239	0.245	0.98 0.72 1.49
	C 0.131	0.118	0.062	0.062	1.11 2.11 1.90
	D 0.299	0.056	0.056	0.058	5.33 5.15 0.96
				t A [ms]
	A	5.9	5.8	5.8	5.6	1.02 1.05 1.00
	B 10.1	10.0	10.0	9.5	1.01 1.06 1.05
	C 24.3	24.9	24.9	24.2	0.97 1.00 1.03
	D 39.3	33.5	33.5	34.0	1.17 1.15 0.98
				t o [ms]
	A	5.7	5.7	5.7	5.8	1.00 0.98 0.98
	B	7.7	7.8	7.8	7.5	0.99 1.03 1.04
	C 13.9	10.4	10.4	10.2	1.34 1.36 1.02
	D 16.9	11.0	11.0	11.4	1.54 1.48 0.96

The same case study has finally been implemented in ABAQUS in order to make a comparison also with the results given by the model CONWEP. The comparison, as already said, will serve to understand the quality of the predictions given by CONWEP and to calibrate its use.


Maytheewat Aramrattana
email: maytheewat.aramrattana@hh.se

Raj Haresh Patel
email: patel@eurecom.fr

Cristofer Englund
email: cristofer.englund@ri.se

Jérôme Härri
email: haerri@eurecom.fr

Jonas Jansson
email: jonas.jansson@vti.se

Christian Bonnet
email: bonnet@eurecom.fr

Evaluating Model Mismatch Impacting CACC Controllers in Mixed Traffic using a Driving Simulator

At early market penetration, automated vehicles will share the road with legacy vehicles. For a safe transportation system, automated vehicle controllers therefore need to estimate the behavior of the legacy vehicles. However, mismatches between the estimated and real human behaviors can lead to inefficient control inputs, and even collisions in the worst case. In this paper, we propose a framework for evaluating the impact of model mismatch by interfacing a controller under test with a driving simulator. As a proofof-concept, an algorithm based on Model Predictive Control (MPC) is evaluated in a braking scenario. We show how model mismatch between estimated and real human behavior can lead to a decrease in avoided collisions by almost 46%, and an increase in discomfort by almost 91%. Model mismatch is therefore non-negligible and the proposed framework is a unique method to evaluate them.

I. INTRODUCTION

Commuters today often face issues such as traffic jams, stop-and-go scenarios, collision avoidance maneuvers, etc. In such circumstances, the vehicles need to frequently accelerate and brake to follow the traffic rhythm, which not only increases fuel consumption, but also increases driving discomfort [START_REF] Rios-Torres | Online optimal control of connected vehicles for efficient traffic flow at merging roads[END_REF]. In the context of connected and automated driving, Cooperative Adaptive Cruise Control (CACC) is posed as a solution to assist the drivers, and ease them from a lot of these issues.

Due to the presence of on-board communication modules that support vehicle-to-vehicle (V2V) and/or vehicleto-infrastructure (V2I) communication, CACC vehicles have the ability to coordinate control actions [START_REF] Shladover | Cooperative Adaptive Cruise Control: Definitions and Operating Concepts[END_REF] with other CACC vehicles. There are many applications of coordinated mobility, including, but not limited to cooperative intersection management [START_REF] Chen | Cooperative intersection management: A survey[END_REF] and merging on highway [START_REF] Rios-Torres | A survey on the coordination of connected and automated vehicles at intersections and merging at highway on-ramps[END_REF]. Nonetheless, introduction of CACC also comes with challenges such as perception and localization errors [START_REF] Patel | Impact of localization errors on automated vehicle control strategies[END_REF], communication delays and packet losses [START_REF] Van Nunen | Robust model predictive cooperative adaptive cruise control subject to V2V impaiements[END_REF], control actuator delays and imperfections [START_REF] Liu | Effects of communication delay on string stability in vehicle platoons[END_REF], etc.

Most algorithms robust to the above issues are developed for a stream of identical CACC vehicles. However, in early deployment phase of future Cooperative Intelligent Transportation Systems (C-ITS), CACC vehicles will need to share roads with other vehicles with little or no automation like manually driven vehicles (MDVs). Unless dedicated lanes are allocated to CACC vehicles, CACC vehicles will need to properly adapt to the behavior of MDVs and vice versa. In order to adapt aptly, controllers of CACC vehicles (on board or remote) will have to make predictions or assumptions on MDVs' behavior. These assumed behavior will not always match the actual behavior of human drivers. In this paper such conflicts are referred to as model mismatch.

Mixed vehicle scenarios have previously been analyzed in [START_REF] Schildbach | A collision avoidance system at intersections using robust model predictive control[END_REF], [START_REF] Jiang | Eco approaching at an isolated signalized intersection under partially connected and automated vehicles environment[END_REF], but model mismatch related issues is yet to receive more attention. To counter model uncertainties, [START_REF] Chen | A robust longitudinal control strategy of platoons under model uncertainties and time delays[END_REF] proposes the use of a robust MPC algorithm. Assumed and actual model of MDVs are represented by IDM+ (IDM+ [START_REF] Schakel | Effects of cooperative adaptive cruise control on traffic flow stability[END_REF] is an extension of Intelligent Driver Model (IDM) [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF]) but with different parameter values to illustrate a model mismatch. Simulations in [START_REF] Chen | A robust longitudinal control strategy of platoons under model uncertainties and time delays[END_REF] do not involve human participation and the performance of the algorithm is yet to be tested in reallife conditions.

Because human drivers have diverse driving behavior under different conditions, it is difficult to represent their behavior mathematically. It is thus necessary to evaluate and validate novel algorithms and innovations by involving real human drivers. At the same time, safety of the human participants also needs to be ensured. Thus, we use a driving simulator to obtain real-time control behavior of human drivers under different driving conditions, and ensure that they are not exposed to the risk of real-life traffic accidents. For instance, a driving simulator was integrated with Plexe (Platooning Extension for Veins) [START_REF] Segata | Plexe: A platooning extension for Veins[END_REF], to study the behavior of a CACC vehicle platoon when a MDV cuts in [START_REF] Aramrattana | Simulation of cut-in by manually driven vehicles in platooning scenarios[END_REF]; however, model mismatch was not considered.

In this paper, we focus on a braking scenario, where a centralized controller intervenes and controls CACC vehicles to avoid collisions. The main contribution is: the evaluation of a centralized controller in presence of model mismatch with real human drivers using the driving simulator from the Swedish National Road and Transport Research Institute (VTI) in a mixed traffic scenario. Moreover, we show the versatility of the controller by simulating braking coordination in a CACC only traffic scenario and compare its performance with two existing CACC controllers.

The remainder of the paper is organized as follows. Section II presents our proposed centralized algorithm. The simulation scenario is described in Section III, with the results presented in Section IV. Section V concludes the paper with suggestions for future work.

II. CENTRALIZED COLLISION AVOIDANCE CONTROLLER

A. System Model

We focus on a mixed traffic scenario consisting of vehicles with different levels of automation. We assume all vehicles to be connected via wireless communication technologies to a centralized controller as illustrated in Fig. 1. All vehicles send their state parameters (position, velocity, acceleration) to the centralized controller over the "uplink". We assume that the centralized controller can control CACC vehicles, which are assumed to have the capability to perform the control action as soon as they receive them. MDVs can not be controlled by the centralized controller. Computations take place on this centralized server, considering vehicles' latest state parameters. Our proposal is to take into account non-automated vehicles such as MDVs, and compute control values for CACC vehicles at each time slot n over a simulation horizon N (n = 1...N ). The computed controls are then transmitted to CACC vehicles in the "downlink".

In this work, the centralized controller is optimized for the braking scenario, where combinations of CACC vehicles and MDVs are driving towards an obstacle. An example of this scenario is described in Section III. The centralized controller computes control inputs for CACC vehicles, such that they avoid front and rear-end collisions.

B. CACC vehicle model

CACC vehicles are controlled by the centralized controller, and are modeled in a way that they avoid both front and rearend collisions. CACC vehicle i (i ∈ Z) ensures distance d i,i-1 between itself (i) and the vehicle in front (i -1) and distance d i+1,i between the vehicle following it (i + 1) and itself are positive. 1 Z is the set of all CACC vehicles.

d i,i-1 > 0 d i+1,i > 0 (1)
where

d i,i-1 (n) = p i (n) -p i-1 (n) -l i-1 n v ≥ i > 1 d i+1,i (n) = p i+1 (n) -p i (n) -l i 1 ≤ i < n v (2)
p i represents the front end position of the vehicle i, and l i represents the length of vehicle i. n v is the total number of vehicles. Note that vehicles i -1 and i + 1 need not be CACC vehicles. We assume perfect localization capability in this work.

C. MDV braking model assumed by centralized controller

At each time slot n in the simulation horizon, when a computation takes place, the centralized controller uses a model to predict braking control values for each MDV over the entire prediction horizon. Different models can be used, we choose a simple but realistic prediction model:

• Before perception response time (Eq. 4 where ζ = c • t i,1 ): the centralized controller assumes that the human driver in vehicle i will start braking after a certain assumed perception response time t i,1 . The driver then increases the braking strength until it reaches a maximum. At maximum braking strength, the vehicle continues to brake until it comes to a halt. • After perception response time:

-If braking magnitude is zero (Eq. 4 where ζ = 0)if vehicle has not started braking, controller assumes that the vehicle will start braking now and continue to increase its braking strength until the vehicle attains maximum braking strength. At maximum braking strength, the vehicle will continue to brake until it comes to a halt. -If braking magnitude is increasing (Eq. 5)the centralized controller assumes that the vehicle will continue to increase its braking strength until the vehicle attains maximum braking strength. At maximum braking strength, the vehicle will continue to brake until it comes to a halt. -If braking magnitude is decreasing or constant (Eq. 6)the model assumes, the vehicle will continue to brake at the previous braking magnitude until halt. As soon as the velocity reaches zero each vehicle stops braking, regardless of its braking strength. We let prediction horizon and the simulation horizon be of equal duration. Let η represent the time slot in the prediction horizon, η = 1...N . The above mentioned braking model can be presented mathematically as:

u i,n =          ψ 1 | ζ=c•ti,1 n ≤ c • t i,1 ψ 1 | ζ=0 u(n) = 0 & n > c • t i,1 ψ 2 ∆u(n) < 0 & n > c • t i,1 ψ 3 ∆u(n) ≥ 0 & n > c • t i,1 ∀ i ∈ Z c
(3) where ψ 1 , ψ 2 , ψ 3 are predicted acceleration values, they are explained in Eq. ( 4), [START_REF] Patel | Impact of localization errors on automated vehicle control strategies[END_REF], and (6) respectively.

ψ 1 =      η • ∆u min i ζ < η ≤ c • t u min i u min i c • t u min i < η < c • t vi=0 0 otherwise (4)
ψ 2 =      u i (n -1) + η • ∆u(n) 1 ≤ η ≤ c • t u min i u min i c • t u min i < η < c • t vi=0 0 otherwise (5)
ψ 3 = u i (n -1) 1 ≤ η ≤ c • t vi=0 0 η > c • t vi=0 (6)
where complementary set of Z which has all MDVs is Z c . For vehicle i at time slot n, u i,n is the predicted acceleration vector (as predicted horizon is equal to simulation horizon, u i,n would have N values). u i (n) is the actual value of acceleration; change in acceleration between two time slots

∆u(n) = u(n) -u(n -1)
is the value of jerk; ∆u min i is the maximum permitted decrease in acceleration between two time slots. t i,1 is the perception response time in seconds. t vi=0 is the time at which velocity reaches zero; t u min i is the time when the acceleration of the vehicle reaches maximum braking magnitude. If n = 1, u i (n -1) = 0; Values in seconds are multiplied with constant c = 10 and converted to time slots (1 second = 10 time slots).

D. Model Mismatch

Model mismatch arises when the actual model is different from the assumed model (described in Section II-C).

In this work, we used two actual models of the MDV, which are derived from either: i) the driving simulator; or ii) a modified version of intelligent driver model (IDM) [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF]. This modified version of IDM involves MDVs implementing controls based on IDM after a particular perception response time. Within the perception response time, acceleration of MDV is zero. This perception response time is added to imitate response delay of a human driver. 2 Model mismatches can only be eliminated if the assumed and the actual model the same. If there is no model mismatch (and no source of error) then the control inputs computed once should be valid and ensure a collision free maneuver. If predicted control values computed at any particular time slot is used in time slots other than the one in which it was computed it might result into collisions. This is because the error introduced by model mismatch would keep accumulating. Thus re-computations of controls is necessary which is the basis of a MPC system.

E. Controller Model

In this work, we propose the use of a centralized controller based on a Model Predictive Control (MPC) system. A MPC controller is able to consider predicted models of different vehicles, different state parameters, various constraints and is able to generate control inputs. In an ideal scenario, at every time slot, fresh state parameters are used to generate fresh control values and are transmitted back to the CACC vehicles. These computations thus take place in a receding horizon method.

The benefits of MPC are as follows: first, at each time slot when control computation takes place, error between predicted model and the actual model can be rectified. Second, MPC based controller generates control inputs over a horizon. In case of computational infeasibility, control inputs generated the previous time slot are to be used. 2 perception response time values are drawn from a normal distribution N (1.33, (0.27) 2 ) and are capped between 0.8 s and 1.8 s. Please refer to [START_REF] Patel | A collision mitigation strategy for intelligent vehicles to compensate for human factors affecting manually driven vehicles[END_REF] for details on computation of perception response time for MDVs.

Next, we formulate a generic version of the above described centralized controller. The state variable x i of a vehicle i (i ∈ 1...n v ) is defined as the position p i , velocity v i tuple in Eq. [START_REF] Liu | Effects of communication delay on string stability in vehicle platoons[END_REF].

x i = [p i v i ] T (7)
The relation between position, velocity and acceleration is given by Eq. ( 8). We assume acceleration to remain constant between two time slots.

∆u i (n + 1) = u i (n + 1) -u i (n) v i (n + 1) = v i (n) + u i (n)∆t p i (n + 1) = p i (n) + v i (n)∆t + 0.5 * u i (n)(∆t) 2 (8)
A discrete time linear control system represented by Eq. ( 9) is used, where values for constants are given by Eq. ( 10),

x i (n + 1) = Ax i (n) + Bu i (n) (9)
A = 1 ∆t 0 1 B = (∆t) 2 /2 ∆t ( 10
)
where ∆t is the time between two consecutive time slots n and n+1. Vehicle and road constraints in terms of minimum and maximum values of position, velocity, acceleration are accounted for in Eq. (11a), and Eq. (11b),

p min i v min i ≤ x i (n) ≤ p max i v max i (11a)
u min i ≤ u i (n) ≤ u max i (11b)
Maximum acceleration and maximum braking is represented by u max i and u min i respectively. To ensure a smooth braking, jerk ∆u(n) is bounded. Note: MDVs are controlled by a human driver, their control values can not be optimized and jerk can not be controlled

∆u min i ≤ ∆u i (n) ≤ ∆u max i (12)
To ensure collision avoidance, the system ensures the distance between vehicles is always positive (Eq. ( 1)). Eq. ( 13) ensures the terminal velocity of CACC vehicles reach zero.

v i (N ) = 0 (13)
Due to the presence of CACCs among MDVs, the behavior of one influences the other. MDV model in this work does not consider the impact of CACC vehicles on MDV. If the actual interdependence can be modeled, the true potential of a centralized controller with CACC vehicles can be evaluated. The cost function (J) is set to penalize strong deviations in acceleration. We can thus set the MPC based system as:

minimize J = nv i=1 N n=2 (u i (n) -u i (n -1)) 2 (14)
subject to Eq. (1), ( 3), ( 7), ( 8), ( 9), [START_REF] Chen | A robust longitudinal control strategy of platoons under model uncertainties and time delays[END_REF], (11a), (11b), ( 12), [START_REF] Segata | Plexe: A platooning extension for Veins[END_REF] The above quadratic optimization problem is to be solved on a centralized controller. We use QUADPROG toolbox on MATLAB to solve Eq. ( 14).

F. Control Buffer

The vector of control inputs over the horizon N generated from Eq. ( 14) is transmitted to CACC vehicles in the downlink. This vector is stored in a buffer. The first value from this buffer is applied by each CACC vehicle.

These computations ideally, should take place with a particular frequency known as the controller update rate. But computations might take longer and controls may not be generated at the desired update rate. Alternately, computations might have terminated but not successful (the optimization problem was infeasible). In such scenario where control inputs for a certain time slot are not computed (due to processing time or infeasibility), we assume CACC vehicles continue to apply next control inputs from the previously received control values stored in the buffer. 3 This concept of the buffer implementation with the control algorithm is explained in Fig. 2.

III. SIMULATION FRAMEWORK

Although a generalized multi-vehicle braking system is formulated, in this work we consider a two vehicle braking scenario where a MDV is following a CACC vehicle. We assume the vehicles are on a single lane road, lane change is prohibited. The CACC vehicle implements control inputs generated by Eq. 14, whereas the MDV is driven by either a human driver, or a modified version of IDM. Following are the parameters used with IDM (parameters follow standard notation as in [START_REF] Treiber | Congested traffic states in empirical observations and microscopic simulations[END_REF]): v 0 , s 0 , T , a, δ, b are 25 m/s, 3 m, 1 s, 1 m/s 2 , 4 and -2 m/s 2 respectively. Length of all vehicles l is set to 4 m. Initial distance between vehicles is 3 m. The frequency of control computation is defined by the controller's update frequency (set to 10 Hz in this paper; ∆t = 0.1s). u min = 5.88 m/s 2 , ∆u min , ∆u max are set to -0.25 and 0.25 respectively. Control and prediction horizon is equal to simulation horizon N . N = 100 unless specified.

Both MDV and CACC vehicles are at halt initially. The first vehicle starts to accelerate with a fixed acceleration of 1 m/s 2 , from 800 meters away from the obstacle. Once the desired speed of 90 km/h (25 m/s) is reached, it cruises until the "notification distance". The notification distance is defined as a distance from the obstacle, where the CACC vehicle detects the obstacle using its sensors, or is notified about the obstacle over V2V/V2I communications. Overview of the scenario is illustrated in Fig. 3. Three different notification distances are used in our study: i) 95.9 meters4 ; ii) 120 meters; and iii) 150 meters. We assume the To involve human drivers in this work, VTI's driving simulation software is used in combination with MATLAB, as illustrated in Fig. 4. All three software blocks in Fig. 4 are running on the same desktop computer. The first MATLAB block, Centralized Controller, is executing the control strategy presented in the Section II. The middle block, Synchronization, is for synchronizing data exchanges between the driving simulator and the centralized controller over a TCP connection. The data exchange happens every 0.1 second, which is also the simulation time step in MATLAB. The last block, VTI's driving simulator, handles the task of fetching inputs from the human driver and displaying positions of the CACC vehicle controlled by MATLAB. The human driver uses a gaming steering wheel and pedals to control the MDV. In the scenario where the MDV was controlled by an actual human, we have six participants (2 women and 4 men) involved in the experiments. For each notification distance, each driver is driving the scenario at least two times, in no particular order. The driver was given an instruction to follow the CACC vehicle in front without making any lane changes.

IV. RESULTS AND ANALYSIS

Results from the simulation scenarios described in Section III are presented in this section. The proposed algorithm is evaluated with the following two metrics: a) percentage of collision avoidance (whether vehicles halted before the obstacle without collisions); and b) discomfort (2 norm of change in acceleration per time slot).

A. Mixed Traffic

The mixed vehicle scenario consists of a MDV following a CACC vehicle, as shown in Fig. 3. We analyze the following model mismatch in the braking scenario under two cases:

• Case A: the assumed driving model is the model introduced in section II-C and the actual driving is based on a modified version of IDM • Case B: the assumed driving model is the model introduced in section II-C and the actual driving inputs are obtained from the driving simulator Results of the above simulation scenarios are plotted in Fig. 5. We observe in A, acceleration and velocity profiles are notably smoother compared to the scenario where the real human is driving the MDV. As one of the evaluation metric is to evaluate discomfort while braking, we compute discomfort values after vehicles are notified at the notification distance. The average discomfort value of Case B is bigger than that of Case A. Moreover, collisions take place when an actual human is driving compared to no collisions when a mathematical model is used. These results are summarized in Table I.

As expected, an increase in notification distance results in higher percentage of collision avoidance. For Case B, the percentage of collisions avoided at 95.9 m, 120 m and 150 m is 35.72%, 57.89% and 60.87% respectively. Overall, 53.57% of collisions were avoided. Whereas in case A, when IDM controls MDV, 100% collisions were avoided. There are several reasons for collisions, due to the differences between assumed value in the centralized controller and the actual value in the driving simulator, e.g. the difference in perception response time, the difference in braking capacity and maximum value of jerk sustainability. In such cases, control optimization computations are infeasible for the assumed values of jerk and the braking capacity. But it may actually be feasible for MDV to brake using the driving simulator.

Furthermore, we observe that the optimization computations takes more time than expected and thus control optimization is not computed at the controller frequency. Computations ideally should take less than 0.1 second to keep real-time performance. However, in our experiments, it takes maximum up to 2 seconds (20 times more) to compute. On an average, we have an optimization computation completed every 0.3 seconds. The previously computed control inputs in the buffer are old and may be from a few seconds ago. These control inputs from the buffer are sometimes not able to avoid collisions.

B. Automated vehicles traffic

Although the proposed controller is designed for mixed vehicle scenarios, it can perform in a CACC-only scenario as well. In this section, the scenario illustrated in Fig. 3, where both vehicles are CACC vehicles, is simulated. To evaluate the performance of our proposed controller, we compare its performance with existing CACC controllers in Plexe such as the one proposed by "Rajamani" [17, Chapter 7], and "Ploeg" [START_REF] Ploeg | Design and experimental evaluation of cooperative adaptive cruise control[END_REF].

Simulation results show that all controllers are able to avoid the collision with the standstill obstacle and come to a halt. Figure 6 plots the acceleration profile of different controllers, when the notification distance is 120 meters. From the plot, we observe that Rajamani controller tends to apply brakes harder than the Ploeg controller. Table II summarizes the discomfort values. The discomfort value of the proposed controller is lower (better) than that of the Rajamani controller and is similar to the value of the Ploeg controller. CACC vehicles based on the Ploeg controller take much longer to come to a halt compared to the proposed controller, with approximately the same cost (discomfort). Acceleration profiles can be influenced by changing the simulation horizon N (refer Fig. 6), and the value of discomfort can be reduced (Table II).

V. CONCLUSIONS AND FUTURE WORK

Difference in the predicted and the actual behaviour of manually driven vehicles (MDVs) gives rise to a model mismatch. In this work, we evaluate the impact of model mismatch on a centralized model predictive control algorithm for braking in mixed traffic. In order to do the evaluation, we propose a simulation framework, which interfaces MATLAB with the VTI's driving simulator.

To simulate a mixed traffic scenario with a CACC and a MDV, first, a mathematical model is used to control the MDV. Second, the controller is interfaced with a driving simulator and a human driver controls the MDV. In the second case, 46% more collisions were recorded and compared to the first case. Moreover, the discomfort of the driver increases by 91%. Based on the significant difference in the results, we conclude that it is necessary to use humans for experimental validation of algorithms. Furthermore, the proposed controller is also compared against Rajamani and Ploeg controllers in a CACC-only scenario. Each controller is able to avoid collisions although the discomfort values are different. The proposed controller thus adapts to both homogeneous and heterogeneous traffic.

An indepth study of the impact of model mismatch can be carried out by using differnet assumed and actual models. As future work, tools like the one proposed in [START_REF] Aramrattana | A simulation framework for cooperative intelligent transport systems testing and evaluation[END_REF] can be used to further add other factors such as communication disturbances and extend the scenario to a more diverse mixed traffic scenario. Also, simulations with more than 2 vehicles shall be carried out in the near future. Moreover, the algorithm needs to be made more robust to failures in mixed traffic scenario. Last but not least, to ensure real time simulation of more than two vehicles, the algorithm needs to be improved and optimized.

Fig. 1 :

 1 Fig. 1: System model: CACC vehicles have both uplink and downlink; MDVs only have uplink.

Fig. 2 :

 2 Fig. 2: Flow chart of the control buffer implementation.Initial value of count is one.

  The CACC vehicle starts 800 meters away from the obstacle at zero speed. MDV obstacle CACC notification distance (b) The CACC vehicle has reached notification distance.MDV obstacle CACC (c) After the notification distance, the CACC vehicle starts braking to stop before the obstacle.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Overview of the simulated scenario.

Fig. 5 :

 5 Fig. 5: Mixed vehicle scenario results: top plots show results from case A; bottom plots show results from case B.

Fig. 6 :

 6 Fig. 6: Acceleration profile of the first vehicle in the automated vehicles traffic scenario, at the notification distance of 120 meters.

TABLE I :

 I Mixed traffic simulation results

		Case A	Case B
	Discomfort	6.66	12.76
	Collisions avoided (%)	100	53.57

TABLE II :

 II Average value of discomfort (2-norm of change in acceleration computed at 10Hz).

		Notification distance (m)
		95.9	120	150
	Rajamani [17]	6.18	6.18	6.18
	Ploeg [18]	1.36	0.90	0.80
	Proposed (N = 100)	1.25	1.15	1.15
	Proposed (N = 150)	1.24	0.99	0.85

vehicle i -1 is followed by i which is followed by i + 1, and so on. vehicles are moving towards the origin where obstacle is located.

If the control computation is infeasible and the control buffer is empty, CACC vehicles applies brakes with maximum jerk added to previous acceleration value.

This is the distance at which at least one DSRC/ITS-G5 safety message would be received with 99.5% probability[START_REF] An | VANET: Is 95% probability of packet reception safe?[END_REF]

ACKNOWLEDGEMENT Raj Haresh Patel is a recipient of a PhD Grant from the Graduate School of the University Pierre Marie Curie (UPMC), Paris. EURECOM acknowledges the support of its industrial members, namely BMW Group, IABG, Monaco Telecom, Orange, SAP and Symantec.