Optimal adjustment of FE model of concrete slab exposed to impact loading

Numerical approach using FEM has been used to describe the behaviour of concrete slab exposed to impact loading. 3D parametrical numerical model has been created, and the influence of various parameters values on model response is being investigated. The analyses have been conducted using explicit numerical solver of commercially available software LS-Dyna. The optimal adjustment of the model has been determined.


Introduction
Concrete is widely used material to create the structures of the modern world. In many cases, the requirements for the constructions is not only to withstand the standard situation loads, but also to retain the resistance in severe extreme loading. For example after exposition to various explosions [1] or impacts [2].In such circumstances, more advanced methods in modelling and non-linear material concrete models needs to be used while conducting the analyses, as e.g. by Neuberger et al [3] or by Králik [4].Vast number of different material models are available to mathematically describe the concrete behaviour. These models include complex constitutive relations. An example of such relations have been presented, evaluated and tested e.g. by Král et al [5]. Material models are defined by various parameters, many of those are complicated to obtain directly from experimental data. In order to obtain the optimal values of material model parameters of concrete exposed to specific loading, a sensitivity study and subsequent optimization process, as a tools of the reverse analysis are often helpful to investigate the optimal adjustment of the numerical model.
The purpose of this study is to investigate the influence of several Johnson-Holmquist-Cook (JHC) [6] constitutive material model parameters on the numerical model response of the concrete structure in the matter of deflections and strains compared to real experimental model data. Crack patterns are discussed as well.
Several authors have devoted their research to concrete constructions exposed to blast loading, e.g. Tai et al [7], Zhao and Chen [1,8] or Thiagarajan et al [9]. This paper focuses on parameter identification of JHC material model suitable for concrete class C30.The original JHC model have been applied for concrete with the compressive strength of 48 MPa [6]. Studies of JHC model parameters suitable for high-strength concrete (67.  MPa) have been conducted by Ren et al [10]. Many researchers have been working on determination of the JHC model parameters also for more standardly used concrete classes, e.g. Xiong et al [11][12], however the language of these useful contributions might cause several obstacles to an average European engineer. Reddy [13] has consider different values of some parameters, based on Thacker's Report [14]. Therefore it is still suitable to study this issue closely.

Physical model
Experimental data are obtained from the physical test of concrete slab with dimensions of 500×500×60 mm, material class C30. In the distance of 9 cm from the concrete slab upper surface (in the middle form the top view), 75 g of TNT explosive has been situated (see Fig. 1). One specimen of this considered distance has been tested.
Concrete specimen has been reinforced in both directions by a steel bars of Ø6 á 100 mm. The steel reinforcement has been situated approximately in the mid-height of the slab. Steel plate with thickness of 15 mm has been modelled by shell elements with Belytschko-Tsay formulation and 2 integration points through thickness. Reinforcement steel bars of the slab and the anchorage bolts of the steel plate have been modelled by beam elements with Hughes-Liu cross section integration, and 2×2 Gauss quadrature rule beam integration. The mesh size of the hexahedron solid elements of the concrete slab is 3 mm.
Segment based formulations of symmetrical contacts (pinball algorithms) have been used. The contacts are defined between the concrete slab and the supportive steel plate. The rubber between the steel plate and the concrete base is in contact with both these parts. Mesh size of the supportive steel plate in the contact area is the same as the mesh size of the concrete slab, and is getting coarser towards the plate edges -see Fig. 2. The maximal sizes of the rubber elements and the concrete base elements are 25 mm and 50 mm respectively.

JHC constitutive concrete material model
To describe the mechanical behaviour of concrete structures exposed to large strains, high strain rates and high pressures, it is possible to use the Johnson-Holmquist-Cook (JHC) [6], which has been implemented in the LS-Dyna [15]. The material model is defined by 19 various parameters to capture the mechanical properties of the concrete, strength features and equation of state (EOS), material damage and the strain rate effects.
Material model yield surface ( Fig. 3) with consideration of damage and strain rate effects is determined by: where σ * = σ / fc and P * = P / fc are normalized equivalent stress and pressure respectively; σ is the actual equivalent stress and P is the actual pressure fc is the unconfined uniaxial compressive strength. The dimensionless strain rate is defined as ̇* = * / 0 * , where * and 0 * are actual and reference (EPS0) strain rates, respectively. Normalized tensile strength T * is defined as T * = T / fc, where T is the unconfined uniaxial tensile strength. The maximal normalized strength of concrete is given by Smax. Material strength constants A, B, N represents: normalized cohesive strength, normalized pressure hardening coefficient and the pressure hardening exponent. Strain rate coefficient is C, and accumulated damage is D (0.0 ≤ D ≤ 1.0).

Fig. 3. Yield surface equation.
Damage D (Fig. 4 left) is accumulated from equivalent plastic strain and plastic volumetric strain. ∆ is the effective plastic strain increment and ∆ is the plastic volumetric strain during a cycle of integration. Total plastic strain + under a constant pressure until fracture is defined as: where D1 and D2 are the damage constants, and EFMIN is a material constant which supress the fracture caused by weak tensile waves. Equation of station (EOS) defines the pressure-compaction response of the concrete, and is divided into three phases (Fig. 4 right), with the first phase (0A) linear elastic between the negative pressure cut-off -T (1-D) and the elastic limit value Pcrush. The second phase (AB) corresponds with transitional region, where the air voids are compressed gradually out of the material and the plastic volumetric strain produces the compaction damage until reaching the value of μlock. The third phase (BC) is the compact region where all the air voids have been removed from the material, and the behaviour is defined by material constants K1, K2 and K3.

JHC material model parameter values
The values of JHC concrete material model parameters are summarized in the table 1, based on various researches [6,10,11,13]. In addition to the parameters for concrete class C30, Xiong et al. [11] have also determined values of strength parameters for concrete class C50 as: A = 0.285; B = 1.51 and N = 0.9, and for low strength concrete (10 MPa) as: A = 0.22; B = 1.96 and N = 0.82. Considered ranges of parameters for the sensitivity study are summarized in the last column of the Table 1.
Young's modulus E has been considered in the range of 25-30 GPa; reference strain rate (EPS0) 0.0001-0.001 Hz. Parameter to distinguish the failure type (FS) of the JHC concrete [15] has been considered as: -1 = fail if damage strength < 0; 0 = fail if P * + T * < 0 (tensile failure); or fail if the effective plastic strain > FS > 0. The resolutions of all the parameters have been considered as continuous, in case of FS parameter as discrete by 7 values with equal probability (-1; 0; 0.1; 0.2; 0.3; 0.4 and 0.5).

Applied loads
The simplified blast model with usage of the pure Lagrangian approach of FEM has been involved. The blast wave is considered as a time dependent pressure load applied at the top surface of the slab. Pressure applied in a certain location is based on the empirical blast loading function by Randers-Pehrson and Bannister [16]. The blast loading equation is: where θ is the incidence angle, ( ) and ( ) are reflected and incident (overpressure) pressures respectively, both dependent on time t, and both defined by Friedlander equation (Eq.5) [16]. In case of ( ), the function is stated as: where Prefl.max is the peak reflected pressure, is the waveform number, and + is the positive phase duration. Parameters are defined in dependence on distance from the blast epicentre R, and the equivalent TNT mass W, in a different way for reflected and incident pressure. Damping has not been involved. Self-weight of the structure has been considered by defining the gravitational acceleration of 9.81 ms −2 .

Sensitivity and optimization analysis
In order to achieve the best match between the numerical analysis and the experimentally measured data (bottom surface strain in time), a sensitivity analysis and subsequent optimization process have been conducted [18]. The sensitivity analysis outcome is the identification of the most important model parameters, which, if changed, have the highest influence on the objective of the analysis. In this study, the objective of this process has been defined as the minimization of the difference between the maximums of the straintime curves of the measurement and the numerical analyses. Space filling Latin hypercube sampling (LHS) method has been used to create 300 various combinations of parameters (design points) to conduct the FE analyses, with automatic adaptation of an initial sampling set (using AMOP -adaptive metamodel of optimal prognosis, with 100 samples in each of 3 iterations).
The subsequent optimization process has been carried on using the regressive metamodel (MOP), and considering the most significant parameters (5 in count). Evolutionary algorithm (EA) has been used, considering the full starting population of 300 calculated samples, with min. 30 and max. 35 number of generations, and 20 parents in each generation.

Results
As a result of the sensitivity analysis, 5 significant parameters have been determined, with the highest influence of the normalized cohesive strength A (Fig. 5). The importance is evaluated by CoP value [18]. The values of these 5 most significant parameters (in order to get the best match between the strain maximums) for the best design point #249 (random realization of parameters during the sensitivity analysis) is depicted in the Fig. 3 (3 rd from left). After FE numerical validation analysis of the optimum design #723, yielded from the optimization process (EA algorithm using the regression model), different values of those less significant parameters have been obtained (Fig. 5 most right). The value of the most significant parameter A is not so rapidly different. The values of the rest parameters which differed in realizations #249 and #723 are summarized in the Table 2 below. Time dependence of mid-span bottom surface strain (in x and y direction) is depicted in the Fig. 6 left. The strains are approximately same for both global directions. Timedeflection dependences of the bottom and upper mid-span elements are in the right part of the Fig. 6, Fig. 7 depicts the crack patterns of the experimental specimen and the selected finite element numerical simulations.

Discussion
The maximal strain of case #249 is closer to the maximal value of 7.27 ‰obtained from the experiments then the maximal value of the optimization analysis result #723 (cca. 8.1‰) (Fig. 6 left). The reson could be running the optimization analysis on the regressive model which was generated with certain error (inperfection). The random design point calculated during the sensitivity analysis #249 was so close to the objective of the process, that the optimization analysis was actually useless. Less realizations of parameters set (design points) would be reccomended during the sensitivity analysis to make the optimization worthfull.
The optimization objective has considered only the difference between the maximal values of strains. Neither have been the deflection in time dependence, or the crack patterns of the slab (which is muchharder to evaluate by a single scalar value). The optimization objective was a single one in this study. In further processes it could be reccomended to change the objective. For example to evaluate the strain-time dependences by more scalar entities then only the difference between the maximal curve values (e.g. slope of the curves at certain points could be considered). To define the Euclidean norm of the global difference between the curves, as Král et al have done [19], has not been considered yet. Perhaps creation of a real concrete test specimens geometry approach by Hušek [20] might help. The time-deflection dependence is however influenced also by the stiffness of the other parts of the model (rubber, steel plate, concrete base, soil etc.), therefore in this case such objective would be slightly disguising not knowing the exact model settings.
The value of parameter A (normalized cohesive strength of the concrete)of the case #249 (Fig. 3) of 0.2215 is quite close to the value of 0.27 considered for C30 by Xiong et al [11]. The crack patterns of the bottom slab surface (Fig. 7) for this case are also in a nice match with the experiment results (Fig. 7 left). Case #723 does not seem to be in such a good match considering the crack patterns.

Conclusion
In order to get the best match considering minimization of the difference between maximal values of strain curves from the numerical analysis and the experimental data, normalized cohesive strength of the concrete (parameter A) is the most important for proper evaluation. It was observed, that value around 0.22 is reasonable for concrete class of C30, which is also in match with foreign research.