Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 51, 4, pp. 1053-1065, Warsaw 2013 THE APPLICATION OF THE SMOOTHED PARTICLE HYDRODYNAMICS (SPH) METHOD AND THE EXPERIMENTAL VERIFICATION OF CUTTING OF SHEET METAL BUNDLES USING A GUILLOTINE Damian Gąsiorek Silesian University of Technology, Department of Theoretical and Applied Mechanics, Gliwice, Poland e-mail: damian.gasiorek@polsl.pl Thearticlepresents the applicationof theSPHmethod for themodeling of cutting composite bundles usingaguillotine.Formanyyears, theSPHmethodhasbeenapplied in themodeling of dynamic phenomena in which rapid strain and large deformations occur. The discrete model has been performed by combining the FEMand the SPH. The SPHmethod has only been applied for to the fragment inwhich thematerial separationduring cutting appears.To verify the numerical results, experiments have been performed using a special test stand for the confirmation of the correctness of the developedmodel and the numerical simulations. Key words: SPH, cutting, guillotine 1. Introduction The problem of cutting of sheetmetal analyzed in this work is very complex, multi-dimensional and unsteady. Due to that fact, computermechanicsmethods are usually used in the simulation processes.Thesemethods include a series of numerical procedures allowing one to solve problems with many variables. The basic numerical methods used in continuum analysis include: the Euler method, Lagrange method, the SPH method (Smoothed Particle Hydrodynamics), the large particle method, the FEM method (Finite Elements Method), the rigid finite elements method, the finite volume method, the boundary element method and many other (Atkins, 2009; Gąsiorek, 2013; Hallquist, 1998; Jach, 2001; Kleiber, 1985). These methods differ mainly in the way of describing the analyzed arrangements and the method of model discretization. One of the challenges related to the analysis of the research problem in question is the scale of the described process. On one hand, there was a necessity to describe a large-scale test stand comprised of a guillotine with a measurement system, and on the other hand the necessity to develop a numerical model of the sheet metal bundle cutting process. This is presently one of the more complex problems and challenges, and concerns the so-called multiscality. The developed discretemodel gives consideration to the arrangement of sheet layers (3×30 layers) as well as the arrangement of 3 cardboard layers separating the 30-layered bundles of aluminum sheets. To solve this problem, the FEM+SPHmodel has been adopted. 2. Theoretical fundamentals of the SPH method The approach named the Smoothed Particle Hydrodynamics (SPH) is a part of the so-called gridless methods group. It has been developed to avoid considerable deterioration of the FEM grid quality in large deformations. Due to the above, the creation of an effective tool allowing for stable and precise solving of integral and differential equations should be regarded themain goal of the development of this method. The SPH, the oldest of the methods, has initially been used by astrophysicists to construct and analyze behavior models of objects in the cosmic 1054 D. Gąsiorek space (Monaghan and Lattanzio, 1985). In this method, the examined field is described by a set of particles with assigned properties and interacting with each other, behaving similarly to particles of the fluid. Their interaction is described with a smoothing function (also called a tapering function) and is stronger with the decrease of the distance between given particles. This work shall present the outline of the method, necessary to understand its principle and basic functional properties. TheSPHmethod consists of twobasic elements. Thefirst is the so-called smoothing function, which describes the parameters of state occurring in a given particle. The second is the set of discrete elements (particles), arranged in the examined space, to which the aforementioned smoothing function is assigned. In the analysis, the state parameters of a given particle are approximated by summation of the values of the smoothing functions of the nearest neighboring particles,with consideration to their distance fromtheparticle in concern (the larger thedistance of the “neighbor” is from the particle in concern, the smaller theweight is). The second element is called particle approximation.Themaximumdistance atwhich there is an interaction between particles is called the smoothing length. Thebasicdependence related to thismethod, the smoothing function, is as follows (Hallquist, 1998; Liu et al., 2003; Liu and Liu, 2010) f(x)= ∫ Ω f(x′)δ(x−x′) dx′ (2.1) where: f(x) is the searched value of the function at a given x point, Ω – examined area, δ – Dirac’s delta, f(x′) – value of the function at points other than x. The above is an interpolation function, allowing for the approximation of the searched values and their derivatives in the Ω examined area/space. Equation (2.1) imposes limitations on the function f,whichmustbeadjusted so as to be integrable in the analyzed, described Ω area.Due to theoccurrenceof theDiracdelta,which isnon-continuousandnon-differentiable, this equation may not be directly used in numericalmodels. Due to the above, the δ function is replacedwith the function W which smoothes the former by introducing an additional dependence of the h variable, which is a finite distance from the particle in concern (Hallquist, 1998; Liu et al., 2003; Liu and Liu, 2010). The replacement of the δ function with the function W is called the nucleus approximation, and it is usually marked as 〈·〉 f(x)≃〈f〉(x)= ∫ Ω f(x′)W(x−x′,h) dx′ (2.2) The function W must fulfill a series of conditions. A simplified and abbreviated description of the smoothing function conditions is often presented in scientific literature. Themost commonly provided conditions are (Liu et al., 2003; Liu and Liu, 2010; Monaghan and Lattanzio, 1985): • the normalization condition: ∫ ΩW(x−x ′,h) dx′ =1 • approaching a value equal to the Dirac delta value at h (smoothing length) approaching zero: limh→0W(x−x ′,h)= δ(x−x′,h) • assuming positive values: W(x−x′)­ 0 • parity, thus ensuring equal interaction of a given particle and other particles at the same distance, yet different directions • highest possible regularity. Any function which fulfills the above conditions may be a smoothing function. In practice, however, the B-spline function, first used byMonaghan and Lattanzio (1985), next byHallquist (1998) and Liu et al. (2003), is applied most often The application of the smoothed particle hydrodynamics ... 1055 W(R,h)=αd    2 3 −R2+ 1 2 R3 for 0¬R< 1 1 6 (2−R)3 for 1¬R< 2 0 for R­ 2 (2.3) where: αd =1/h, 15/(7πh 2), 3/(2πh3) respectively for one-, two-, and three-dimensional space. The second step necessary to develop a SPHmodel is the division of the Ω area/space into particles to which the smoothing function is assigned as provided in a diagram in Fig. 1. The value of the f function (2.3) may then be expressed as 〈f〉(x)= ∫ Ω f(x′)W(x−x′,h) dx′ ≃ N∑ q=1 f(xq)W(x−xq,h)∆Vq (2.4) where: N is the number of particles present within the limits of the support of the smoothing function, ∆Vq – part of the area covered by the smoothing function. Fig. 1. The division of the Ω area into qi particles Assuming that a particle is a finite, concentrated volume (sub-area), it may be presented as a ratio of the particle mass to its density: mi/ρi. In the conducted tests, Lagrange’s description of particlemotion has been assumed. Following the process of the discretization of the examined Ω area/space with a finite number of particles, equation (2.2) may be presented as follows (Hallquist, 1998; Liu et al., 2003; Liu and Liu, 2010) 〈f〉(x)= N∑ q=1 mq ρq f(xq)W(x−xq,h)∆Vq (2.5) where: m is the particle mass, ρ – density of the particle, N – number of particles within the limits of support of the smoothing function. The above approach using the SPHmethod requires also the application of the principles of conservation of mass, momentum and energy. These principles assume the following form (Liu et al., 2003; Liu and Liu, 2010) ∂ρp ∂t = N∑ q=1 mqvpq ∂Wpq ∂xp ∂vp ∂t = N∑ q=1 mq (σp ρ2p + σq ρ2q )∂Wpq ∂xp ∂up ∂t = N∑ q=1 mq σpσq ρpρq vpq ∂Wpq ∂xp (2.6) where: vi = dxi/dt, σ – stress, v – speed, ij – indices indicating the components, W – internal energy. 1056 D. Gąsiorek Due to such an approach, both the approximation of solids and fluids is possible, for the description of which such a method has been developed and implemented into computer codes. The coupling of the SPHandFEMareasmaybedone in twoways: byusinghybrid elements, in which the SPH andFEMdescriptions are coupled or by using a typical contact algorithm. In this work, a penalty contact algorithm has been applied to describe the interaction between the SPH area which undergoes the damaging process (loss of mass integrity), and the area/space described by FEM. 2.1. Constitutive model adapted for numerical calculations In order to describe the behavior of an aluminumsheetmetal, the Johnson-Cook constitutive model has been applied. Themodel allows one to describe plastic-elastic media with considera- tion given to the Huber-Mises-Hencky (H-M-H) plasticity conditions and to relate them to the speed of deformation. The simplified parameters for this model are presented in Table 1. Table 1. Parameter values of the Johnson-Cook model for aluminum sheet metal (Iqbal and Gupta, 2008) Density Kirchhoff’s A [MPa] B [MPa] Material ρ modulus C n [kg/m3] G [GPa] Aluminum 2787.0 28.0 156.0 345.51 0.001 0.18 To take into consideration the damage process – the loss of material integrity, the use of an extended damage accumulation model has been proposed, which is based upon the so-called equivalent value of the εpl eff failure strain. The value gives consideration to the triaxial state of tensions in the material, the effect of speed of deformations and the effect of temperature. The form of the equation is a representation of the Johnson-Cookmodel (Atkins, 2009) ε pl eff = [ D1+D2exp ( D3 σm σ̄ )]( 1+D4 ln ε̇ ε̇ ) (1+D5T̂) (2.7) where: D1-D5 are material constants and σm – mean stress. Table 2.Parameter values for the aluminum failure model (Atkins, 2009) Material D1 D2 D3 D4 D5 m Aluminium 0.07 1.24 −1.14 0.14 0.0 0.85 A plastic-elastic model of the material (the so-called bilinear model) has been applied to describe the behavior of the cardboard layers. The value of stress in the material is determined in line with the equation (Hallquist, 1998) σy =σo+β 3 2 Epε p eff (2.8) where: σo – yield point, β – strengthening parameter assumming values from 0 to 1 range, Ep – strengtheningmodulus and ε p eff = t∫ 0 √ 2 3 ε̇Pijε̇ P ij dt where ε̇p is the speed of the plastic deformations. The application of the smoothed particle hydrodynamics ... 1057 Table 3. Parameter values for the plastic-elastic model for a layer of cardboard (Giampieri et al., 2011; Murayama et al., 2005) Material E [MPa] σy [MPa] Cardboard 3626.0 26.4 2.2. Description of the solution method of the dynamic equilibrium equation A variational approach has been adopted in the work to solve the analyzed problem of numerical representation of the sheet metal cutting (the so-called weak formulation). The weak form of the equilibrium equation is derived using Hamilton’s variational principle, according to which the state of an elastic body may be described as a sum of kinetic energy, deformation energy and potential energy of the external forces acting on the body in concern (Dacko et al., 1994; Kleiber, 1985) Π =Ek−Ew+Ep = 1 2 ∫ Ω ρu̇Tu̇ dΩ− 1 2 ∫ Ω ε T σ dΩ+ ∫ Ω uTb dΩ+ ∫ Γ uTf dΓ (2.9) where: Ω denotes the volume of the body in concern, Γ – area of the bodyatwhich the external forces are acting, f – vector of the external forces (not body forces) acting on the body. The determination of the approximated solution of the boundary value problem with uni- form conditions for a differential equation consists in searching the minimum of the stationary functional ∓, which is a necessary condition for optimality. With the application of Hamilton’s principle, we have (Małachowski and Łazowski, 2007; Rakowski and Kacprzyk, 2005) δΠ =0 (2.10) which leads to the following dependence ∫ Ω δεTσ dΩ− ∫ Ω δuTb dΩ− ∫ Γ δuTf dΓ + ∫ Ω ρδuTü dΩ =0 (2.11) Physically, this conditionmeans that the total energyof the systemmustbe stationary (reach the minimum) at any displacement variation. This equation may be solved using various methods, such as the FEM. If the Ω body in concern is divided into smaller volumes (elements/ mass points) with predefined shapes describedwith simple polynomial functions, it will be possible to replace the stationarity conditionwith a sumof stationarity conditions, also in smaller sub-areas. This allows us to construct, respectively, the matrices of mass and rigidity (Rakowski and Kacprzyk, 2005) M= ∫ Ω ρNTN dΩ K= ∫ Ω BTCB dΩ (2.12) where: N is the matrix of interpolation functions/shape functions, M – mass matrix, K – rigidity matrix. By fulfilling condition (2.10) in the individual sub-areas andby solving system (2.11) (that is, by determining the U displacement vector) and with consideration to the boundary conditions and remembering that the problem in concern is changing in time, we finally reach the known equation of equilibrium (Dacko et al., 1994; Rakowski and Kacprzyk, 2005) MÜ+KU=P (2.13) Thephenomenonofdamping is an extremely importantandcomplexproblem in theexamination of the dynamic response of structural systems. When the displacement of a body changes over 1058 D. Gąsiorek time, two groups of additional forces occur (Hallquist, 1998). The first group is the result of the viscous drag of themedium, while the second is the result of the internal friction. The damping forces arising in this way are of two different natures (Hallquist, 1998). Assuming the above, the element describing damping forces should be included in equation (2.13), which then results in the complete equation of dynamic equilibrium (Dacko et al., 1994) MÜ+CU̇+KU=P (2.14) where: C = α1M+α2K is the damping matrix, α1, α2 – factors determined based upon the modal determinations (Rakowski and Kacprzyk, 2005). Theuniversalmethodwhich solves dynamic equilibriumequation (2.14) is a direct numerical integration method. The equation of motion is integrated relative to time (the variable t) using the step method in consecutive moments. A solution in the range [t0, tn) is searched. The complete solution is achieved by assuming its division into n number of intervals of ∆t and thus, n∆t= tn−t0. In themethod of numerical integration of the equation ofmotion, the result is searched in form of a recurrence formula. The state of displacement is determined (speed, acceleration) based on the displacement state (speed, acceleration) in the previous moments. The result is obtained by determining the searched values in consecutive moments, using a step-by-step method. It is assumed that the displacements, speeds and accelerations of the system are known for the initialmoment, that is t= t0 and these are equal to x0, ẋ0, ẍ0.Naturally, the acceleration ẍ0 is not an initial condition and is calculated based on the x0, ẋ0 conditions (e.g. by the equation of motion). To express the speed ẋ and the acceleration ẍ in the displacement x function, a differential approximation of time derivatives is applied (Kleiber, 1985; Małachowski and Łazowski, 2007; Rakowski and Kacprzyk, 2005). The following is assumed in the central differences method (Hallquist, 1998; Małachowski and Łazowski, 2007; Rakowski and Kacprzyk, 2005): — speed ẋ n+1 2 = 1 ∆t n+1 2 (xn+1−xn) (2.15) — acceleration ẍn = 1 ∆tn ( x n+1 n −x n−1 n ) (2.16) while the equationofmotion for the timemoment tn for thenon-linear case assumes the following form Mẍ=Pn−F int n −Cẋn (2.17) where: ẋn is an unknown speed value at the tn moment. Simultaneously, an asynchronous damping effect is assumed in the following form Cẋn−1 n ≈Cẋn (2.18) The updated value of acceleration is described by the equation ẍn =M −1 ( Pn−F int n −Cẋn−1 n ) (2.19) The updated speed and displacement formulas are as follows x n+ 1 n = ẋ n− 1 n +∆tnẍn xn+1 =xn+∆tn+1 2 ẋ n+ 1 2 (2.20) The application of the smoothed particle hydrodynamics ... 1059 The above procedure is called the explicit integration method, and to find the xt+∆t displa- cement value, a condition of dynamic equilibrium in the tmoment is used. All analyses may be conducted at the level of subsequently processed elements, which increases the speed of calcula- tions and fundamentally decreases the required internalmemory of the computer (Kleiber, 1985; Małachowski and Łazowski, 2007; Rakowski and Kacprzyk, 2005). This allows for an effective solution of systems with a large number of degrees of freedom on relatively small machines. The disadvantage of this method, however, is the lack of non-conditional stability of the solution. The stability of the above system of equations is dependent upon the necessity of the step length selection after ∆t, so as it was lower or equal to ∆tcr critical time, which is dependent upon the properties of the entire examined system (Kleiber, 1985; Małachowski and Łazowski, 2007; Rakowski and Kacprzyk, 2005). Such stability may be achieved through the fulfillment of the Courant-Friedrichs-Levy criterion (Kleiber, 1985; Małachowski and Łazowski, 2007) ∆ts = Le Q+ √ Q2+ c2 (2.21) where Q is the function of the viscous coefficients C0 and C1 formulated as follows Q= { C1c+C0Le|ε̇kk| for ε̇kk < 0 0 for ε̇kk ­ 0 (2.22) where: Le = υe/Aemax is the characteristic length of the element, υe – volume of the element, Aemax – area of the largest side of the finite element, c – the local speed of sound. The above analysis is conducted for all bodies of a solid medium described in a discrete manner inLagrange’s coordinate system, and theminimal value is searched for. Such an analysis concerns the bodies described using the FEM, the SPHmedium and with consideration to the contact rigidity. Based upon the performed analyses for the three aforementioned groups (FEM, SPH and contact rigidity), the minimum value is searched for in the current time step, in line with the following dependence ∆tn+1 = amin(∆tMES,∆tSPH,∆tcontact) (2.23) Additionally, to improve the stability, a scaling factor in form of a constant a is introduced. Its value formechanical problems is 0.9 (Hallquist, 1998). In other cases, there is a high probability that the stability of the result will not be achieved, and thus it will be impossible to achieve the result for the final time. 2.3. Numerical model in the FEM+SPH method In the suggested approach to the modeling of the cutting process, a coupling of the FEM model and themodel based on hydrodynamic particles (the SPHmethod) has been proposed. It has been decided tomodify the presentedFEMmodel by separating an area/space inwhich the- re is a direct interaction between the guillotine cutter and the sheet metal, and by describing it with the SPHmethod.Themain goal that has been assumed is to performamore precise repre- sentation of the integrity loss process of the aluminum sheetmaterial under the influence of the moving cutter.Due to the above, andwhile understanding the advantages of the aforementioned method – especially in the description of large deformations and structure defragmentation – the guillotine blade contact area and the decohesion area has been described with a regular “grid” of 18000 SPH particles (Fig. 2). The remaining area/space has been described, as in the first model, with a grid of 8-node finite elements. Both spaces, that is the the discrete FEM space and the gridless SPH space, have been coupled using the kinematic constraints method. 1060 D. Gąsiorek Fig. 2. FEM+SPH coupledmodel for the performance of numerical tests of the aluminum sheets cutting process 2.4. Numerical tests of the cutting process as described by Lagrange from the perspective of the SPH method While observing the process of aluminum plate deformation (Fig. 3), it may be noticed that the plates are damaged mostly locally and no raising of the plates occurs as a result of the elastic response. The only observable effects are high local accumulations of deformations. In the process, the layers move vertically as a result of the cutter movement towards the plates. In the area of direct interaction with the cutter, the process of material decohesion occurs, in the area of which the phenomenon of strain concentration is present. This value reaches the level of approximately 500MPa (Fig. 4) at the face of the cutter blade in sheet metal material and is caused by viscous processes occurring in the matter. Fig. 3. Deformation of sheet metal obtained in the process of numerical cutting for the FEM+SPH model in selected moments of time Fig. 4. Reduced stress according to H-M-H in the contact area of the guillotine blade and the layers of the cut aluminum plate in consecutivemoments of time The application of the smoothed particle hydrodynamics ... 1061 Fig. 5. Behavior of the separated particles of the cut material in the SPHmethod In the analyses, for realization of thedecohesionprocess, thedeformation failure criterion has been applied. In the further analysis, the decohesion process of the elements of the equilibrium equation does not occur, and – as presented in Fig. 6 – after reaching the critical value, only the process of integrity loss occurs between individual particles. In the cutting area, material concentration takes place. The performed analysis of deformation speed has also confirmed the occurrence of the aforementioned viscous effects in the plate material. The value, however, is nearly three times lower and amounts to approximately 1 ·102 s−1. Fig. 6. Characteristics of the plasticization process for a selected SPH particle describing the aluminum plate in the cutting area, up to the moment of decohesion (particle deactivation) Fig. 7. Diagram of the deformation speed for a selected SPH particle describing the aluminum plate in the cutting area, up to the moment of decohesion (particle deactivation) The performed energy analysis has confirmed (Fig. 8) the facts mentioned above. The total increase of energy was nearly completely related to the increase of deformations in the integrity loss area. Theminimal increase of the kinetic energy in the final stage of the simulation results from the process of horizontal movement of a larger number of aluminum plates. 1062 D. Gąsiorek Fig. 8. Energy balance characteristics [10−6kJ] The analysis of the resultant force (Fig. 9) in the base of the stand shows no visible effects in the first stage (especially the typical drops to zero value). This is caused by amore accurate representation of layers with a larger number of SPH particles, which, in consequence, causes that the cutter contacts next particles immediately after cutting a layer. It is also related to the nature of the method, in which, as presented before, the radius of interaction and the shape of the smoothing function are decisive for the influence of particle interaction. Fig. 9. Characteristics of the force measured in the indicated point, resulting from the sheet metal cutting process for the FEM+SPHmodel This shape of results is considerablymore resembling of the results obtained at the test stand in form of a registered resultant force in the base. 3. Experiments To verify the numerical model, tests have been performed using an original test stand. Similar tests have been described in works by Kaczmarczyk et al. (2007), Mężyk et al. (2010). The difference between the stands consists, most importantly, in the placement of a series of sensors, which unequivocally map the initial boundary conditions assumed in the cutting process. The cutting tests have been performed in line with the numerical model description, for a bundle of printing plates (3×30pcs), separated every 30th plate with a cardboard separator. 3.1. Experiment results As a result of the tests, waveforms of forces under the test stand in time have been obtained, which are presented in Fig. 10. The application of the smoothed particle hydrodynamics ... 1063 Fig. 10. Dependence between the cutting force and the blade displacement In the case of cutting a bundle of plates, the visible amplitudes correspond to the cutting force registered at an actuator placed at the bottom of the test stand. During the cutting of a cardboard, a significant decrease of the force necessary to cut the bundle is visible. In a further step, a comparison of cutting forces obtained by means of the SPH numerical method and using the test stand was presented (Fig. 11). Fig. 11. A comparison of the resultant force of the experiment and the FEM+SPHmethod Fig. 12. Characteristics of the force measured in the indicated point – comparison of the results 4. Summary of the numerical research Based upon the performed numerical tests in the field of numerical cutting of the sheet metal simulation process using the FEM+SPHmethods, the following conclusions may be reached: • the developed model is capable of estimating a mean reaction value resulting from the cutting process, 1064 D. Gąsiorek • due to amore precise description of the deformation field in the contact area of the cutter blade and the material, it is possible to create a numerical representation of the local decohesion process and accompanying phenomena by the coupling of the global FEM model with the local SPH description, • the application of advanced constitutive models which give consideration to the plastic, elastic and viscous traits with an implemented criterion of material deformation accumu- lation, allows one to consider many effects accompanying the cutting process, including the integrity loss processes of the structure, • the developed model allows one to conduct further sensitivity-oriented analyses of the developedmethodsand the search for optimalparameters for cutting sheetmetal, including the stand parameters. To recapitulate, it must be stated that the application of an advanced numerical apparatus, presently allowing one to couple numerical methods (e.g. in the case of the FEM analysis and theSPHmethod), enables one to obtain a faithful representation ofmanyphenomenaand search for results similar to the results obtained using the test stand on that basis. Theproposedmodel is basedupon the so-called Lagrange description,whichdoes not findan application in the description of bodies characterized by low rigidity. This, however, opens a new area of work in further research, in which bodies and phenomena/processesmay be described in Euler’s system.This is significant due to the possibility of numerical description of bodieswhich undergo high deformations in dynamic reaction processes. Such works are presently the latest direction of research and numerical analysis in this dynamically developing branch, namely the computer mechanics. Acknowledgement The work has been accomplished under the research project No. N N503 047040 financed by the Scientific Research Committee. References 1. AtkinsT., 2009,The Science and Engineering of Cutting, Butterworth-Heinemann,GreatBritain 2. Bajer C., 1997,Numerical Modelling of Dynamic Space-time Contact Problems (in Polish),Wy- dawnictwo IPPTPAN,Warszawa 3. Belytschko T., Liu W.K., Moran B., 2000, Nonlinear Finite Elements for Continua and Structures, JohnWiley & Sons, England 4. Dacko M., Borkowski W., Dobrociński S., Niezgoda T., Wieczorek M., 1994, Finite Elements Method in Construction Mechanics (in Polish), Arkady,Warszawa 5. Gąsiorek D., 2013, Modelling and Experimental Investigation of Dynamic Processes Occurring During Cutting Lithograpic Plates Using Guillotines (in Polish), Wydawnictwo Politechniki Ślą- skiej, Gliwice 6. Giampieri A., Perego U., Borsari R., 2011,A constitutivemodel for themechanical response of the folding of creased paperboard, International Journal of Solids and Structures,48, 2275-2287 7. Hallquist J.O., 1998, LS-Dyna. Theoretical Manual, California Livermore Software Technology Corporation 8. Iqbal M.A., Gupta N.K., 2008, Energy absorption characteristics of aluminum plates subjected to projectile impact, Latin American Journal of Solids and Structures, 5, 259-287 9. JachK., 2001,ComputerModelling ofDynamic Interactions ofBodiesUsingFreeParticlesMethod (in Polish),WydawnictwoNaukowe PWN,Warszawa The application of the smoothed particle hydrodynamics ... 1065 10. Kaczmarczyk J., Gąsiorek D., Mężyk A., 2007, Numerical analysis of defect formation in guillotine plate cutting process (in Polish),Modelowanie Inżynierskie, 3, 34, 61-66 11. KleiberM., 1985,Finite ElementsMethod inNon-linearContinuumMechanics (inPolish),PWN Warszawa-Poznań 12. Liu M.B., Liu G.R., 2010, Smoothed particle hydrodynamics (SPH): an overview and recent developments,Archives of Computational Methods in Engineering, 17, 25-76 13. Liu M.B., Liu G.R., Lam K.Y., 2003, Constructing smoothing functions in smoothed particle hydrodynamics with applications, Journal of Computational and Applied Mathematics, 155, 2, 263-284 14. Małachowski J., Łazowski J., 2007, Numerical aspects of barrel-bullet modeling (in Polish), Przegląd Mechaniczny, 6, 33-38 15. Mężyk A., Gąsiorek D., Kaczmarczyk J., Rak Z., Skibniewski A., 2010, Experimental research on guillotine cutting of plate stacks process (in Polish),Przegląd Mechaniczny, 10, 5 16. Monaghan J.J., Lattanzio J.C., 1985, A refined particle method for astrophysical problems, Astronomy and Astrophysics, 149, 1, 135-143 17. MurayamaM.,NagasawaS.,FukuzawaY.,YamaguchiT.,Katayama I., 2005,Orthotropic effect and strain dependency of paperboard on load characteristic of center bevel cutter indented on paperboard, Journal of Materials Processing Technology, 159, 199-205 18. Rakowski G., Kacprzyk Z., 2005,Finite Elements Method in Construction Mechanics (in Po- lish),Wydawnictwo PolitechnikiWarszawskiej Manuscript received February 11, 2013; accepted for print June 18, 2013