Numerical Simulation of Metal Machining Process with Eulerian and Total Lagrangian SPH

TL;DR

This study compares Eulerian and Total Lagrangian SPH methods for simulating metal machining, demonstrating that TLSPH offers improved stability and accuracy in modeling strain localization and chip formation.

Contribution

The paper introduces the application of TLSPH to metal machining, highlighting its advantages over ESPH in stability and simulation fidelity.

Findings

Both SPH methods capture key metal processing phenomena.

TLSPH provides better strain localization and chip morphology.

TLSPH is more stable and efficient for machining simulations.

Abstract

This paper presents numerical simulations of metal machining processes with Eulerian and Total Lagrangian Smoothed Particle Hydrodynamics (SPH). Being a mesh-free method, SPH can conveniently handle large deformation and material separation. However, the Eulerian SPH (ESPH) in which the kernel functions are computed based on the current particle positions suffers from the tensile instability. The Total Lagrangian SPH (TLSPH) is free of this instability as the kernel functions are calculated in the reference configurations. In this work, the metals are modelled using the Johnson-Cook constitutive model, which can capture strain hardening and thermal softening in metals. The processing/cutting tools are modelled as rigid bodies, while the metal-tool contact forces are considered using the standard SPH interaction and the particle-particle pinball contact in ESPH and TLSPH, respectively. The two methods are employed to model several cases with impact, pressing, and cutting; the results are compared with reference experimental and numerical results. It is found that both the two SPH methods can capture the salient phenomena in metal processing, e.g. strain localisation, large deformation, and material separation. However, the TLSPH approach provides a better simulation of strain localisation and chip morphology. This work shows that the TLSPH method has the potential to model the metal machining processes efficiently without any numerical instabilities.

Tables

Table 1. Table 1 : Johnson-Cook (J-C) Material parameters of low strength 4340 steel

Mechanical properties	J-C strength parameters	J-C thermal parameters
$\rho =7830$ kg/m3	$A$= 792 MPa	$T_r$= 293 K
$\mu =82.9$ GPa	$B$= 510 MPa	$T_m$= 1293 K
$K=169.1$ GPa	$n$= 0.26	$m$= 1.03
	$C$= 0.014	$C_p=460$ J/(kg$\cdot$K)

Table 2. Table 2 : Computational parameters for Taylor impact test

Inter-particle	Smoothing	Time step	Artificial viscosity
spacing ($\Delta p$)	length ($h$)	($\Delta t$)	parameters (${\beta }_1,{\beta }_2$)
0.5 mm	0.75 mm	$1\times {10}^{-8}$ s	(1.5, 1.5)

Table 3. Table 3 : Johnson-Cook (J-C) Material parameters for AISI4340 steel workpiece

Mechanical properties	J-C strength parameters	J-C thermal parameters
$\rho =7830$ kg/m3	$A$= 792 MPa	$T_r$= 298 K
$E=207$ GPa	$B$= 510 MPa	$T_m$= 1573 K
$\nu =0.3$	$n$= 0.26	$m$= 1.03
	$C$= 0.014	$C_p=480$ J/(kg$\cdot$K)

Table 4. Table 4 : Computational parameters for metal pressing

Inter-particle	Smoothing	Time step	Artificial viscosity
spacing ($\Delta p$)	length ($h$)	($\Delta t$)	parameters (${\beta }_1,{\beta }_2$)
0.00667 mm	0.01 mm	$1\times {10}^{-10}$ s	(0.5, 0.5)

Table 5. Table 5 : Computational parameters for metal cutting

Inter-particle	Smoothing	Time-step	Artificial viscosity
spacing ($\Delta p$)	length ($h$)	($\Delta t$)	parameters (${\beta }_1,{\beta }_2$)
0.01 mm	0.018 mm	$1\times {10}^{-10}$ s	(0.5, 0.5)

Equations

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}=-\rho \nabla \cdot \mathbf{v}$$

$$\rho =J^{-1}{\rho }_0$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=\frac{1}{\rho }\nabla \mathbf{\sigma }$$

$$\frac{\mathrm{d}v}{\mathrm{d}t}=\frac{1}{{\rho }_0}{\nabla }_0\cdot \mathbf{P}$$

$$\mathbf{F}=\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}\mathbf{X}}=\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}\mathbf{X}}+\mathbf{I}$$

$$\dot{\mathbf{F}}=\frac{\mathrm{d}\dot{\mathbf{x}}}{\mathrm{d}\mathbf{X}}=\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}\mathbf{X}}$$

