Using the Abaqus CDP model in impact simulations

The understanding and assessment of accidental crash scenarios on reinforced concrete structures is of high interest for safety issues. Although experimental research on this topic has already a long and successful history, there are many issues waiting to be solved as far as numerical simulations are concerned. In this article some particularities of the Abaqus Concrete Damaged Plasticity (CDP) model are investigated with the purpose of using efficiently the CDP model in impact loaded reinforced concrete structure simulations. In particular, the sensitivity of the simulation response with respect to model parameters and element size is studied. The simulation response is compared to measurements from benchmark impact tests on reinforced concrete plates.


Introduction
On the role of simulations in structural integrity assessments In the life cycle of critical infrastructures, such as nuclear power plants, large civil engineering structures and industrial complexes it is often the case that the original life span needs to be extended and the original design criteria need to be upgraded to meet the current requirements. These requirements may change due to changes in the society, climate or geology or due to increased risk awareness and public opinion. For example, in the context of Finnish nuclear industry, the nuclear safety and security regulations (YVL), [1], the criteria are subdivided into Design Basis Criteria (DBC), which include design for normal operation and anticipated accidents, and Design Extension Criteria (DEC), which include rare and extreme external hazards. New built power plants must meet both DBC and DEC, but there is also strong pressure for old plants, which seek for long term operation extension, to prove that design extension criteria are met, as well.
From the simulation point of view, the design basis criteria can mostly be addressed with linear static and dynamic structural analysis, and there is a strong confidence within the engineering community that the linear methods produce valid results within their domain of application. On the contrary, design extension criteria often involve non-linear phenomena and simulations that are run beyond the failure point are needed to assess structural integrity. Airplane crash simulations on nuclear power plants can be stated for example, as illustrated in [2,3]. The problem with full scale structural simulations involving non-linear material models, contact algorithms and possibly runtime changes in mesh topology, is the lack of confidence due to scarcity of validated benchmarks at various structural scales.
The present work aims at remedying this situation by providing validation data on hard missile impacts on reinforced concrete slabs. The material model to be validated is the Abaqus Concrete Damaged Plasticity (CDP) model, [4], customized with the user field extension and element removal algorithm originally proposed in [5,6] and later enhanced in [7]. The material model is incorporated in several benchmark impact simulations that are compared to the respective experimental results obtained from the VTT international share benefit IMPACT programme, phases I, II and III and the Finnish national SAFIR 2 programme. In particular, the tests A11 and A12 from IMPACT I and II, respectively, and E4 and E5 from SAFIR ERNEST 3 are used, [8,9,10].

Concrete behavior under dynamic loading
Dynamic behavior of concrete, whether caused by monotonic or cyclic loading, is different from the quasi static behavior. In aircraft impact, which is the main application domain of this study, strain rates in the concrete varying from 1 to 100 s −1 , are typically observed. However, strain rates up to 10 6 s −1 can be observed in ballistic impact and blast loading. In severe earthquake loading the strain rates are lower, but the loading is cyclic. Studies on the dynamic increase factor (DIF) of plain concrete, which are performed using the Split Hopkinson Pressure Bar (SHPB) both in compression, [11,12], and in tension, [13], show a noticeable increase in the DIF for strain rates higher than 1 s −1 . A dynamic increase in the fracture energy, too, can be observed in SHPB tests, [14]. Dynamic increase for lower strain rates has also been observed in direct tensile tests, [15], with strain rates varying from 10 −6 to 10 −1 1/s. Figure 1 shows a DIF increase both in compression and in tension.
(a) Compressive dynamic strength increase [11] (b) Tensile dynamic strength increase [16] [6]. From the Abaqus CDP model field dependency implementation point of view, there is one particular physical phenomena that is of primary interest. In uniaxial compression, the triaxial stress state resulting from a confinement due to the inertia effects has been acknowledged as the main reason behind the dynamic increase of the peak stress, [17]. Confined uniaxial peak stress and corresponding strain data has been collected from two sources, [18,19], and plotted in Figure 2. Notice a linear dependency of the peak stress, f cm , on the confinement ratio and a quadratic dependency of the corresponding strain, c 1 . The Eurocode proposes a piecewise linear dependency for the peak stress dependency on the confinement ratio, [20], with a threshold at 5% confinement ratio. (b) Concrete strain at peak ( c 1 ) Figure 2: Concrete compressive strength and strain dependence on confinement stress according to Sfer, D. et Al. [18] and Vu, X.D. [19] The confinement stress dependency of concrete hardening in a general stress state is a necessary ingredient in any hard missile impact simulation. Without this extension of the standard CDP material model, the concrete in front of the missile is simply too weak. The simulations described in [6] show that the concrete plug that is formed in front of the missile head has a very high confinement ratio due to inertia effects.

Element removal in hard missile impacts
In the context of finite elements, an important topic of discussion that frequently comes forth in publications on impact simulations is element removal, which is used to materialize macroscopic concrete fracture, the formation of the shear cone as well as spalling and scabbing at the front and back surfaces. In many simulations, brittle behavior of the material and large rigid body motion of the shear cone is not possible to simulate without proper element removal algorithm. Also steel reinforcement tensile breaking is materialized using element removal, but this type of element removal algorithm is based on different algorithms, for example the Johnson-Cook damage initiation and evolution, [4]. Although there are other numerical solutions for modeling crack formation and growth in the finite element context, such as XFEM, adaptive re-meshing to follow crack propagation or using bond interfaces between elements, these more sophisticated tools are too costly to be used in full scale impact simulations. Discrete Element Method (DEM) based approaches, [21], as well as meshless finite element method approaches are also used, but then the question of adequate reinforcement and steel-concrete bond modeling becomes an issue.
The element deletion algorithm essentially consists in removing elements in the correct place and at the correct time. After the element is removed, the contact algorithm updates the contact surfaces to include also facets from the elements adjacent to the removed one. In [22], the equivalent plastic strain in compression (PEEQ) is monitored for each element with the CDP model. When the arbitrarily chosen cutoff value of 0.2 is reached at a given integration point, the corresponding element is removed from the simulation. The justification of this choice is primarily numerical, since the explicit central difference time integration algorithm proposed by Abaqus, fails if an element undergoes too large deformations in a time increment. 4 In contrast, [23] proposes an element deletion criterion where the cutoff value of the equivalent plastic strain depends on the confining pressure. This seems to be a natural choice, since the hardening evolution of the cohesion stress in compression depends on the confinement rate. In [6], the proposed element removal criterion is based on a logical "or" operation of two conditions: when, in any integration point, the cutoff value of the equivalent plastic strain in compression (PEEQ) is reached, or, when the cutoff value of the equivalent plastic strain in tension (PEEQT) is reached. The cutoff value of the equivalent plastic strain in compression is a function of the confinement rate, whereas the cutoff value of the equivalent plastic strain in tension is a function of the maximum principal strain rate.
One of the drawbacks of all of these element removal criteria is that the element removal does not depend on the actual stress or strain state, only on the value of the internal hardening parameters. Therefore, in [7], an element removal algorithm based on actual shear strain values was proposed. One of the targets in this work is to test the proposed shear strain based element removal algorithm against a larger collection of benchmark tests.

Description of the user extended CDP model
The built-in Abaqus CDP model is an elastic-plastic model with optional scalar stiffness degradation accumulation. The yield surface is given by the so called "Barcelona" yield function, shown in Equation (1), and isotropic strain hardening is defined with respect to two independent internal hardening variables, describing tensile and compressive behavior, respectively. Historically, the "Barcelona" yield function was first proposed in [24], but the isotropic hardening behavior that is based on compressive and tensile equivalent plastic strain was introduced in [25,26]. The yield surface can be described as a modified Drucker-Prager cone with a Rankine tension cutoff. The yield function shape is therefore determined by both the compressive cohesive stress,σ c and the tensile cohesive stress, σ t . In the effective stress space 5 , the yield function can be written as follows in terms of the effective Mises stress,q, the effective hydrostatic stress,p and the effective maximum principal stressσ max : The material parameter α is determined from compressive equibiaxial tests. The material parameter γ = 3(1 − K c )/(2K c − 1) controls the magnitude of the corrective term, as shown in Figure 3. The effect K c can be seen in compressive stress states, in particular in the confined uniaxial compression state. The flow rule used in the CDP model is non-associative and the direction of the flow is obtained from the outward pointing normal of a convex flow potential, which is a smooth Drucker-Prager hyperboloid (i.e. a hyperbola of revolution around the hydrostatic axis). The eccentricity of the hyperbola, eσ t 0 , controls the flow behavior on the positive (tension) side of the hydrostatic axis. The angle of the hyperbola (i.e. the angle between the hydrostatic axis and the hyperbola asymptote) is called the angle of dilation, φ, and it controls the amount of plastic volumetric dilation that a material point undergoes when subject to shear deformation.
The compressive cohesive stress,σ c and the tensile cohesive stress,σ t , evolve with respect to the scalar internal hardening variables, the equivalent plastic strain in compression, p c , and the equivalent plastic strain in tension, p t , respectively 6 . The rate of the equivalent plastic strain in compression is proportional to the minimum principal plastic strain rate, whereas the rate of the equivalent plastic strain in tension is proportional to the maximum principal plastic strain rate. The coefficients of proportionality depend on the current stress state at the material point.
Although the built-in Abaqus CDP model is not able to realistically describe impact simulations on concrete, it is possible to enhance the built-in model by user-defined field variable dependencies. In this study, field variable dependencies based on the confinement stress dependency of the compression hardening and strain rate dependency of tension softening is proposed.

Compression hardening with confinement dependency
By substituting a uniaxial confined stress stateσ max =σ mid = −σ cnf and |σ min | ≤σ cnf in the yield function, (1), and considering the additive elastic-plastic strain rate decomposition the following expressions for axial effective stress and total strain are given: The parametric plot obtained from Equations 2a and 2b is the one that has to be matched with the confined uniaxial experimental test results on concrete cylinders. Further, consider a specific hardening evolution in compression, for example a polynomial of exponentials, as suggested in [26]: By fixing the initial effective yield stress, as suggested in the Eurocode, toσ c 0 = 0.4 f cm , the material parameters a c and b c can be computed as a function of the confinement ratio, CR =σ cnf /f cm . The mean concrete cylinder strength at current confinement can be computed as the product of the unconfined mean concrete cylinder strength multiplied by a confinement increase factor. The Eurocode, [20], suggests that the confinement increase factor on the mean concrete cylinder strength, ζ, is a piecewise linear function of the confinement ratio: ζ(CR) = 1.0 + 5.0 CR for CR ≤ 0.05 and ζ(CR) = 1.125 + 2.5 CR for CR ≥ 0.05. The confinement increase factor on the strain is equal to the confinement increase factor on the mean concrete cylinder strength to the second power, c 1 (CR) = ζ 2 (CR) c 1 . Then, by substituting the Eurocode 2 values for mean concrete cylinder strength and the corresponding total strain at current confinement ratio in Equations 2a and 2b, the expression of the equivalent plastic strain when the peak stress is reached is given as a function of the current confinement ratio: For shorthand, consider the notation µ(CR) =σ c p The material parameters a c (CR) and b c (CR) are then given by Equations (5) and (6).
Equation (5) leads to the constraint µ ≥ 1, which gives the theoretical upper limit for the confinement ratio, CR UL . This upper limit for the confinement ratio is computed from 2.5 ζ(CR UL ) − 1+2α+γ 1−α CR UL = 1. The value of upper limit, CR UL , depends on the parameters α and γ, which themselves depend on the ratios σ b 0 /σ c 0 and K c , respectively. For example, for parameter values σ b 0 /σ c 0 = 1.15 and K c = 0.8, the upper limit CR UL tends to infinity, which means, in this specific example case, that the expressions given in (5) and (6) are valid for any confinement ration value. Figure 4 shows the evolution of the  The confinement ratio field variable is defined in the USDFLD/VUSDFLD user subroutine. It is assumed here, that the confinement stress is a weighted average between the current value of the maximum principal stress and the maximum value of the maximum principal stress with respect to time. In other words, we compute the confinement ratio at the beginning of each time increment as per Equation (7) in the user subroutine. Here, the weight coefficients 0.4 and 0.6 are arbitrarily chosen.

Tension softening with strain rate dependency
It has been shown in many experiments that the behavior of concrete in uniaxial tension is sensitive to the size of the specimen, [27,28]. Following the procedure indicated by Hillerborg and his co-authors, [29,30,31], it is assumed that ahead of the (micro)crack tip in a material RVE loaded in tension, there is a fracture process zone. In this zone, the gradual coalescence of micro voids and other defects in the concrete matrix induce propagation of the actual crack tip when tension load increases. On the other hand, the fracture process zone is capable of transferring small stresses, which diminish as the crack opening increases and the crack tip moves forward. The established way to describe the fracture related phenomena in continuum mechanics is material softening. The CDP model enables both isotropic strain softening and scalar stiffness degradation as means to introduce softening in the model, but in this context only isotropic strain softening is used. It is assumed that the fracture energy needed to transform a fracture process zone per unit of surface into an actual crack is a material constant, denoted G F , and has the units of N/m. Therefore, the fracture energy is defined as the area under the stress -displacement curve in the fracture process zone. Since the built-in Abaqus CDP model is a local continuum formulation, the only way to attenuate the effect of mesh dependency in numerical computations is to introduce the fracture energy as a material parameter. Although this approach is theoretically unsatisfactory, since it does not solve the issue of ill-posed initial value and boundary value problem, it nevertheless provides a practical solution to circumvent mesh size dependency. Values for the fracture energy can be computed for various concrete grades, as shown in [32]. Lee and Fenves, [26], suggest that the shape of the softening stress strain curve is exponential.
whereσ t 0 and b t are material constants, which may depend on the strain rate. Introduction of the dynamic increase factor, DIF f , which multiplies the Eurocode concrete mean tensile stress at quasi-static rate, f ctm , is one way of expressing strain rate dependency. Hence, one may define the strain rate dependent peak stress in uniaxial tension asσ t 0 = DIF f f ctm . As shown in [14], the fracture energy is a rate dependent quantity. It can therefore be computed from the expression Computation of the value of the integral leads to the following value for the parameter b t : In Abaqus the user may enter cohesive tensile stress input data as a function of plastic displacements instead of plastic strains. Then, by defining the plastic displacement w p t = l ch p t , the cohesive stresses are computed from the relationσ Figure 5 shows the evolution of the cohesion stress in tension as a function of the equivalent plastic strain/displacement in tension for various strain rates (Figures 5a and 5c) and the corresponding parametric stress-strain/displacement plots (Figure 5b and 5d).
The dynamic increase factors DIF f and DIF g are assumed to be functions of the relative tensile strain rate, SR =˙ max /˙ QS max , where˙ QS max = 10 −6 1/s denotes the quasistatic strain rate. Consider, for instance, piecewise linear functions on the logarithmic scale. In these computations we have used DIF f = 1 + 1 12 log(SR) for log(SR) < 6 and DIF f = −15 + 11 4 log(SR) for log(SR) ≥ 6. Likewise, we have used DIF g = 1 + 5 6 log(SR) for log(SR) < 6 and DIF g = −60 + 11 log(SR) for log(SR) ≥ 6. The strain rate field variable is also defined in the USDFLD/VUSDFLD user subroutine. It is assumed here, that the strain rate is a weighted average between the current value of the strain rate and the maximum value of the strain rate with respect to time. In other words, we compute the strain rate at the beginning of each time increment as per Equation (10) in the user subroutine. Here, the weight coefficients 0.4 and 0.6 are arbitrarily chosen.
The element removal algorithm From a physical point of view, element removal in impact simulations is a means of materializing macroscopic fracture of concrete in the continuum finite element framework. From a numerical point of view, element removal enables necessary topological changes in the finite element mesh so that small body deformations but large rigid body motion can be successfully computed by the explicit time integration. In the studies done in [23,22], the element removal algorithm is based on the evolution of the internal hardening variable p c . This choice is clearly the simplest one. It is motivated by the fact that in pure shear behavior the equivalent plastic strain in tension, p t , rapidly plateaus to relatively |σ max | (MPa)ǫ max = 1 · 10 −6 1/ṡ ǫmax = 1 · 10 −3 1/ṡ ǫmax = 1 · 10 +0 1/ṡ ǫmax = 2 · 10 +0 1/ṡ ǫmax = 1 · 10 +1 1/s (b) Axial stress -total stain |σ max | (MPa)ǫ max = 1 · 10 −6 1/ṡ ǫmax = 1 · 10 −3 1/ṡ ǫmax = 1 · 10 +0 1/ṡ ǫmax = 2 · 10 +0 1/ṡ ǫmax = 1 · 10 +1 1/s (d) Axial stress -total displacement Figure 5: Stress-strain and Stress-displacement plots in uniaxial tensile state (C60 concrete, solid: l ch = 100mm, dashed: l ch = 10mm) small values whereas the equivalent plastic strain in compression continues to grow along with the shear, as shown in [7]. However, as it was shown in the same study, a more adequate element removal algorithm can be obtained by basing the element removal on shear behavior.
In this shear behavior based removal algorithm, it is assumed that in a concrete shear band, on the boundary of the shear cone, the strain state is close to the pure shear strain state. Therefore the octahedral strain (multiplied by √ 2) is used as a measure of the shear strain. This choice is motivated by the fact that in a pure shear strain state = diag(1, − 1 2 , − 1 2 ) , the value of the octahedral strain is oct = / √ 2. To rule out strain states that are not close to the pure shear state, but still have large octahedral strain values, the measure that gives the pure shear content of an arbitrary strain state is given by a combination of the following coefficients: The pure shear strain content of an arbitrary strain state can thus be defined as follows: Hence, by setting the condition psh ≥ cutoff max whenever the the material point is not properly confined, i.e. when CR ≤ Γ cnf CR UL (K c = 0.75), we can assume the element at that material point can be removed from the finite element analysis. The "safety coefficients" Γ psh and Γ cnf are provided here to prevent unduly early element deletion. The coefficient Γ psh can be chosen by trial and error method in such a way that the element removal algorithm produces global responses that are mesh size independent. In this study we have chosen this coefficient to be where the characteristic length l ch is given in meters. Recall also, that the cutoff strain, cutoff max , depends on the value of the dilation angle, as shown in [7]. The coefficient Γ cnf is obtained by observing the size of the high confinement zone in the beginning of the simulation. These observations lead to the choice of the value Γ cnf = 0.75.

Description of the benchmark experimental tests
The benchmark experiments that the simulations are validated against are hard missile impacts on reinforced concrete slabs. The missile is a stainless steel tube filled with nonreinforced light weight concrete and it is shot against the target slab using pressurized air. The target slab is held in place by a steel frame, which can be assumed stiff compared to the stiffness of the concrete slab. The free span of the target slab is 2m by 2m, and the steel frame connection to the slab can be assumed to be a hinged one. Detailed descriptions of the experimental setup can be found in [9].  As shown in Figure 6, four distinct tests, the E4 and E5 tests 7 from SAFIR ERNEST project, the A11 and A12 tests from VTT IMPACT phases I and II project are considered here as candidates for model validation purposes. The expected main failure mechanism is local shear occurring at the boundary of the plug. In addition, front surface spalling and back surface scabbing of concrete occurs together with concrete tearing at back surface caused by the movement of the reinforcement mesh. The missile and slab parameters for each test can be found in Table 1. The main test output values are summarized in Table  2.

Description of the benchmark simulations
The benchmark simulations are carried out using Abaqus with explicit central difference time integration. A quarter of the model is considered, as shown in Figure 7. The concrete slab is bounded by steel channels in order to withstand the contact pressure from the steel rods on both sides of the slab. These steel rods are attached to the frame that holds the slab in place. In these simulation models, the steel frame is not modeled at all. The steel rods, which are assumed to simulate hinged boundary conditions, are modeled as rigid bodies and they are rigidly fixed to the frame of reference.  The concrete of the slab is modeled using eight-node reduced integration elements with Abaqus default distortion and hourglass control (C3D8R in Abaqus nomenclature). Mesh sensitivity analysis is conducted with two mesh sizes: one with characteristic element length of 10mm (25 elements over the thickness of the slab) and another characteristic element length of 5mm (50 elements over the thickness of the slab). The steel reinforcement is modeled using two-node linear beam elements (B31 in Abaqus nomenclature) that are tied node-to-node to the C3D8R elements. Table 3 summarizes the essential input data used in the simulations.

