Concrete Targets with Heterogeneities under Impact Loading

: The main aim of the paper is to present an elegant procedure for introducing heterogeneity to simplified computational models. A computational model is considered here to be the-oretically any numerical model that includes material with structural strength. Heterogeneity then means the randomness in material parameters, which allows a simplified model to behave like the comprehensively described structure of a material. The whole procedure is explained using the example of a material model of concrete in the form of a concrete block (target) which is exposed to high-speed impact loading. Even though the described procedure for the introduction of heterogeneity to computational models can be applied during the use of practically all numerical methods, the Smoothed Particle Hydrodynamics (SPH) method will be used due to the type of loading involved. The creation of the presented procedure was prompted by the constantly increasing number of input parameters used with material models. Elementary material models with a minimum of inputs can be used when applying the procedure. The end of the paper deals with the improvement in use that enables heterogeneity (i.e. randomness of material structure) to be created in such a way that a simplified model reflects as much as possible the real material structure, not only the variations of mathematical functions.


Introduction
The original idea in the very foundations of structure designing and related domains consisted in creating a simple and robust concept. With the progress of time, simple outlines gradually became more complicated, eventually reaching a stage where a mere pencil with a sheet of paper had ceased to suffice. Interestingly, at present we are often 108 M. Hušek, P. Král, J. Kala, P. Hradil and P. Maňas faced with efforts to design structures as simple as reasonably possible; after being completed, however, such structures may not necessarily be the simplest ones in terms of their behaviour and the presumptions of structural mechanics [1][2][3].
In concrete structures, the material heterogeneity offers a new dimension of complexity. Concrete, as the mixture of a concrete additive (or aggregate for short), a cement binder, and water, does not constitute a regular structure from any perspective or at any scale; with respect to this fact, not even the behaviour of the given structures can be assumed to remain identical in all cases. The introduced randomness is further multiplied by the type of load, which is invariably comprised in structural mechanics analyses [4][5][6].
In real-world conditions, countless loading tests can be carried out, and the rates of occurrence of an effect enable us to presume the most probable result of the series of tests to follow. However, numerical analyses and simulations, using common computational techniques such as the Finite Element Method (FEM), Finite Difference Method (FDM), or Smoothed Particle Hydrodynamics method (SPH), will not provide divergent results even with an infinite number of computations unless special add-ons are applied or changes in input parameters are made.
In simulations where concrete-based materials are involved, it is always a problem related to what material model should be used [7,8]. The results obtained from varied models and their modified versions can differ fundamentally. It is not advisable to introduce the heterogeneity directly into the description of the material model either, mainly due to the increasing complexity of the operations and the lack of knowledge concerning the sensitivities of the model's individual parameters before running the actual computation [9,10]. For example, in the case of the FEM, a material model would have to be generated for each finite element to ensure that the material heterogeneity is contained in the numerical model.
A very scarcely employed option consists in introducing the heterogeneity directly into the source code of the numerical method; however, there still remain the questions of whether this approach can be used to define the random behaviour of the monitored area, or the heterogeneity, and whether such numerical treatment can be further interconnected and explained using a real concrete structure. Numerical attenuation is utilized in fracture mechanics simulations, where the FEM often finds effective application. As the material model is the same in the entire monitored area, the finite element geometry can be modified in such a way that any failure is immediately localizable at a pre-selected point (for example, via notching on a concrete beam). In high-speed loading, however, the FEM does not offer a suitable solution [11,12].
Conversely, the attractiveness of the SPH approach intensifies with the increasing loading speed; yet it also needs to be noted that, in the case of the SPH, there is no physical mesh to interconnect the individual particles, and a mere change in the configuration of the particles may not provide the desired results (in fact, this type of action may even cause adverse effects, including the formation of numerical cracks) [13,14].
The paper is divided into a total of three parts in such a manner that they are linked to one another as logical units. The first part describes the procedure, which can generally be used to introduce heterogeneities (i.e. oscillations of material parameters) into the source code of the SPH method. The second part shows a practical application of this procedure in which a concrete block (target) incorporating the mentioned heterogeneities (material oscillations) is exposed to impact loading. The third part explains how the application of heterogeneities can be exploited further to enable a simplified model to reproduce the real structure of a modelled material as much as possible.

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 f (x) used in the SPH method starts from the following identity where f is a function of the three-dimensional position vector x, and δ (x -x′) is the Dirac delta function given by In Eq. (1), Ω is the integration region of the integral that contains x. Eq. (1) implies that a function can be represented in an integral form. Since the Dirac delta function is used, the integral representation in Eq. (1) is exact or rigorous as long as f (x) is defined and continuous in Ω [13]. If the Delta function δ (x -x′) is replaced by a smoothing function W (x -x′, h), the integral representation of f (x) is given by where W is the so-called smoothing function and h is the smoothing length defining the influence area of the smoothing function W. Note that as long as W is not the Dirac delta function, the integral representation in Eq. (3) can only be an approximation [13]. The continuous integral representations concerning the SPH integral approximation in Eq.
(3) can be converted into discretized forms of summation over all the particles in the support domain shown in Fig. 1. The corresponding discretized process of summation over the particles is commonly known as particle approximation.
If the infinitesimal volume dx′ in Eq. (3) at the location of particle j is replaced by the finite volume of the particle ∆Vj that is related to the mass of the particles mj by where ρj is the density of particle j (= 1, 2,…, N) in which N is the number of particles within the support domain of particle j, then the continuous SPH integral representation for f (x) can be written in the following form of discretized particle approximation [13] as Eq. (6) states that the value of a function at particle i is approximated using the average of those values of the function at all the particles in the support domain of particle i weighted by the smoothing function shown in Fig. 1. Fig. 1 Particle approximation using particles within the support domain of the smoothing function W for particle i

