Abstract
The subject of the paper is a description of a simple test from the field of terminal ballistics and the handling of issues arising during its simulation using the numerical techniques of the finite element method. With regard to the possible excessive reshaping of the finite element mesh there is a danger that problems will arise such as the locking of elements or the appearance of negative volumes. It is often necessary to introduce numerical extensions so that the simulations can be carried out at all. When examining local damage to structures, such as the penetration of the outer shell or its perforation, it is almost essential to introduce the numerical erosion of elements into the simulations. However, when using numerical erosion, the dissipation of matter and energy from the computational model occurs in the mathematical background to the calculation. It is a phenomenon which can reveal itself in the final result when a discrepancy appears between the simulations and the experiments. This issue can be solved by transforming the eroded elements into smoothed particle hydrodynamics particles. These newly created particles can then assume the characteristics of the original elements and preserve the matter and energy of the numerical model.
1. Introduction
When a projectile flying at high speed collides with a concrete surface, its kinetic energy drops to zero but the internal energy of the system increases sharply. The response of the concrete surface to this type of load often takes the form of irreversible (plastic) deformations. If the kinetic energy of the projectile is sufficiently great and the body of the projectile is sufficiently stiff, penetration of the concrete surface occurs. This type of failure is also accompanied by the chipping off of the concrete and the development of dynamically propagating cracks. The successful execution of numerical simulations of this phenomenon is very difficult, however, particularly when the finite element method (FEM) is used. In order for the results of FEM simulations to correspond with the results of experiments, it is necessary to combine suitable material models with numerical model failure techniques, but it is also essential to avoid numerical problems which often negatively affect the results of the simulations.
Today, thanks to constantly ongoing research into concrete structures, extensive concrete material model databases are available and often implemented in commercial programs such as LS-DYNA [1]. The options and conditions for the use of a given selected material model are, however, often open to debate [2]. On top of that, before the execution of the simulation an assumption often has to be made regarding the type of failure that will be decisive so it is possible to select a material model at all [3]. The response of the concrete also depends on the character of the load [4–7], which makes the choice even more complicated. Generally, in the case of high-speed loading it is necessary to use a material model which takes strain rate into consideration [8–10]. In such cases, equations of state (EOS) are often used to enable the successful description of the material model due to the fact that bodies behave in a similar manner to fluids when under high-speed loading. As far as the software is concerned, an algorithm has to be accessible which will interpret such information appropriately for the computational process. This is enabled, for example, by the previously mentioned LS-DYNA program [1] and also AUTODYN [11]. The HJC model [12, 13] (named after its authors Holmquist, Johnson, and Cook) can then be a suitable material model, as it has the above-mentioned properties; it considers the strain rate and utilizes EOS for description.
However, simulations also need to include a numerical technique which will enable the simulation of continuum failure; in the case of the FEM method, this takes the form of the failure of the finite element mesh. If this technique was not used, the simulation of, for example, the chipping off of material or the perforation of the loaded structure would not be possible. Additionally, numerical problems (the locking of elements, negative volumes, etc.) could arise due to the influence of excessive deformations (distortion of elements). The numerical erosion of finite elements can be such a technique. Even though numerical erosion is primarily used to filter out problem elements from the calculation, it can also be suitably used to aid in the simulation of cracks, penetration or perforation, and also the fragmentation of matter. The technique of combining the FEM and element erosion is often used in high-speed simulations, particularly in the simulation of penetration and perforation [14–19]. However, the criterion of element erosion is not unambiguous [20] and can affect the results of the simulation [11]. For example, the damage parameter is selected as an erosion parameter in [14], while in [15] it is maximum tensile stress, in [16, 18] it is geometric strain, in [17] it is fracture strain, and in [19] the erosion parameter is a combination of the damage parameter value and the maximum principal strain. Despite the use of advanced numerical techniques, the results of simulations are still being compared, most frequently with the values from analytical relationships obtained from an extensive amount of experiments [21, 22]. In the case of specific simulations, results are compared directly with experiments such as [23]. Unfortunately, as a result of the heterogeneity of the structure of concrete, agreement between simulations [24] and experiments cannot be guaranteed despite the use of complex modelling techniques.
The dissipation of matter from calculations as a result of the deletion of elements can be a significant problem in cases where the numerical erosion technique is used. In the majority of cases, this problem is omitted and left unsolved [14–19]. So, why use the FEM method for high-speed load simulations when so many problems and complications arise during calculations? With regard to the existence of mesh-free methods such as smoothed particle hydrodynamics (SPH), it is possible to realize high-speed load simulations without needing to include numerical erosion or deal with problems concerning the dissipation of matter [25, 26]. There are several reasons for preferring the FEM method, the main one being the low computational requirements compared with the aforementioned SPH method. Even though the SPH method can be a good choice for high-speed simulations of dynamic events, even in this case many numerical difficulties arise (tensile instability, zero-energy modes, etc.) which need to be resolved; see also [27–29]. However, if the strong points of the FEM and SPH methods were combined, a very useful apparatus for dealing with high-speed simulations could be created.
With regard to the facts mentioned above, the aim of the paper is to propose a procedure for the execution of simulations of the high-speed loading of concrete by steel projectiles. This approach will combine the FEM method with the numerical erosion of elements, and its primary focus will be on solving the problem with matter dissipation. This dissipation can be prevented via the transformation of eroded elements into SPH particles. In addition, current procedures for dealing with such simulations which do not introduce either the erosion of elements or the transformation of eroded matter will be compared. The concept of combining the FEM and SPH techniques with the inclusion of element erosion aims mainly at the improvement of simulation techniques in cases of high-speed loading in such a way that possible numerical issues are minimized and the results of simulations correspond better to the results of experiments.
2. The Erosion Problem
The element erosion function, while not a material property or physics-based phenomenon, provides a useful means of simulating the spalling of concrete and provides a more realistic graphical representation of actual impact events. Erosion is characterized by the physical separation of the eroded solid element from the rest of the mesh [30]. Though element removal (erosion) associated with total element failure has the appearance of physical material erosion, it is, in fact, a numerical technique used to permit the extension of the computation. Without numerical erosion, severely crushed elements in Lagrangian calculations would lead to a very small time step, resulting in the use of many computational cycles with a negligible advance in the simulation time. Moreover, Lagrangian elements which have become very distorted have a tendency to “lock up,” thereby inducing unrealistic distortions in the computational mesh [31]. The erosion function allows the removal of such Lagrangian cells from the calculation if a predefined criterion is reached. When a cell is removed from the calculation process, the mass within the cell can either be discarded or be distributed to the corner nodes of the cell. If the mass is retained, the conservation and spatial continuity of inertia are maintained.
However, the compressive strength and internal energy of the material within the cell are lost whether the mass is retained or not. Even though the filtering out of unsuitable (unneeded) elements is more a matter for numerical simulations, it can be connected (to a certain degree) with the physical matter of the material model.
2.1. Residual Compressive Strength with SPH
The moment at which an element of the Lagrangian mesh erodes is in conflict with what happens in real life, however. In reality, the material does not cease to exist but is only crushed and flakes off; see Figure 1. Even though one cannot speak about the strength of the material as such, the particles and wedge-shaped fragments that fly off can create secondary or residual strength. In order to better approximate reality, it is advantageous to maintain the presence of even such particles in the simulation. This can be done via the transformation of eroded elements into SPH particles. Subsequently, these freely moving particles with the characteristics of the original materials will interact with the rest of the computational model and thus better reflect reality.
(a) Penetration and spalling
(b) Scabbing
(c) Perforation
3. The SPH Method
Smoothed particle hydrodynamics (SPH) is a mesh-free, Lagrangian particle method, originally developed for modelling fluid flows [32]. In 1982 Zukas first used the SPH method for modelling high-speed impact events [33]. Since high-speed impacts and penetration processes usually involve large deformations, their simulation is generally difficult for traditional grid-based numerical methods such as the FEM. The SPH method is, on the other hand, a particularly suitable candidate for these types of problems.
3.1. Essential Formulation of the SPH Method
The formulation of the SPH method is often divided into two key steps. The first step is the integral representation of field functions and the second is particle approximation. The concept of the integral representation of a function used in the SPH method starts from the following identity:where is a function of the three-dimensional position vector and is the Dirac delta function given byIn (1), is the volume of the integral that contains . Equation (1) implies that a function can be represented in an integral form. Since the Dirac delta function is used, the integral representation in (1) is exact or rigorous as long as is defined and continuous in [32]. If the Delta function is replaced by a smoothing function , the integral representation of is given bywhere is the so-called smoothing function and is the smoothing length defining the influence area of the smoothing function . Note that as long as is not the Dirac delta function, the integral representation in (3) can only be an approximation [32].
The continuous integral representations concerning the SPH integral approximation in (3) can be converted into discretized forms of summation over all the particles in the support domain shown in Figure 2. The corresponding discretized process of summation over the particles is commonly known as particle approximation.
If the infinitesimal volume in (3) at the location of particle is replaced by the finite volume of the particle that is related to the mass of the particles bywhere is the density of particle in which is the number of particles within the support domain of particle , then the continuous SPH integral representation for can be written in the following form of discretized particle approximation [32]:Equation (5) states that the value of a function at particle is approximated using the average of those values of the function at all the particles in the support domain of particle weighted by the smoothing function, shown in Figure 2.
In Figure 2, is a constant related to the smoothing function for particle and defines the effective (nonzero) area of the smoothing function. This effective area is called the support domain for the smoothing function of particle .
The adaptability of SPH is achieved by performing particle approximation at each time step based on particles arbitrarily distributed in the current support domain. Because of this adaptive SPH approximation performed at the very early stage of field variable approximation, the formulation of SPH is not affected by the arbitrariness of the particle distribution as it changes with time. It can therefore naturally handle problems involving extremely large deformation [32]. More information about the construction of the smoothing function can be found in [34].
4. The HJC Material Model
Under shock wave compression, plastic deformation and crack-induced damage should be taken into account in the modelling of concrete behaviour. A constitutive law combining pressure-dependent plastic hardening, damage-softening, and the strain rate effect, especially suited for the prediction of concrete response under dynamic loading such as blasts and impacts [35], was developed by Holmquist et al. [12]. In this constitutive law the normalized equivalent stress is defined aswhere and are the normalized equivalent stress and pressure, with and being the actual equivalent stress and the quasi-static uniaxial compressive strength, respectively. Scalar damage is a value from 0 to 1 that describes the accumulation of damage as a percentage of the full cohesive strength that the material possesses. When , the material is undamaged and exhibits its full strength, but at the material is fully damaged and retains the least confined shear strength. is the dimensionless strain rate, where and are the actual and reference strain rates, respectively. The material constants are , , , , and , where is the normalized cohesive strength, is the normalized pressure hardening coefficient, is the strain rate coefficient, is the pressure hardening exponent, and is the normalized maximum strength that can be developed.
The model accumulates damage from both equivalent plastic strain and plastic volumetric strain . The compression damage is expressed aswhere and are material constants; is the normalized maximum tensile hydrostatic pressure that the material can withstand; and is the maximum tensile hydrostatic pressure.
The equation of state (EOS linear elastic region, transition region, and compact region) is expressed as follows [12, 13] for loading:where is the standard volumetric strain; is the elastic bulk modulus; and are the current and initial densities, respectively; and are the crushing and locking volumetric strains, respectively; and are the crushing and locking pressures; and , , and are pressure constants; see Figure 3. The modified volumetric strain is defined as The tensile pressure is limited to during numerical calculations.
4.1. Implementation of Erosion
The HJC material model offers the implementation of element erosion on the basis of the size of the damage strength, . The damage strength, , is defined in compression when asor in tension if , asThe element erodes at the moment when its damage strength falls below 0. From the physical point of view the implementation of such an algorithm is justified [1].
4.2. Transformation of Eroded Elements
If the above-mentioned problem with eroded matter and its dissipation were not solved, in extreme cases results might be gained which appear not to make sense. An example of this can be found in Figure 4, where a concrete specimen impacts a rigid base at high speed. At the moment when the sample comes into contact with the base, it starts to be crushed. The kinetic energy of the specimen transforms into plastic deformations of its body, and the specimen thus starts to slow down. Figure 4 compares the progression of the impact of the specimen when its numerical description does not enable erosion of the elements (on the left) and when it does enable such erosion (in the middle). It is obvious that the results are significantly different in the final stage. The top surface of the specimen with element erosion almost falls onto the base, while the top part of the specimen without such erosion remains at half its original height.
Even though the assumption was introduced that the only elements which will be eliminated from the calculation are those that would be in the stage of full failure from the perspective of the material model and therefore would not influence the further course of the simulation significantly, it is obvious from the results that such a simple implementation of numerical erosion is inadequate, even if it is linked to the material model.
The third variant of the sample in Figure 4 (on the right) includes the adaptive transformation of eroded elements into SPH particles. These particles have assumed the material properties of the eroded elements and their weight and speed, but also their stress state and so forth. More information on SPH can be found in [32].
It can be concluded from the results that the described residual strength appears thanks to SPH particles. This strength is represented by the interaction of SPH particles with the rest of the numerical model and simultaneously with one another. The variant of the model with SPH particles then corresponds well with the original model where erosion was not included. Mass of the specimen through the simulation is shown in Figure 5.
5. Combination of FEM and SPH
The motivation for the combination of the FEM and SPH methods is mainly the need to filter out numerical problems in cases involving large deformations (FEM), but it is also due to the need to lower the computational requirements of the simulations (SPH). As was stated in Section 3.1, the SPH method can deal with problems concerning large deformations while avoiding numerical complications thanks to its adaptive nature. However, the issues with computational requirements cannot be solved absolutely conclusively, particularly in the case of the SPH method.
It is true for both methods that the computational requirements (time requirements) increase as the number of elements or particles grows. These requirements can be reduced to a certain degree.
In the case of the FEM, the quality and computational requirements of the solution are closely connected with the number of integration points in the element. If numerical complications of the hourglass type are solved, reduced element integration can be used, one integration point per element. Moreover, numerical integration takes place most frequently in the reference element, which is an advantage [36].
In the case of the SPH method, the integration points also represent particles (one integration point per particle), while the numerical integration takes place in the form (5). This means that, in the case of both the FEM and the SPH methods, the computational requirements do not necessarily have to increase above all limits as a result of the increasing number of elements or SPH particles.
However, there is a problem with the SPH method: it is necessary to search for neighbouring particles within the support domain and perform other associated operations in each time step of the explicit algorithm for the solution of a nonstationary task. The size of the support domain can also be different in every step (and actually is in the majority of cases). A small domain size can result in a low-accuracy solution, while a large support domain can cause the smoothing of local properties [32]. If, however, there was no search for nearby particles in each time step, the adaptability of the SPH method would be lost. This would mean losing the capability to successfully investigate excessive deformations.
With a combination of the FEM and SPH method a good balance is attained between computational requirements and the quality of the achieved results. Figure 6 shows a comparison of the calculation lengths required by simulations of the experiment described in Section 6, for a 500 ms−1 projectile impact speed and different numerical model variants.
5.1. Mesh Density and SPH Regularity
In the case of the FEM [37] and SPH [32] methods, the accuracy of the result generally increases with the increasing density of the finite element mesh or SPH particle distribution; see the diagram in Figure 7.
In the case of a combination of material damage models and FEM, localized damage can occur if such dependencies are not filtered out [38]. In the majority of cases such damage concentrates in the smallest element, or the element with the lowest stiffness. With increasing FEM mesh density, failure types can occur which do not correspond with reality. This issue can be dealt with in several ways. If there are links to the fracture energy of the material in the description of the material model, the dependence between the size of the finite element and the shape of the stress-strain diagram [38] can be used; the stress-strain diagram will be altered for various element sizes. If the material model does not include fracture energy, a solution can be achieved via the connection of a nonlocal model [39]. The nonlocal treatment basically attempts to average out the failure/damage values of neighbouring elements in order to minimize the mesh dependency of the results [40]. This technique was also applied to the HJC material model.
The SPH method is also influenced by the density of the particle distribution. However, the regularity of this distribution in the initial stage of the simulation is a greater problem [1]. In the case of a bad distribution, particle clusters can be created. Such locations in the support domain would then contain large quantities of particles. In contrast, very sparsely populated areas can appear. As a result of such clusters, false cracks may occur. They may also have an effect on the behaviour of real cracks, which have the tendency to bypass places where clusters occur. This problem with the SPH method still has not been resolved. Due to this, a regular finite element mesh distribution was maintained in the simulation so that no clusters occur in the SPH particles as they appear.
In the simulation of the experiment from Section 6, studies examining the dependence of the result on the density of the finite element mesh or the regularity of such distribution were not carried out. Nevertheless, such studies are planned for future research. The FEM elements in the area around the point of impact were approximately 1/20 of the projectile’s diameter in size.
5.2. From FEM to SPH
An SPH particle is created at the moment when an FEM element erodes. The conditions for element erosion during simulations are shown in (10) and (11). The SPH particles are included in the calculation from the beginning so that their position and other characteristics can be determined as they appear during the simulation. The position most frequently selected for them is in the centre of gravity of the FEM elements. However, the particles are inactive so that there is no increase in computational requirements. Inactive particles are marked grey in Figure 8.
After the erosion of FEM elements, there are two ways in which interaction can take place between the SPH particles and the undamaged FEM model. The first option is via a contact formulation, most frequently the penalty method. This is drawn schematically in Figure 9. In Figure 9, already active SPH particles are shown in green, while inactive SPH particles are coloured grey. The layer of FEM elements for which the contact algorithm is active is marked in red.
The second option is the creation of a transition layer on the interface between the FEM elements and SPH particles. This transition layer then contains FEM elements as well as active SPH particles, which are coupled. This is drawn schematically in Figure 10. Again, active SPH particles are green, inactive SPH particles are grey, and SPH particles that are coupled with FEM elements to create a transition layer are marked blue.
Transition layers were used in the simulation of the experiment from Section 6. It is expected that future studies will be conducted where the variable will be the amount of SPH particles placed in the element.
6. The Experiment and Simulation
In 2002 Buchar et al. performed an experiment with projectiles; see Figure 11, defined by the NATO STANAG 2920 standard, which they shot at speeds of 300–1400 ms−1 into concrete targets with the dimensions 0.1 × 0.1 × 0.1 m and a strength in compression of 43.1 MPa [41]. The mass of the projectiles was 1.102 g. A good correlation was then found between the performed experiments and the theory applied according to [42].
6.1. Setting Up the Simulation
The most important simulation settings are described in the tables, where Tables 1 and 2 describe material properties, Table 3 describes the setting of the SPH formulation, and Table 4 describes the setting of the FEM mesh. The simulations were carried out in the LS-DYNA program [1], in which robust explicit method solvers are implemented. LS-DYNA also enables the execution of simulations using the FEM as well as SPH methods, or their combination. The chosen size of the coefficient of friction between the concrete block and the projectile was between 0.1 and 0.45. However, studies have shown that the amount of friction does not have a significant influence on penetration depth. This is because a ballistic cavity is created in the surroundings of the projectile as it impacts and subsequently penetrates the concrete slab at high speed. As a result of the creation of this cavity, the influence of side friction between the projectile and the concrete slab is eliminated. However, this can be just a temporary extension, a temporary cavity. In the majority of cases, however, the surrounding area chips off.
6.2. Results of the Simulation
The results of the numerical simulation of the experiment are presented in graphical form in Figure 12. The image shows the stages of penetration of a projectile with an impact speed of 500 ms−1. Once again, on the left, there is the numerical model without element erosion, in the middle there is the model with erosion, and on the right there is the model which also contained the transformation of eroded elements into SPH particles.
The size of the crater in the final stage of penetration (i.e., the quantity of eroded elements) increased exponentially with the growing impact speed of the projectile. This was echoed by the considerably greater depth of penetration for model with erosion in contrast with the model without erosion and the model with SPH particles.
However, this is not in accordance with the experiment. The area where the data were measured is marked in blue in Figure 13. The coloured area in Figure 13 represents acceptable results according to the theory [42], expanded by the standard deviation measured in [41]. Even though other parameters were studied in [41], for example, the shape of the crater or the constants of deceleration forces (see [43]), they were not needed for the evaluation of the functionality of the adaptive transformation algorithm of SPH particles in the simulations that were carried out. In the future, it is expected that studies will be performed that will also include these aspects.
At higher speeds it is thus unacceptable to choose the permitted numerical erosion function without further investigation with regard to the fact that the penetration depth increases excessively with increasing projectile impact speed. The model without the erosion of elements approximates the measured data well. However, due to the large degree to which the elements deform and the very small time steps resulting from this, the time required for the calculations grew uncontrollably. The model in which the elements can transform into SPH particles also approximates the measured data well. Thanks to the erosion, it eliminates deformed elements and additionally provides a great deal of information about the behaviour of fragmenting particles of the model.
Figure 14 shows the area where failure occurs as the projectile impacts the specimen at a speed of 500 ms−1. It shows the damage parameter , which has the values 0-1. When comparing the area of failure for all variants it is obvious that the numerical models with and without erosion are almost identical. Another essential fact is that their areas of failure are regular, oval shaped. However, this behaviour does not correspond to that seen in experiments [43]. In contrast, the model with erosion and SPH particles produces an area of failure which is larger and also irregular. Even though this behaviour can be explained by the irregular interaction of the SPH particles and the interaction between them and the rest of the structure, the failure results correspond a lot more with those obtained from experiments [43].
Figure 15 shows a comparison of the speed of the projectile from the moment of impact up to when it comes to rest. As has already been mentioned several times, the model without erosion of elements and the model with erosion of elements and SPH are very similar from the aspect of the conservation of mass (see also Figure 4). With regard to the fact that the penetration depths of both variants are also similar, it was assumed that the deceleration force and thus the course of deceleration would also be similar [43]. The results in Figure 15 also prove this.
As already mentioned the SPH particles are present in the computational model right from the beginning of the simulation; they are inside the FEM elements, but inactive. The activation of the particles takes place at the moment when the erosion of the element commences, that is, at .
Figure 16 shows the stages of penetration by a projectile with an impact speed of 500 ms−1. The picture also shows the gradual activation of the SPH particles. The interaction between the SPH particles and the still noneroded elements is mediated by SPH particles inside the elements which are on the outer surface [1].
The number of active particles during the simulation with the projectile impact speed of 500 ms−1 is shown in Figure 17.
7. Conclusion
The contribution describes a simple experiment from the scientific field of terminal ballistics. It lists possible negative aspects of numerical simulations and simultaneously proposes a solution to these problems which consists in the numerical erosion of FEM elements and the subsequent preservation of their matter in the form of its adaptive transformation into SPH particles.
The implementation of numerical erosion in simulations can provide a useful tool for the removal of excessively deformed elements which can cause not only the extreme prolongation of the calculation but also its undesired locking, known as element locking. As matter and internal energy can be lost simply by implementing this numerical erosion, which can lead to the production of incorrect results, such dissipations have to be prevented. A suitable solution is therefore to use the adaptive transformation of eroded elements into SPH particles. These new particles assume the properties of elements which have eroded and subsequently create a certain residual strength due to their interaction with the surroundings, but in the sense of already damaged (i.e., loose) material.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This contribution was created with the financial aid of Project GACR 14-25320S, Aspects of the Use of Complex Nonlinear Material Models provided by the Czech Science Foundation and under Project no. LO1408 AdMaS UP, Advanced Materials, Structures and Technologies, supported by Ministry of Education, Youth and Sports under the “National Sustainability Programme I.”