Blast load estimation using Finite Volume Method and linear heat transfer

From the point of view of people and building security one of the main destroying factor is the blast load. Rational estimating of its results should be preceded with knowledge of complex wave field distribution in time and space. As a result one can estimate the blast load distribution in time. In considered conditions, the values of blast load are estimating using the empirical functions of overpressure distribution in time (Δp(t)). The Δp(t) functions are monotonic and are the approximation of reality. The distributions of these functions are often linearized due to simplifying of estimating the blast reaction of elements. The article presents a method of numerical analysis of the phenomenon of the air shock wave propagation. The main scope of this paper is getting the ability to make more realistic the Δp(t) functions. An explicit own solution using Finite Volume Method was used. This method considers changes in energy due to heat transfer with conservation of linear heat transfer. For validation, the results of numerical analysis were compared with the literature reports. Values of impulse, pressure, and its duration were studied.


Introduction
The need to accurately quantify the blast overpressure loadings is important because detonations represent a common threat for the security design of building (i.eg.progressive collapse [1]) as a result of big pressures and also fire [2].
For several years great effort has been devoted to the study of blast response of structural elements such as columns [3], slabs [4], RC walls [5], brick walls [6] or prestressed elements [7].In those papers and in related references, however, the blast load overpressure is estimated with simplified triangular distribution in time.It is well-known that exponential function approximates the values of the overpressure distribution in time more closely [8,9].As mentioned in [10], this indicates changes in the pressure-impulse diagrams (Δp(I)) proposed in [3,4].
The paper presents a method of numerical analysis of the phenomenon of the air shock wave propagation that can be used to make the Δp(t) functions more realistic.This method describes an explicit own solution using Finite Volume Method (FVM) and considers changes in energy due to linear heat transfer.For validation, the results of numerical analysis were compared with the literature reports.The free field explosion and also onedimensional flow (an explosion in the pipe) and three-dimensional flow (explosion within the compartment) of the shock wave were analyzed.Values of impulse, pressure, and its duration were studied.Finally, an overall good convergence of numerical results with experiments was achieved.

