Open access peer-reviewed chapter

Hot Compression Tests Using Total Lagrangian SPH Formulation in Energy-Based Framework

Written By

Kadiata Ba

Submitted: November 7th, 2018 Reviewed: March 18th, 2019 Published: April 10th, 2019

DOI: 10.5772/intechopen.85930

From the Edited Volume

Progress in Relativity

Edited by Calin Gheorghe Buzea, Maricel Agop and Leo Butler

Chapter metrics overview

606 Chapter Downloads

View Full Metrics


Limitations of the finite element method (FEM) in some cases involving large deformations as in forging or high compression tests are overcome nowadays by meshless methods such as the smoothed particle hydrodynamic (SPH) method. This paper presents a corrected total Lagrangian SPH formulation for problems encountering large deformations in solid mechanics. The continuum is modeled as a Hamiltonian system of particles (energy-based framework). The total Lagrangian formulation proposed overcomes some problems faced by standard SPH in simulating solid mechanic problems such as tensile instability. Numerical applications compared with experimental results are presented to show the capabilities of this novel formulation.


  • SPH
  • Hamiltonian system
  • total Lagrangian
  • thermomechanical
  • solid mechanics

1. Introduction

The use of the smoothed particle hydrodynamic (SPH) method [1, 2, 3, 4, 5] (Figure 1) in solid mechanics is quite recent (in the 1990s) compared to the SPH fluid formulation. Libersky and Petschek [6] and Libersky et al. [4] are cited as the first to use SPH in solid mechanics, for impacts modeling at high speeds and phenomena of rupture, perforation, and fragmentation. As SPH is a meshless method, there is no mesh to distort; it can efficiently handle large deformations. SPH is an efficient numerical method for applications in forging processes [7], machining [8, 9, 10], and welding [11]. Classical approach is widely used to describe SPH equations but faces drawbacks such as lack of completeness and tensile instability (numerical fragmentation). Total Lagrangian corrected SPH formulation is then used to fix the cited problems. In this paper, a Hamiltonian formulation is proposed for dynamic and steady-state problem simulation focusing on numerical efficiency such as accuracy and simulation time. The governing equations are derived following a Lagrangian variational principle leading to a Hamiltonian system of particles [12, 13, 14]. With the Hamiltonian SPH formulation, local conservation of momentum between particles is respected, and they remain locally ordered during the process as wanted in solid mechanic problems.

Figure 1.

(a) Schematic representation of the discretization of the domain Ω by a set of points i [15] and (b) seen in the space of a B-spline [16].

Total Lagrangian formulation reduces also the computational time by avoiding the search for neighboring particles for the construction of the kernel function at each time step. Through axial and lateral compression tests, the accuracy of the new formulation is shown. Results are compared to a classical formulation based on differential equations for solid mechanic applications.


2. Discrete equations of motion from energy-based formulation

The governing equations are derived following a Lagrangian variational principle leading to a Hamiltonian system of particles (energy-based) [12, 17, 18, 19] where the motion of each particle is given by the classical Lagrange equations. Therefore, as explained by Bonet et al. [18], constants of the motion such as linear and angular momentum are conserved.

For each particle, the physical quantities are calculated through interpolation over neighbor particles. Every particle is considered as a moving thermodynamic subsystem [12]. The volume of each particle is given by


where mi is the mass of the particle and ρi its density.

To proceed with a variational formulation of the equations of motion of the continuum, the kinetic K, internal int, and external energy ext of the system need to be defined.

With dissipative effects such as plasticity, the equations of motion of the system of particles representing the continuum can be evaluated following the classical Lagrangian formalism. For more details, readers can refer to Ba and Gakwaya [12]:


where xi and vi are the spatial position and velocity of the particle i, Πdis is the dissipative energy, and L is the Lagrangian given by


By substituting Eq. (3) into Eq. (2), it leads to


The total kinetic energy of the system can be approximated as the sum of the kinetic energy of each particle:


For a common case where the external forces result from a gravitational field g, the total external energy is


