Impact Resistance of Concrete Structures

We present simulations on the resistance of concrete structures due to impact loading. Therefore, a mesh-free method is exploited in combination with a viscous damage model that accounts for strain rate effects and high pressures that commonly occur in those events. In contrast to many other studies in the literature, we account for uncertainties in the material parameters by subjecting them to a probability distribution function. Furthermore, we also consider the geometric correlation inside the concrete and show that this geometric correlation has a minor influence on the predicted results. All numerical results are compared to the impact experiments performed by Hanchak and coworkers as well as our own impact experiments.


Introduction
Predicting the impact resistance of concrete and reinforced concrete structures has been of interest for many years.Constitutive models that account for strain rate effects, high pressure, and the associated dynamic strength increase are the key to successful predictions.Constitutive models accounting for strain rate effects can be classified into damage models, plasticity models, and damage-plasticity models [1][2][3][4][5][6].Damage models for concrete accounting for high strain rates were developed, for instance, by [7][8][9][10].A popular plasticity model for concrete was proposed by [11] as an extension of the Johnson-Cook [12] model originally developed for metals.In this model, the dynamic strength increase was related to the current strain rate.However, memory effects were not accounted for.Furthermore, discontinuities in the stress-strain response are allowed which are physically not meaningful.In [13], the dynamic buckling strength was therefore dependent on the strain rate history and an average strain rate.A viscous damage-plasticity model accounting for high strain rates and high pressures was proposed by [14].Therefore, a dynamic damage variable was introduced that decayed the static damage evolution.The decaying function depends on previous damage-rates or damage increments.The scalar damage was later on extended to a damage tensor accounting for the anisotropy for concrete under tensile loading [15] while a scalar damage was maintained in compression.
Most of the studies on the impact resistance of concrete were done with the finite element method [16][17][18][19][20][21][22][23][24].Commonly, element-deletion techniques [25,26] were applied to allow for large deformations and perforations.A strong competitor to finite elements is mesh-free methods [27][28][29][30][31][32][33][34][35][36][37] that allow for a more robust and accurate prediction of dynamic fracture and fragmentation; see [38] for a recent overview of mesh-free methods and their applications.Reference [39], for instance, was able to predict the impact resistance of the experiments carried out by [40] with mesh-free methods quite accurately, while the FE simulations done in [40] led always to an underestimation of the impact resistance.Also, penetration was not predicted correctly in those simulations.Dynamic fracture in early mesh-free methods was treated by simple separation of particles [41][42][43][44][45][46][47].However, [30,48] showed that care needs to be taken as the fracture often is caused by numerical artefacts.More advanced methods such as the extended finite element method [49,50], extended mesh-free methods [34,[51][52][53][54][55][56][57][58][59], the phantom node method [60][61][62][63][64], or efficient remeshing techniques [65][66][67][68][69][70] have been applied to dynamic fracture of concrete as well [71,72].However, they do not seem suitable to model the impact resistance of concrete with large deformations and an enormous number of cracks.A simpler and more efficient pathway is the Cracking Particles Method (CPM) [73,74] which models the crack as set of crack segments.Similar methods were used by [75][76][77][78][79][80].Other researchers try to combine the advantages 2 Mathematical Problems in Engineering of mesh-free methods and finite elements.Therefore, only the area around the perforation was modelled by a meshfree method whereas a finite element discretization was employed elsewhere reducing computational cost [31,80].Most simulations were based on pure deterministic models that can lead to unrealistic crack patterns as the ones shown in [39].In these simulations, the radial main cracks appear too close to each other.In order to achieve unsymmetrical crack patterns, [81] introduced randomness in the strength to avoid the unrealistic symmetric patterns they modelled before in [82][83][84].However, all other parameters remained deterministic.
Several benchmark experiments have also been conducted to test the influence of the thickness of the specimen or the speed of the impactor on the impact resistance of concrete or reinforced concrete structures [85][86][87][88][89][90][91].One classical experiment was done by [92] who subjected concrete structures to impact with different speeds of the impactor.They measured the residual velocity of the impactor.The recent experiments by [40] also report penetration depth.
This paper presents one of the first studies on stochastic modelling of dynamic fracture of concrete.Our constitutive model is based on the model described in [14] and combined with a Lagrangian mesh-free method for the computational studies.Fracture will be modelled by a simple particle splitting as described in [77].This methodology will be described in the first sections before we present computational results that will be compared to the experiments of Hanchak et al. [92].We also carried out our own experiments that will serve as validation.

Governing Equations and Discretization
The governing equation to be solved is the balance of linear momentum that is given in weak form of the updated description of motion by where  int ,  ext , and  kin denote the internal, external, and kinetic energy, respectively.The domain is denoted by Ω whereas Γ  ∪ Γ  = Γ with Γ  ∩ Γ  = 0 denotes the Lipschitz continuous boundary;   and   indicate the strain tensor and Cauchy stress tensor,   and   are the body forces and tractions, respectively,   is the displacement field, and  is the density while  denotes the first variation and the superimposed dote material time derivatives which are identical to the partial derivatives in an updated Lagrangian framework.
We have employed a stabilised [48] nodally regularised [30] element-free Galerkin method [93].Therefore, the discretization of the displacement field is given by where u  denote the nodal parameters, which are not equal to the physical displacement value at that point, that is, u(x  ) ̸ = u  , and   (x) indicate the shape functions.They are derived from minimising a discrete weighted norm as explained in detail, for example, in [93]: where ) is called the moment-matrix and the matrix B(x) is obtained by Substituting the discretization of the displacement field and virtual displacement field into the weak form yields the well-known equation of motion that is solved for the nodal parameter u  () that is stored in the global vector D: with the mass matrix M = ∫ Ω N  N Ω and the vector of external and internal forces where N and B are matrices containing the mesh-free shape functions and their spatial derivatives, respectively.The integrals are evaluated by stabilised nodal integration which ensures bending exactness in a Galerkin method.The use of an updated Lagrangian kernel formulation as suggested in [39] ensures the stability of the method while maintaining the applicability to extremely large deformations needed for dynamic fracture and fragmentation.

