Craters in concrete slabs due to detonation – drawbacks of material models with a Mohr-Coulomb yield surface

Numerical simulations have been performed with a commercial distributed explicit FE-solver and the results have been compared with experiments. High explosive was placed in front of different concrete slabs with the dimension 100 × 100 × 16 cm. Some of the results of the simulations, in particular the profile of the craters, are not in agreement with the test results. Therefore the key characteristics of the constitutive equation based on Mohr-Coulomb yield surfaces and a damage evolution linked to the plastic strain has been reviewed.


Introduction
This paper shows drawbacks of finite element simulations of concrete walls: Five different composite concrete walls following to Fig. 1 have been investigated and compared to experiments.
The amounts of explosive have been chosen in order to create a hole through the wall thickness (called "big" quantity hereafter) or to have undamaged rear surface (called "small" quantity hereafter).Since the corresponding amount was not known in advance, a huge number of tests have been performed and compared to numerical results.
The constitutive relation with parameters from in the literature and already implemented in the material libraries of commercially distributed finite element programs yield in general good results -even with isotropic material behavior!The amount of explosive or the kinetic energy of a penetrator [1] needed for perforation of concrete is predicted well.However, very simple empirical formulas based on experiments yield the same results.Unfortunately, the standard models built in commercial finite element programs did not correctly predict the cross section of a crater inside concrete.The diameters of the crater on the front and rear surface of the concrete can be calculated precisely, but the diameter of the hole inside the concrete wall was always much bigger than in reality.
This work gives a short introduction into constitutive relations which are used in commercial available finite element codes.The material models are based on Mohr-Coulomb yield surfaces and an isotropic damage parameter which is linked to the effective plastic strain.Emphasis was put on the key characteristic, in particular on the additive decomposition of the stress tensor into pressure and stress deviators as well as the dependence of the model on the fineness of the finite element.It will be shown, that the sensitivity of material parameters onto the dimension of the crater is completely different for each a Corresponding author: markus.conrad@saabgroup.comtype of composite concrete slab shown in Fig. 1.The introduction of an orthotropic damage parameter which is able to reproduce the orthotropic behaviour of a crack only improved the qualitative issues of a crater.Since the orthotropic damage is based on the principal strain, the result is still dependent on the fineness of the finite element.There is no convergence of the results by mesh refinement.On the contrary the material parameters have been optimized for element size of 5-10 mm.

Constitutive relations
Our FE-Programs uses an additive decomposition of the stress and strain tensor into dilation and deviation.The stress deviator is linked to the strain deviator by the shear modulus.Beyond the yield surface the total strain tensor is decomposed into an elastic strain and plastic strain.The compression is linked to the pressure by the bulk modulus.The equation of states allows for crushing.However, there is no cross link between deviator and dilation in the constitutive relation.There is no volumetric expansion or compression linked to pure shear loads.

Equation of State
The equation of state links the pressure to the compression.After a tiny elastic zone with reversible compression the concrete starts to crush (irreversible compaction) up to full compaction.During crushing, the bulk modulus varies from the initial to the fully compacted value.Autodyn's RHT-Concrete model uses the p-a equation of state [1], some concretes in LS-Dyna [2] (e.g.*mat concrete damage rel3) use look up table.In other application the bulk modulus is the first constant of a cubic equation of state.