The total internal energy can be expressed as the sum of the products of particle masses by the amount of energy accumulated per unit mass π that depends on the deformation, density, or other constitutive parameters:


The dissipative energy can also be expressed as


where πdisp is the dissipative energy per unit mass and d is the rate of deformation tensor.


with the velocity gradient given by


where Wxixjh is the SPH kernel function and h is the smoothing length.


3. Corrected total Lagrangian SPH formulation for solid mechanics

Total Lagrangian formulation [20, 21] is well suited for solid mechanic problems as the SPH particles change less often their neighbors than in fluid mechanics [12]. The SPH kernels and their gradients are then expressed in the initial configuration (material coordinates X are used). The proposed corrected kernel is to address the lack of completeness and interpolation consistency; the smoothing length h is considered as a functional variable in the calculation of the gradient of the kernel function W [19].

Lagrangian and spatial coordinates are connected through the gradient of deformation tensor F:


where u is the displacement of a material point.

The expression of the corrected gradient of deformation tensor F, in total Lagrangian formulation, is given by


where Xj is the gradient with respect to a material point X, Vj0 is the initial volume of particle j, h0 is the initial smoothing length, and I is the identity matrix.

B is the expression of the correction of the gradient expressed as [20]


The corrected mass conservation equation for particle i is


where J and ρ0 are the Jacobian and the initial density.

The corrected momentum equation for a particle i is


where a, W˜, and fi are the acceleration, the normalized smoothing function, and the body force.

P is the first Piola-Kirchhoff stress.


where σ is the Cauchy stress tensor and FT is the inverse of the transpose of the gradient of deformation tensor.

The corrected energy conservation equation for particle i is


where e˙ is the energy rate, rpl is the mechanical contribution (heat generated by the plastic dissipation), k is the conductivity, and T is the temperature of the particle.

The equation of motion and the equation of the thermal energy of each particle can be put after discretization and evaluation of all the interactions in the following forms:


where finti and fexti are the internal and external forces, hintik and hextik are the internal and the external heat flux, and C is the capacitance matrix.

The expression for the internal force for a given particle can be expressed by differentiating the internal energy per unit mass with respect to the nodal positions as


where G is the gradient function and contains the corrected kernel gradients W at the initial reference configuration.


Internal heat flux can be expressed as


where k is the heat conductivity matrix and T the vector of nodal temperatures.

Explicit finite difference method is used to solve numerically the differential equation presented in this section through explicit dynamic algorithm to update the velocity, position, and temperature of each SPH particle.


4. Temporal integration scheme

A typical integration scheme used for integrating SPH equations is the leapfrog algorithm (Figure 2), an extension of the Verlet algorithm with low storage memory during computation.

Figure 2.

SPH code structure [23].

The heat transfer equations are integrated using the explicit forward-difference time integration rule [22].


T˙t is computed at the beginning of the increment by


The stability time is given by


where Δrmin is the smallest interparticle distance and α is the diffusivity of the material.


where k is the conductivity and s is the specific heat.

For the mechanical part, an explicit central-difference integration rule is used to integrate the equation of motion. The nodal accelerations u¨at time t is given by


where M, Pt, and It represent the mass matrix and the external and internal forces.

The integration leads to the nodal velocity u˙ (Eq. 28) and the nodal displacement u (Eq. 29).


The stable time is calculated as follows:


where Le and cd are, respectively, the characteristic length of the element and the dilatational wave speed of the material.

The whole thermomechanical problem is solved by explicit coupling; both the forward-difference (for the thermal problem) and central-difference (for the mechanical problem) integrations are explicit.

The structure of the SPH code is described below (Figure 2). The time integration routine is the main subroutine. It calculates the new variables (density, acceleration, external force, internal force).


5. Material behavior

Johnson-Cook model [24, 25, 26] is used in this work, and the flow stress is expressed as follows:


