Modeling of shock wave propagation over the obstacles using supercomputers

A computational technology has been developed for modeling the propagation of shock waves in volumes with obstacles. It was shown that the proposed approach allows to simulate real situations of shock wave propagation with an accuracy of several percent (for pressure wave amplitude). This accuracy is achieved through parallelization of the computational process and the use of restructured meshes. An example of calculation is given.


Introduction
The problem of the shock wave propagation in a cluttered space is of great interest from a practical point of view. In particular, to assess the consequences of the explosion of fuel-air clouds. Such clouds can form in the atmosphere during an emergency release of fire and explosive liquids and gases.
To solve this problem, it is necessary to describe a number of rather complex processes:  release of gas, liquid or multiphase flow from equipment;  propagation of the substance in the environment (dispersion in the air, spreading over the surface);  combustion/explosion of clouds of the fuel-air mixture, generation and propagation of shock waves;  action of shock waves and other damaging factors on various objects. Each of these stages is described in terms of its own models, which in turn can introduce various errors in the solution.
It should be noted that potentially the largest errors are introduced at those stages of the calculation, where there may be abrupt changes in physical quantities that will determine the level of damage. In this case, for a correct assessment of the negative impact, it is necessary to resolve the space-time changes with the required accuracy.
To such abruptly changing parameters, first of all, one can attribute pressure Indeed, in the combustion/explosion of fuel-air clouds, the pressure sharply increases in the generated shock waves / compression waves, and then decreases in the rarefaction waves propagating after them. As a result, the pressure profile looks as shown in Figure 1 by the red line (a). In numerical simulations, if the spatial resolution is not sufficient the gradients in the pressure profiles are "smeared out", the pressure maxima are "cut off," and the picture of the "fine" wave propagation structure is lost. As a result, the wave profile (see Figure 1, the red line (a)), is reproduced approximately (see Figure 1, blue line (b)). The difference in the peak pressures between the "exact" shock wave and the "approximate" numerical solution (namely, this peak pressures determine the nature of the consequences) can reach several tens or even hundreds of percent. Figure 1. Shock wave pressure profiles: highaccuracy calculation (a) and calculation performed using a "coarse" mesh (b). So it is clear that there is a need to develop methods and software for high-accuracy modeling of shock wave propagation processes in 3D volume with obstacles.
Until recently, the necessary accuracy of the calculation could be achieved only for onedimensional [1] - [3] or two-dimensional [4] - [7] problems. But now calculations in a threedimensional setting can be made with the appropriate accuracy. There is a relatively small number of software codes capable to adequately describe the processes of generation and propagation of shock waves; FLACS [8][9][10] and KFX [11,12] are the most important codes for simulating industrial accidents. These codes are, in fact, "black boxes" and have a fairly high cost, so the development of their analogs, based on published models and methods, is an urgent task.
In this paper, we propose approaches to mathematical modeling of explosive processes in various emergency situations. In this case, the proposed approaches use such elements as:  technology of parallel calculations;  technology of restructured mesh.
Using the high-accuracy numerical methods, the proposed approach allows in many cases to achieve the reasonable accuracy necessary for predicting the parameters of shock waves with an error of several percent.

Basic equations and method of solution
To calculate the propagation of shock waves in a gas mixture, the following standard equations for the transfer of mass, momentum, energy, and individual components of the gas mixture are used (without taking into account chemical reactions, wall friction, turbulence and buoyancy): The equation of mass conservation: The equation of momentum conservation: The equations of conservation of individual gas mixture components: here x ispatial coordinates, ttime, ρdensity, u i -velocity components, ppressure, Ethe specific total energy equal to the sum of the thermal component, chemical component and kinetic energy, Y imass fraction of the i-th gas component of the mixture.
The system of equations (1) -(4) is closed by the equation of state for a mixture of ideal gases. Equations (1) -(4) were integrated using the second-order Godunov-Kolgan method [13]. The Godunov-Kolgan method is a conservative monotone method. This method allows to resolve the discontinuous solutions (shock waves and contact discontinuities) that are smeared over the mesh. This smearing for the shock wave is usually 3-4 cells and does not depend on the wave intensity. This is acceptable with the available processing power.