Constitutive Model
We have used a modified version of the constitutive model presented in [14].Instead of using a viscous damage-plasticity model, we removed the plasticity part of this model as it adds material parameters that are difficult to calibrate.Furthermore, preliminary studies indicated that the plastic response is of minor importance as the results in the next section will show.The basic idea of the damage model is given subsequently.For more detailed information, the reader is referred to [14].

Results
Two types of experiments are employed in order to validate the proposed computational framework: (i) the impact experiments carried out by Hanchak [92], (ii) contact explosion experiments done at the laboratory of our university.
Hanchak tested the impact resistance of concrete structures of different concrete strength.The main goal of the experiments of Hanchak studies was to test their high performance concrete compared to standard concrete.The high performance concrete had a compressive strength of 140 MPA while the standard concrete's compressive strength was nearly 50 MPA.Therefore, concrete specimens were subjected to impacts of different projectile velocities.The experimental set-up is illustrated in Figure 1.As our focus is not on high performance concrete but rather on validation of our framework, we present only the results of the impact simulations on standard concrete.In the results, all plates were perforated and the residual impact velocity was measured.We will focus on experiments with penetration later.
A full three-dimensional discretization has been utilized.Preliminary studies showed substantial differences in results obtained from 2D simulations based on plane strain assumptions.Around 2 million nodes were needed in order to get convergence in the residual impact velocity and simultaneously in the damage pattern.The (deterministic) parameters of the constitutive model are listed in Table 1.They are classified in different categories as indicated by the different lines: (1) bulk parameters, that is, Young's modulus, Poisson's ratio, and density; (2) damage parameters  0 ,   , and   determining the "static" damage evolution: the key parameters are   and   controlling the ductility and maximum strength of the concrete and the parameters determining the dynamic part of the damage responses  1 and  2 : they determine how the dynamic damage is decaying over time;  5) the parameters of the Hugoniot-curve accounting for the pressure dependence of concrete under high strain rates and high volumetric strains.All parameters, except the ones associated with the damage surface, are assumed to be stochastic in our computations and are subjected to a uniform probability distribution function with a scatter of 5% and 10%, respectively.In order to identify the influence of each physical mechanism, we modified only parameters associated with this physical mechanism according to the probability distribution function while keeping all other parameters constant.We are aware that a global sensitivity analysis as proposed, for example, by [94][95][96][97] in the context of studies on polymeric composites or by [98] in the context of optimization will provide more quantitative results.However, such studies are out of the scope of this paper and left for the future.
Before discussing the results presented in Table 2, let us discuss the different physical phenomena and the related input parameters: (1) The parameters related to the first line of 1 govern the initial elastic response of the material in the bulk.
(2) The parameters  0 ,   , and   determine the stressstrain response under static conditions.The first parameter determines the transition to nonlinear material behaviour while the other two parameters determine the tensile/compressive strength and failure strain.
(3) The parameters  1 ,  2 govern the dynamic damage evolution and are responsible for the strength increase under high strain rates.
(4) The parameters   ,   are responsible for the amount of compressive damage caused under tensile loading and vice versa.
(5) The last set of parameters in Table 1 account for the pressure dependence of the concrete.
The results in Table 2 show the lower and upper bounds of the residual impact velocity for the different input parameters for different impact velocities.The values of the deterministic simulation are depicted in Table 2 as well.As can be seen, uncertainties in the input substantially influence the residual velocity of the impactor.While parameter sets 1 (parameters related to the initial elastic response) and 4 barely have an effect on the residual velocity (the scatter in the output is much smaller than the scatter in the input), the static damage evolution (parameter set 2) and pressure dependence (parameter set 5) have a substantial influence in the order of the scatter in the input.This is observed consistently for both 5% and 10% scatter.The biggest influence is achieved by the parameters related to the dynamic strength increase which is more pronounced for increasing scatter.Finally, we remark that we also varied two sets of input parameters while keeping the other parameters constant.Since parameter sets 3 and 5 had the biggest influence on the output, we varied those parameters.However, the upper and lower bounds of the residual velocity barely changed.The scatter in the output is slightly higher than the scatter in the input.Note that the effect gets more pronounced with increasing impact velocity.Overall, one can see an excellent agreement between computational and experimental results.
The second test example is a concrete slab under contact explosion as depicted in Figure 2. The explosion was modelled by a (JWL) John-Wilkonson-Lee equation of state.The explosive is also modelled with a mesh-free Lagrangian method, which is one of the advantages of the proposed method: the fluid-structure interaction remains fairly simple compared to a coupled Lagrangian-Eulerian coupling such as (ALE) which is required for finite elements, for instance.
In the experiment, a penetration was observed and the penetration depth as well as the diameter of the crater at the top surface was measured.We used the same concrete as in the experiments performed by Hanchak and therefore the same material parameters of the constitutive model were used.Figure 3 shows the damage shape towards the end of the simulation; the damage is indicated by the red color.The full deformation has not been reached, yet.However, the damage evolution is fully completed already.Table 3 summarizes the key results of the simulation.We tested again the influence of uncertain input parameters for the same probability distribution function and scatter as discussed previously.The observations are somewhat similar compared to the observations for the Hanchak experiments.However, the uncertainties in the input parameters seem to be less pronounced here.The largest effect is attributed to parameter set 3 related to the dynamic damage evolution followed by parameter sets 2 and 5. Parameter sets 1 and 4 have nearly a negligible influence on the penetration depth and damage diameter.