Explicit solution using Finite Volume Method
The unknown Δp(t) function reflects the thermodynamic state in each point of disturbed gaseous medium, and also reflects the influence of the boundary conditions.The most advanced models available nowadays for Computational Fluid Dynamics are based on Finite Volume formulations, in which the governing equations for the fluid domain (equations for a compressible inviscid fluid, expressing the conservation of mass, momentum and energy) are formulated and solved in conservative form.The Finite Volume Method (FVM) is one of the best methods of solving issues related to gas flow [11].Thus conservation of mass, momentum and energy can be written as follows [12]: When analyzing the state of energy, the change of kinetic energy in heat energy was considered, taking into account the linear heat transfer [13].It was assumed that the thickness of the shock wave front is equal to three cell sizes.If the particle velocity in the first cell is equal u, then the particle velocities in second and third are equal to 0.67u and 0.33u respectively.The temperature decrease was assumed in the same way.Before arrival of the air shock wave the particle velocity is equal to 0 and the temperature is equal to the ambient gas temperature [14].
For further analysis finite difference scheme was applied.In the fluid domain a nodecentered Finite Volume scheme is adopted and uses the classical explicit Finite Difference time integration scheme.Finite Volume Method on a full unstructured grid is introduced to treat the fluid domain.Equations (1) to (3) have the following difference solution in each of the three orthogonal directions (to compare the 1D case see [11]):  The energy of the system is conditioned by the state of the kinetic energy (velocity dependent) and heat.For the purposes of this paper the simplest model of heat transfer was applied (namely the linear temperature decrease) which is expressed as follows [13]: where: c p -thermal conductivity coefficient of gaseous medium.
This strategy requires careful consideration and proper synchronization of the time integration schemes used in the two subdomains.The error introduced should be negligible given the shortness of time increments in explicit time stepping.This paper presents also an extended version of equation of the maximum time step in the 3D as given: (10) where: C FL -non dimensional parameter at values from 0 to 1, a 1 -sound propagation velocity in gaseous medium.
The system of equations ( 4) to ( 9) expresses the gaseous medium flow in the free field region.Graphical interpretation is presented in Fig. 1.The considered volume of gaseous medium is outlined by thick lines and the adjacent volumes by dotted lines.It is assumed that the two parallel sides are being displaced with velocities u n and u n+1 due to changes in energy.This results in displacement of the mass (M) of gaseous medium to a finite volume, which is highlighted by thick line, and hence the change of density and mass of this volume.Weight increase is associated with a pressure change (p).Consequently, there is also a change in energy.Next, a loop is performed over all finite volumes for the following time steps in order to compute the internal forces.Assume that a complete solution, i.eg.all discrete quantities related to the gas state, are known at time t n and one wants to find the solution at t n+1 .First, the velocity from equation ( 9) is computed.Therefore, the new mass is evaluated via equation (4).Then, the internal forces are computed via equations ( 5) to (7).Divide values of computed external forces by surface areas to give the pressure values, then subtract the value of ambient air pressure from the external forces to compute the overpressure values (Δp(t)).The system of equations must be supplemented by the boundary conditions at the interface of building compartments with adjacent finite volumes and at the point of detonation of condensed explosive.When assuming boundary conditions one can simulate the inhibition of gaseous flux by the building compartments.The assumption of the compartment velocity equal to 0 (boundary condition) in case of high mass of compartments (concrete, RC) is reasonable, because the blast loading is completed before the compartment deformation started [15].When considering light compartments (made of steel or glass) a suitable coupling strategy must be chosen.One of the best is the strategy based upon suitable kinematic constraints on the velocities of the fluid and of the structure along the fluid-structure interface [16].In the point, when the shock wave region starts, some boundary condition should be known.These are the particle velocity u 1/2 1 , the shock wave pressure p 1/2 1 , the density of post-explosion gases ρ 1 1 and the temperature of postexplosion gases.They can be computed based on the rules presented in [17] and as follows below.
This study examines also numerical modeling of detonations of spherical and cubical high explosives and characterizes their effects in the near field, called the region of postexplosion gases.This is the region defined by a distance equal to approximately from 10 to 15 times the charge radius, within which the shock wave is affected by local phenomena, including expansion of the detonation products.In case of spherical charge gasses has the shape of a sphere, and in case of cubical charge the shape of an octahedron [18].
Then the procedure based on [17] is applied.In the first step the volume of the postexplosion gases should be known.When considering the discrete space the quantity of cells is the main factor (s).In case of a sphere, the finite volume at coordinates (l; k; m) is inside the region of post-explosion gases in n th time step if the following inequality is satisfied: (11) and in case of an octahedron (the plane equation): . n m k l d (12) Note that in the time-stepping procedure the gas density is obtained first.Then, the particle velocity is solved on the current configuration.Finally, the gas pressure is obtained as the last result as follows: where: ρ fdensity of post-explosion gases, mmass of explosive charge, s -quantity of cells, v fparticle velocity of post-explosion gases, ρ 1ambient air density, ρ fambient air density, Ddetonation velocity of considered explosive material, p fpressure of postexplosion gases, p 1ambient air pressure.
Then the density of post-explosion gases, computed from equation ( 13) (assume that the cell size is equal to the charge dimension) is equal to the air density in 7 th time step in case of spherical charge (2 x 6,5 = 13 -the medium between 10 and 15) and till this point the shock wave begins to propagate.In case of a cubical charge the region of post-explosion gases finishes in 9 th time step.This difference is due to bigger volume of a sphere than an octahedron entered into this sphere.Discrete regions of post-explosion gases in case of spherical charge (7 th time step, squares) and cubical charge (9 th time step, dots) are presented in Fig. 2. Fig. 2. Discrete spheres (drawn using squares) and octahedrons (drawn using dots) in the following time steps When calculating the overpressure values in the shock wave region, assume that some essential boundary conditions are imposed, which can be expressed by the initial velocity (u 1/2 1), the shock wave pressure (p 1/2 1) and the density (ρ 1 1 ) equal to the particle velocity (v f ), the pressure of post-explosion gases (p f ) and the gas density (ρ f ) respectively.The first test problem is the classical free field detonation of 1 kg TNT at a standoff distance of 5 m, for which an empirical solution is available.Nowadays scientists do not carry on research in such a simple configuration (without compartments reflecting the shock wave), so the empirical equations based on [17] were applied.The relative overpressure proposed by Sadowski (Δp f1 /p 1 ) or the relative overpressure proposed by Brode (Δp f2 /p 1 ), the overpressure duration (τ in ms), the exponential parameter (n) and the overpressure distribution in time (Δp) are given by: °°°°® ' ' ' where: r -standoff distance between the considered point and the detonation point.The overpressure diagram is presented in Fig. 3a.
The second case considers a steel cube with insusceptible walls.One cube wall was joined to a tube of 16.8 cm in diameter and 128 cm in length, which was in turn connected to another tube at some distance (see Fig. 3b) [19].A TNT explosive charge with mass of 18.5 g was detonated at the joint of the cube with the tube.Graphs of the overpressure recorded by sensors located in points 1 and 2 provide the best illustration of the assumptions for the 1D model.The overpressure diagram is presented in Fig. 3c.
To validate correctness of three-dimensional solution, blast pressure distribution from a detonation of 1 kg TNT charge in the center point of a vented room was examined [20].The room was a composite steel and concrete structure of a horizontal square projection with a side of 2.9 m and height 2.7 m (internal dimensions) with a hole in the roof of 1.20 m in diameter (see Fig. 3d).Nine pressure gauges were installed on one of the walls (G1 to G9).The overpressure diagram is presented in Fig. 3e.
Table 1 presents initial and boundary conditions.The cell size should not be bigger than the charge dimension.The compartment velocity is equal to 0. The time step value (Δt) is equal to Δx/7000.This is the time to reach the distance of first cell by the shock wave (the TNT detonation velocity -7000 m/s [21]).