$$f({\mathbf{x}}_i)=\sum _{j}f({\mathbf{x}}_j)W({\mathbf{x}}_{ij})\frac{m_j}{{\rho }_j}$$

$$\nabla f({\mathbf{x}}_i)=-\sum _j\left[f({\mathbf{x}}_i)-f({\mathbf{x}}_j)\right]{\nabla }_{i}W({\mathbf{x}}_{ij})\frac{m_j}{{\rho }_j}$$

$$W(q,h)={\alpha }_d\{\begin{array}{ll}1-\frac{3}{2}{q}^2+\frac{3}{4}{q}^3,& \text{if }0\le q\le 1\\ \frac{1}{4}(2-q{)}^3,& \text{if }1\le q\le 2\\ 0& \text{otherwise}\end{array}$$

$$\frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}=\sum _{j=1}{m}_j\left({\mathbf{v}}_i-{\mathbf{v}}_j\right)\cdot {\nabla }_{i}W({\mathbf{x}}_{ij})$$

$$\frac{\mathrm{d}{\mathbf{v}}_i}{dt}=\sum _{j=1}{m}_j\left(\frac{{\mathbf{\sigma }}_i}{{\rho }_i^2}+\frac{{\mathbf{\sigma }}_j}{{\rho }_j^2}-{\pi }_{ij}\mathbf{I}-p_{ij}^a\mathbf{I}\right){\nabla }_{i}W({\mathbf{x}}_{ij})$$

