A study on the effects of an explosion in the Pantheon of Rome

The response of an emblematic monumental structure, the Pantheon in Rome, to an internal blast is addressed. The analysis is a coupled solid-fluid numerical simulation done using JWL equation for the simulation of the blast event and considering the pre-existent cracks in the dome as well as the material nonlinearities of low-tensile strength concrete aggregates. We identify two main phenomena; a focalization of shock waves inside dome-vaulted like structures due to their tendency to concentrate the blast energy and the role played by the pre-existent cracks in the evolution of the structural damage


A study on the effects of an explosion in the Pantheon of Rome
Research on the effects of an explosion on monumental architectures are hence interesting for assessing the potential effects of a blast onto a monument structure and also for helping in the design of reinforcements or any other possible passive protection device.
This domain is still almost unexplored; in fact, the most part of papers concerning the effects of an explosion on a civil structure regard modern reinforced concrete or steel D R A F T structures with simple geometries, normally squared buildings, [Remennikov, 2003], [Ngo et al., 2007], [Koccaz et al., 2008], [Draganic and Sigmund, 2012].
This is not the case of monumental structures that have often a complex geometry, sometimes very articulated, which renders the assessment of the blast loads strongly case dependent and affects the type of simulation. This point is tackled in Sect. 3, where a short account of the state of the art for what concerns the simulation of blast loads is given and, on its base, the choice of the method used for the case of Pantheon is justified.
On the other hand, monumental structures are either masonry-like or timber structures, or both of them at the same time. In particular, monuments constituted by masonry-like materials have a structural response strongly affected by the no-tension behavior of the material. This is a key point, considered in Sect. 4.
A particular attention must hence be paid to the procedure to be used for numerical simulations, that must account, on one hand, for the hypervelocity of the phenomenon and, on the other hand, for the peculiar constitutive law of the material, that must be able to describe the non-linear phenomenon of damage, i.e. of cracks propagating into the body of the structure as a consequence of the blast actions. Considerations about the procedure used for the numerical simulations are given in Sect. 5 and the results are presented in Sect. 6. But let us start by briefly introducing, in Sect. 2, the object of this study, the Pantheon of Rome.