Size of the Support Domain
The extent of the support domain is defined according to Fig. 1 as the size of the generally variable parameter h, which is called the smoothing length. Parameter h can also be multiplied by constant κ. Particles which are inside the support domain attributable to particle i are called neighbouring particles. If the resultant value of the product κh in each time step of the numerical simulation is the same, there can be the decrease in the number of neighbouring particles and thus also the decrease in the accuracy of the solution due the effect of excessive deformations (i.e. during the mutual divergence of the SPH particles). It is advisable to change the size of the support domain during the calculation in such a way that the number of neighbouring particles is constant. There are many ways to dynamically develop h so that the number of neighbouring particles remains relatively constant. In 1989, Benz [14] suggested a method of developing the smoothing length. This method uses the time derivative of the smoothing function in terms of the continuity equation where d is the number of dimensions and ∇ ⋅ v is the divergence of the flow. This means that the smoothing length increases when particles separate from each other and reduces when the concentration of particles is significant. It varies in order to keep the same number of particles in the neighbourhood. Eq. (7) can be discretized using SPH approximations and calculated with other differential equations in parallel [13].

Eulerian and Lagrangian Kernels
The approach in Eq. (7) is applicable when the integral representation of field functions is formulated in spatial coordinates (Eulerian kernel). With an Eulerian kernel, the smoothing length of a particle changes through the calculation. As a consequence, the neighbourhood of each particle needs to be updated at each time step [15]. However, nothing exists to prevent the number of neighbouring particles changing. Despite the implementation of Eq. (7), the particles can enter and leave the support domain and thus tensile instability can occur. In other words, this means that the possibility of simulating ductile failure during the excessive divergence of particles from one another disappears. The behaviour of an Eulerian kernel during a calculation can be seen in Fig. 2.

Fig. 2 Eulerian kernel -the amount of particles within the support domain might not be constant
However, when the integral representation of field functions is formulated in the material coordinates (Lagrangian kernel), the neighbours' list of each SPH particle is defined in the initial configuration and remains constant throughout the whole calculation. It means that the support domain of a particle follows material deformation in order to always keep the same neighbours. It provides a way of solving tensile instabilities [13]. The behaviour of a Lagrangian kernel during a calculation can be seen in Fig. 3. To eliminate the tensile instabilities in the simulations the Lagrangian kernel was selected. The value of parameter κ = 1.2 was used. In other words, the extent of the support domain was 20 % bigger than using Eq. (7).

