Finite element analysis of concrete slab exposed to high velocity pressure wave – simplified vs. smoothed-particle hydrodynamics (SPH) method

05005 Finite element analysis of concrete slab exposed to high velocity pressure wave – simplified vs. smoothed-particle hydrodynamics (SPH) method Daniel Jindra1*, Petr Hradil1 and Jiří Kala1 1Brno University of Technology, Faculty of Civil Engineering, Institute of Structural Mechanics, Veveří 331/95 60200 Brno, Czech Republic Abstract. Many structures are required to sustain the structural resistance also under extreme loading conditions, for example impacts of high-velocity objects (airplane crash into nuclear power plant), impacts of projectiles (defence structures), or while exposed to high velocity pressure wave caused e.g. by explosion of various chemicals in industry, nature gas, or also conventional weapons. Numerical analyses of these phenomena are feasible while utilizing explicit approach of the finite element method (FEM), available in commercially accessible software LS-Dyna. In order to predict the behaviour of the structure properly, advanced nonlinear material models are required to be considered, which are often mathematically described by numerous input parameters. Several approaches to model the exposure to blast load exist, from simplified, where the blast wave is considered as timedependent pressure based on empirical equations, to more advanced ones, where the propagation of the pressure wave itself through the surrounding environment is being modelled Arbitrary Lagrangian Eulerian, (ALE method), or so called smoothed-particle hydrodynamics (SPH) method, which might be used to model the blast itself. In this paper, FEM analyses of a simply supported concrete slab with basalt fibre reinforced polymer (BFRP) bars exposed to close range explosion of TNT charge are presented. 3D numerical models are analysed utilizing explicit solver of LS-Dyna. Karagozian and Case (K&C) nonlinear material model for concrete is used, which is suitable when high strain rates are present in the quasi brittle materials. Two variants of the blast loads modelling are compared. The simplified empirical approach, which is less demanding on computational power, and feasible for utilization in case of simple structure geometry, and more demanding method using SPH method to model the TNT detonation and interaction with the exposed concrete slab. The results of these numerical analyses are compared with experimental data based on available literature, and properly discussed. 1 Introduction Numerous structures are required to resist extreme loadings, for example nuclear power plant needs to withstand an airplane crash into the place where reactor is located, as modelled by * Corresponding author: ji[email protected] © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (https://creativecommons.org/licenses/by/4.0/). MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 Králik [1]. Bridge pylons or water dams needs to retain structural resistance after crashes of heavy traffic, or in case of military defense structures after impacts of projectiles. Another example of impact load is high velocity pressure blast wave from explosion of various chemicals in industry, natural gas or conventional weapons. With utilization of modern calculation technology [2], it is possible to numerically model and analyze the effects of the blast load on various structures, hence to determine the structural resistance, or improve the design so the structure can withstand certain load. Several approaches of modelling the blast load are available. The effects of the blast wave might be considered as a pressure load, which is varying in time, and is applied on the exposed surface of the structure, which is modelled by Lagrangian mesh of finite elements. The time dependency of the load is based on empirical equations [3], [4] and numerous experimental data to define the parameters for the inputs [5]. This simplified approach is also implemented in LS-Dyna [2] explicit solver, where it is referred under “load blast enhanced”, abbreviation “LBE”. In case of more complex structural geometries, or when interference of several blast pressure waves from more epicentres of explosion needs to be considered, it is possible to model the propagation of the pressure waves in the surrounding environment (air, water). “Arbitrary Lagrangian Eulerian” (ALE) method is suitable to be used, where the exposed structure is modelled by Lagrangian mesh, but for the surrounding domain, multi material ALE mesh (MMALE) [2] needs to be utilized. This method has significantly larger requirements for CPU time, but it is more robust. For explosive material, equation of state (EOS) is required to be incorporated [6], as well as the material parameters for the explosive itself [7], [8]. It is possible to combine these two approaches (LBE and ALE), and the ALE mesh might be modelled only in the closest surrounding of the structure. The exterior surface of this domain (air) facing the blast epicentre is then covered by special ambient elements [2], which are loaded by the pressure-time functions (LBE). The purpose of these ambient elements is to determine the thermodynamic state data (based on the load pressure data) for the surrounding ALE air domain. Rankine-Hugoniot relations [9] are used to determine the particle velocity and density. The combination of these two methods allows to analyze the impact of the pressure wave on a geometrically more complex surface, while saving some computational time not modelling the propagation of the wave all the way from the blast epicentre which might be located in further distance. Several studies comparing these approaches are available, e.g. by Tabatabaei et al. [10], or Slavik [11]. Another approach is to model the blast itself using the smoothed particle hydrodynamics (SPH) method, which was presented by Monaghan [12] and Lucy [13] independently. Although this method was originally introduced to simulate astrophysical phenomena, the approach has been utilized in numerous tasks, e.g. in simulations of high velocity impacts by Libersky [14], analyses of rock and soil exposed to blast pressure by Pramanik [15] and Chen [16]. Schwer et al. [17] and Trajkovski [18] also modelled air blasts by SPH and compared performance with previous methods. The material response is also dependent on the strain rates – hence, the same material behaves differently under different velocities of the loading. Higher tensile strength of the concrete material has been observed when exposed to high velocity impact loads [19]. Concrete structures exposed to blast loading has been studied e.g. by Tai et al. [20], Zhao and Chen [21], [22], Thiagarajan et al. [23] and Dubec, et al. [24]. The phenomena of structural response under high velocity impact loads is still object of continuous research. This paper presents numerical finite element (FEM) analysis of recent physical experiments by Gao et al. [25], where the performances of sea-sand sea-water concrete slabs reinforced by basalt fibre reinforced polymer (BFRP) are compared with ordinary concrete BFRP slabs, while exposed to close-range TNT blasts. In this study, the analyses of ordinary concrete slabs are 2 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 presented and discussed. In order to properly describe the material response when high strain rates are involved, Karagozian and Case (K&C) [26], material model is utilized. The strain rate effects of concrete are defined in accordance with research by Malvar et al. [19]. The strain rate effects of the reinforcement are neglected, as will be explained later. Otherwise if applicable, it would be possible to consider these effects based on the review of static and dynamic properties of steel bars by Malvar et al. [27]. 2 Physical model and experiments In detail the physical experiments are described by Gao et al. [25]. In this study, the slab noted by Gao et al. as BRPS1 [25] – the ordinary (plain) concrete slab is numerically analyzed. The concrete class is C40 with the average compressive strength of 49.34 MPa, based on tests of 6 concrete cubes (150 mm) cured for 28 days at room temperature [25]. Concrete tensile strength (at static strain rates) is considered as 3 MPa, Poisson’s ratio as 0.2. The dimensions of concrete slab are: 1100 mm in length, 50 mm in height (thickness) and 500 mm in width. Slab is reinforced by Ø 6 mm BFRP bars at the bottom surface (10 mm concrete cover) in both directions in regular grid with span of 100 mm. The effective depth of the slab is 37 mm. Material response of the BFRP reinforcing bars has been tested [25], the average tensile strength is 1.53 GPa (see Fig. 1 a, based on [25], [28]) and the Young’s elastic modulus is 57.68 GPa. Fig. 1 b depicts the experimental setup. Slab is simply supported on a steel frame, structural span is 1000 mm. At the edges, the slab is secured by the frame from the top side to avoid the post-blast uplift. The TNT charge is located above the mid-span in the stand-off distance of 1000 mm. In this study, the variant with 0.4 kg of TNT is analyzed (by Gao et al., also 0.8 kg and 2.0 kg was tested [25]). Strain gauge has been installed directly on the BFRP bar in the mid-span before casting of the concrete. On the top surface of the slab, air pressure sensors are located, and one free-field air pressure sensor is in 3 m distance from the charge (Fig. 1 b, c). (a) (b) (c) Fig. 1. a) tensile tests of the BFRP bars, scheme based on graph by Gao, Feng et al. [25],[28]; b), c) Experiment set-up 3D perspective schemes based on figures in Gao, Feng et al. [25]. 3 Numerical analyses 3.1 Simplified blast model – method #LBE In the simplified approach known as “load blast enhanced” (LBE) in LS-Dyna [2], the blast pressure wave is considered as time-dependent surface pressure load, which is determined by the empirical blast loading function defined by Randers-Pehrson and Bannister [3]: ),coscos)(t(Pcos)t(P)t(P sr  21 22  (1) 3 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 where θ is the angle of incidence, Ps(t) and Pr(t) [Pa] are time dependent incident (free-air) and reflected overpressures respectively, both defined in accordance with Friedlander equation [4], which is in case of free-air overpressure Ps(t) determined as: ,e)t/t(P)t(P o t/tb osos  1 (2) where Pso is the peak incident (free-air) overpressure [Pa], b, is waveform decay coefficient [-], and to is the duration of the positive phase [s] (afterwards the overpressure wave, there is a wave of lower pressure). These parameters are determined based on scaled distance Z introduced by Hopkinson [29] and Cranz [30]: , W R Z /31  (3) where W [kg] is equivalent mass of TNT, and R [m] is the distance from blast epicentre. Values of these parameters in SI units might be obtained from JRC technical report [5]. The analyzed case of 0.4 kg of TNT in the 1 m stand-of distance results in scaled distance of 1.357 m·kg−1/3, and the contact time of the blast wave with the structure surface ta is expected to be approximately 700 μs. Experimental reference time-pressure data depicted in Fig. 2 a) are based on study by Gao et al. [25]. Position of the sensor A is in the mid-span right beneath the TNT charge and sensor B is at the slab edge (Fig. 1 b). The experimental time-pressure curves are in a rather nice match with the empirical estimations (equations 1, 2), and the values of the pressure impulses i [Pa·s] obtained by integration of these pressure-time curves approximately corresponds (280 ≈ 256 for sensor A; 270 ≈ 214 for sensor B, what is considered as a nice match). (a) (b) Fig. 2. a) Pressure-time experimental references and empirical estimations; b) estimation of the pressure impulse based on empirical approach – symmetric 1/4 of the slab. 3.2 Smoothed particle hydrodynamics – method #SPH The SPH method is a mesh-free (Lagrangian) solver originally developed for hydrodynamics. The governing equations of fluid dynamics are in a form of partial differential equations, and are determined by interpolation from the particles. Discretized formulation derivation of the SPH method is divided into two steps. In the first step called kernel approximation, the arbitrary function and its gradient are introduced in dependence on smoothing length h and smoothing kernel function W (adopted as cubic B-spline [2]):    'dx)h,'xx(W)'x(f)x(f (4)   'dx)h,'xx(W)'x(f)x(f (5) 4 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 In the second step, which is termed as particle approximation, the integral forms of the function and the function gradient are approximated by the summarization of the nearest particle values:   N jj j ijji m W)x(f)x(f 1  (6)    N jj j ijji m W)x(f)x(f 1  (7) where mj and ρj are mass and density respectively. Wij = W (xi − xj,h). For the explosive TNT material, Jones Wilkins Lee (JWL) Equation of State (EOS) is used [6], where the pressure is defined as a function of the internal energy per volume E, and relative volume V: V E e VR Be VR AP VRVR                     21 21 11 (8) where EOS parameters A, B, R1, R2, ω and also TNT material parameters as ChapmanJouguet Pressure (PCJ), detonation velocity, density and internal energy per reference volume E0 [2] are defined in accordance with [7],[8]. 3.3 Nonlinear material model for concrete - Karagozian and Case (K&C) K&C material model [26] is suitable for numerical analysis of concrete structures exposed to high strain rates. The first release of the material model is dated to year 1994 [31]. In the second release, some additional aspects as shear dilation has been implemented [32]. The third release [33] has introduced automatic parameter generation based on uniaxial compressive strength (used in this study). The dependence on finite element mesh geometry due to strain-softening of K&C has been reduced in 2010 [34]. More details are available is in K&C report [35] which describes the use and validation of this material model. K&C model is constitutive model defined by three invariants with using of three shear failure surfaces: initial, maximal and residual. These surfaces are mutually independent, and the general definition is mathematically described by equation (8) as: paa p a)p(F ii ii 21 0  , (9) where y, m and r is substituted for index i to describe initial yield strength surface, maximum shear failure surface and residual failure surface respectively. Values of the parameters 𝑎𝑎𝑗𝑗𝑖𝑖 (j = 0,1,2; i = y, m, r) are required to be calibrated based on experimental data. Pressure 𝑝𝑝 [Pa] is dependent on the first invariant of stress tensor: 3 1 I p (10) The resultant failure surface is interpolated between the maximal and initial (Equation 11) or between the maximal and residual (Equation 12) failure surface based on provided formulas:     myym )p(F)p(F)p(F)()J(r)J,J,I(F   3321 (11)     mrrm )p(F)p(F)p(F)()J(r)J,J,I(F   3321 (12) where 𝐼𝐼1 is the first invariant of the stress tensor, 𝐽𝐽2 and 𝐽𝐽3 are the second and the third invariants of the deviatoric stress tensor. 𝜆𝜆 express modified effective plastic strain, 𝜂𝜂(𝜆𝜆) is a function of internal damage dependent on 𝜆𝜆, with values: 𝜂𝜂(0) = 0 , 𝜂𝜂(𝜆𝜆𝑚𝑚) = 1, 𝜂𝜂(𝜆𝜆 > 𝜆𝜆𝑚𝑚) = 0. This means, that the failure surface begins at the initial yield strength surface, and reaches 5 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 the maximum shear surface as 𝜆𝜆 is increasing to 𝜆𝜆𝑚𝑚. Afterwards the failure surface decreases to the residual surface for further increasing 𝜆𝜆 up to the value of 𝜆𝜆𝑚𝑚. The relations between 𝜆𝜆, 𝜆𝜆𝑚𝑚ax and 𝜂𝜂(𝜆𝜆) are calibrated based on experimental data. 𝑟𝑟(𝐽𝐽3) is a scaling factor in a form of equation by William Warnke [36], which express the dependence on 𝐽𝐽3 in a way that the transition between brittle and ductile (under higher confinement) response is well described. 3.4 Strain rates material dependency The material strength parameters increase rapidly under high strain rates   . During the exposure to blast loads strain rates are in range from 10 to 1000 s−1, and the increase is about 100% for concrete compressive strength, and 600% for the concrete tensile strength [19]. The effect of strain rate on concrete parameters is expressed by dynamic increase factors DIF, noted TDIF for tensile strength increase and CDIF for the compressive strength increase, and depends on strain rate   [19]: 1 0161 30            sCDIF f f . cscs c        ; 1 31 30           sCDIF / cs        (13) cu f.  7505 1  ; 21566  .)log( ; ccu f.f 2051 (14) 1 1            sTDIF f f tsts t        ; 1 31 1          sTDIF / ts        (15) cocf/f 81 1  ; 26   )log( ; MPaf co 10 (16) Where fc [Pa] is compressive strength at the dynamic strain rate   in range 3·10−5 – 300 s−1, fcs [Pa] is compressive strength at the static loading strain rate cs   = 3·10−5 s−1, and the relation is determined by parameters δ, β and concrete cubic strength fcu. Analogically for tension, ft [Pa] is tensile strength at dynamic the dynamic strain rate   in range 1·10−6 – 300 s−1, fts [Pa] is tensile strength at static loading strain rate ts   = 1·10−6 s−1, and α γ parameters are involved. For reinforcing material (BFRP bars), the increase in yield strength at higher strain rates is not considered, due to the fact that the maximal strain experimentally detected in a bar was 1.5% [25], what corresponds to stress of 865 MPa (E = 57.68 GPa), smaller than the tensile strength of 1.53 GPa with a brittle failure (Fig. 1 a), so there is no need to incorporate the dynamic increase factors for the BFRP bars. The material model in the numerical analyses is considered as bilinear, with negligible hardening (linear elastic, ideal plastic), with the yield stress of 1.53 GPa. The strains are monitored carefully that are below 2.65%, what is limit value at the tensile strength. 3.5 Numerical finite element models In this study, performances of two different approaches of blast loads modelling are compared, the simplified approach (LBE), and the smooth particle hydrodynamics (SPH). Numerical finite element models are depicted in the Fig. 3. In both cases, the concrete slab is being modelled the same way (Fig. 3 a), by hexahedral 8 nodal solid elements (with 3 translational degrees of freedom per each node). Near the line supports at the edges, the concrete solid elements are cubes (either 10 or 5 mm), and in the ¼ of structural span, there is a mesh refinement area, where two dimensions of the cubes are refined into half size, resulting in prisms (either 10×5×5 mm or 5×2.5×2.5 mm). The element size along the slab 6 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 width remains unchanged. These two mesh sizes are further noted as 10×5 mm, and 5×2.5 mm for the cases of 10 mm and 5 mm edge sizes respectively. Solid elements with constant stress formulation along with hourglass control have been used. Boundary conditions are applied to a layer of contact material elements, and linear elastic concrete is defined for solid elements near these contact areas (Fig. 3 a). These boundary conditions are applied also from upper surface in order to simulate supports of the steel frame (Fig. 1 b). The reinforcement bars are modelled by beam elements, with Hughes-Liu formulation with cross section integration [2]. This option is more robust, suitable also for plastic materials, but demands more CPU time then formulation by Belytschko-Schwer, which would be suitable only for linear materials [2]. Mesh size of the reinforcement is the same as the refined mesh in the mid-span of the slab (hence either 5 or 2.5 mm). The history of axial strain is being monitored in a longitudinal reinforcing bar in the mid-span of the slab. While utilizing the SPH approach, in addition to concrete slab, it is required to discretize also the TNT charge (Fig. 3 b). The SPH nodes of the explosive are required to be aligned in a regular cubic grid, radial grid is not suitable [2]. The exact dimensions of the TNT block, as referred by Gao et al. [25] was unknown, so it is assumed that the 0.4 kg TNT charge is in a shape of cylinder of height 0.1 m. Based on TNT density, which is considered as 1650 kg m–3 [7], [8], the diameter has been determined. Two grid sizes of the SPH nodes are compared in this study, 2 mm and 1 mm. The masses of each SPH node (for certain case of the grid size) are the same, and are determined from the total number of modelled SPH nodes and the total weight of the TNT charge (0.4 kg). Two variants of blast epicentre are compared – in the cylinder centre of gravity (#mid), and at the upper side of the cylinder (#up). (a) (b) Fig. 3. Finite element model for: a) the simplified blast load approach (LBE); b) smooth particle hydrodynamics (SPH) During the explicit analysis, the time is required to be discretized into small steps, where the minimal critical time step is determined based on mesh size and the wave propagation velocity through the material continuum (depends e.g. on density, elastic modulus) [2]. The smaller the mesh, the smaller is the critical time step. The automatic default time step might be altered by the ts parameter (time step scale factor). Several variants of this parameter have also been considered and the results are compared. 4 Analyses results The analyses were calculated using Intel Xeon E5 1620 0 (Sandy Bridge-EP) 3.6 GHz, with total number of 4 CPU's, using 4 SMP threads. Summary of all the analysed cases is in the table 1. The name of the case contains information about the slab mesh density (either 10×5 or 5×2.5), SPH grid (1 or 2 mm), time step scale factor (ts) and the position of the initial detonation epicentre (#mid or #up), as closely described in previous the chapters. In the second column, there is actual discretization time (for mesh cases 5×2.5 these are average 7 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 values). If the time step is not fine enough, penetration of SPH particles through the slab might occur – as noted in the last column (see also Fig. 4 for insight, what is subjectively considered as “large penetration” by the authors of this study). Table 1. Summary of analyzed cases. Case Analysis time step [s] Physical calculation time [h:m:s] SPH penetration LBE 5×2.5 2.76E-07 00:37:32 – 10×5, SPH1, ts 0.05, #mid 1.53E-08 22:25:55 None 10×5, SPH1, ts 0.05, #up 1.53E-08 21:11:49 None 10×5, SPH1, ts 0.3, #up 9.19E-08 3:28:41 Large 10×5, SPH1, ts 0.6, #up 1.84E-07 1:44:11 Very large 5×2.5, SPH1, ts 0.3, #mid 5.67E-08 20:32:20 Large 5×2.5, SPH1, ts 0.3, #up 4.46E-08 25:34:16 Large 5×2.5, SPH1, ts 0.05, #up 1.12E-08 104:04:09 None 10×5, SPH2, ts 0.3, #up 9.19E-08 1:43:37 Small 10×5, SPH2, ts 0.1, #up 3.06E-08 4:50:01 None 10×5, SPH2, ts 0.1, #mid 3.06E-08 4:53:21 None In Fig. 4, there are plots of SPH resultant velocities at time 450 μs after detonation. The difference between the SPH clouds based on detonation point location (#mid or #up) is also nicely visible (Fig. 4 a and b). The color scale is the same for all the cases in the Fig. 4. Resultant velocity is smaller for larger more massive particles (Fig. 4 c). Time histories of the axial strain are depicted in the Fig. 5. Fig. 5 a summarizes the cases with 1 mm SPH grid and 10×5 mesh size. Coarser 2 mm SPH grid with the same mesh is in Fig. 5 b. The finest mesh 5×2.5 with 1 mm SPH grid is in Fig. 5 c. Analogically, the midspan displacements of the slab are depicted in Fig. 6. (a) 10×5, SPH1, ts 0.6, #up (b) 5×2.5, SPH1, ts 0.3, #mid (c) 10×5, SPH2, ts 0.1, #up (d) scale Fig. 4. Resultant velocity of SPH particles at time 450 μs after detonation, front ant bottom views. 8 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023 (a) (b) (c) Fig. 5. Axial stress of reinforcement in time: reference data compared with numerical analyses. (a) (b) (c) Fig. 6. Mid-span displacement in time: reference data compared with numerical analyses. 5 Discussion Sufficiently small time step is required in order to avoid the penetration of the SPH particles with the exposed structure. To achieve this, the default critical time step of the explicit analysis [2] needs to be decreased by the time step scale factor. This value is dependent on mesh and SPH grid size. The detonation velocity of the TNT explosive material considered in these analyses is 6930 m s–1 [7], [8], and no velocity of any SPH particle is larger then this in the moment after detonation (and before reaching the contact surface). Although, the velocity values of the particles which have penetrated through the slab have been numerically increased (Fig. 4 a). No velocity is increased when all the particles bounce at the contact surface. The time a particle with the TNT detonation velocity would travel the distance of 2.5 mm (finer concrete slab mesh) is 0.36 μs. Based on Table 1 data, to achieve no penetration, time step of the analysis needs to be smaller than 12.36 % of these 0.36 μs (5×2.5, SPH1, ts 0.3), 3.1 % might already be unnecessarily too small (5×2.5, SPH1, ts 0.05). For 5 mm distance, this time would be 0.72 μs, and sufficiently small time step to avoid penetration was 4.24 % of this value (case 10×5, SPH2, ts 0.1). Small penetration occurred with time step size equal to 12.7 % of 0.72 μs (10×5, SPH2, ts 0.3). Although the exact limit value of the time step where no penetration of the SPH particles occurs is not determined, it appears that time step scale factor ts of 0.1 might be sufficiently small. This step size corresponds with circa 5 % of the time the SPH particle at the TNT detonation velocity needs to travel the distance of the exposed mesh elements. Comparing the performances of the simplified LBE method (case LBE 5×2.5) with the SPH, for the geometry of a simple slab the utilization of the SPH is rather ineffective (see calculation times in table 1). In order to achieve comparable performance in matter of displacement and reinforcement axial strain time histories (Fig. 5, Fig. 6), approximately 36 – 160 times more physical time is required (considering cases 10×5, SPH1, ts 0.05; and 5×2.5, SPH1, ts 0.05). 9 MATEC Web of Conferences 396, 05005 (2024) https://doi.org/10.1051/matecconf/202439605005 WMCAUS 2023