Obstacle accounting
The correct formulation of the problem requires that equations (1) -(4) must be supplemented by boundary and initial conditions.
In the simplest case, the boundary conditions must be specified only at the boundaries of the computational domain, which for rectangular decomposition is a parallelepiped. Accordingly, the boundary conditions must be defined on each of the six faces of the original parallelepiped.
In the presence of obstacles, a part of the space within the calculation region is filled with objects that are bounded by impenetrable rigid surfaces (examples of such objects are, for example, buildings, equipment, fences, etc.). In this case, the computational domain can be treated as an initial parallelepiped, from which obstacles are removed. The surfaces on which it is necessary to set the boundary conditions are now not only the six surfaces of the original parallelepiped, but the entire set of surfaces that restricted obstacles. Accordingly, the definition of boundary conditions on these surfaces will allow modeling of the actual presence of obstacles. In practice, this is achieved by combining the boundaries of the obstacle with the boundaries between the difference cells (with the some error that naturally arises in this case) and setting conditions on the boundaries of the difference cells such as a "rigid wall" without sticking, which corresponds to setting zero normal velocity and tangential velocities from the center of the nearest cell from the calculation region.

Use of restructured meshes
As already noted above, the solution of most problems arisen in calculation of shock wave the propagation faces the problem of obtaining the required calculation accuracy.
The most obvious way to increase the accuracy of the calculation (for the given numerical method) is to reduce the size of the mesh cell. However, this method is associated with a significant increase in the cost of computational resources: since explicit schemes are usually used to calculate the propagation of shock waves, a double decrease in the space step size leads to 2 4 increase in the number of operations performed, and, consequently, to the same growth in the need for time of calculation.
So the question arises how to increase the accuracy of the calculation without increasing (or significantly not increasing) the costs (primarily computer memory and time of calculation).
The most effective way to improve accuracy is to use non-uniform meshes. Irregularity of the mesh suggests that the difference cells in different parts of space and at different times have different sizes. In the case of 3D problems with parallelepiped computational region there are three main approaches for constructing an irregular mesh:  concentration/rarefaction of spatial steps along one direction; such a construction of the mesh can lead to a decrease in the accuracy of the numerical method due to a decrease in the accuracy of the approximation of spatial derivatives;  nesting of meshes ("patches"/mesh refinment) is organized by dividing the large cells into several smaller ones; such construction of a mesh requires the use of rather complex algorithms for hierarchical data organization; also, for such a mesh, complex algorithms of moving / moving in the flow of a region with a nested mesh may be required;  the restructuring of the mesh is extremely effective when the process evolves over space over time, covering more and more new regions; in this case, there is no need to immediately perform a calculation on the entire calculation region, since the parameters may not change in some regions.. In this paper, the problem of propagation of shock waves is considered -a kind of a flow that progressively develops into newer and new regions, therefore, the reorganization of the mesh can be used: at the beginning of the calculation, a relatively small calculation region is chosen and calculation is carried out only in it with a small step size; when the flow develops to such an extent that it is about to go beyond the boundaries of the calculation region, the sizes of the calculation region and the sizes of the mesh are increased by a predetermined value to allow the further flow evolution. An example of this approach is shown in Figure 2.