Results and discussion
Fig. 3a, 3c and 3e present graphs with an overpressure versus time in 1 st case, 2 nd case and 3 rd case respectively.As can be seen, numerical results reflect well the time to reach the shock wave and the duration of the shock wave.It can also be observed that when the overpressure obtained in the tests increases, the overpressure obtained numerically also increases.The same applies to the overpressure decrease.The values of maximum overpressures and impulses (area under the overpressure graph) obtained numerically and those from literature reports are similar.As mentioned in [10] this results in the change of pressure-impulse diagram location.In consequence, the point at coordinates (I,p max ), which expresses the detonation parameters, could be beneath Δp(I) diagram location.This means the possibility of surviving the element during the detonation, which normally should be destroyed.(see Fig. 4).
The method can also be used for analyzing internal explosion, which tends to be difficult to determine in the case of using empirical equations [9].

2 2( 3 )F
where : t -time, ρ -gas density, V -finite volume, S -surface area, v n -flow rate of the gaseous medium by the surface S, o v -velocity vector consisting of components [u,v,w] in each of the three orthogonal directions, p n -pressure by the surface S, o m -unit vector of internal forces, c v -specific heat of the gaseous medium at constant volume, Ttemperature, x m q -heat flux density related to the unit mass of the gas, x n q -surface density of the heat flux.
volume mass in cell l and in time step n, Δt -time step, Δx,Δy,Δz -cell size in each of the three orthogonal directions, in each of the three orthogonal directions in cell l(k; m) and in time step n, on the surface of the finite volume, in each of the three orthogonal directions Δqheat flux.

Fig. 1 .
Fig. 1.Diagram of gas flow in Finite Volume Method.

Fig. 4 .
Fig. 4. Influence of change of Δp(I) diagram location on assessment of bearing capacity of element.predict incident overpressures and impulses and provide guidance on the use of reflecting boundaries.The algorithm can be used to analyze the internal explosions, where the wave field is more complex.The paper is the result of research tasks carried out under the research RMN No. 800/2016, implemented in the Faculty of Civil Engineering and Geodesy in Jaroslaw Dabrowski Military University of Technology. )

Table 1 .
Initial and boundary conditions.