The Pantheon of Rome
The Pantheon of Rome, one of the most admired and studied monuments ever, was built upon the rests of previous temples by the Emperor Hadrian, since A.D. 118 to about 128, or later, perhaps until 140, under Emperor Antoninus Pius, and probably it is the joint work of Hadrian and of the Nabatean great architect Apollodorus of Damascus (see [Marder and Jones, 2015]).
The main body of the Pantheon, the so-called Rotonda, is composed by a cylinder whose inscribed sphere is coincident, for its upper part, with the dome, while its bottom touches the ground, Fig. 2. The coffering is sculpted in the intrados of the dome offering a high aesthetic value, but also reducing the dome's weight.
With a diameter of 43.30 m, according to the measurements of de Fine Licht, [K. de Fine Licht, 1968], [Mark, 1990], or of 150 Roman feet, i.e. 44.55 m, according to Wilson-Jones, [Wilson-Jones, 2000], [Como, 2013], the Pantheon's dome is, still today, the largest dome in the world, apart the modern realizations in reinforced concrete. In fact, it is larger than the dome of Saint Peter in Vatican, whose diameter is 42 m, and also of the octagonal dome of F. Brunelleschi, in Santa Maria del Fiore at Florence, whose base is circumscribed to a circle of 41.57 m diameter, [Cowan, 1977]. However, unlike these two famous domes, and also of other ones, made by bricks, the Pantheon's one is made of concrete, a technique already mastered by Apollodorus in other previous works (e.g. the vault of the Great Hall in the Trajan's markets in Rome, see [Perucchio and Brune, 2008], [Perucchio and Brune, 2009]). D R A F T Figure 1: The Pantheon of Rome (from [Pulvirenti, 2014]).  [Pulvirenti, 2014]).
If the intrados of the dome is sculpted by a coffering, the external lowest part of the dome is modeled by stepped rings, whose function has been the object of different investigations, see Fig. 2. In an early application of finite element analysis to Roman structures, Mark and Hutchinson [Mark and Hutchinson, 1986], [Mark, 1990] used a simplified two-dimensional model of the Rotonda, exploiting the axial symmetry of the structure. According to this study, and contrary to what commonly thought, the step-rings do not contrast the hoop tension in the lowest part of the dome, because actually the dome is spread of cracks that have an almost meridional direction and stop at a latitude of approximately 57 • , where the hoop stress begins to be compressive.
The distribution of these cracks, Fig. 3, was detailed in 1934 by A. Terenzio, the Superintendent of the Monuments of Latium, who had carried on a series of inspections on the dome of the Rotonda after that some fragments had fallen down, [Terenzio, 1934].

D R A F T
We cannot establish here what is the true origin of these cracks, but certainly they are to be modeled to have a more realistic response of the structure to gravity and blast loads.
In fact, on one hand, the cracks undoubtedly change the stress regime of the Rotonda, on the other hand, they considerably affect the response of the structure to a blast by a special local mechanism, as we will see in Sect. 6.2.  Mark and Hutchinson 1986 pg. 31. 29 It may also be rooted in the ambition to make sweeping, generalized statements; e.g., "Roman structural development was nowhere near so radical as that of the late nineteenth century, when the introduction of the new industrial materials brought forth a true revolution in building design." And yet, somehow, "A far better analogy is provided by the late twelfth-century development of the High Gothic cathedral…" (Ibid. pg. 33).

Modeling the blast actions
An explosion is an extremely rapid and exothermal chemical reaction, which lasts just few milliseconds. During detonation, hot gases, pertaining to the chemical process, expand quickly and, for the hot temperatures produced almost instantaneously, the air around the blast expands. 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. We describe below this phenomenon, more details can be found in [Vannucci et al., 2017a].
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 s : overpressure due to the blast; it is the pressure in the air relative to P o .
• P r : reflected overpressure, i.e., the pressure, relative to P o , acting at a point q of a solid surface when hit orthogonally by a shock-wave.
• t A : arrival time, i.e., the instant at which the shock-wave peak arrives at q.
• t o : positive phase duration, defined below.
• t o− : negative phase duration, defined below. Figure 4 represents an ideal profile of the overpressure P s (q, t) produced by a blast. When the shock wave arrives at q, after t A , the pressure instantaneously increases from the ambient pressure P o to a peak for P s : a strong discontinuity, indeed.

D R A F T
For t > t A the overpressure decreases exponentially until time t A +t o , when P s = P o , which marks the end of the so-called positive phase. After t A + t o , we have the negative phase: the pressure decreases with respect to P o and comes back to P o after t o− > t o . During the negative phase the decrement of the pressure is in modulus much lower than the peak pressure of the positive phase. Consequently, the negative phase can be neglected for structural analyses, though it can be important in some special cases, due to its duration, always much longer than the positive phase. Such a behavior is idealized: perturbations can occur, due to different circumstances.
The shock wave is the main mechanical effect of a blast on a structure, but not the only one: hot gases, expanding, produce the so-called dynamic pressure, least in value with respect to the shock wave and propagating at a lesser speed, while the impinging shock wave can be reflected by solid surfaces and acts again on other surfaces as reflected shock wave.
The overpressure P s at a point q decreases with increasing time t and distance R. Generally, the time rate decrement is much greater than the space rate one: the blast overpressure is like a localized pressure wave propagating at high speed and decreasing intensity in the distance.
P r is the pressure acting on a surface impinged by the incident overpressure P s . The peak of P r is normally much greater than the one of P s , measured at the same point and assuming the absence of any surface.
The simulation of a blast can be conducted by using different approaches. Here we refer to three phenomenological approaches: JWL, CONWEP, and TM5-1300 models. JWL stands for Jones, Wilkins and Lee, the authors of this model, [Jones and Miller, 1948], [Wilkins, 1964], [Lee et al., 1968]. The JWL model is physically based: the laws of thermodynamics are used to recover a phenomenological description of a chemical blast. It allows a detailed description of a blast phenomenon, including the propagation of the shock-wave in a medium, e.g. air, its reflection on solid surfaces and the expansion of the hot gases, i.e. the dynamic pressure.
The JWL model gives the overpressure P s as function of different parameters:

D R A F T
A, B, R 1 , R 2 and ω are parameters depending upon the explosive type, along with ρ 0 , its density, while ρ is the density of the detonation products and E m is the internal energy per unit mass. The detonation velocity v D and the Chapman-Jouguet pressure p cj need to be specified too. All their values are selected to fit experimental results on the cylinder expansion test (see Table 1).  [Lee et al., 1973]. Though the detailed modeling of a blast, JWL model needs of discretizing, finely, the charge and the fluid domain, that can be very large, besides the structure if a coupled structural analysis is to be done.
Multi-physics transient problems, with a strong fluid-structure coupling, lead to numerical simulations that can be heavy. Empirical methods are often used for their reduced computational costs. The two most commonly used empirical models rest upon different but related studies of the U.S. Army Corps of Engineers (USACE): the document [USACE, 1986] and the Technical Manual TM5-1300 [USACE, 1990]. They contain the model CONWEP, completed by successive documents [USACE, 2008]. The Joint Research Center of the European Union has produced in 2013 a Technical Report [Karlos and Solomos, 2013] substantially referring to these two last USACE documents and to another Technical Report of the U.S. Army [Kingery and Bulmash, 1984].
Such empirical models are less precise in predictions than JWL, especially because they neglect reflected waves, which may have, in contrast, prominent effects, especially for internal blasts. Depending upon geometry, the concentration of the reflected waves can give rise to local effects that can be greater than the original shock wave. In the case of vaulted structures limited laterally by walls, a localized shock wave produced by the reflected waves can hit the vault with an overpressure far greater than that produced directly by the original impinging shock wave [Vannucci et al., 2017a], [Masi, 2017].
The case of the Rotonda is hence particularly interesting: by its geometry, focalizing effects can happen and be determinant. That is why we have chosen to study the problem of a blast inside the Pantheon using the JWL model, that allows for modeling completely the explosion. This needs the discretization not only of the structure, to perform the coupled structural analysis, but also of the air volume inside and partly outside the Pantheon. This is a step forward with respect to a previous study of the authors, where the CONWEP model was used, [Vannucci et al., 2017b]. Finally, we have a very huge numerical model, needing a considerable computational effort; this aspect is detailed in Sect. 5.

Modeling the material behavior
The Pantheon is made by different materials: bricks, mortar and concrete, for the cylindrical part of the Rotonda, granite for the columns inside and concrete for the dome. The D R A F T cylindrical part, whose thickness varies from ∼ 2.4 m to ∼ 6 m, is actually composed also by pillars and arches, having a structural role and merged in the wall. Concerning the dome, it is more correct to say that it is composed of concretes, because the Romans used different types of concrete, from the heaviest one in the lower part, whose thickness is ∼ 5.9 m, to the lightest one in the upper part of the dome, where the thickness decreases to ∼ 1.5 m, see Fig. 5.

Characterisation of the materials
In this study, we have used the same mechanical model for the constitutive law of the concretes composing the dome or the foundations, the masonry of the cylindrical part and the granite of the columns. We precise that we have modeled the masonry as an equivalent concrete, as it is almost impossible, for such a large structure, to model precisely the actual disposition of the bricks composing the pillars and numerous arches included in the masonry of the Rotonda, see Fig. 5.
The materials are assumed to be linearly elastic in compression. It is commonly considered that in monumental structures compression stresses are very low, usually an order of magnitude below the compressive strength of the material, cf. [Heyman, 1995, Stefanou et al., 2015, so that there is no possibility of crushing. However, in the present study where the main actions are produced by a blast, this assumption may be questionable: high compressive stresses produced by the shock wave cannot be excluded a priori in some parts of the structure. That is why in our approach the consistency of the subsequent calculations is checked by monitoring that the compressive stresses do not exceed the compressive strength, assumed to be, for all the concretes, equal to ∼ 5 MPa, according to the experimental results by [Samuelli-Ferretti, 1996] and [Giavarini et al., 2006] and the models proposed in  and [Ivancic et al., 2014]. The verification is given in the Appendix.

D R A F T
Masonry and concrete show a low strength to tensile stresses, so small that often they are modeled as no-tension materials. In our study, a small, but not null, tensile strength f t is considered for the material. When f t is exceeded, damage occurs and the following tensile softening is characterised, as proposed in [Hillerborg et al., 1976], in terms of the fracture energy G f , i.e., the energy required to open a unit area of crack in Mode I.
Different material characteristics must be determined in order to adequately define their mechanical response: the density ρ, the Young's modulus E, the Poisson's ratio ν, the tensile strength f t , the fracture energy G f and the softening behavior.
We have considered five zones for the materials of the Rotonda, corresponding to those indicated in Fig. 5: foundations, brick-faced concrete of the cylindrical wall, lower, intermediate and upper zone of the dome. In addition, we have considered apart the material of the columns in the interior of the Rotonda, made of granite.
The concrete densities have been fixed upon the indications given in [Mark and Hutchinson, 1986]. About the other quantities, the choice is much more problematic. Though we know today rather well the behavior of light concretes, according to historians and researchers in mechanics of ancient monuments, the physical properties of concretes and mortars of the past have quite different characteristics. Moreover, the data that can be found in the literature about ancient concretes are often quite different and most of all fragmentary. The only reliable data about the Pantheon were the densities. Therefore we needed a rationale linking in some way the other material parameters to densities in such a way that, starting from few known data, the other ones could be determined.
For all the concretes, we have put ν = 0.2, a mean value often used for the concrete Poisson's coefficient in ancient monuments, see e.g. [Perucchio and Brune, 2009].
About the other material parameters, the question was more delicate. First of all, it is worth noting that experimental tests show that the tensile behavior of Roman concrete follows that of mortar (see [Samuelli-Ferretti, 1997, Giavarini et al., 2006, Jackson et al., 2009). It has to be noted that in experiments the cracks propagate either in the mortar volume itself, in the aggregates, or at the interfacial adhesive zones . Therefore, it seems that the contrast of the tensile strength between the Roman concrete's aggregates and the mortar it self is rather low [Jackson et al., 2009], giving a less heterogeneous material that behaves more as a mortar and less as a concrete, at least in the traction regime, which is of interest here. Being aware of this, Brune, in his PhD thesis , conducted experimetal tests on a reproduction of ancient mortar with different uring time, see Fig. 6. Brune introduces a bilinear approximation of the experimental softening curve as shown in Fig. 7. Following his assumption and referring to Fig. 8, we can compute the fracture energy for normal tensile stresses as where σ is the maximum principal stress; w is the displacement normal to the crack surface; w k is the normal displacement relative to the kink point, with stress σ k ; and w f the one corresponding to the complete loss of strength. Applying the bilinear approximation D R A F T Figure 6: Average characteristic response for different curing times, according to .
Figure 7: Best-fit traction-separation descriptions for each curing duration, according to .
of  in terms of the following parameters, the two integrals in eq. (2) become so that This expression allows to compute w f if G f , f t , Ψ and ζ are known:   ) . Regarding the Young's modulus, after comparing the value of E given by Brune with other tests, see e.g. [Jackson et al., 2009], [Jackson et al., 2014], , [Samuelli-Ferretti, 1996], it is likely to consider the value E = 3.37 GPa, more or less, as that of a concrete whose density is 1520 kg/m 3 .
Another value that can be extrapolated from the literature, is that concerning E for the cylindrical wall; it is composed by concrete and masonry and we have found a reliable value for ancient masonry in [Como, 2013]: E = 6.6 GPa. Consequently, rounding the above values, we have put E = 3 GPa for concrete 5, the lightest one, and E = 7 GPa for concrete 2. The other values can be determined upon a rule stating the dependence of E on ρ. We did not find it in the literature about ancient concretes and the scarcity of data about these materials did not allow us to extrapolate such a rule.
That is why we have considered modern technical recommandations about lightweight concretes, [ACI, 1999], [FIB, 2000]. The following law is proposed in [Sanpaolesi and Formichi, 2009] about the Young's modulus: D R A F T with E in GPa, ρ in kg/m 3 and f c , the mean compressive strength, in MPa. What is apparent from this relation is that E is a quadratic function of ρ. We have hence looked for a law representing E satisfying this rule, giving a minimum of E for ρ = 1350 kg/m 3 and the two values given above for concrete 5 and concrete 2. Some simple passages give the relation E = 2.5 × 10 −5 ρ 2 − 0.0675ρ + 48.5625, with E in GPa and ρ in kg/m 3 . This allows to have an estimation of E also for the other concretes.
About the tension strength f t , the only reliable datum is that of Brune, see Tab. 2: f t = 0.55 MPa. We attribute this value to concrete 5, the lightest one. Like for E, we need hence to put in relation f t with ρ. We still use a relationship proposed for modern lightweight concrete, [Sanpaolesi and Formichi, 2009]: where f ck is the characteristic value of the statistical distribution of f c . It is likely that the compressive strength of the concretes in Tab. 3 do not change considerably. Then, making the assumption that f ck is practically constant, we get the following relation between the values of f t for two different concretes: Then, putting f t1 = 0.55 MPa and ρ 1 = 1350 kg/m 3 , we get the following linear relation for f t as function of ρ: f t = 3.25 × 10 −4 (880 + 0.6ρ), with f t in MPa and ρ in kg/m 3 . This relation allows us to have an estimation of f t also for the other concretes of the Rotonda.
Concerning the fracture energy G f , once more we have attributed the value of 55 J/m 2 of Tab. 2 to concrete 5. The data concerning G f in the literature are even more questionable than those of the other material parameters: they are considerably affected by the experimental or numerical procedure used for its evaluation besides the intrinsic properties of the material, like the size of the aggregates, the curing process among others. As a consequence, the spread of data is considerable, see for instance [Hillerborg, 1985], [Weerheijm and Vegt, 2010]. Hence, once more we have looked for a relation between G f and another known material parameter of concrete.
In [Dehn, 2004] the author gives a linear relation between G f and the tensile strength f t for lightweight concretes with natural sand, which is the case of the concretes of the Rotonda: with G f in J/m 2 and f t in MPa. Assuming a linear relation between G f and f t is a strong assumption, however is justified by [Dehn, 2004] and it is frequently the case in D R A F T geomaterials such mortars and masonry structures [van der Pluijm, 1999]. Therefore, for two different concretes, we get and putting, as indicated above, G f 1 = 55 J/m 2 and f t1 = 0.55 MPa, data relative to concrete 5 in Tab. 3, we get the linear relation allowing us to obtain an estimation of G f for all the other concretes.
For what concerns the ratios Ψ and ζ we have put Ψ = 0.35 and ζ = 0.15, as suggested by Brune, see Tab. 2.
For the granite of the columns, we have taken the values found in [Buyukozturk, 1993] and we have put Ψ = 0.25 and ζ = 0.50.
Using eqs. (3) and (6), we get also all the other parameters defining the tensile softening represented in Fig. 8. All the material parameters so found are shown in Tab. 3. In Fig. 9 we show the diagrams of the constitutive law σ − w. The diagram of the granite has been reduced by a factor 4, for graphical reasons.
It is worth noting that the values of the different material parameters are close to experimental data for ancient Roman concretes, see [Lamprecht, 1984], [Jackson et al., 2009], ,  and [Jackson et al., 2014]. They represent hence a plausible variation of such parameters with the concrete density.

The brittle cracking model
The nonlinear brittle cracking model proposed in [Hillerborg et al., 1976] is implemented in Abaqus using a smeared cracked approach, [ABAQUS, 2016]. According to the model, when the maximum principal stress at a Gauss point reaches the material tensile strength f t , a fictitious crack is formed in the plane orthogonal to the direction of the principal tensile stress. The orientation of the planar Mode I crack in the three-dimensional space remains constant. This is, of course, a modeling assumption but a reasonable one for the material and loads considered. Once the peak tensile stress is reached, a linear softening branch is followed based on the fracture energy G f . This allows to minimize mesh dependency due to the softening behavior and to dissipate adequately the energy.

The numerical procedure
The use of the model JWL for the blast simulation needs, as already mentioned, a meshing not only of the Pantheon's structure, but also of the air volume interested by the explosion, in this case the internal volume of the Rotonda and also a volume around the two openings, the entrance door and the oculus, the round opening at the dome's top, with a neat diameter of ∼ 7 m, [K. de Fine Licht, 1968]. The details on the models of the air volume and of the structure and the procedure used for their numerical validation are given below, separately for the two parts. Using the symmetry of the structure, see below Fig. 14, we have modeled just one half of the structure and of the air domain. This allows, for the same computing effort, to obtain a more detailed model, i.e. a finer discretization.

Validation of the fluid domain finite volumes model
The model of the air domain is shown in Fig. 10. The air domain of half the structure has been decomposed into 2.068 × 10 6 8-nodes hexahedral volume elements, for 10.34 × 10 6 degrees of freedom (DOF) on the whole (5 DOF for each volume element: the 3 displacement components and the pressure of the volume element centroid plus the the Eulerian Volume Fraction (EVF), [ABAQUS, 2016]).
In Coupled-Eulerian-Lagrangian (CEL) simulations, ABAQUS/Explicit takes into account for the Eulerian fluid domain through the so-called volume-of-fluid method: the material flowing through the mesh is tracked by the definition of an additional variable within each element, the EVF. In our case, this allows to compute not only the propagation of the shock waves, but also the diffusion of the explosive material inside the air domain. The coupling between the Lagrangian (solid) domain and the Eulerian (fluid) one, i.e. the fluid-structure interaction, is achieved by a general contact algorithm, with a null interface friction coefficient and a penalty method, [ABAQUS, 2016].

D R A F T
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.
We have modeled an Eulerian domain as deep as 5 meters above the oculus and a volume of 2 meters of air on the lateral sides (see Fig. 10). The choice of the fineness of the mesh has been done upon a convergence analysis: four different meshes for the air domain have been considered: M1, with 1.48 × 10 6 DOF, M2 with 6.3 × 10 6 DOF, M3, the chosen mesh detailed above, and a reference mesh M4, the finest one, with 16.94 × 10 6 DOF.
For each mesh, we have computed the blast overpressure due to a TNT equivalent charge of 30 kg in correspondence of three different observation points: P1, P2 and P3, indicated in Fig. 10, respectively close to the base of the cylindrical wall, to its top and to the oculus. The time histories of the overpressure, for the three different points and for each mesh, are plotted in Fig. 11. In Fig. 12 we show, for points P1 to P3, the relative error evaluated for meshes M1 to M3 with respect to mesh M4; such an error, in percent, is computed as where p M i is the peak of the overpressure, in P1, P2 or P3, evaluated for mesh Mi, while p M 4 is the same for the reference mesh M4.
Considering the results shown in Fig. 11 and 12, we can see that the choice of the mesh M3 is a fair compromise: it guarantees a good quality of the result (the maximum relative error is 3.8%, for point P3) with a problem size that is not extreme.
We have also compared the results given by JWL with those given by CONWEP for the reflected pressure P r on the solid points corresponding to P1, P2 and P3, see Fig. 13. The comparison clearly shows how much the use of CONWEP in a case like this one can lead to erroneous results: the time history of the overpressure is completely different, mainly due to the reflected shock waves.

Validation of the Pantheon finite element model
The entire structure of the Pantheon has been modeled using CAD programs, see Fig. 14.
The model comprehends the coffering of the dome intrados, the step-rings, the interior columns, absides and cavities present in the cylindrical wall and the foundation ring.
The CAD model has been successively meshed to obtain a finite element (FE) model. The FE model comprehends only the Rotonda, because the structural role of the pronaos is undoubtedly little, for the blast actions of an explosion in the interior of the Pantheon; so, we have neglected it.
Just like for the mesh of the air domain, we have performed a convergence analysis for validating the choice of the Rotonda mesh, composed of four-nodes tetrahedral elements, supported by both the standard and explicit solvers of ABAQUS. Also in this case, the fineness of the mesh has been chosen after a convergence analysis, performed in order to obtain a reliable degree of accuracy with an acceptable size of the whole numerical model (air domain plus structure).
The convergence analysis has been done using the standard (implicit) scheme for eleven different structural meshes MSi, i=1,...,11, whose characteristics are given in Tab. 4. The standard analysis made for each mesh MSi is divided into two parts: a modal analysis and a static one, for the only gravity load. We have monitored the eigenfrequencies of the first twenty vibration modes of the structure and, for the static analysis, the vertical displacement of point P3 in Fig. 11. The twenty eigenfrequencies f j , j = 1, 2, . . . , 20, and the vertical displacement u of the point P3, are calculated for each structural mesh MS. The convergence has been evaluated like for the air domain, calculating for each mesh MSi=1,...,10 the errors ∆f j of the frequencies f j and ∆u of the displacement u, relatively to the same quantities f r j and u r , calculated for the reference mesh MS11, the D R A F T finest one: In Fig. 15 we show the diagrams of ∆f 1 , ∆f 10 and ∆f 20 along with ∆u as functions of the DOF number.  Figure 16 shows the relative errors ∆f 1 , ∆f 10 , ∆f 20 and ∆u versus the normalized CPU time necessary for th static and modal analyses; for each mesh MSi, i=1,..., 11, this is the ratio between the CPU time corresponding to MSi and that of MS1. After evaluation of the results in Fig. 15 and 16, we have selected the mesh MS11, the finest one, for the computations. In fact, we have estimated that the convergence has been reached (the diagrams in Fig. 15 and 16 show an asymptotic behavior) while the increase of the computation time can still be considered as tolerable. A detail of the structural mesh MS11 is shown in Fig. 14.

D R A F T
To take into account for the existing cracks in the dome, Fig. 3, we have modeled 14 cracks, propagating along the meridians up to 57 • above the plane passing through the sphere center in Fig. 2; they are indicated in Fig. 14 c). Each crack is modeled as an interface, with a hard normal contact behavior (no penetration between interacting parts of the model) and a Coulomb friction coefficient of 0.5 (cf. [Brune and Perucchio, 2012]). A penalty algorithm was used [ABAQUS, 2016] for the tangential components of the contact stress.
Parallel computing is used to reduce significantly the duration of the analyses. All the simulations have been performed using a 24-cores workstation; the entire model is, thus, divided into 24 geometric domains, taking advantage of the 24 processors of the machine.
The energy balance of the coupled fluid-solid numerical simulation is verified by monitoring the total energy, E total , of the whole model defined as (the superscripts E and L refer to the Eulerian and Lagrangian domains, respectively) where E v is the viscous energy dissipated by bulk viscosity; E ke is the kinetic energy; E w is the work done by gravity load and E i is the internal energy, which is given by the sum of elastic strain energy E e , artificial strain energy E a (that is the energy stored in hourglass resistances) and energy dissipated through damage E dmd , namely for the Eulerian and Langrangian domains, respectively. The consistency of the numerical simulation is checked by computing the variation of the total energy with respect to its initial value (E total | t=0 , when t = 0) during the mechanical calculations, i.e.
The total energy is found to be constant within an error ∆E total < 1 %, that is acceptable for an explicit scheme.

Numerical simulations
We have simulated the effects of an explosion in the center of the Rotonda. It is worth noting that the quantity of explosive has been selected by searching the minimum quantity of TNT able to cause the partial collapse of the structure. Figures 17 to 19 represent the pressure field at different (normalized) times in the interior of the Rotonda and, starting from the moment where the first cracks appear in the dome, the effects on the structure, where the propagation of cracks is clearly visible. We comment separately the pressure field and the effects on the structure. D R A F T

Evolution of the pressure field
The pressure field represented in Fig. 17 is comprehensive of all the blast phenomena: shock wave, reflected shock waves and dynamic pressure. The first two sketches in Fig. 17 show the very initial phases of the blast, with the shock wave propagating hemispherically until it touches the cylindrical wall. At this moment, the reflected shock waves add to the principal one, creating a complex pressure field and, though the principal shock wave decreases in intensity, important local concentrations of pressure are possible due to the interaction between the principal shock wave and the reflected ones.
Local concentration of the pressure are clearly visible in the interior of the niches and in the parts of the coffering looking downward, but the most important concentration happens exactly in the central axis of the Rotonda, of course as an effect of its cylindrical symmetry. An important zone of high pressure, in the form of a butterfly, forms and progresses towards the oculus. This last represents an important escape way for the pressure and the gases, as clearly visible in the last pictures of Fig. 17. Also the entrance door is an escape way, but apparently it has a least effect of the oculus, because the focalized shock waves move toward the dome's top.
D R A F T Figure 17: The evolution of the (normalized) overpressure field inside the Rotonda. Dark blue regions represent negative or zero overpressure.
Such pressure concentrations can lead to local value of the overpressure far greater than that produced by the original shock wave. This is clearly visible, e.g., in Fig. 11, where D R A F T the pressure for point P3, close to the oculus, clearly shows secondary peaks far higher than the first one, due to the impinging shock wave. These secondary peaks mark the passage, through the oculus, of the "butterfly" pressure wave. Another important fact, still visible in Fig. 11, is the fact that the overpressure lasts much more than the single original shock wave. Of course, this has a strong effect on the damage of the structure, because it is much longer exposed to high impact pressures.
Rather interesting is the fact that the coffers act like dampers of the pressure impulse. The sharp-cornered surfaces of the coffers are nothing but obstacles to the propagation of the internal blast wave: the waves reflected by the coffering array looking downward contrast the rising incident wave through the particular combination of incident and reflection angles, as clearly shown in Fig. 17. The arrangement of the coffers themselves results in a step-by-step attenuation, mitigating the incident pressure as getting closer to the upper part of the dome and the intense, butterfly-shaped, concentration that rises in the center of the Rotonda.

Evolution of the structural damage
It is exactly the "butterfly" zone of focalized pressure, result of the interaction of the principal and reflected shock waves, that produces the principal damages to the Rotonda, and in particular it is likely to produce the destruction of the dome. This is apparent looking at Fig. 18: the first cracks appear in the compression ring of the oculus when the focalized pressure approaches the dome and when it passes through the oculus, the existing cracks begin to propagate toward the dome's top and eventually they reach it. The static regime of the dome is completely changed, because a shell behavior is no more possible and, each part of the dome behaving like a wedged cantilever, a circular crack appears in the upper part of the dome as a consequence of the bending tensile stresses produced by gravity so that, in the end, the upper part of the dome collapses.
In the same time, in the the cylindrical wall the existing cracks propagate through the entire wall thickness, as depicted in Figs. 18 and 19. Moreover it is shown that the granite columns are fractured by the blast.
It is worth mentioning that tension is the main damage mechanism. The developped compressive stresses are much lower than the typical compressive strength of the material involved, even under blast actions. This is discussed further in the Appendix.
The failure and collapse of the dome is hence the consequence of different causes: the blast, the propagation of the existing cracks and the selfweight of the dome. The real extent of the existing cracks in the Rotonda's dome is uncertain. Therefore, for the sake of completeness, we performed the same simulation described above and on exactly the same model of the Pantheon without the existing cracks. All the other parameters were left unchanged.
The evolution of the structural damage of the structure in this case is depicted in Fig.s  20 and 21. It is evident that the case without cracks has a different failure mechanism, i.e. the upper part of the dome is damaged, but to a smaller extent with respect to the case with cracks. This is because in the case of the cracks, as mentioned above, the blast effect is a propagation of these toward the top of the dome and the creation of the cantilever effect, that destroy a large part of the dome under the action of the gravity. In the case without cracks, the upper part of the dome is destroyed mainly by the cracks produced at the intrados by the blown of the blast. Moreover, while in the case with cracks these propagates also downward in the cylindrical wall of the Rotonda, this last is practically unaffected by the explosion. As a conclusion, we can say that also if the existing cracks are ignored in the simulation, the dome is destroyed by the blast, though to a less extent, while the cylindrical wall remains intact (see Fig. 22). The existing cracks constitute hence a weakness of the whole structure of the Pantheon for what concerns the effects of an explosion.

Considerations about the numerical simulations
The results presented hereon have been produced by a numerical simulation; though there is always, in such a kind of problems, a certain ineradicable difference between reality and simulation, nevertheless they certainly grasp the main response of the structure to a blast.
For what concerns the simulation of the blast itself, it clearly shows how much it is important to take into account for reflected shock waves, a fundamental aspect of the phenomenon for internal explosions. This means that a study of such phenomena done using empirical models like CONWEP, not taking into account for reflected shock waves, can lead to results really far from reality and principally underestimate the effects of the blast on the structure.
The response of the structure deserves a commentary. The results shown hereon have been D R A F T obtained for the least quantity of explosive able to destroy, at least in part, the Pantheon. Without entering in detail, it is thought that bringing such a quantity of explosive in the interior of the monument would be impossible, if not with some mechanical means.
In other words, the Pantheon is intrinsically safe with respect to explosions. This is due to different reasons. On one hand, its dimensions: it is a large building and the structure itself is really massive: the cylindrical wall thickness varies from ∼ 2.4 m to ∼ 6 m, while that of the dome of ∼ 1.5 m at its top, as already mentioned. Then, its form: the cylindrical form is by itself more likely to resist to impact pressures than a flat one, more exposed to bending stresses. Finally, and this point is really surprising, the material. As already highlighted in Sec. 4.1, Roman concrete has mechanical properties quite different from modern light concretes: its Young's modulus (∼ 3 to ∼ 13 GPa, , [Lamprecht, 1984]), good tensile strength, ∼ 0.6 MPa, and fracture energy, ∼ 55 ÷ 60 J/m 2 , confer to Roman concrete very good properties of resilience and shock absorbing. Though, of course, ancient Romans did not conceive this material for absorbing a blast effects, nevertheless it reveals once more to be a key of the structural success of the Pantheon.