Parallel calculations
The numerical method chosen for solving the equations of gas dynamics uses a three-dimensional structured rectangular mesh. The complexity of the computations in each cell is similar, since the calculations are based on the same formulas of the numerical method. Thus, if each of the parallel processes performs calculations with the same number of cells, then the load on the processors will be the same, and they will not expect each other to exchange data. Therefore, in this case it is more efficient to parallelize the calculated algorithm not at the level of tasks that are the same for each cell of the difference mesh, but use the so-called "geometric parallelization".
With "geometric parallelization", each compute node will show approximately equal performance. In this case, the efficiency of the calculation depends on the time spent on exchanging information between the processors. The time spent on exchanging information is determined by the number of exchanges that each process must perform and the amount of information sent in one exchange.
Each process must, independently of other processes, perform calculations for one step in time in the cells that are part of its zone of responsibility. For this, the process must have information not only about cells in its parallelepiped, but also about cells that are in the zone of responsibility of other processes. When we use the Godunov-Kolgan method, these will be cells that are in two layers surrounding the parallelepiped (see Figure 3). Formally, these two layers of cells provide the setting of "boundary conditions" for the area of responsibility of the processor. The process does not make calculations for these "boundary cells", they are only needed to perform calculations in the cells lying inside the parallelepiped (zone of responsibility).
Exchange of information on the parameters of "boundary cells" occurs between the processes of the area of responsibility, which have a common surface. For each process, whose parallelepiped does not have a common border with the boundary of the entire calculation area, it is necessary to perform up to 6 such exchanges, that is, 6 times to send data about its cells and 6 times to receive data from neighboring processes. Figure 3. "Cross-linking" of computational domains, which are computed on different computational nodes/processors.. In this case, it is necessary to coordinate the sending and receiving of information between the processes. At the beginning, each process calculates the process numbers from which it should receive data, and to which it must send the data. After that, the exchange processes occur according to the following algorithm. First, the processes exchange information on the surfaces perpendicular to the X axis, then the Y axes and after that the Z axis. On the processed axis, the process sends information to the process that serves the "left" area from it, and receives data from this process After that, similar operations are carried out with the process, which is "on the right."