Discussions
We presented a computational framework for predicting the impact resistance of reinforced concrete structure.This framework is based on a mesh-free Lagrangian method which guarantees the applicability to large deformations as they occur under high velocity impact, explosion, dynamic fracture, and fragmentation.A viscous dynamic damage model is exploited which accounts for the strength increase of concrete at high strain rates and the pressure dependent response of concrete.Furthermore, we consider the influence of stochastic input parameters.
The framework is validated by comparison of numerical predictions to experimental data: (1) the famous Hanchak impact experiments which measure the residual velocity and (2) our own explosion experiment comparing the penetration depth and damage diameter at the top surface.Our results show excellent results with the experimental results.Furthermore, we account for stochastic influence.Therefore, we classify the parameters of the constitutive model into 5 categories related to (1) the initial mechanical properties in the bulk, (2) static damage evolution, (3) dynamic damage evolution, (4) reduction factor of compressive damage at tensile load and vice versa, and (5) pressure dependence.By varying each set of parameters while keeping the others constant, we show that the key factor influencing our output is related to parameter set 3. We also reveal that the scatter in the output is not small.

( 4 )
(X − X  ) denoting the kernel function and p  (X) the polynomial basis.To ensure Galilean invariance, we opted to choose a linear basis p  (x) = [1  ].The exponential weighting function has been employed in all our simulations.

Table 1 :
Parameters of the constitutive model.the parameters of the damage surface  1 to  4 ; (4) the reduction factors   and   controlling how much compressive damage can be done by tensile loading and vice versa; and (

Table 2 :
Residual velocity ([m/s]) of the impactor in dependence of the impact velocity ([m/s]) and the scatter in the input parameters: the first line of the same impact velocity refers to a 5% scatter while the second line refers to the 10%.

Table 3 :
Penetration depth and damage diameter for the explosion experiment: the first line refers to a scatter of 5% while the second line refers to the 10% scatter in the input.