Numerical Heterogeneity
Each SPH particle can be considered a Lagrange element with regard to the fact that the mass mj allocated to a particle moves together with the particle during the simulation. Also, it can be stated that mass mj acts in Eq. (6) as a weight coefficient. The higher the mj value, the more particle j is going to influence its surroundings. This information can be utilized very simply to create numerical heterogeneity [16].
Numerical heterogeneity can be considered to be an adaptation of the computational model in the sense of the modification of its numerical code or of the numerical method with which the simulation is calculated. Combination is also possible. The following approach can be considered a combination of both methods. It modifies the computational model as well as the numerical method. The modification consists in the introduction of the oscillation of the masses of SPH particles. However, this oscillation will directly precede the compilation of Eq. (6) for each i particle. In other words, the oscillation of the weight coefficients occurs directly. Other impacts may include the creation of virtual geometry in the background of the calculation; more information can be found in [16]; also see Fig. 4.
With regard to that, it can be stated that the introduction of the material heterogeneity into the tested models was performed via modifying the weight coefficient value of the individual SPH particles; in other words, the particle mass mj included in Eq. (6) was modified. Further, the same density value was preserved in all the particles, meaning that -with respect to the validity of Eq. (4) -the volume assigned to the individual particles was virtually modified (enlarged / reduced).
However, due to the necessity to maintain the best possible particle distribution regularity in the tested models (an identical smoothing length for all the particles), the initial particle distributions were generated for the state where all the particles are assigned the same volume ∆Vj. This operation is graphically represented in Fig. 4.

Mass Distribution and Experiment Set Up
As the individual SPH particles had to be assigned various masses or weight coefficients, it was necessary to build the entire distribution procedure upon the precondition of a constant weight of the resulting body -concrete block (target) in the tested scenario.
The procedure can nevertheless satisfy the conditions of any distribution function (including, for example, the normal, uniform, or Poisson distribution variants). For the actual testing, the uniform distribution of the occurrence of masses was selected. The reason for choosing uniform distribution consisted in testing the computational stability. The properties of the smoothing function enable the SPH method to suppress or highlight the heterogeneities contained in a computational procedure [13]. By another definition, the sensitivity of the computational stability in the presence of particles with almost zero mass (such as air pores in concrete) was tested.

Fig. 5 Diagram of the numerical experiment
A brief description of the simulated experiment is presented in Fig. 5. The scheme shows a rigid projectile striking a concrete block at the speed of 500 ms -1 . Rigid projectile was simulated with the Finite Element Method (FEM) and the concrete block with the SPH method. The concrete block was fixed on the edges. A penalty based contact algorithm was used for the interaction between FEM and SPH. As one of the preconditions for the functionality of the introduced procedure consists in regular initial particle distribution, the block was discretized with 10 particles in the depth and 90 particles in the height and width. The particles were arranged to form a regular, grid-based field, as it is also indicated in Fig. 6. The simulations involved 81 000 SPH particles in total and were performed via the LS-DYNA program [17]. The Continuous Surface Cap Model (CSCM) was chosen as the material model of concrete to be used [18,19]. Tab. 1 shows the parameters employed in the simulations. Fig. 6 demonstrates tested models with the representation of the masses. Before running a computation, different fields within the concrete block were selected where randomness consisting in uniform distribution of masses was taken into account.
Based on the type of field (size, shape, volume or algorithm) and a distribution function, the random masses are generated to create modified models according to the procedure in Fig. 4. More about the topic can be found in [20,21].

Numerical Simulations Results
The light-grey portions of the images invariably represent a crackless material, whereas the dark-red sectors stand for complete failures (main cracks mostly). All the results were captured at the same moment. Fig. 7 shows deformed models into which material heterogeneity (i.e. oscillation of material parameters) has been introduced. It is clear from Fig. 7 that when heterogeneity is introduced into models, they become prone to random failures. This demonstrates (inter alia) the functionality of the introduced procedure. If the loaded targets are examined in a greater detail, several conclusions can be drawn.

