Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 51, 3, pp. 615-625, Warsaw 2013 COMPARISON OF NUMERICAL TESTING METHODS IN TERMS OF IMPULSE LOADING APPLIED TO STRUCTURAL ELEMENTS Łukasz Mazurkiewicz, Jerzy Małachowski, Paweł Baranowski, Krzysztof Damaziak Military University of Technology, Department of Mechanics and Applied Computer Science, Warsaw, Poland e-mail: lmazurkiewicz@wat.edu.pl; jerzy.malachowski@wat.edu.pl; pbaranowski@wat.edu.pl; kdamaziak@wat.edu.pl This paper presents comparison of numerical testing methods of an impulse loading which comes fromadetonationprocess, i.e. blastwavepropagation in a gasmedium. Investigations were carried out using an analytical and numerical model based on the Finite Element Method. In order to reduce computational time, the substitute analytical model with one degree of freedomwas implemented, which replaced a chosen actual system (I-section steel column).For structuremodelling, the constitutivemodelwasused,which included the strain rate effect. From the performed analyses, an acceptable similaritywas noticed, although the discretemodel due to greater forces gave inflated results. Nevertheless, it should be pointed out that simplifiedmethods do not take anywave and flowaround effects into consideration, which have an influence on the dynamical response of the structure and are possible to implement in the gas medium coupling. Key words: blast wave, FEM, Eulerianmesh 1. Introduction Within mechanics, some analyses are carried out which enclose the effects of pressure wave propagation and their impact on encountered structures. Hence, two main groups of methods for numerical implementation of pressure wave propagation can be distinguished, i.e. methods based on analytical models andmethods which use discrete models. The first ones are methods which substitute the problem into a system with one degree of freedom (SDOF) or multiple degrees of freedom (MDOF) (US ArmyCorps of Engineers, 2008). The others are based on the Finite Element Method (FEM), which use several algorithms for determination of dynamical loading of the structure, e.g. ConWep blast wave function, Smoothed Particle Hydrodynamics (SPH) method and Arbitrary Lagrangian-Eulerian formulation (ALE). A detonation process and pressure wave propagation and its interaction with structures were thoroughly investigated by other authors (Lu et al., 2003;Małachowski, 2010a,b;Mazurkiewicz andMałachowski, 2011). The presented paper encloses comparison analysis of numerical methods based on discrete models, where the detonation process was simulated using ALE formulation, ConWep function and incorporate analytical model with one degree of freedom (SDOF). The first two methods are available in LS-Dyna code, whereas the systemwith SDOF is implemented in “Single degree of freedom Blast Effects Design Spreadsheet” (SBEDS) software developed by US ArmyCorps of Engineers (USArmyCorps of Engineers, 2008). Themain objective of the presented paper is to assess effectiveness and reliability of the mentioned methods on the example of a structural element subjected to a dynamic impulse. 616 Ł. Mazurkiewicz et al. 2. Theoretical basis of implemented methods As stated before, there are two methods of analysis of a structure subjected to an impulse loading: methods based on analytical models (beam elements) and methods which use discrete models reflecting full three-dimensional actual geometry of the structure. 2.1. Analytical models The simplest models for dynamical responses of structural elements are models which use systems with one degree of freedom which replaces the actual structural element with lumped stiffness, inertia and damping (Fig. 1). Fig. 1. Loaded beammodel replaced by a systemwith SDOF (USArmy Corps of Engineers, 2008) Consequently, the following dynamic equation of dynamic equilibrium is formulated (US ArmyCorps of Engineers, 2008) Meü(t)+Ceu̇(t)+Re(u)u(t)=Fe(t) (2.1) where Me is the effective mass of the equivalent SDOF system, Ce – effective viscous dam- ping coefficient (constant), Re(u) – effective stiffness, ü(t) – acceleration, u̇(t) – velocity, u(t) – displacement of the mass, Fe(t) – effective load acting on the equivalent SDOF system, t – time. Theequivalent SDOFsystem isdefinedas a systemthathas the sameenergy in termsofwork energy, strain energy and kinetic energy as the real beamelement subjected to the pressurewave from a detonation. This is achieved through relating the applicable values of the blast loaded component and the corresponding effective parameters of the equivalent SDOF systemwith the relevant transformation coefficients, as shown in (2.2) (US ArmyCorps of Engineers, 2008) KLmMcü(t)+Ccu̇(t)+Rc(u)u(t)=Fc(t) (2.2) where Mc is the mass of the blast-loaded component, Cc – viscous damping constant, Rc(u) – stiffness resistance, Fc(t) – history of load acting on the blast-loaded component, KLm – load-mass factor equal to Km/KL, Km –mass transformation factor, KL – load trans- formation factor. The mass transformation factor Km is calculated by comparing the work of external forces for the component and SDOF system. Consequently, for the component, this work is obtained by (US ArmyCorps of Engineers, 2008) WE(t)element = L ∫ 0 P(x,t)U(x,t)dx (2.3) where L is the length of the blast-loaded component, U(x,t) – its displacement, P(x,t) – dy- namic load acting on the blast-loaded component. Comparison of numerical testing methods in terms of impulse ... 617 Assuming that U(x,t) = Umax(t)φ(x) and P(x,t) = Pmax(t)p(x), which gives: φ(x)=U(x,t)/Umax(t) and p(x)=P(x,t)/Pmax(t), the following equation is obtained WE(t)element =Pmax(t)Umax(t) L ∫ 0 p(x)φ(x)dx (2.4) Then, for the SDOF system, the work of external forces can be described by WE(t)SDOF =Fe(t)u(t)=KLFc(t)u(t) (2.5) Assuming that u(t)=Umax(t) and Fc(t)=Pmax(t) ∫L 0 p(x)dx it gives WE(t)SDOF =KLPmax(t)Umax(t) L ∫ 0 p(x)dx (2.6) By comparing equations WE(t)element =WE(t)SDOF, after some transformations, KL co- efficient is obtained KL = L ∫ 0 p(x)φ(x)dx L ∫ 0 p(x)dx (2.7) Analogously, by comparing kinetic energy for the component and SDOF, the system Km coefficient is obtained Km = L ∫ 0 m(x)φ2(x)dx L ∫ 0 m(x)dx (2.8) The formulated equivalent system is then loaded with a pressure impulse. 2.2. FE analyses Numerical simulationswere carried outusinganLS-Dynaexplicit codededicated fordynamic simulations using the Finite Element Method (FEM) central difference scheme with modified time integration of the equation ofmotion. The basis of thismethodwas frequently described in the literature (Belytschko et al., 2000; Hallquist, 2003). Themain assumption of it is to replace the continuous componentbyadiscrete onewhich consists of finite elements connected bynodes. For the definedmodel, a well-known matrix equation of motion is defined by (Hallquist, 2003) Mü(t)+Cu̇(t)+K(u)u(t)=F(t) (2.9) where M is themassmatrix, C – dampingmatrix, K(u) – stiffnessmatrix, ü(t) – acceleration vector, u̇(t) – velocity vector, u(t) – displacement vector, F(t) – vector of external forces, t – time. 618 Ł. Mazurkiewicz et al. 3. Blast wave loading modelling 3.1. SBEDS pressure function Both in the analytical and discrete model, the dynamic pressure loading is described by a mathematical function and an idealized shape of the pressure pulse at a point caused by the shock wave from a high explosive detonation, see Fig. 2. Fig. 2. Typical pressure-time history of the pressure wave from an explosion (US Army Corps of Engineers, 2008); i s i s – positive and negative impulse, respectively, t o t o – positive and negative phase duration, respectively The shape of the positive phase blast load from an open air explosion in Fig. 2 can be represented mathematically with Eq. (3.1) (US Army Corps of Engineers, 2008) Ps(t)=Pso [ 1− (t− tA to )] exp ( − t− tA θ ) tA ¬ t¬ tA+ to (3.1) where Ps(t) is the shockoverpressureasa functionof time [MPa], Pso –peakside-onoverpressure [MPa], t – detonation time [ms], tA – arrival time of the initial shock front [ms], to – positive phase duration [ms], θ – shape constant of the pressure waveform. Fig. 3. Positive phase parameters of the blast wave for a spherical TNT explosion in air at the sea level (US ArmyCorps of Engineers, 2008) In order to obtain peak pressure values, the characteristics obtained from empirical tests are usedwith the reduced distance of the charge (Z parameter) which is described byEq. (3.2) (US Army Corps of Engineers, 2008). Characteristics of examplary explosion parameters with the side-on (subscript so) and reflected (subscript r) pressure wave are presented in Fig. 3 R1 R2 = 3 √ W1 W2 → Z = R 3 √ W (3.2) Comparison of numerical testing methods in terms of impulse ... 619 where Wi ismass of the explosive [kg] for the case i, Ri – standoff distance from the detonation point of Wi, Zi – scaled standoff for the case i. 3.2. ConWep function pressure impulse loading This method is similar to the one mentioned before and used in this paper, i.e. the SBEDS system. It was developed by Randers-Person and Bannister in 1997. In this case, the load blast is described by the ConWep function which substitutes the wave propagation with a pressure function. In order to simulate the blast loading, it is necessary to determine the equivalent of TNTmass, the type of explosion (on the surface or in the air), charge localization and surfaces segments interactingwith theblastwave.Themaindifference between those twomethods is that theConWep function applies the loading separately on every finite element of amodel. Thus the effect of pressure wave propagation within the loaded structure can be observed from elements which are closest to the explosive. ConWepapproach allows generation of the equivalent pressure value (Hyde, 1988; Hallquist, 2003) P(t)=Pref cos 2θ+Pin(1+cos 2θ−2cosθ) (3.3) where θ is the incidence angle of the wave, Pref – overpressure value of the reflected pressure wave, Pin – overpressure value of the incident pressure wave. 3.3. Arbitrary Lagrangian-Eulerian formulation (ALE) with ConWep implementation Another way to represent the blast wave interaction with a structure is the Arbitrary Lagrangian-Eulerian formulation, where it is necessary to define an Eulerian air domain, in which the explosive pressure wave propagates. In the presented case, the pressure impulse was applied using the ConWep model, Eq. (3.3), to the Eulerian air domain surface. The air is considered as a simple ideal gas with a linear polynomial equation of state (Hallquist, 2003) p=(C4+C5µ)E (3.4) where µ= ρ/ρ0 is density, ρ0 is the initial density, C4 and C5 are polynomial equation coeffi- cients, E is the internal energy. TheALEprocedure consists of twomajor steps: the classical Lagrangian step and the advec- tionEulerian step.The advection step is carried outwith the assumption that the displacements of nodes are very small in comparison to characteristics of the elements surrounding these nodes, e.g. dimensions. Moreover, in this procedure, a constant topology of the mesh is provided. The governing equations for the fluid domain (Euler domain) describe the conservation of mass, momentum and energy (Hallquist, 2003) dM dt = d dt ∫ V (t) ρvdV = ∮ S(t) ρ(w−v) ·nd dQ dt = d dt ∫ V (t) ρvdV = ∮ S(t) ρv[(w−v) ·n]dS− ∫ V (t) ∇pdV + ∫ V (t) νgdV dE dt = d dt ∫ V (t) ρedV = ∮ S(t) ρe(w−v) ·ndS− ∫ S(t) pv ·ndS+ ∫ V (t) ρg ·vdV (3.5) where ρ is fluidmass density, p – pressure, g – acceleration of gravity, e – total specific energy. The quantities M,Q and E are the total mass, total momentum and total energy, respectively, of the control volume V (t) bounded by a surface S which moves in the fluid (gas-air) with an arbitraryvelocity wwhichmaybe zero inEulearian coordinates or v inLagrangian coordinates. The vector n is the outwards normal to the surface S. 620 Ł. Mazurkiewicz et al. 4. Object of investigations The research object in the present paper is an I-section steel column 203× 203× 86 (British standards) made of A992 steel subjected to a blast load generated by 250kg of TNT at 10m distance from the component. Also, static axial loading 300000N is applied and fixed-simple boundary conditions are assumed (Fig. 4) Fig. 4. I-section steel column used in the analysis 4.1. FE model Adiscretemodel of thementioned column (Fig. 5) has been developed usingBelytschko-Lin- -Tsay shell elements (Belytschko et al., 2000). These elements are effective in simulating several numerical problems due to a reduced number of important mathematical operations that have to be performed (Hallquist, 2003) Fig. 5. Numerical model of the column with the applied axial loading 4.2. Material model The load due to outbreak causes very rapid changes of the strain field in the structure. Therefore, an elastic-plastic material model with isotropic hardening and strain rate effects has been applied in order to describe properties of the beam material. Thus, simplified Johnson- -Cook (J-C) model was implemented, which provides a prediction of the flow stress σflow for large strains and high strain rates, when its dependence on the strain rate is linear in a semi- logarithmic scale. Material parameters presented in Table 1 for the simplified J-Cmaterial were Comparison of numerical testing methods in terms of impulse ... 621 taken from the literature (Małachowski, 2010a,b). The mathematical formula which describes this model is as follows (Hallquist, 2003) σflow = [A+B(ε p)n](1+C ln ε̇p ∗ ) (4.1) where A,B,C,n arematerial constants, εp – effective plastic strain, ε̇ p ∗ – effective plastic strain rate. Table 1. Steel parameters for the simplified Johnson-Cookmodel used in the analysis ρ [kg/m3] A [MPa] B [MPa] C n 7850 365 510 0.0936 0.9 5. Numerical analyses 5.1. Pre-stress (static) simulation Simulation of two numerical models (ConWep and by ALE+ConWep cases) was preceded by the so called pre-stress stage. In this stage, the I-beam structure was subjected to nominal load Pn, whichwas equal to the load carried by the pillar during its normal service. Incremental static analysis was performed using the full Newton-Rapshon algorithm. The equation solved at this stage has the following form (Hallquist, 2003) Ki∆xi−1 =∆Qi (5.1) where Ki is the stiffness matrix, ∆xi−1 – displacement vector, ∆Qi – external force vector. Convergence of the solution was controlled by two criteria: displacement and energy relative convergence tolerance. In the next stages, the results from static analysis were taken into account as the pre-stress field present in the column. It was obtained using a dynamic relaxation procedure. Dynamic relaxation allows one to quickly reach the preloaded state by the linear ramping nodal displa- cement field to prescribed values over 100 time steps. It should be noted that the comparison of the stress field taken from the static analysis and the stress field generated by the dynamic relaxation procedure showed small differences introduced by the latter. On the other hand, the procedure allows for application of the predefined stress filed on the selected part of the FE model in a very effective manner. 5.2. ConWep function loading – conditions of analysis As mentioned before, the ConWep function simulates the pressure wave loading from an explosive by generating an applicable pressure value P(X,t),X =X(x,y,z). Consequently, the pressure is then distributed on chosen surface segmentswith time and element localization taken into consideration (Fig. 6). 5.3. Lagrangian-Eulerian with ConWep – conditions of analysis The most advanced case of the proposed analysis is Arbitrary-Lagrangian-Eulerian formu- lation where Euler elements define surrounding air. The interaction between the Lagrangian structure, which in this case is the I-section column, and air is obtained through numerical coupling of these two media (Małachowski, 2010a,b; Mazurkiewicz, 2011). The pressure wave was applied analogously to the previous case with the ConWep function, but with one major difference: it was applied to the surrounding air surface, not directly to the column surface. Thus, the pressure P(X,t) was a particular boundary condition on one of 622 Ł. Mazurkiewicz et al. Fig. 6. Pressure impulse loading – schematic presentation the outer surfaces of the gas medium described by Eulerian coordinates (Fig. 7). Additionally, on the outer walls of it, a non-reflecting option was applied, which considers the flow of the pressure outside the domain. The blast wave propagated in the air domain and interacted with the structure, which can be noticed in Fig. 8. Fig. 7. Numerical model of the column inside the Euler air medium Fig. 8. Pressure wave flow effect in the air (Eulerian) medium Comparison of numerical testing methods in terms of impulse ... 623 6. Results and conclusions From the performed analysis, dynamic responses and deformations of the structural element were obtained. The column structure deflection versus time characteristics for all three cases are presented in Fig. 9, which also shows the dynamic response of these components. It should be pointed out that during the analysis, the column structure did not exceed the yield stress. From the carried out numerical simulations, it can be concluded that both mixed ALE-ConWep and ConWep models give similar results. Fig. 9. Column deflection versus time characteristics for all studied cases For the SBEDS system, the obtained dynamic response is higher due to a greater pressure impulse energy (longer duration time) generated by the algorithm. It was confirmed through further analyses, where the pressure impulse wasmeasured in all methods (Fig. 10). Due to the wave effects implemented in ALE + ConWep case, the incident and reflected pressure can be noticed, and also the pressurewave characteristic is not ideal as in the othermethods.The initial value of pressure was caused by the atmospheric pressure, which was defined by the polynomial EOS for the air Euler medium. Taking a closer look at Fig. 10, it can be concluded that both mixed ALE-ConWep and SBEDSmodels give similar results (pmax ≈ 2.2MPa), whereas in the ConWep model pmax ≈ 1.3MPa. Fig. 10. Pressure impulse versus time – comparison of the methods 624 Ł. Mazurkiewicz et al. It is worth tomention the effectiveness of implementedmethods. TheALE formulation, due to the greater number of elementsmodeling the airmediumand its complexmathematical back- ground, is much more computationally expensive – the analyses were carried out for ∼ 33min. For the ConWep model, the simulation ended after 10min, whereas for the SBEDS system, it was approximately 1min (Fig. 11). It is also noteworthy that in order to perform simulations using discrete models, there is a need to develop geometry of the structure and, consequently, a discrete model with finite elements, apply boundary initial conditions and implement a consti- tutivematerial model suitable for a specific type of the simulated problem. In contrast with the SBEDS system, the overall model is chosen from the available software database. Fig. 11. Comparison of CPU time for all cases The presented above comparison results show great efficiency of the analytical models. Ne- vertheless, the scope of application of this method is very limited and only allows for analysing basic structural elements. For the preliminary analyses of the response of components subjected to the pressure impulse from an explosion, the ConWep function is adequate enough. It gives muchmore information about the structure deformation (local deformation of the element), al- lows formore complexmaterial models (with strain rate taking into account) and does not have so many limitations of its usage. The most advanced Arbitrary-Lagrangian-Eulerian formula- tion, where Euler elements define the surrounding air, is muchmore computationally expensive than the other two. Although, with wave and flow around effects taken into consideration, it has an influence on analysis of more complex structures, which was confirmed bymany authors (Włodarczyk, 1995; Lu et al., 2003; Małachowski, 2010a,b; Mazurkiewicz and Małachowski, 2011). The carried out studies allow for choosing themost effective and reliable method for finding the responseof a structure subjected to dynamic loading.Thiswill give thepossibility for further simulations at the developing stage of safety structures. Acknowledgements The research was carried out under a research grant from Polish Ministry of Science and Higher Education No. 0097/R/T00/2010/12.This support is gratefully acknowledged. References 1. Belytschko T., Liu W.K., Moran B., 2000, Nonlinear Finite Elements for Continua and Structures, JohnWiley & Sons Ltd 2. Hallquist J.O., 2003, Theorymanual, California Livermore Software Technology Corporation 3. HydeD., 1988,User’s Guide forMicrocomputer ProgramsConWep and FUNPROApplications of TM 5- 855-1, Fundamentals of Protective Design for Conventional Weapons, US Army Engineers 4. LuY.,XuK.,LimH.S., 2003,Numerical simulationof concretebreak-upunder explosive loading, Design andAnalysis ofProtective Structures against Impact/Impulsive/Shock Loads, Proc.DAPSIL Comparison of numerical testing methods in terms of impulse ... 625 5. Małachowski J., 2010a, Influence of HE location on elastic-plastic tube response under blast loading, Shell Structures Theory and Applications, 2, 179-182 6. Małachowski J., 2010b, Modelowanie i badania interakcji ciało stałe-gaz przy oddziaływaniu impulsu ciśnienia na elementy konstrukcji rurociągu, BEL Studio,Warszawa 7. Mazurkiewicz Ł.,Małachowski J., 2011,Elastic-plastic half cylindrical surface responseunder blast loading,CMM-2011 – Computer Methods in Mechanics, Warsaw, Poland 8. US Army Corps of Engineers, 2008, User Guide for the Single degree of freedom Blast Effects Design Spreadsheet (SBEDS) 9. Włodarczyk E., 1995,Podstawy detonacji, WAT,Warszawa Manuscript received September 6, 2012; accepted for print October 30, 2012