$${\pi }_{ij}=\{\begin{array}{ll}\frac{-{\beta }_1{\bar{C}}_{ij}{\mu }_{ij}+{\beta }_2{\mu }_{ij}^2}{{\bar{\rho }}_{ij}},& \text{if }{\mathbf{v}}_{ij}\cdot {\mathbf{x}}_{ij}\le 0\\ 0,& \text{otherwise}\end{array}$$

$${\mu }_{ij}=\frac{h{\mathbf{v}}_{ij}\cdot {\mathbf{x}}_{ij}}{{\mathbf{x}}_{ij}^2+ϵh^2}$$

$$p_{ij}^a=\gamma \left(\frac{|p_i|}{{\rho }_i^2}+\frac{|p_j|}{{\rho }_j^2}\right){\left[\frac{W({\mathbf{x}}_{ij})}{W(\Delta p)}\right]}^n$$

$$f({\mathbf{X}}_i)=\sum _{j}f({\mathbf{X}}_j)W({\mathbf{X}}_{ij})\frac{m_j}{{\rho }_{0j}}$$

$$\nabla f({\mathbf{X}}_i)=-\sum _j\left[f({\mathbf{X}}_i)-f({\mathbf{X}}_j)\right]{\nabla }_{i}W({\mathbf{X}}_{ij})\frac{m_j}{{\rho }_{0j}}$$

$$\frac{\mathrm{d}{\mathbf{v}}_i}{\mathrm{d}t}=\sum _j{m}_j\left(\frac{{\mathbf{P}}_i}{{\rho }_{0i}^2}+\frac{{\mathbf{P}}_j}{{\rho }_{0j}^2}-{\mathbf{\Pi }}_{ij}\right){\nabla }_{i}W({\mathbf{X}}_{ij})$$

$$\begin{gathered} {\mathbf{F}}_i=-\sum _j\left({\mathbf{u}}_i-{\mathbf{u}}_j\right)\otimes {\nabla }_{i}W({\mathbf{X}}_{ij})\frac{m_j}{{\rho }_{0j}}+\mathbf{I} \\ =-\sum _j\left({\mathbf{x}}_i-{\mathbf{x}}_j\right)\otimes {\nabla }_{i}W({\mathbf{X}}_{ij})\frac{m_j}{{\rho }_{0j}} \end{gathered}$$

$${\dot{\mathbf{F}}}_i=-\sum _j\left({\mathbf{v}}_i-{\mathbf{v}}_j\right)\otimes {\nabla }_{i}W({\mathbf{X}}_{ij})\frac{m_j}{{\rho }_{0j}}$$

$$\dot{\mathbf{s}}=2\mu \left(\dot{\mathbf{ϵ}}-\frac{1}{3}\dot{\epsilon }\mathbf{I}\right)+\mathbf{s}\mathbf{\omega }-\mathbf{\omega }\mathbf{s}$$

$$\dot{\mathbf{ϵ}}=\frac{1}{2}\left(\mathbf{l}+{\mathbf{l}}^{\mathrm{T}}\right)$$

$$\mathbf{\omega }=\frac{1}{2}\left(\mathbf{l}-{\mathbf{l}}^{\mathrm{T}}\right)$$

$$\mathbf{l}=-\sum _j({\mathbf{v}}_i-{\mathbf{v}}_j)\cdot {\nabla }_{i}W({\mathbf{x}}_{ij})\frac{m_j}{{\rho }_j}$$

$$p=K\left(\frac{\rho }{{\rho }_0}-1\right)$$

$${\sigma }_y=[A+B{ϵ}_{pl}^n][1+C\log ({\dot{ϵ}}_{pl}^{\ast })][1-T{{}^{\ast }}^m]$$

$${\dot{ϵ}}_{pl}^{\ast }=\frac{{\dot{ϵ}}_{pl}}{{\dot{ϵ}}_0}$$

$$T^{\ast }=\frac{T-T_0}{T_m-T_0}$$

$$\frac{\mathrm{d}T}{\mathrm{d}t}=\chi \frac{W_p}{\rho C_p}$$

$$\Delta {\mathbf{ϵ}}_{pl}=\frac{1-c_f}{2\mu }\mathbf{s}$$

$$\Delta {ϵ}_{pl}=\sqrt{\frac{2}{3}\Delta {\mathbf{ϵ}}_{pl}:\Delta {\mathbf{ϵ}}_{pl}}=\frac{1-c_f}{3\mu }\sqrt{\frac{3}{2}\mathbf{s}:\mathbf{s}}$$

Full text

††journal: Computational Particle Mechanics

1 Introduction

Metal machining processes such as cutting, pressing, forming, drilling and milling are one of the most important and widely-used procedures in the manufacturing industry. The optimisation of tools, specimens and processing flows are essential for achieving higher performance in manufacturing. However, the nature of the operation makes it a challenging task for numerical modelling. During processing, the specimen/workpiece undergoes finite deformation and the behaviour of the material changes based on the accumulated strains and strain rates. The machining procedure also involves heat generation; hence, the thermal softening behaviour needs to be considered. This entire process becomes even more complicated with the possibility of strain localisation, cracking, and material separation. Experimental and empirical observations [1, 2, 3] are used to understand the metal cutting. Simple analytical models [4, 5] are also used. However, the cost and complexity of the experiment limit its practical application, whereas the analytical models are less reliable due to simplifications. Therefore, the numerical models are used as an alternative to model the metal machining process.

As the most widely-used numerical method, the finite element method (FEM) was employed to model orthogonal metal cutting in [6]. However, the FEM has inherent difficulties in modelling problems involving finite deformations and material separations due to the utilisation of mesh. Mesh distortion and discontinuities require a special technique to handle. Several studies employed the Arbitrary Lagrangian-Eulerian (ALE) method [7, 8] to treat the excessive element distortion and material separation, in which the Johnson-Cook material model [9] is used to account for the influence of plastic strain, strain rate and thermal softening. In ALE, the issues of element distortion, entanglement, and material separation are handled through continuous re-meshing. However, the re-meshing is not only computationally intensive but also difficult to perform in three-dimensional complex cases. Moreover, it is well-known that FEM suffers from mesh-dependency in modelling strain localisation, a dominant mechanism in metal processing. For instance, the influence of mesh size and orientation in metal cutting are investigated in [10]. It is found that the mesh influences the formation of shear bands, chips and the values of the computed cutting force. It is, therefore, challenging to obtain mesh-objective results with the ALE approach.

Lagrangian particle-based methods are naturally capable of capturing large material deformation due to the absence of meshes. The issues related to mesh, e.g. distortion, entanglement are not encountered. Therefore, some of the meshless methods are employed to model metal processing. For instance, the Material Point Method (MPM) is used to simulate the metal pressing/cutting process in [11]. Another popular meshless method for metal cutting modelling is Smoothed Particle Hydrodynamics (SPH) [12, 13, 14, 15, 16, 17, 18, 19]. SPH is a truly meshless method as it does not require any background mesh. It was first developed for the astronomical problems [20, 21] and later gradually used for the simulation of fluid flows and solid deformations [22, 23]. In SPH, the computational domain is discretised into Lagrangian particles carrying field variables and moving with the material. The meshless and Lagrangian nature of SPH makes it a perfect candidate for metal processing modelling. It is shown in previous works that shear bands and formation of chips are captured naturally without any requirement of special treatment such as nodal enrichment or re-meshing. SPH based simulations in LS-DYNA are performed [13] to predict the cutting forces, plastic strains and cutting planes with reasonable accuracy. A coupled FEM and SPH framework is also developed in [16].

In conventional SPH, particle positions are updated after each computational step. The kernel functions are computed based on the updated particle positions; therefore, the traditional SPH (ESPH) kernel function is termed as Eulerian kernel function in the sense that particles can move in and out of its influence domain. Although widely used, the Eulerian kernel function leads to the well-known tensile instability [24], i.e. a local clustering of particles forming unphysical numerical fracture. A popular correction is the artificial pressure/stress [25, 26], which can effectively alleviate the tensile instability. However, this approach requires a proper tuning of the parameter to avoid any unphysical effect in the simulation. This tuning process is case-dependent thus highly undesirable. Studies reveal that tensile instability is associated with the use of Eulerian kernel function [27]. It is found that this instability can be avoided if the reference particle position is used to calculate the kernel functions [28, 29]. The reference configuration-based kernel function is called Lagrangian kernel, and the corresponding SPH formulation is termed Total Lagrangian SPH (TLSPH). The TLSPH is entirely free of tensile instability [30, 27, 31]; therefore, it is an ideal option for metal processing modelling.

The present work aims to simulate the metal pressing/cutting processes in the conventional/ Eulerian SPH (ESPH) and TLSPH frameworks and provide a comparison between these two methods. In our study, the metal materials are modelled using the Johnson-Cook model. The processing tools are considered as rigid bodies in contact with the specimens. In the ESPH and TLSPH, the tool-metal contact forces are modelled using the standard SPH interaction and the particle-particle contact algorithm, respectively. The presented methods are first validated using a Taylor impact case. Then metal pressing and cutting problems are modelled. The results from the two formulations are compared. Numerical results from the literature are also considered as a reference to evaluate the simulations.

The remaining part of the paper is organised as follows. The formulations of the ESPH and TLSPH are introduced in section 2. The validation of the two SPH frameworks and the numerical simulations of the metal pressing and cutting are presented in section 3. Some conclusions are drawn in section 4.

2 Computational Frameworks

In this section, the ESPH and TLSPH are introduced. In SPH, the computational domain is represented by Lagrangian particles [32, 33, 34, 35]. The particles interact with each other through a kernel function. In ESPH, the kernel function is calculated based on the current configuration, which is updated in every computational step. On the other hand, the Lagrangian kernel, which is based on the undeformed reference configuration, is used in the TLSPH.

The governing equations in the current [35] and reference [36] configurations are

Mass conservation:

[TABLE]

[TABLE]

Momentum conservation:

[TABLE]

[TABLE]

where $\rho$ is the density, $\bm{v}$ is the velocity vector in the current configuration $\bm{x}$. The subscript [math] indicates the values are computed at the reference configuration $\bm{X}$. The current $\bm{x}$ and reference $\bm{X}$ descriptions are related through a mapping $\bm{\phi}$ as $\bm{x}=\bm{\phi}\left(\bm{X},t\right)$. $J$ is the Jacobian computed based on the deformation gradient matrix $\bm{F}$ as $J=\mathrm{det}(\bm{F})$. $\bm{\sigma}$ is the Cauchy stress tensor, $\bm{P}$ the first Piola Kirchhoff stress tensor. $\nabla$ and $\otimes$ denote the divergence/gradient operator, and the outer product between two vectors, respectively.

The deformation gradient $\bm{F}$ is defined as

[TABLE]

where $\bm{u}$ is the displacement calculated as $\bm{u}=\bm{x}-\bm{X}$. $\bm{I}$ represents the identity tensor. Similarly, the rate of deformation gradient is

[TABLE]

where $\bm{v}$ is the velocity.

2.1 Discretised forms in Eulerian SPH (ESPH)

In ESPH, any field variable $f(\bm{x}_{i})$ and its derivative $\nabla f(\bm{x}_{i})$ are approximated based on the current configuration as follows

[TABLE]

[TABLE]

where $f(\bm{x}_{j})$ is the field variable value at the $j$-th particle, $\bm{x}_{ij}=\bm{x}_{i}-\bm{x}_{j}$ is the vector from particle $j$ to particle $i$. $W(\bm{x}_{ij})$ is the kernel function defined in the current configuration. $m_{j}/\rho_{j}$ represents the volume of the $j$-th particle in the current configuration. The following cubic B-spline function is used as kernel function in this work.

[TABLE]

where $\alpha_{d}$ is the normalisation coefficient taken as $10/(7\pi h^{2})$ in 2D and $\alpha_{d}=1/(\pi h^{3})$ in 3D, $h$ is the smoothing length and $q=|\bm{x}_{i}-\bm{x}_{j}|/h$ is the normalised distance associated with the particle pair.

The governing Eqs. (1) and (3) in the current configuration can be discretised [37] as,

[TABLE]

[TABLE]

where $\nabla_{i}W_{ij}$ is the kernel gradient evaluated at the current position of particle $i$. $\pi_{ij}$ is the artificial viscosity used to stabilise the computation in the presence of shock in field variables. The Monaghan type artificial viscosity [38] with controlling parameters $\beta_{1}$, $\beta_{2}$ is used.

[TABLE]

where

[TABLE]

with $\bm{v}_{ij}=\bm{v}_{i}-\bm{v}_{j}$ the relative velocity, $\epsilon=0.01$ the coefficient used to prevent singularity when two particles are too close. $C=\sqrt{E/\rho}$ is the sound speed in the medium with $E$ Young’s modulus. A bar above a quantity indicates that it is the averaged value over the two particles.

In the ESPH, the tensile instability is common in the simulations of solids. The artificial pressure [25, 26] is used to treat the tensile instability. The form of artificial pressure $p^{a}_{ij}$ writes

[TABLE]

where $p=-(\sigma_{xx}+\sigma_{yy}+\sigma_{zz})/3$ is the pressure. $\gamma$ is a case-dependent tuning parameter in the magnitude of 0.1. It is zero if $p_{i}>0$ and $p_{j}>0$. $n$ is a constant taken as $4$, $\Delta p$ is the initial particle spacing.

2.2 Discretised forms in Total Lagrangian SPH (TLSPH)

In the TLSPH description, a function $f(\bm{X}_{i})$ and its derivative $\nabla f(\bm{X}_{i})$ are approximated in the reference configuration as

[TABLE]

[TABLE]

where $\nabla_{i}W(\bm{X}_{ij})$ is the kernel gradient, $\bm{X}_{ij}=\bm{X}_{i}-\bm{X}_{j}$ and $m_{j}/\rho_{0j}$ represents the distance vector and the volume of the $j$-th particle in the reference configuration.

In the TLSPH formulation, the mass conservation equation is trivially satisfied; therefore, it is not needed to be solved numerically. The discretised form of the momentum conservation Eq. (4) writes

[TABLE]

where $\bm{P}=J\bm{F}^{-1}\bm{\sigma}$ is the first Piola Korchhoff stress. $\bm{\Pi}_{ij}=J\bm{F}^{-1}\pi_{ij}$ is the artificial viscosity. The TLSPH is free from any tensile instability due to the use of reference configuration in the computation. Therefore, there is no need of any numerical stabilization term.

Additionally, the discretised forms of the deformation gradient and its rate are

[TABLE]

[TABLE]

2.3 Constitutive equations

The Cauchy stress tensor $\bm{\sigma}$ can be expressed in terms of deviatoric stress $\bm{s}$ and hydrostatic pressure $p$ as $\bm{\sigma}=\bm{s}-p\bm{I}$. The deviatoric stress rate is calculated from the material frame invariant Jaumann stress rate as

[TABLE]

where $\mu$ is the shear modulus, $\dot{\bm{\epsilon}}$ is the strain rate tensor and $\bm{\omega}$ is the spin tensor. $\dot{\varepsilon}=\bm{\dot{}}{\epsilon}_{xx}+\dot{\epsilon}_{yy}+\dot{\epsilon}_{zz}$ is the trace of the strain rate tensor $\dot{\bm{\epsilon}}$. The strain rate $\dot{\bm{\epsilon}}$ and spin tensors $\bm{\omega}$ are calculated from the velocity gradient tensor $\bm{l}$ as follows

[TABLE]

[TABLE]

In ESPH, the velocity gradient tensor $\bm{l}$ is obtained as

[TABLE]

whereas in TLSPH, the velocity gradient tensor is computed as $\bm{l}=\bm{\dot{F}F^{-1}}$.

In this work, the hydrostatic pressure $p$ is obtained from a linear equation of state

[TABLE]

where bulk modulus $K=E/[3(1-2\nu)]$ with $\nu$ the poison’s ratio.

2.3.1 Viscoplastic material model

Metals are modelled using the Johnson-Cook [9] (J-C) model, which can consider salient metal behaviours such as plastic hardening, rate dependency, and thermal softening. In the J-C model, the von Mises yield stress $\sigma_{y}$ is calculated as

[TABLE]

where $A$ is the initial yield stress of the material; $B$, $n$, and $C$ are the hardening parameters. $\epsilon_{pl}$ is the accumulated effective plastic strain. $\dot{\epsilon}_{pl}^{\ast}$ denotes the dimensionless plastic strain rate defined as

[TABLE]

where $\dot{\epsilon}_{pl}$ and $\dot{{\epsilon}}_{0}$ are the current effective plastic strain rate and the reference strain rate. The homologous temperature $T^{\ast}$ is defined as

[TABLE]

where $T$, $T_{0}$ and $T_{m}$ are the current temperature, the room temperature, and the melting temperature of the metal, respectively. It is assumed that the temperature increases due to some fraction of the plastic work being converted into heat.

[TABLE]

where $W_{p}$, $C_{p}$ and $\chi$ are the accumulated plastic work density, the specific heat density at constant pressure, and the Taylor-Quinney empirical constant. $\chi=0.9$ denotes the percentage of plastic work under the adiabatic condition that is dissipated as heat.

The von Mises yield criterion $y_{f}=\sqrt{J_{2}}-\sigma_{y}/\sqrt{3}$ is used to account for metal plasticity, where $\sigma_{y}$ is the yield stress obtained using Eq. (25), and $J_{2}=\bm{s}:\bm{s}/2$ is the second invariant of the deviatoric stress tensor. The Wilkins criterion $\bm{s}_{n}=c_{f}\bm{s}$ is used for return mapping when the trial elastic stress state is beyond the yield surface, with $c_{f}=\min\left(\sigma_{y}/\sqrt{3J_{2}},1\right)$ and $\bm{s}_{n}$ the corrected deviatoric stress tensor. Finally, the following equations are used to compute the increment of plastic strain, the increment of effective plastic strain, and the accumulated plastic work density

[TABLE]

[TABLE]

[TABLE]

2.4 Contact Force

The interaction between different bodies in ESPH can be modelled naturally with the conventional SPH particle interaction. However, the contact force cannot be modelled using the SPH particle interaction in TLSPH due to the use of the reference configuration. Consequently, an explicit algorithm for contact force is essential for the TLSPH for modelling multibody interaction. The pin-ball contact [39] is applied in this work. Here, the contact algorithm is discussed briefly. For more details, the readers may refer to [40].

In the pin-ball contact, the particles are assumed to have a predefined finite influence zone, as shown in Figure 1. The influence zone is represented by a sphere and a circle in 3D and 2D, respectively. The radius of the circle or the sphere is $kh$, where $h$ and $k$ are the smoothing length and a constant factor, respectively. For one particle, the contact force is computed if another particle from a different body comes within its influence radius $kh$.

The contact force is a function of material parameters, penetration depth, and rate of penetration. First, the magnitude of penetration $p_{d}$ of the particle pair is computed as

[TABLE]

where $R_{i}$ and $R_{j}$ are the influence radius of $i$-th and $j$-th particle, respectively. $p_{d}>0$ indicates material penetration, which results in the contact force [41]

[TABLE]

where

[TABLE]

where $\Delta t$, $\dot{p_{d}}=|\bm{v}_{i}-\bm{v}_{j}|$, and $K_{p}$ are the time step, the rate of penetration, and the scaling factor generally chosen through numerical experiment. The momentum Eq. (17) in the TLSPH is modified by including the contact force as

[TABLE]

Note that the contact force is only employed in simulations using TLSPH.

3 Numerical simulation

In this section, the ESPH and TLSPH methods are first validated using a Taylor impact test [42]. The numerical simulations of metal pressing and cutting are then performed. The predictor-corrector time integration is used to integrate the discretised governing equations. The Courant–Friedrichs–Lewy (CFL) condition is applied to determine the stable time step. In all the simulations, the tool-metal contact is assumed to be frictionless. In all simulations with the ESPH, the artificial pressure coefficient is taken as 0.3. In TLSPH, the reference configuration is updated to the current configuration when the deviation between the current and reference frames is considerable [43]. The particle volumes and the material coordinates are updated as $V\leftarrow JV$ and $\mathbf{X}\leftarrow\mathbf{x}$ respectively. The updates preserve the stress tensors. The updates of the reference frames are essential for problems involving a large plastic flow of materials to capture the changes in particle interaction [43].

3.1 Taylor impact test

A low strength 4340 steel cylinder of 37.97 mm length and 7.595 mm diameter strikes a rigid wall with an initial velocity of 181.0 m/s [42] as shown in Figure 2. The material and computational parameters are given in Table 1 and 2, respectively.

The numerical results of the final length of the cylinder after impact and the final diameter at the mushrooming end are compared with the experimental results [44] and the numerical results obtained with LS-DYNA [42]. The final diameters of the mushrooming end are found to be 9.80 mm and 9.83 mm with the ESPH and TLSPH, respectively. The experimental value is 9.5 mm, and the computed value from LS-DYNA is 10.6 mm. The predictions of the final diameter of the cylinder at the mushrooming end are very close to the experimental observation. The current predictions deviate from the experimental observation by 3.16$\%$ and 3.47$\%$ with the ESPH and TLSPH, respectively. However, the deviation of the computed diameter by LS-DYNA is 11.58$\%$ to the experimental result. The final lengths of the deformed cylinder are also calculated and are found to be 34.19 mm and 34.18 mm with the ESPH and TLSPH, respectively. The experimental and computed values with LS-DYNA are found to be 34.6 mm and 34.4 mm. The present predictions are close to experimental and numerical observation. The differences between our results and the experimental results are 1.18% and 1.21% for the ESPH and TLSPH, respectively.

The change in the length of the cylinder is also compared with the LS-DYNA results in Figure 3. The present predictions are in good agreement with the LS-DYNA results. The contours of the accumulated effective plastic strain and temperature are shown in Figure 4 and 5 respectively. It can be readily observed that the distributions of the plastic strain and temperature are captured quite well and are similar to the LS-DYNA results. However, the contour plot from the TLSPH provides a better agreement in the distribution of plastic strain and temperature with the LS-DYNA. Furthermore, the particle distributions in the TLSPH results are more regular, whereas, in the results obtained from the ESPH, the irregular particle distribution caused by tensile instability can be observed even with the artificial pressure.

3.2 Metal pressing

In this section, the pressing of a metal block of AISI 4340 steel by a rigid tool is modelled, as shown in Figure 6. The tool is pressed from the free right side, while the left side of the metal block is fixed. The workpiece is of 1 mm $\times$ 0.5 mm $\times$ 0.1 mm, and the pressing tool is of 0.1 mm $\times$ 0.7 mm $\times$ 0.1 mm. The pressing tool moves with a constant velocity of 200 m/s. A three-dimensional simulation is performed, using 184,624 particles for the workpiece and 25,199 particles for the pressing tool. The response of the workpiece is modelled by the Johnson-Cook material model and the pressing tool is modelled as a rigid body. The constitutive parameters are listed in Table 3. The other relevant computational parameters are given in Table 4.

The numerical results from the ESPH and TLSPH are compared with those from the Material Point Method (MPM) and Finite Element Method (FEM) simulations [11]. The final shape and distribution of equivalent plastic strain are shown in Figure 7, with results from all the four approaches. It is observed that the results from the four simulations have similar final shape and distribution of equivalent plastic strain, indicating the both the ESPH and TLSPH can correctly capture the large deformation and strain localisation in the pressing. Compared with MPM, the final shapes obtained from SPH simulations are closer to the FEM results, as the MPM results show large deformation along the secondary group of shear bands, leading to somewhat different final shape. Also, the equivalent plastic strains in the MPM results are diffusive and lack distinct shear bands. On the other hand, the results from the two SPH methods show clear shear bands, like those from FEM. However, the ESPH results show slightly more diffused distribution of equivalent plastic strain. This is due to the application of the updated particle position in the computation.

The accumulated equivalent plastic strain across the metal block is plotted in Figure 8. The two methods give very similar results except for the most left part. The figure shows a consistent trend that the right part of the metal block is subjected to more plastic deformation. The increase in temperature is shown in Figure 9. It is observed that the area and magnitude of temperature change are consistent with the equivalent plastic strain. This increase in temperature can further reduce the yield stress in shear bands, making the sling of material in these regions more easily.

To analyse the effect of initial particle spacing, four simulations with different resolution ($\Delta p=0.0167$ mm, 0.0133 mm, 0.01 mm, 0.0067 mm) are performed. In the four simulations, the ratio between the smoothing length and initial particle spacing is kept constant as 1.5. The results are compared with the results obtained from the particle spacing 0.0067 mm in Figure 10. It can be observed that the plastic strain bands become more distinct as the particle resolution is refined. Moreover, results from TLSPH show more regular particle distribution than those from ESPH.

3.3 Metal cutting

Formation of shear bands and accumulation of plastic shear strain play an essential role in metal cutting processes. As the cut materials often undergo considerable plastic deformation, high strain localisation, and material separation, mesh-based methods such as FEM have difficulties because of mesh distortion and entanglement. The discontinuities in metal cutting processes are also problematic for mesh-based methods. As SPH is meshless, it is free from mesh deficiencies and naturally capable of simulating large deformation. Moreover, material separations can be simulated naturally in SPH.

The setup for the orthogonal metal cutting process is shown in Figure 11. The workpiece (3 mm $\times$ 1 mm) is made of AISI 4340 steel. The JC material model is used to simulate the material behaviour. The JC parameters are shown in Table 3. In the metal cutting process, a rigid cutting tool (0.5 mm $\times$ 1 mm) moves with a constant velocity of 50 m/s [11]. In the present work, it is further assumed that there is no friction between the workpiece and the cutting tool. Other computational parameters are listed in Table 5. The pinball contact model is used for the computation of the contact forces in the TLSPH.

3.3.1 Results with feeding length of 0.3 mm

The plastic strain and temperature distribution are important output variables in the metal cutting process. The results of accumulated equivalent plastic strain from the ESPH and TLSPH are shown in Figure 12 and 13, respectively. It can be observed that both approaches produce continuous chip morphologies in the metal cutting process. As the tool cuts further into the workpiece, plastic strain accumulates in front of the toe of the tool, and then the strain localisation propagates towards the surface in an inclination approximately equal to 45∘. This numerical observation is inconsistent with the plastic theory. Once a shear band is fully developed, the cut material above the newly-formed shear band slides along it, leading to a new chip. With this mechanism, results from both the ESPH and TLSPH show a number of chips separated by shear bands. The increase of temperature is also shown in Figure 14 and 15. The increase in temperature is consistent with the distribution of equivalent plastic strain.

The formation of chips in the metal cutting process is compared with results from MPM and FEM in Figure 16. Numerical results of chip formation and distribution of shear bands from ESPH and TLSPH are in good agreement with the MPM and FEM results. However, the formation of chips and the deformation of the cut stripe captured using TLSPH are closer to the MPM and FEM results. The reason may lie in the formulation of TLSPH in which the number of instabilities is less than the ESPH. Moreover, the bending of the cut stripe in the TLSPH is more than the ESPH simulation, which is closer to the results from MPM and FEM. Furthermore, the shear bands are more distinct in results from the TLSPH. More importantly, the simulations using TLSPH does not need to tune the coefficient of artificial pressure like in ESPH.

3.3.2 Effect of feeding length on shear band

The shear bands and the accumulation of plastic strains play an important role in the metal cutting process. The increase of the feeding length increases the cutting force in the simulation. This, in turn, increases the amount of plastic deformation and heat dissipation in the workpiece. The depth of feeding is increased to 0.5 mm from 0.3 mm used in the previous metal cutting simulation. This influences the metal cutting process and changes the chip morphology. The deformation of the workpiece with 0.5 mm feeding length are shown in Figure 17 and 18 by the ESPH and TLSPH respectively. It is observed that the increase in the cutting forces leads to the increment of the size of the chips, i.e. the number of shear bands decreases. This is due to the increase of the localisation of material plastic strains and leads to larger chips between shear bands as shown in Figure 17 and 18. The computed results are compared with the results from MPM [11]. It can be observed in Figure 19 that the present results are in good agreement with the results from [11].

4 Conclusion

The numerical predictions of metal machining process such as metal pressing and cutting are a complex problem. These processes involve localisation of plastic strain, shear band formation, material separation and temperature increment leading to material softening. The mesh-based methods such as FEM are not a suitable candidate due to mesh distortion, entanglement etc. In this paper, the ESPH and TLSPH frameworks are used to model the metal machining process. Due to the particle nature, SPH is inherently capable of capturing the large plastic deformation of materials. However, the ESPH suffers from tensile instability leading to numerical fracture of material. The TLSPH is free from this due to the computation of kernel functions in the reference configuration. The ESPH and TLSPH frameworks are first verified using the Taylor impact test. The computed final length, diameter and the change in length over time are in good agreement with the experimental and numerical observations.

The ESPH and TLSPH are used to simulate the metal pressing and cutting process. The accumulated effective plastic strains and the distribution of temperature are compared with the numerical result available in the literature. The presented results are consistent with the observations from the literature. It is further observed that the TLSPH provides a more realistic chip morphology, deformation of the workpiece and generation of shear bands. The change in feeding length influences the material deformation. With the increase in the feeding length, the cutting force increases. This leads to a reduction in the number of shear bands and an increase in the width between two adjacent shear bands. In the current work, the effect of friction in the metal cutting process is neglected. However, the influence of frictional force and material failure may have impacts on metal machine processing. In future, more advanced contact algorithm will be considered for the frictional contact.

Acknowledgement

This work has received funding from the European Unions Horizon 2020 research and innovations program under grant agreement No. 778627 and from FFG the Austrian Funding Agency under the project No. 865963.