Fig. 7 Deformation of the heterogeneous models after direct impact of the rigid projectile
Even though the variant target a) was created with the idea of incorporating the oscillation of the masses of all SPH particles (hence the name completely random), symmetricity of failure can be observed in Fig. 7. This phenomenon occurs most frequently in homogeneous models -i.e. in models with no introduced heterogeneity. One explanation is that it was caused by the use of a uniform distribution of masses. As far as the functional aspect is concerned, once again a regular structure was actually createdheterogeneous when viewed from close up but homogeneous as a whole.
The variant targets b) and c) do not show any signs of symmetricity. This is mainly due to the fact that the oscillation of masses itself was only performed at several positions on the model (target) which were also random -generated with noise functions, see also [20,21]. Again, distinctions can be made between the individual types of failure. The random fields of generated variant b) were created in regions extending depthwise through the target. On the other hand, the random fields of generated variant c) were created in regions extending laminarly along the surface of the target. As a result, target b) has a tendency to split, with the formation of many major cracks, while in contrast, target c) tends to break at the interface of the transition of the laminar field of heterogeneities -in the same manner as lithospheric plates.
Even though the graphs of the speed of a rigid projectile in Fig. 8 are almost identical, which is mainly due to the use of small oscillation of masses, a large support domain size and the selected smoothing function, it can be concluded that certain procedures are not suitable. Using the procedure for the creation of target a), a random failure can be obtained which is symmetrical and nearly corresponds to the model without any introduced oscillation of masses -generally, it is not suitable. Using the procedure for the creation of target c), a random failure can be obtained, but unfortunately, laminar fields of heterogeneity only appear in large structures where they are created via a poorly executed compaction process -they are not suitable for small structures or models. Using the procedure for the creation of target b), a random failure can be obtained and simultaneously the shape of the heterogeneity fields can be justified as being the result of prefabricated production processes -it can generally be considered suitable.
From the practical point of view, the presented procedure is numerically stable, mathematically transparent and is able to simply produce various results. More about the topic can be found in [16,20].

Transition Layer
The presented procedure for the introduction of heterogeneity was based on the idea of the oscillation of masses or weight coefficients, which were applied to the whole model, or only in certain fields of the model. However, these fields were selected on the basis of mathematical functions, and not with regard to the structure of the material.
The whole procedure can only be improved by the application of the oscillation of masses at the interface between aggregate and cement binder -the transition layer, in the case of concrete. In other words, the process of the uneven solidification and hardening of cement binder on grains of aggregate can be simulated. The creation of the transition layer is very simple; a thorough description of this algorithm can be found in [20,21]. By way of illustration, Fig. 9 shows several examples of the structures of the transition layer of concrete material for various aggregate grains.

Conclusions
The paper discusses a simple procedure for introducing numerical heterogeneity into SPH-based models. Interestingly, such numerical heterogeneity originates in a real property of concrete-based materials. The principles of the relevant procedure and the conditions that, if satisfied, enable to avoid numerical cracks, which are unrelated to the numerical heterogeneity are also characterized. The whole procedure for the heterogeneity introduction of a material is based on the modification of the weight coefficient of individual SPH particles. In other words, it is not necessary to generate a great number of material models as is often the case with other methods -on the contrary, only one is sufficient. This greatly decreases the requirements for the initial stages of the calculation. Moreover, it is very easy to introduce various heterogeneity fields, which are related to any distribution function via the direct modification of the SPH source code. It can be said that another advantage of the presented method is its simplicity and transparency.
The entire principle is presented using a simple fracture mechanics test where a rigid projectile strikes a concrete block (target); the results clearly point to the functionality of the procedure. In subsequent research, the authors would like to investigate the generation of random fields with predefined distribution functions, which are variable in time. This would make it possible to consider phenomena related to the ageing of materials, potentially even with their cyclic loading.