where ε is the plastic strain, ε̇s1 is the plastic strain rate, ε̇0s1 is the reference plastic strain rate, Tmis the melting temperature, Tr is the reference temperature, T is the current temperature, A is the yield stress, B is the coefficient of strain hardening, C is the coefficient of strain rate hardening, n is the strain hardening exponent, and m is the thermal softening exponent.

The material used for the simulations (see Section 6) is an Al-Zn-Mg-Cu aluminum alloy. The Johnson-Cook material parameters are shown in Table 1.

A (MPa)B (MPa)cnmε̇0 (s−1)Tm (°C)Tr (°C)

Table 1.

Johnson-Cook material parameters [12].


6. Applications

6.1 Axial compression test

A cylindrical sample (diameter, 25.4 mm; length, 25.4 mm) was subjected to the uniaxial compression test at constant velocity (2.54 mm s−1) and high temperature (400°C). Both experimental and numerical tests were performed (Figure 3). The aim of this test is to demonstrate the efficiency of the proposed total Lagrangian SPH formulation. We compared the numerical stress-strain curve with the experimental ones to verify the accuracy and the stability of the code (see Figure 4). To confirm the validity of the experimental result, the tests were repeated three times.

Figure 3.

Axial compression test setup in SPH.

Figure 4.

Stress-strain curves: experimental vs. SPH.

6.2 Lateral compression test

Lateral compression test is performed in this section (Figure 5). Eulerian and Lagrangian simulation time are compared, and tensile instability is verified. The test is carried out at 30 mm s−1, and cylindrical sample (diameter 25.4 mm, length 25.4 mm) with 5313 particles was used. This is a case of a large deformation test; the initial diameter of the sample was reduced over 50% during the test. See comparison of results at Table 2.

Figure 5.

Lateral compression test: Eulerian vs. Lagrangian.

Simulation timeNumber of particles
Eulerian SPH4 h 04 min5313
Lagrangian SPH1 h 36 min5313

Table 2.

Simulation time comparison.

6.3 Discussion

Figures 4 and 5 and Table 2 gather the tests results. From Figure 4 (axial compression test), we can see that the SPH result is very accurate compared to the experimental ones. Less than 5% of error is noted between the curves. The simulated sample shows no clustered particles, meaning there is no tensile instability.

Figure 5 and Table 2 show the results of the lateral compression test and confirm the previous result. Even in very large deformation test, particles keep their initial neighbors and do not suffer from tensile instability. In addition, the simulation time is very interesting compared to classical SPH formulation; simulation time is reduced drastically (from 4 h 04 min to 1 h 36 min); a good numerical efficiency is reached.


7. Conclusion

A corrected SPH particle approximation in energy-based framework is presented. Stability (no tensile instability), accuracy, and fast result production are shown leading to the conclusion that the total Lagrangian SPH formulation is very well suited to simulate solid mechanic problems. This is particularly interesting in simulating large deformation problems with physical fragmentation where the numerical fragmentation (tensile instability) will not corrupt the results.



The author wishes to acknowledge Augustin Gakwaya for his appreciated help.


  1. 1. Fries TP, Matthies HG. Classification and Overview of Meshfree Methods. Brunswick, Germany: Technical University Braunschweig; 2004
  2. 2. Li S, Liu WK. Meshfree and particle methods and their applications. Applied Mechanics Reviews. 2002;55(1):1-34
  3. 3. Liu GR, Liu MB. Smoothed Particle Hydrodynamics, a Meshfree Particle Method. Livre: Springer; 2005
  4. 4. Libersky L, Petschek A, Carney T, Hipp J, Allahdadi F. High strain Lagrangian hydrodynamics. Journal of Computational Physics. 1993;109:67-75
  5. 5. Wolf S. Méthode Sans Maillage. Laboratoire de mathématiques appliquées Aux systèmes; 2007
  6. 6. Libersky L, Petschek A. Smooth particle hydrodynamics with strength of materials. In: Advances in the Free-Lagrange Method Including Contributions on Adaptive Gridding and the Smooth Particle Hydrodynamics Method. Springer; 1991. pp. 248-257
  7. 7. Cleary PW, Prakash M, Das R, Ha J. Modelling of metal forging using SPH. Applied Mathematical Modelling. 2012;36:3836-3855
  8. 8. Limido J, Espinosa C. Modélisation numérique de la Coupe Orthogonale en Ugv. National Conference Proceedings; 2005
  9. 9. Limido J, Espinosa C, Salaun M, Lacome JL. SPH method applied to high speed cutting modelling. International Journal of Mechanical Sciences. 2007;49(7):898-908
  10. 10. Limido J. Étude de l’effet de l’usinage grande vitesse sur la tenue en fatigue de pièces aéronautiques [Thèse de doctorat]. Université de Toulouse; 2008
  11. 11. Timesli A. Simulation du soudage par friction et malaxage à l’aide de méthodes sans maillage [PhD dissertation]. Université de Lorraine; 2013
  12. 12. Ba K, Gakwaya A. Thermomechanical total Lagrangian SPH formulation for solid mechanics in large deformation problems. Computer Methods in Applied Mechanics and Engineering. 2018;342(1):458-473
  13. 13. Grenier N, Antuono M, Colagrossi A, Le Touzé D, Alessandrini B. An Hamiltonian interface SPH formulation for multi-fluid and free surface flows. Journal of Computational Physics. 2009;228:8380-8393
  14. 14. Hu XY, Adams NA. A multi-phase SPH method for macroscopic and mesoscopic flows. Journal of Computational Physics. 2006;213:844-861
  15. 15. Caleyron F, Combescure A, Faucher V, Potapov S. Une méthode sans maillage pour la modélisation des interactions fluide-structure: Application à la rupture d’un réservoir sous impact. 10e Colloque National en Calcul des Structures, 9-13 Mai 2011, Presqu'île de Giens (Var); 2011
  16. 16. Maurel B. Modélisation par la méthode SPH de l’impact d’un réservoir rempli de fluide [thèse de doctorat]. Institut National des Sciences Appliquées de Lyon; 2008
  17. 17. Bonet J, Rogríguez-Paz MX. Hamiltonian formulation of the variable-h SPH equations. Journal of Computational Physics. 2005;209:541-558
  18. 18. Bonet J, Kulasegaram S, Rodriguez-Paz MX, Profit M. Variational formulation for the smooth particle hydrodynamics (SPH) simulation of fluid and solid problems. Computer Methods in Applied Mechanics and Engineering. 2004;193:1245-1256
  19. 19. Lavoie MA, Gakwaya A, Nejad Ensan M. Variable-h and energy conserving SPH formulation with application in aerospace engineering. Journal Mathematics in Engineering, Science and Aerospace. 2010;1(1):27-70
  20. 20. Reveles J. Development of a total lagrangian SPH code for the simulation of solids under dynamic loading [PhD dissertation]. Cranfield University; 2007
  21. 21. Vignjevic R. Review of development of the smooth particle hydrodynamics (SPH) method. In: Predictive Modeling of Dynamic Processes. UK: Cranfield University; 2009. pp. 367-396
  22. 22. Abaqus documentation v6.14
  23. 23. Loosveldt C, Watrin D. Alternative Numerical Methods in Continuum Mechanics, Smoothed Particle Hydrodynamics. Université de Liège; 2010
  24. 24. Johnson GR, Cook WH. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In: Seventh International Symposium on Ballistics; The Hague, The Netherlands. 1983. pp. 541-547
  25. 25. Hor A. Simulation physique des conditions thermomécanique de forgeage et d’usinage-caractérisation et modélisation de la rhéologie et de l’endommagement [Thèse de doctorat]. École Nationale Supérieure d’Arts et Métiers; 2011
  26. 26. Lin YC, Li L-T, Fu Y-X, Jiang Y-Q. Hot compressive deformation behavior of 7075 Al alloy under elevated temperature. Journal of Materials Science. 2012;47(3):1306-1318

Written By

Kadiata Ba

Submitted: November 7th, 2018 Reviewed: March 18th, 2019 Published: April 10th, 2019