Strength
The strength of the Material is described by three surfaces: Elastic Limit, failure surface and residual surface.Beyond the elastic limit, there are plastic deformations with strain hardening effect up to the failure surface.Hereafter the material strength softens up to fully damage where the residual surface is applied.The three yield surfaces are shown in Fig. 3 as a function to the pressure and the principal stresses.The residual surface is relatively close to the failure surface except in region with low pressure.However, exactly this low pressure region and the expanding regime is in our focus, i.e. −0.1 < p/ f c < 0.2.
The use of yield surfaces is established for ductile materials such as metals.The surfaces in Fig. 3 are based on Mohr-Coulomb's, Johnson-Cook's, or a similar law and take the brittleness of concrete in a threefold manner into account: (1) The residual surface is zero for negative pressure, i.e. broken concrete is expanded without any forces.(2) The elastic and failure surfaces allow for reaction forces in expanding finite elements.
(3) The contours of the surfaces are not circles but take a triaxial loading into account using a function of the principal stresses following to Willam and Warnke [1,2].The absolute value of the yielding is reduced for tensile stresses.
Malvar et al. [2] shows how the surfaces have to be adapted to static tests with compressive strength f c and tensile strength f t .For example, the absolute value σ m of the compressive failure surface is given by with r the relation of the compressive strength of a "new" concrete to the compressive strength of a concrete for which the material parameters a 0 , a 1 , and a 2 are validated.
Finally, the strength increases with high strain rate.