Results
The results of the simulations for the marker quantities (penetration depth, residual velocity 8 , maximum deflection at mid-slab, shear cone angle 9 , tunneling depth and weight of spalled concrete) are presented in Table 4. These quantities are all evaluated at a time frame sufficiently far from the impact time frame to be considered as "permanent" values. Since the shear cone and tunneling depth values are measured geometrically from the deformed element meshes, there is an important part of subjectivity in the readings. The marker quantities can be used to investigate the sensitivity of the simulations with respect to a change in a given parameter. On the other hand, the marker quantities from simulations can also be compared to the marker quantities measured in the experiments (see Table 2) in order to assess the validity of the proposed model. In addition to the marker quantities, the missile residual plots are shown in Figures  8, 9, 10, 11 for tests E4, E5, A11 and A12, respectively. Field output plots at selected time frames are shown in Figures 12 and 13 for the test E4, in Figures 14 and 15 for the test E5, in Figures 16 and 17 for the test A11 and in Figures 18 and 19 for the test A12. These field output plots give an indication on the progress of the failure mode.

Discussion
This study has two objectives. First, it aims at assessing the sensitivity of the simulations with respect to a given parameter variation. The parameters considered here are the angle of dilation and the element characteristic length. Second, the study aims at validating the proposed concrete material model by comparing the simulation results to the experimental 8 The residual velocity is the velocity of the missile after penetration of the slab 9 The shear cone angle is defined as the angle between the axis of the shear cone and its generatrix test results. Both the sensitivity analysis and the model validation are carried out by considering each marker quantity separately.

Mesh size sensitivity
The chosen element mesh sizes are motivated as follows. For the coarse mesh, l ch = 10mm is a maximum element size above which various local physical phenomena can not be described accurately. For the dense mesh, l ch = 5mm is a minimum element size from computational cost point of view. An inspection of the missile residual speed plots (Figures 8, 9, 10, 11) tells that the missile residual speed is not sensitive to the mesh size for the given simulations. The same conclusion can be drawn by considering the relative differences shown in Table 5, where relative differences are computed by considering the quantities computed with the dense mesh as the reference quantities.

Sensitivity with respect to the angle of dilation
Two values of the angle of dilation are considered. The angle of 30 degrees represents a low value and the angle of 35 degrees represents a high value of the angle of dilation. A comparison of shear cone angles as the angle of dilation increases from 30 to 35 degrees indicates that on average there is an increase both in the angle of the shear cone and the weight of spalled concrete, as shown in Table 6, where relative differences are computed by considering the quantities computed with the low dilation angle as the reference quantities.
On the other hand, if one compares the residual speeds, there is, on average, a decrease of the residual speed as the angle of dilation increases form 30 to 35 degrees. To summarize, the simulations are sensitive to the value of the dilation angle, and therefore particular care should be paid in the simulations when choosing the value of the dilation angle. Validation against experimental test results Table 7 shows the relative error of the simulated marker quantities when experimental marker quantities are considered as reference values. Based on these error estimates, the following statements can be made. The proposed user extension to the CDP concrete model and element deletion algorithm estimates rather accurately the residual velocity (except the A12 benchmark test 10 ). The estimates of spalled concrete, shear cone angle tensile to compressive meridians slope ratio

Conclusions
The proposed user expanded Concrete Damaged Plasticity (CDP) model together with a shear strain based element removal algorithm proposal has been investigated by simulating four benchmark impact tests. For each benchmark test the corresponding simulations have been carried out with various material parameter values and mesh sizes. The first conclusion is that the proposed material model and the element removal algorithm is mesh size independent. The second conclusion is that the material model and the element removal algorithm is sensitive to the variation of a particular material parameter, namely the angle of dilation. Therefore particular care should be paid in impact simulations to the choice of that parameter. Validation of the proposed material model and element removal algorithm against experimental results yields rather satisfactory results for the selected benchmark tests. Nevertheless, it must be acknowledged that the range of benchmark tests should be widened in order to generalize these conclusions to arbitrary impact simulations.