Calculation example
To test the proposed model and the developed program (including a parallel version), the following task was solved: on a flat surface there are several rectangular parallelepipeds, each of which is an object with absolutely rigid walls. The layout of the parallelepipeds is shown in Figures 4a (top view) and 4b (three-dimensional plan). From the scale segment in Figure 4a it can be seen that the length of the parallelepipeds is of the order of a several tens of meters. And their width is several meters. Such scales are typical for typical industrial buildings Usually on industrial sites there are buildings in 2-3 floors, therefore as height of considered parallelepipeds (buildings) was chosen to be 8 m. To generate the shock wave, a sphere containing compressed air was specified (in Figure 4, this sphere is marked in red). The expansion of such a sphere will allow to simulate realistic shock waves.
Next we will use the following agreement on directions in Figure 4a: the direction to the top will be considered as the north, the direction to the right to the east, and so on.
The lower point of the sphere with compressed gas touches the ground surface. The radius of the sphere is equal to 15 m. The pressure in the sphere is equal to 10 atm. This situation corresponds to the value of internal energy in the cloud in the amount of 456000 MJ, which corresponds to the combustion of approximately 9.9 tons of conventional fuel with a combustion heat of 46 MJ/kg. Combustion of this amount of fuel is possible with sufficiently large releases.
As already noted above, on the ground and on all surfaces of the "buildings" the "rigid wall" condition is specified without sticking.
Calculations using the approach outlined above were carried out in four versions:  "coarse mesh" with a step of 0.96 m on a mesh of 520х520х400 cells; the calculation time for 256 processors was 4 hours 56 minutes; during this time the process was calculated for 1 s;  "medium mesh" with a step of 0.4 m on the mesh of 744x744x500 cells; the calculation time for 320 processors was 11 hours 31 minutes; during this time, 0.33 s have been calculated; during this time the process was calculated for 0.33 s  "fine mesh" with a step of 0.2 m on the mesh of 1000х1000х500 cells; the calculation time for 320 processors was 15 hours 43 minutes during this time the process was calculated for 0.19 s  calculation on a mesh with a large step (2 m) on a single-processor computer; such a calculation can be considered typical from the point of view of obtaining information about the propagation of shock waves in a obstacle volume using a single-processor computer in an acceptable time frame (several hours). Figure 5 shows the loading of the north side of building 2 (in the middle of the wall) at points located above the ground at altitudes of 0 m (Figure 5a), 4 m ( Figure 5b) and 8m (Figure 5c), the last point is already on the edge of the roof. From this figure one can see that the loading along the lower part of the building (at heights of 0-4 m) is noticeably larger than on the upper part of the north wall. This is due to the fact that after a sudden increase in the pressure after the arrival of the shock wave, its decrease immediately begins, and this decrease takes place due to both the rarefaction wave propagating behind the incoming shock wave and the rarefaction waves generated at the boundaries of the buildings when they are flowing along the incident wave. And, if the decrease in pressure from a rarefaction wave "sitting" on an incident shock wave is the same for almost all points of the northern facade, the rarefaction formed at the edge of the "northern wall-roof" of building 2 begins to affect the upper part of the wall earlier and only after some time this rarefaction reaches the bottom. A similar loading takes place in the sections of the northern wall adjacent to the western and eastern walls of building 2.
The loading at the lower point of the northern wall of building 2 is 237 kPa, the loading at 4 m is 233 kPa, and the loading at the roof level (8 m) -182 kPa (calculation data for a "fine mesh").
From the data presented, one can see that using the proposed approach, it is possible to significantly increase the number of cells (in tens or even hundreds of times), expand the calculation region and duration in comparison with the calculation on a single-processor system where it is possible to carry out calculations on meshes with only a few million cells..
As one can see from Figure 5 calculation on a single-processor computer in 2-2.5 times underestimates the peak amplitude of shock waves.
From Figure 5 it is clearly seen that a decrease in the size of the space step leads to a significant improvement in the accuracy of the solution. From Figure 5 it follows that at a step of 0.2 m accuracy is attained a value which can be considered acceptable: the difference between calculations on a "fine mesh" (0.2 m) and "medium mesh" (0.4 m) is several per cent -5-10%; this value can be regarded as the accuracy of the calculation. This means that with the next step size decrease, the improvement in accuracy can reach a value of 3-5%. Accordingly, the size of the spatial step of not more than 0.2 m can be considered as the step, the calculations with which can be considered as reliable. So the value In Figure 5, in addition to the first peak of pressure (the effect of the primary wave formed in the atmosphere with the motion of the resting air), a second pressure peak is also observed. It is a secondary shock wave, arising upon repeated expansion of the cloud after it's collapse, which in turn follows the first expansion of the cloud (when the first shock wave was generated). This wave comes to the building 2 at 0.4-0.45 s and is calculated with some "smearing", which means a decrease in the gradient at the front and smoothing its maximum.
In addition to the primary and secondary waves, pressure dependencies on time can also be demonstrated waves reflected from neighboring buildings. This situation is observed for building 1, where the building itself is relatively small, and the nearby building 2, because of its larger size, is capable to reflect waves of greater intensity, since these waves will be weakly attenuated at distances comparable to the size of building 2. Figure 6 shows the pressure versus time dependences for a point located in the middle of the west wall of building 1 (at a height of 4 m from the ground). Due to the location of the buildings 1 and 2 and the small distance between them, a region of high pressure is created (see Figure 4), when shock wave arrives. Unloading of this corner space occurs through the gap between buildings 1 and 2. Accordingly, at the exit from the corner of buildings 1 and 2, a shock wave is formed, which creates additional loading on the western wall of building 1. This is clearly seen in Figure 6, where the primary shock wave hasn't a ordinary triangular profile, but has a trapezoidal profile with two local maxima. And the second pressure peak at this trapezium is higher than the first one. This second pressure peak was formed by the interaction of waves that rounded building 1, and above all by the wave coming from the northern wall of building 2. The results of the calculation using a restructured mesh are shown in Figure 6 also. Calculations on 256 processors took 8 hours, about the same as on a "coarse" mesh. However, in terms of accuracy, the solution on a restructured mesh significantly exceeds the solution on a "coarse" mesh in calculating the fine structure of the wave.
Actually, the restructured mesh initially had a size of 60x60x40 m with the initial spatial step of 0.111x0.111x0.1125 m (520x520x400 cells -the size same as for the "coarse" mesh). During the calculation of a process with a duration of 1 s, 5 rearrangements of the meshes were made. At the same time, each time the spatial step of the next mesh increased by 1.553092 times, reaching 1,039156497 m by the end of the calculation.
Similar loading patterns are observed for other buildings. In conclusion of this paper, it is interesting to consider not only the behavior of pressure at individual points, but also the spatiotemporal variation of pressure. The pressure fields at ground level are shown in Figure 7  From Figure 7a for the time moment 0.113 c. one can see how the shock wave that came to building 2 (and also to the building located north-east of the site of the explosion) is reflected from the northern and eastern walls of the building. In this case, the pressure significantly increases at the point of reflection (this pressure corresponds to the first peak at the gage dependences). It is also clearly seen from Figure 7a that a region of low pressure forms at the point of the "explosion". It is also clear that the incident shock wave has an axisymmetric shape.
In Figure 7b for the moment of time 0.196 s one can see further reflection of the shock wave from two buildings in the north and south-east (areas of intense red color). These reflected waves, when approaching the edge of the building, exit into already unconfined volume, which somewhat amplifies the waves that go further through the air. This can also be seen in the case of a building located northwest of the "explosion" point (see Figure 7b and 7c)).
From Figure 7b one can see how the reflection occurs in the corner formed by the buildings 1 and 2. Due to the partial confinement of space in this corner, an area of high pressure exists for some time (see Figure 7b). And in this region, one can clearly see how the high-pressure region formed in the corner is unloaded through the gap between buildings 1 and 2. As a result of which there is a rather strong wave (area of high pressure) behind the building 1 on its western side, which increases the loading of building 1 on the western wall (this loading in Figure 6 can be related to the presence of a second peak at the pressure-time dependence) From Figure 7b it is also clear that loads on long buildings located so that the shock wave "slides" along the building (for example, a building in the east of the "explosion" point) are much less effective.
It should be noted one more feature, clearly visible in Figure 7b, by the time 0.196 s the expansion of the cloud in the place of the "explosion" is replaced by its collapse, as a result, in the epicenter of the "explosion" the pressure increase begins.. In Figure 7b, the zone of increasing pressure at the place of the "explosion" just begins to form, and in Figure 7c it has already formed in place of the initial cloud (red circle in Figure 7c. This area will subsequently expand and generate a secondary shock wave that is clearly visible on the dependences of most gauges (see Figure 5-6).
It should be noted one more feature of shock wave propagation in a space with obstacles: in the initially spherical wave, after its interaction with obstacles, there are sections of the front with different pressures on it. Thus, when the shock wave propagates in a volume with obstacles, the inhomogeneity of the parameters along the front of the shock wave remains after the release of the shock wave from the clutter. Thus, after the propagation of the shock wave in a volume with obstacles, even when the shock wave leaves this volume, the inhomogeneity of the parameters along the front of the shock wave is presented. The magnitude of pressure variation along the front can reach up to several tens of percent. Figure 7d shows the situation when the primary wave has interacted, practically with all the buildings. This figure also clearly shows the secondary shock wave propagating from the place of the "explosion" as a result of repeated (after the collapse) expansion of the cloud in the place of the "explosion". When the secondary wave propagates, the same effects are observed as in the propagation of the primary wave: reflection from the walls, the appearance of zones of increased pressure in partially confined spaces (for example, in the corner arrangement of buildings), generation of shock waves in the case of unloading of high pressure regions through the canyons between buildings, weaker loading of buildings when shock waves slide along them, etc.

Conclusions
A computational technology has been developed for calculating the propagation of shock waves in a space with obstacles that allows us to take into account not only the action of the primary wave with an acceptable accuracy, but also the effect of waves reflected from other obstacles, secondary shock waves, and also to see with good resolution the nature of the load of various parts and elements of buildings, which increases the reliability of assessing the impact of shock waves on them.
It was shown that  for realistic estimation of shock wave parameters in modeling typical industrial accidents, a spatial resolution of the mesh of 5-10 cm is required (for explosions up to several fuel tons);  the modeling of the shock wave propagation with a spatial step of several meters introduces an error in the amplitude of the shock wave of 200-300%; from this it follows that now simulation of shock wave propagation in a three-dimensional formulation is possible only on clusters with several hundreds of processors;  the use of restructured (expanding) meshes is a promising way to optimize the solution of such problems and allows to resolve pressure peaks with an error of several percent.