Damage
The damage is a parameter in volumetric finite elements which characterises the influence of cracks in brittle materials.The stresses σ are calculated on the bases of the nodal forces and the cross section A 0 of the element.However, an expanding crack reduces this cross section by its surface A R .As a consequence the local stresses in the element σ are increased as The damage is linked to the relation ( As a consequence, a fully damaged material does not anymore resist to expansion or shear.There is either a continuing softening between the failure and residual surfaces or the strength follows to the failure surface until a critical value of the accumulated damage (typically D c = 0.7) is reached.Beyond this value the material strength follows the residual surface.The first technique is called continuously softening and the latter is called spontaneous failure.

Damage evolution
The damage growth follows to an empirical formula.Each of the authors in the reference [1][2][3] proposes its own formula.They agree that the damage is linked to the plastic strain and total expansion (or compression).Figure 4 shows the damage as a look up table between effective plastic strain and pressure.The damage is given by the user as a function of a modified plastic strain λ with  With constants b 1 , b 2 , and r f which takes the compressive strength of a "new" concrete into account.Malvar et al. [2] put a lot of emphasis on the damage evolution:" If the analysis yields a different localization width than anticipated, this should be corrected and the calculation restarted.These parameters will be of importance when the structure analyzed is lightly reinforced or is tension-or shear-critical." Malvar et al. [2] describe the strain hardening also by the damage parameter: For very small plastic strain, there is a continuous transition from the elastic limit to the failure surface.Since the strain hardening effect is accomplished at really small effective plastic strains a simplified strength algorithm without separate elastic limit and strain hardening gives nearly the same results: The energy accumulated during strain hardening is small compared to the total energy needed for perforation of concrete walls.

Localisation
The energy consumed by the crack growth should be equal to the energy dissipated by the damage and the plastic deformation.Lemaître et al. [3] showed that this has an impact to the correct mesh size: The energy c released by crack growth is proportional to the increase of the crack surface and equal to the energies released by plastic deformation p and damage growth D which are proportional to the volume involved.A characteristic length l 0 results by equating these two energies, i.e.
With ε u the ultimate strain at rupture, G c the fracture energy release rate, f t the tensile strength and D c the critical Damage.Lemaître et al. [3] showed that the characteristic length in concrete is in the order of 10-100 mm, depending on f t , ε u and Young's modulus E.
In case the mesh size is refined, the parameters for the damage evolution law have to be adapted.The simulation results are mesh dependent, unless some average value of the damage over a region of the characteristic length is used.For a fine mesh, the damage is concentrated in few heavily distorted elements.The volume involved is small and thus the total energy dissipated by plastic deformation and damage growth is much less than in reality.Several techniques to overcome this drawback have been developed and published in the literature.Peerlings et al. [4] compare the different averaging techniques with different weight functions with the gradient enhance method.In the present work, the simplest method has been chosen: The parameters for damage evolution have simple been optimized for mesh sizes between 5-10 mm.

Anisotropy
A fracture reduces the strength of a structure only for tensional and shear loads which have components along the normal to the fracture surface.Compressive loading with closed cracks are equal to loadings without any cracks.In addition tensional loads parallel to the crack surface are also not altered by the presence of crack.

Damage evolution
The anisotropic influence of cracks in a structure is taken into account with a second order damage tensor D. Following to Lemaître et al. [3], the damage grows proportional to the strain rate only in the direction with positive principal strain. with and the material constants A and a.The operator ":" is the scalar product of two tensors: x:x = x i j x i j with Einstein's convention of summation over repeated indices.The symbol • is the MacCauley bracket.
Finally the threshold of damage growth does not simply depend on a ignition limit strain ε lim 0 , but also on the progress of the fracture (damage).Damage growth starts if the strain exceeds a critical value ε crit .

Local stress tensor
The local stress σ is increased due to that fact, that the crack reduces the original cross section.With A 0 the cross section normal to the direction t of the traction and A R the surface of the crack projected onto A 0 the local stress tensor is derived from the global stress tensor σ as

EPJ Web of Conferences
The direction of the normal vector of the crack surface is generally not parallel to the direction t of the traction.The term A 0 /(A 0 − A R ) is consequently a second order tensor and in principle the inverse of the damage tensor D. A formal relationship taking also the direction of the traction and the normal vector of the crack surface into account is given by Voyiadjis et al. [5]: The components of the fourth order tensor M consists on terms of the components of the damage tensor only.Details are given in the original reference [5]: Using the Voigt's notation for the stress tensor σ = (σ x , σ y , σ z , τ xy , τ xz , τ zx ) the fourth order tensor M is transformed to a second order tensor as:

Results using models with isotropic damage
The results of the finite element analysis were only roughly precise: The amount of explosive for creation of a hole is predicted quite well even with isotropic material models.
Even the biggest diameter of the craters measured on the front and rear surfaces are close to the reality.However the cross section of the crater in the model was wrong not only in quantity but also in quality.The hole inside the concrete wall is much smaller than the diameter of the crater on the front surface.In case the damage is taken as a measure for the crater, the calculated diameter of the hole is close to biggest diameter of the craters!This is a consequence of the localization and its characteristic length l 0 which is in the order of 10-100 mm.The material parameters and mesh size are linked such that the damaged is smeared over a rather big region.Let us discuss on behalf of Fig. 5 the nonlinear synergetic effect of the mesh size and the solver used for model set-up: The craters of the two figures in the middle have almost the same dimension (Ø 12-16 mm) although the quantity of the explosive is different by factor two. Indeed the model with the fine mesh is "soft".However, the damage zone is expected to increase not proportional to the quantity of the explosive.On the contrary, following facts are evident: The detonation shock pressure inside the explosive is constant and a characteristic of the explosive.The mass of the explosive is proportional to the volume and the momentum applied to the concrete is roughly proportional to the contact surface between explosive and concrete.
The explosive and its detonation have been modelled with an Eulerian mesh fine enough to resolve the over pressure of the shock pressure which is reflected by the concrete all.The concrete of the models on the outer left and outer right in Fig. 5 is modelled with Lagrangian mesh, of course including the reinforcement bars.Autodyn has a built-in algorithm which transfers the momentum from the Eulerian mesh to the Lagrangian mesh and viz.versa.Unfortunately, there is no such implementation between an Eulerian mesh and a SPH part.However there is a work around: The momentum is first transferred from the Eulerian zone to a thin foil modelled with Lagrangian elements.This foil is attached by a contact algorithm to the SPH-part.This twice fold transfer from the detonation shock in the Eulerian mesh to the SPH-solver dissipates numerically more energy than the direct transfer to the Lagrangian mesh.As a consequence, the momentum transferred to a concrete wall modelled with smooth hydrodynamic particles is less than the one transferred to concrete wall modelled with Lagrangian elements.

Parameter study
The influence of three settings in the models has been investigated fully factorial as follows: • Element length: 5 mm or 10 mm • Solver: Lagrangian or Smooth Particle Hydro-dynamics (SPH) • Damage growth rate: "quick" or "slow"; "quick" means, that the element is fully damaged at small plastic strain such that the residual strength is applied already in an early state of the simulation.
The pairs of the variation above have a weaker (5 mm element length, Lagrangian Solver, quick damage growth rate) or a stronger (10 mm element length, SPH Solver, slow damage growth rate) resistance to perforation.Simulations with weak parameters are expected to yield bigger dimensions in craters as simulations with stronger parameters.Unfortunately some simulations yielded the opposite behavior.The fully factorial evaluation showed the influence of the parameters onto the characteristic of the craters was different or even contrary for each configuration shown in Fig. 1.The evaluation of the performance is based on quantities such as diameter of the crater and zone of delamination of the different simulation are within ±15% of the measurements.One example on the basis of the configuration in the middle in Fig. 1: A quick damage growth rate should yield an increased diameter of the crater compared to the one of a simulation with a slow damage rate.The same was anticipated for the mesh refinement.This was true for small quantities of the explosive (called "small" hereafter) which were not enough for the perforation of the wall.In contrast, the anticipation did not correspond for quantities of explosive (called "big" hereafter) needed for perforation.To make matters worse, a simultaneous refinement of the mesh and quick damage growth rate yielded an opposite result, i.e. a smaller diameters of craters for the "small" quantity of explosive.Finally a combination of coarse mesh and slow damage growth rate yielded a smaller diameter in crater for the "big" quantity of explosive, as expected by assuming positive synergetic effects.

Composite normal concrete with UHPFRC
A delamination of a composite concrete slab occurred not exactly between the interface of the normal concrete and the UHPFRC but inside the normal concrete.The failure occurred in typically 1-2 cm below the interface such that the rear plane of reinforcement bars is still covered by normal concrete.The delamination is not planar but there are concentric rings detected as shown in Fig. 6.
The simulation gives some hints how the concentric rings emerge.Even the dimension and location of the delamination of some simulations was in quite good correspondence with the experimental results.

Influence of anisotropy in normal concrete
In the beginning of Sect.3.1 has been shown that profile of the crater of the simulation are even qualitatively wrong.This drawback is mainly observed for the first two types  On the left hand side is shown an example with a Lagrangian mesh and a small quantity of explosive which did not create a hole thru the complete wall.On the right hand side is shown an example with SPH model and a big quantity of explosive far above the quantity needed for perforation.The solid and dotted lines show the profile of the craters of corresponding field test.The simulation shows the damage between 0 (black) and 0.7 (white) which corresponds to the critical damage D c . of concrete walls shown in Fig. 1.The main reason is assumed to be due to the characteristic length l 0 lined to the damage parameter.In this section the improvements due to introduction of an orthotropic damage parameter will be show on the basis of Fig. 8.There is some qualitative improvement of the profile.However a good set of parameters (damage growth rate, critical damage D c for spontaneous failure, mesh size, Mohr-Coulomb yield surfaces) have not been found yet.Since the list of parameters is so long, a set of parameters optimized of the basis of some tests may yield wrong results for another test setup.In addition, a result of the fully factorised parameter 04019-p.5 study is the dependency of the numerical results on the combination of parameters in a way which is not at all intuitively predicable.

EPJ Web of Conferences
The example above shows the significant the difference of a simulation with continuous decrease of the strength from failure to limit yield surface (with increase of the damage from zero to one) compared to a simulation with spontaneous failure with only two distinct state of the strength: The accumulated damage within the "strong" sate is below a critical value and the strength is defined by the residual surface.After the accumulated damage exceeds a critical value the material response is "weak" and the strength follows the residual surface.As a consequence, the first principal stresses, which result from the reflection of the shock wave and are detected close to the rear surface and in a distance above 16 cm from the centre, are higher, because the strength material in the centre is hold on the failure surface for a "quite long" duration.Therefore, the last model in Fig. 8 predicts on the one hand a damaged zone just under the rear surface in a distance beyond 16 cm from the centre.On the other hand, the diameter of the hole in the concrete is closer to the reality.However the shape of the crater in the simulation is "distorted".

Conclusions
It has been shown that the results of the simulations with commercial available programs with standard parameters for the constitutive relation of concrete yields generally good results.However, the carter profile of the simulation does not correspond to the reality.The introduction of orthotropic damage improved the results only slightly.The finite element models are dependent on the meshing because the damage evolution is linked to the plastic strain.Moreover, it is quite cumbersome to define a right set of parameters not only because of numbers of parameters but also due to their synergetic effects such that the relationship with the simulation results is highly nonlinear or even chaotic.
E. Cadoni and G. Riganti [6] are developing a constitutive equation for concrete under highly dynamic loading which is not based on a damage parameter.They reconsider to the definition of the damage parameter and calculate the fracture surface using a crack propagation speed, a defined density of nucleation, and a limit load at which a crack starts to evolve from a nucleation.Spontaneous failure starts if the local stress tensor based on formula (8) exceeds a limit.In principle, such a model does not depend on the fineness of the mesh.Therefore, the results shall converge with mesh refinements.
This study has been financed by the arma suisse and is a part of a long term research project ARAMIS 042-39.The author was always invited to participate to meetings and field tests.There was free access to internal reports, in particular to the evaluation of modified Split Hopkinson Bar tests (6) and field tests (6) sketched in Fig. 1.

Figure 1 .
Figure 1.Five different concrete walls with high explosive attached on the front surface.All walls are 16 cm thick except the ultra high performance fibre reinforced (UHPFRC) wall which is only 6 cm thick.The normal concrete is double reinforced, i.e. a mesh of steel reinforcement bars are placed approximately 3 cm below each the surfaces with dimension 100 cm × 100 cm.

Figure 2 .Figure 3 .
Figure 2. The pressure is a function of compression.The dashed lines -which slopes are the bulk moduli -define the volumetric crushing.

Figure 4 .
Figure 4.The damage is a function and plastic strain ε p .The modified plastic strain λ is a function of the plastic stain and pressure.Depending on λ, the strain hardening (transition from elastic limit to failure surface) and softening (transition from failure surface to residual surface) effect is evaluated.

Figure 5 .
Figure 5. Crater in UHPFRC following to the second configuration of Fig. 1 in the scale 1:5.On the left hand side are shown a Lagrangian model and an SPH-Model with 5 mm element size and a small quantity of explosive which did not perforate the wall.The dotted and solid lines are the crater profile of the field tests.On the right hand side are shown the SPH and Lagrangian model with 10 mm element size and the double amount of explosive needed for perforation.The models shows the damage between 0 (black) and 0.7 (white) which corresponds to the critical damage D c .

Figure 6 .
Figure 6.The delamination is approximately 1-2 cm below the interface between the double reinforced normal concrete and the UHPFRC as shown with the dashed line marked by two arrows on the right hand side.The photograph on the right hand side shows the delamination surface with three concentric rings as hills and valleys.The white dashed circles mark the valleys.

Figure 7 .
Figure 7. Crater and delamination in composite concrete slabs.This is the configuration in the middle of Fig. 1 in the scale 1:5.On the left hand side is shown an example with a Lagrangian mesh and a small quantity of explosive which did not create a hole thru the complete wall.On the right hand side is shown an example with SPH model and a big quantity of explosive far above the quantity needed for perforation.The solid and dotted lines show the profile of the craters of corresponding field test.The simulation shows the damage between 0 (black) and 0.7 (white) which corresponds to the critical damage D c .

Figure 8 .
Figure 8. Crater in double reinforced normal concrete corresponding to the first configuration in Fig. 1 in the scale 1:5.On the left hand side is a Lagrangian model with isotropic damage parameter.The models in the middle and on the right hand side are examples with orthotropic damage parameter.Shown is the first principal damage from 0 (black) to 0.7 (white) which corresponds to the critical damage D c .The curves (solid, dotted, and dashed lines) are the profiles of the craters measured after corresponding field tests.
P xx P yy P zz − D 2 yz P xx − D 2 xz P yy − D 2 xy P zz − 2D xx D yy D zz and D i j the components of the orthotropic damage tensor and P i j is chosen such that D i j + P i j = δ i j (with Kronecker symbol δ i j ).The sub-matrices M 11 , M 12 , M 22 are: