Applied Mathematical Modelling 36 (2012) 3836–3855
Contents lists available at SciVerse ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Modelling of metal forging using SPH Paul W. Cleary a, Mahesh Prakash a, Raj Das b,⇑, Joseph Ha a a b
CSIRO Mathematics, Informatics and Statistics Normanby Road, Clayton, Victoria 3168, Australia Department of Mechanical Engineering, University of Auckland, Auckland 1010, New Zealand
a r t i c l e
i n f o
Article history: Received 1 July 2010 Received in revised form 14 October 2011 Accepted 2 November 2011 Available online 10 November 2011 Keywords: SPH Elastoplastic deformation Forging Metal forming Defect Plasticity
a b s t r a c t In solid metal forming processes, such as forging, large distortions in the material present challenging problems for numerical simulation using grid based methods. Computations invariably fail after some level of mesh distortion is reached unless suitable re-meshing is implemented to cope with the mesh distortion arising from the material deformation. The issue of mesh distortion and the subsequent re-meshing are topics of much research for grid based methods. These problems can be overcome by using a mesh-less numerical framework. In this paper, the application of a mesh-less method called Smoothed Particle Hydrodynamics (SPH) for modelling three-dimensional complex forging processes is demonstrated. It is shown that SPH is a useful simulation method for obtaining insights into the material deformation and flow pattern during forging of realistic industrial components. The effect of process parameters and material properties on the quality of the forged component is evaluated via SPH simulations. This includes the determination of forging force required for adequate die filling which is an important criterion for die designs. Material hardening, controlled by the degree of heat treatment, is found to have a profound effect on the material deformation pattern and the final product. Forging defects such as incomplete die filling, asymmetry in forged components, flashing and lap formation are shown to be predicted by SPH. SPH can thus potentially be used both for assessment of the quality of forged products and evaluation of prototype forging system designs. Ó 2011 Elsevier Inc. All rights reserved.
1. Introduction Forging is one of the most economical and efficient methods for fabricating complex metal components. In this process, metal is subjected to high compressive force, usually via hydraulic press or hammer blows. The metal is plastically deformed and assumes the shape of the dies. The resulting components have refined grain structure and improved mechanical properties, such as strength, toughness, and ductility [1–3]. In manufacturing industries, the choice of suitable forging process parameters and the design of tools and dies are still largely based on trial and error. Numerical simulations can be a useful tool in obtaining understanding and knowledge about the material deformation process in forging. The main reasons for carrying out forging simulations include: Optimisation of number of operations to form the component. Optimisation of the size and shape of the blank. Optimisation of material flow.
⇑ Corresponding author. E-mail address:
[email protected] (R. Das). 0307-904X/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2011.11.019
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3837
Prediction of forging defects such as fold, wrinkle and under fill. The understanding of material deformation and flow pattern during forging provides key inputs to the prototype design of industrial forging systems and guides the choice of process parameters. For example, one important criterion of forging process design is to minimise the generation of defects in the final forged product. Forging defects and flaws in mechanical components reduce their structural integrity and can lead to premature failure of structures, thereby leading to downtime maintenance and repair cost [4,5]. The ability to predict the generation of possible defect types and their locations via numerical modelling will enable assessment of the effect of die design and/or processing conditions on the quality of forged products. This can potentially help in appropriate equipment design and selection of process parameters to reduce forging defects. The numerical modelling work undertaken so far has not been adequate to predict forging defects. The key requirements of modelling forging are handling deforming geometry and implementing history-dependent constitutive behaviour. The difficulties experienced by mesh-based numerical methods in addressing these issues and how Smoothed Particle Hydrodynamics can be effectively used in such situations are discussed below. Forging processes are characterised by rapid deformations under large forces that generate high intensity plastic deformations. When the conventional finite element method (FEM) is used for modelling large strain problems, high deformation of the geometry inevitably leads to mesh distortion. The topology and shape of the elements have profound effect on the quality and accuracy of FE solutions [6]. Mesh distortion leads to spurious numerical errors and mesh sensitivity of the solutions. Severely distorted elements may even lead to a singular stiffness matrix (non-positive Jacobian) and consequent analysis failure [7]. This is further compounded in forging simulations, where very high displacements may produce mesh inter-penetration and entanglement, resulting in the loss of connectivity and topological integrity of the elements. In such cases, re-meshing part or all of the geometry is the only alternative solution. As a result, dynamic and adaptive re-meshing techniques have received a great deal of attention [8–11]. Maillard et al. [8] developed an automated re-meshing technique for two-dimensional shapes and simulated cold forming using an elastoplastic finite element method. In the late 1990s, a number of adaptive three-dimensional re-meshing algorithms emerged for FE simulations of metal forming [9–11]. Intermediate re-meshing is computationally expensive and is difficult to implement for complex geometries, particularly for those undergoing rapid change in shape during forging. Re-meshing in FEM involves redistribution of existing nodes and/or reconstruction of new elements. This can lead to violation of consistency conditions, especially for tool-workpiece contact implementation in material forging [9]. Forging processes inherently rely on flow of plastically deformed solids. Implementation of many elastoplastic constitutive models requires information on the prior material deformation (stress–strain) histories. Accurate calculation of the time variation of state variables for a given volume of material are difficult for FEM and become stronger when used with local or global re-meshing of the geometry. Additional enhancements in the re-meshing algorithms are required to capture the time history of relevant field variables. A multi-mesh algorithm was used in [12] that re-meshes the high deformation regions with a fine mesh and transfers data between previous and newly re-meshed zones. Recently, Chan Chin [13] proposed the use of particle trackers in FEM for history tracking and flow behaviour monitoring in forging simulations. However, progressive data mapping between successive meshes introduces inaccurate history tracking for a given volume of material due to interpolation (mapping) of data from one mesh to another, thus making it difficult to include path-dependent nonlinear material behaviour. It can also affect the accuracy of time integration in each stage of re-meshing. The errors introduced in the solutions due to data mapping between meshes have not been reported in the literature. Many advantages of mesh-less methods have led to their applications in material forming simulations. Three categories of mesh-less methods are commonly used: coupled mesh-less-mesh-based (e.g. FEM-EFG), partial mesh-less methods that use a background mesh for interpolating and/or storing state variables (e.g. EFG, RKPM and PIC methods), and fully mesh-less methods that solely use particles for interpolation, but may use a background grid for neighbour search only (e.g. SPH). Coupling of mesh-based and mesh-less methods for numerical simulation of metal forming is becoming popular to avoid some of the difficulties associated with a mesh-based geometry representation [14–16]. Liu et al. [16] proposed a combined Finite Element and Element Free Galerkin (FE–EFG) method for bulk metal forming simulations. Adaptive conversion between FE and EFG domains was based on local strain histories and error estimates. It was found that this procedure could only be effectively applied for forming of simple shapes. Instabilities at the coupling region were observed for relatively complex three-dimensional geometries. In general, Galerkin mesh-free methods are prone to spatial instabilities [17] and usually require additional stabilization techniques. For example, Chen et al. [17] used a strain smoothing technique to reduce this. So it appears that coupled methods are not easy to implement, and may even have poor accuracy compared to that of full FE simulations (in conjunction with remeshing). Full mesh-less methods, such as Reproducing Kernel Particle Method (RKPM) and Smoothed Particle Hydrodynamics (SPH), offer many advantages for modelling continuous high deformation processes, such as forging. RKPM has been used for modelling of material forming processes [18]. Examples of three-dimensional bulk material forming simulations using RKPM can be found in [19,20]. Specific applications of RKPM also include sheet metal (super-plastic tension) forming [21,22]. Park [23,24] used a rigid-plastic RKPM for metal forming simulations. This approach uses a constitutive flow formulation based on the assumption that elastic effects are insignificant during metal forming operations. However, for certain categories of metal forming, such as forging, elastic spring back plays an important role in determining the quality of the final product. Guang-ru et al. [25] used RKPM for thermo-mechanical analysis of forging processes. Since RKPM uses a background mesh (cells) for interpolating some of the state variables, it also suffers from unsatisfactory performance (although to
3838
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
a lesser extent) for reasons similar to those encountered with FEM. RKPM often requires adaptive generation of background cells at each time step [26], which is also difficult to implement for complex geometries. Finally, RKPM solutions are found to be sensitive to the nature of background cells employed. The type of the background cells (triangular or quadrilateral) influences the prediction of flow pattern, forging load, and strain distribution which is strongly undesirable [27]. A fully mesh-less method, such as SPH, offers many advantages for modelling continuous high deformation processes, such as forging. SPH is a particle based numerical method which is used to obtain solutions to systems of partial differential equations. The problem geometry is discretised into ‘particles’ that represent specific material volumes. The method is based on the use of local interpolation from the surrounding discrete particles to construct continuous field approximations that can be used in the discretisation of field equations. For details of SPH fundamentals, see [28–30]. SPH is an established numerical method and has been traditionally used for modelling fluid flows [28,31,32,30]. In recent years, there has been a growing interest in applying SPH to a wide variety of solid mechanics problems [33–38]. SPH does not use any underlying grid or mesh structure to represent the problem geometry and so eliminates time-consuming effort of mesh generation and refinement. It also avoids the need to address the difficulties associated with traditional mesh-based methods in maintaining the integrity and quality of the mesh under large deformation thereby reducing associated errors. The grid-free nature of SPH makes this method ideally suited to modelling material forging processes that involve large deformations and discontinuities. Moreover, SPH can handle complex free surface deformation which enables prediction of surface contacts during forging and naturally predicts self-contact of surfaces after folding. In SPH, a given particle always represents a specific mass of material and carries that information with it. This is a critical attribute of Lagrangian methods. This means that information on the precise state of each piece of metal is recorded at all times and the history of each piece of metal is built into the particle data. So SPH can predict the thermal and mechanical history of the workpiece without the need for data mapping between evolving (discretised) geometric configurations. This provides significant capability to track properties, such as cumulative plastic strain, damage, metal composition, trapped gas, metallic phase, microstructure, grain flow, and surface oxide. Some or all of these properties can then be used to feed back into the flow dynamics using suitable constitutive models. This also facilitates easy incorporation of rate-dependent deformation behaviour, and so the effect of process parameters, such as axial feed rate, on the thermo-mechanical responses can be included. Another advantage of SPH is that it can easily handle coupled physical processes, such as elastoplastic flow and heat transfer, in a unified analysis framework. The modelling of all components of a multi-physics application using the same numerical framework means that the physical laws, interfacial boundary conditions, and associated evolution of the material properties and parameters can be seamlessly incorporated into an integrated analysis system. SPH can also easily incorporate geometric nonlinearity by using a large strain formulation and material nonlinearity due to its Lagrangian formulation. All these analysis features are essential requirements of modelling material forming processes. Despite many attractive features, application of SPH for metal forging has been quite limited. Bonet and Kulasegaram [39] used a variant of SPH, known as Corrected SPH (CSPH), and demonstrated its potential in metal forming modelling via simple 2D forging examples. Cleary et al. [38] also demonstrated the strengths of SPH for modelling of forging and extrusion via 2D simulations. In this paper, we demonstrate the potential of SPH for simulation of three-dimensional forging processes involving realistic complex shaped industrial components and multiple moving and stationary dies. First, the formulation of the SPH method for elastoplastic deformation is presented. The ability of the SPH method to model severe plastic deformation during forging is then demonstrated via two dimensional examples. It is shown that SPH is a useful simulation method for explicitly predicting several types of forging defects (e.g. incomplete die filling, flashing, and lap formation) that commonly occur in industrial forged components. SPH is then applied to model forging processes of three-dimensional industrial components. For these examples, the material flow is modelled to investigate the plastic deformation pattern and to establish the capabilities of SPH as a useful simulation tool for industrial scale material forming process modelling. 2. Sph formulation for elastodynamics SPH is a particle based numerical method, used to obtain solutions to systems of partial differential equations. The problem geometry is discretised into ‘particles’ that represent specific material volumes. For more comprehensive details, one can refer to Monaghan [29] and Cleary et al. [30]. In recent years, SPH has been developed for and applied to modelling solid deformation problems in a range of applications [33–38]. The interpolated value of a function A at any position r can be calculated from the discrete values Ab on the surrounding particles using SPH interpolation by
AðrÞ ¼
X b
mb
Ab
qb
W ðr rb ; hÞ;
ð1Þ
The sum is over all particles b within a radius 2h of r. Here W(r, h) is a C2 spline based interpolation or smoothing kernel with radius 2h, that approximates the shape of a Gaussian function but has compact support and mb and qb are the mass and density of particle b. The gradient of the function A is given by differentiating the interpolation Eq. (1) to give:
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
rAðrÞ ¼
X
mb
b
Ab
qb
rWðr rb ; hÞ:
3839
ð2Þ
Next we present the SPH equations for transient elastoplastic problems. 2.1. Continuity equation The continuity equation for an elastoplastic solid is given as:
dq ¼ q r v dt
ð3Þ
where q is the density and v is the velocity. From Monaghan [28], the SPH form of the continuity equation is:
dqa X mb ðv a v b Þ rW ab : ¼ dt b
ð4Þ
We denote the position vector from particle b to particle a by rab ¼ ra rb and let W ab ¼ Wðrab ; hÞ be the interpolation kernel with smoothing length h evaluated for the distance jrab j. This form of the continuity equation is Galilean invariant (since the positions and velocities appear only as differences), has good numerical conservation properties and is not affected by free surfaces or density discontinuities. 2.2. Momentum equation The momentum equation used for the elastoplastic deformation of the solids is:
dv 1 ¼ r r þ g; dt q
ð5Þ
where r is the stress tensor, and g is the gravity. The stress tensor can be written as:
r ¼ PI þ S;
ð6Þ
where P is the pressure and S is the deviatoric stress tensor. We use a neo-Hookian elasticity model with a bulk modulus K and a shear modulus G. The SPH form of Eq. (5) expressed in terms of overall stress is:
dv a X ra rb mb 2 þ 2 þ Pab I ra W ab þ ga ; ¼ dt qa qb b
ð7Þ
where ra and rb are the stress tensors of particles a and b, and ga is the body force at particle a. Pab is an artificial viscous stress term with shear and bulk viscosity components and is given by Monaghan [28] as:
( ac
2 ab gab þbgab
Pab ¼
q ab
ifðv a v b Þ ðra rb Þ < 0;
0
ifðv a v b Þ ðra rb Þ P 0;
ð8Þ
where
gab ¼
hðv a v b Þ ðra rb Þ jra rb j2
ab ¼ and q
ðqa þ qb Þ : 2
Here c0 is the sound speed for the particles, a and b are shear and bulk viscosity coefficients. We use a = 1.0 and b = 2.0. The evolution equation for the deviatoric stress S is given by Gray et al. [35] as:
ij dS 1 ¼ 2G e_ ij dij e_ ij þ Sik Xjk þ Xik Skj ; 3 dt
ð9Þ
where the strain rates are:
e_ ij ¼
1 ov i ov j ; þ 2 oxj oxi
ð10Þ
1 ov i ov j 2 oxj oxi
ð11Þ
and
Xij ¼
is the Jaumann rotation tensor that accounts for the large rotational effect.
3840
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
2.3. Equation of state and time stepping The SPH method is compressible with the pressure calculated from the particle density using an equation of state:
P ¼ c2 ðq q0 Þ;
ð12Þ
where P and q are the current pressure and density, q0 is the reference density and c is the speed of sound given by
sffiffiffiffiffiffi K : c¼
ð13Þ
q0
Unlike with SPH for incompressible fluids where an artificially low speed of sound is used to accelerate the computations, here we use the correct material properties and speed of sound so as to correctly resolve stress wave propagation. The Improved Euler scheme given in [29] is used for explicit integration of the SPH equations. For elastic solids, this time step Dts is limited by the Courant condition, and is given as:
Dts ¼ minf0:5h=ca g: a
ð14Þ
where ca is the speed of sound for the material of particle a. 2.4. Plasticity model: Von Mises plasticity with linear isotropic hardening The radial return plasticity model by Wilkins [40] is used to model elastoplastic deformation. A trial deviatoric stress Str is calculated assuming an initial elastic response. The increment of plastic strain is:
Dep ¼ where
rv m ry 3G þ H
;
ð15Þ
rv m is the von Mises stress and ry is the current yield stress. The plastic strain is then incremented:
ep ðt þ DtÞ ¼ ep ðtÞ þ Dep :
ð16Þ
The yield stress increment Dry at each time step is calculated as:
Dry ¼ HDep ;
ð17Þ
where H is the hardening modulus. The deviatoric stress S at the end of a time step (t + Dt) is given by
Sðt þ DtÞ ¼ r s SðtÞ;
ð18Þ
where rs is the radial scale factor given by
rs ¼
ry : rv m
ð19Þ
2.5. SPH solution process The governing partial differential equations for elastoplastic deformations Eqs. (5)–(11), (15)–(19) are converted to ordinary differential equations by approximating the spatial derivatives and functions using the SPH interpolations Eqs. (1) and (2). Any higher order derivative can be approximated by the next lower order derivatives using the finite difference method. Thus, higher order derivatives in the PDEs are ultimately expressed in terms of the first order derivatives, which in turn can be evaluated via the interpolation relation in Eq. (2). So using SPH interpolants one is able to convert parabolic partial differential equations into a coupled set of ordinary differential equations for the motion of the particles and the rates of change of their properties. These ordinary differential equations are integrated using a second order Improved Euler predictor–corrector. Firstly Eqs. (9) and (4) are used to give the current deviatoric stress and density. The density is then used in Eq. (12) to give the pressure. Any limiting of the deviatoric stress by the plasticity model is then evaluated using Eqs. (15)–(19). The final deviatoric stress and the pressure along with the current velocity are then used in Eq. (7) to calculate the acceleration of the particles. This equation is then integrated along with the kinematic equations to give updated particle positions. SPH is robust in that it can operate on a field of disordered particles provided the kernel is differentiable. There is no need for a fixed, structured grid as required in finite-volume or finite element methods. 2.6. SPH boundary conditions Modelling of forging requires contact simulation between dies and metal. The dies are modelled as rigid bodies. The contact between the rigid die surface and the metal is implemented using a ‘boundary normal force’, where a rigid surface is discretised into particles, called ‘boundary particles’. These ‘‘boundary particles’’ are assigned properties similar to the
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3841
moving ‘‘metal’’ particles, e.g. mass, position, density, viscosity, etc. To avoid penetration of metal particles through the rigid surfaces of dies, the boundary particles exert forces on metal particles in the normal direction of the rigid die surface. The boundary forces are typically of Lennard–Jones form [29]:
8 s > < Aðl=aÞ fb ðaÞ ¼ Dða bÞ2 > : 0
if 0 < a < l; ð20Þ
if l < a < b; otherwise;
where fb is the boundary force experienced at a perpendicular distance a from the boundary particle, b is the distance that limits the range of boundary force, i.e. beyond which no boundary force is exerted, s is an exponent that determines the steepness of the boundary force and s = 8 was used here, A is a boundary force factor that determines the magnitude of the boundary force, and is generally 0:01c2 . The value of A can be varied to avoid die surface penetration. The factors D sþ2 2 sb and l are given as: D ¼ A 2b and l ¼ sþ2 . 3. Material flow and deformation modelling including defect prediction for simple forging processes First, we demonstrate the application of SPH for simulating forging processes via simple examples. The ability of SPH to predict material flow and plastic deformation pattern is investigated by evaluating the evolution of flow velocity, plastic strain and stress distribution in the metal during the forging operation. The aim of this section is to establish the utilities of the SPH method to estimate appropriate process parameters (forging force), to explicitly predict various defects due to forging, and to assist in the design of workpiece and forging equipment (e.g. die). To avoid complexities in geometry and reduce computation time, two-dimensional simulations are presented. The material properties used in the forging simulations are those of an Aluminium alloy (A6061) and are given in Table 1. Fig. 1 shows a simple process of an impression die forging where two dies are brought together by a constant applied force on the top die which is free to move in response to the balance of applied force and reaction force from the workpiece. The workpiece is placed between the die halves and is compressed. It first undergoes elastic and then plastic deformation until it takes on the shape of the cavity formed by the dies. The dies are both 1 m wide. The upper die is 0.2 m high and the lower die is 0.4 m. The workpiece is 0.56 m wide and 0.4 m high. The upper moving die is located above the workpiece which is held by the stationary lower die. For this case, the applied force on the upper die is 400 kN. In this SPH simulation, 3027 particles with 10 mm spacing are used to model the dies and workpiece. In Fig. 1, the particles are coloured by the horizontal velocity component with blue indicating 2.5 m/s, green indicating 0.0 m/s and red indicating 2.5 m/s. The force from the upper die starts compressing the workpiece at 10 ms. As the top die moves downwards, the region near the two contact points with the lower die begins moving sideways in opposite directions at a velocity of around 1.5 m/s. After around 40 ms the metal is stretched horizontally in all regions that are in contact with the lower die surface. The horizontal velocity varies from 2.5 to 2.5 m/s at this stage. By 50 ms the workpiece has stretched sufficiently sideways to almost reach the outer edges of the die cavity. The bottom of the die free surface is curved, and as a result the metal moves slower along the sides due to the resistance to flow. The forging process is complete at 60 ms when the metal completely fills the die cavity. In Fig. 2, the particles are coloured by equivalent plastic strain with blue indicating 0.0 and red indicating 1.0. The figure shows the deformation process for the workpiece and the locations of high plastic strain. The force from the upper die is experienced by the workpiece at 10 ms. As the top die moves down, plastic deformation of the workpiece occurs in the regions near the two contact points with the lower die. Due to the horizontal stretching motion of the metal, the deformation is mainly sideways forcing solid metal into the side cavities. The central part of the workpiece moves steadily down into the open region of the bottom die with little deformation. After about 40 ms, high deformations are evident in all regions of the workpiece that are in contact with the lower die surfaces. The central section of the workpiece continues to move unimpeded down into the lower die space. Moderate plastic strains of at least 0.5 are now experienced by all the material above the high strain contacts as this material is progressively forced sideways. By 50 ms, the outer edges of the contact regions between the workpiece and the lower die have moved outward almost up to the opening between the dies. All the metal being forced into these side gaps between the dies has now experienced medium to high levels of plastic deformation. The free surface of the workpiece at the bottom in the middle of the die has now become curved as the resistance to flow along the walls slows the outer material. The forging operation is completed soon after 60 ms and the workpiece takes on the shape of the cavity formed by the upper and lower dies. The high plastic strain is
Table 1 Properties of Al alloy A6061 used for forging. Bulk modulus (GPa) Shear modulus (GPa) Initial yield stress (MPa) Hardening modulus (MPa) Density (kg/m3)
70.0 27.0 55.2 1.67 2700.0
3842
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
Time = 0.01 s
Time = 0.02 s
Time = 0.03 s
Time = 0.04 s
Time = 0.05 s
Time = 0.06 s
Fig. 1. Simple forging example using SPH. The particles are coloured by the horizontal velocity component with blue being 2.5 m/s, green 0.0 m/s and red indicating 2.5 m/s. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
concentrated in the side regions of the forged specimen. High strain is important since it affects material properties and mechanical performance of the die and the creation and migration of grain scale defects. The plastic strain is moderate along the bottom of the forging and low across the top and throughout much of the middle regions. In Fig. 3 the particles are coloured by von Mises stress with blue indicating values up to 1.0 106 Pa and red indicating a value of 5.5 107 Pa and above. The workpiece has a yield stress of 5.5 107 Pa. Regions experiencing stresses consistently above this value will undergo plastic deformation. The corners of the workpiece in contact with the lower part of the die have stress values above the yield stress limit at 10 ms. At this stage the central core of the die also experiences a high level of stress. However the stress in this central section becomes lower at 20 ms. High stress levels need to be maintained for a prolonged period for the metal to deform plastically. Stress values above the yield stress limit are maintained in the region around the contact points of the workpiece with the lower half of the die at 30 and 40 ms, resulting in high levels of plastic
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3843
Fig. 2. Simple forging example using SPH. The material is coloured by the level of plastic strain with blue for 0.0, green for 0.5 and red for 1.0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
deformation in these regions. At the same time consistent high stress levels are also experienced by the metal that is forced into the side cavities of the die at 40 and 50 ms. These are the regions experiencing moderate to high plastic strain levels. At 60 ms when the component is completely forged, high stress levels are observed in the workpiece on the sides and almost in the entire lower half. The high stress level in the lower half of the workpiece is however not sustained for a long duration thereby resulting in insignificant plastic deformation in this region.
3.1. Effect of forging force on die filling pattern Next we study the effect of forging force on the filling pattern of the die. The magnitude of the applied force used to move the top die is an important process parameter in controlling the forging process. If the forging force applied to the top die is low, then it is not sufficient to completely close the dies.
3844
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
Time = 0.01 s
T ime = 0.02 s
Time = 0.03 s
T ime = 0.04 s
Time = 0.05 s
T ime = 0.06 s
Fig. 3. Simple forging example using SPH. The material is coloured by the value of Von Misses stress with blue for values up to 1.0 104 and red for 5.5 105 Pa and above. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 4 shows the effect of forging force on the filling pattern of the dies. The final state of the incomplete forging is shown in Fig. 4(a) for a force of 0.3 KN. The final forged component is shown in Fig. 4(b) when a higher forging force of 0.8 KN is used. Incomplete forging due to insufficient force is one of the simplest forging defects to predict for SPH since its automatic free surface predictions require no additional efforts in modelling. In this case, a force of 0.4 KN is found to be sufficient to produce a well forged component. The use of larger forces than 0.4 KN does not change the final forged state, rather it only decreases the time required to forge the component. In Fig. 5, the motion of the top die with time is shown for different applied forces. The die accelerates due to the applied force leading to a decrease in the height of the dies. As the workpiece deforms and becomes elastically loaded, a resisting force builds up and the top die decelerates and gradually comes to rest. The figure shows that the forging operation took 60 ms for an applied force of 0.4 KN. Doubling the applied forging force reduces the forging time to only 36 ms. The cases
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3845
(a)
(b)
Fig. 4. Effect of different applied force: (a) inadequate applied force of 0.3 kN, (b) sufficient forging force of 0.8 kN. The particles are coloured by plastic strain.
Vertical Distance (m)
0.30
kN Force = 0.2 2E-05 kN Force = 0.3 3E-05
0.25
kN Force = 0.4 4E-05 kN Force = 0.8 8E-05
0.20 0.15 0.10 0.05 0.00 0.00
0.02
0.04
0.06
0.08
0.10
0.12
Time (s) Fig. 5. Position of the top die as a function of time for different applied forging forces.
where the vertical position of the dies does not reach 0.025 m indicate that the die has not been able to fully close and that the force applied was inadequate for complete forging. The magnitude of the applied force needed to produce complete filling of the die depends on the bulk modulus of the material (which controls the stiffness) and particularly on the plasticity properties. The majority of the deformation is plastic flow and this is controlled by the yield stress and the hardening modulus. 3.2. Prediction of forging defects One advantage of such simulations for forging is to reduce material waste and under-filled dies by predicting the correct size of the initial workpiece needed to exactly produce the required forging. Fig. 6 shows the results of using the incorrect size for the workpiece. In Fig. 6(a), the workpiece is too small leading to incomplete filling of the die, with long narrow gaps at the bottom of the die left unfilled. Incomplete filling is a key defect that requires the part to be discarded. Fig. 6(b) shows
3846
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
(a)
(b)
Fig. 6. Forging defects due to improper size of the workpiece: (a) incomplete filling due to a workpiece that is too small, and (b) flashing of excess material between the top and bottom dies. Particles are coloured by plastic strain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
the flashing produced by using an overly large workpiece. Flashing occurs when metal is forced, in a thin layer, out between the two dies by the very high compressive pressures developed in the die when the workpiece volume is too large. If the flashing is thin, this can be removed by trimming, but this adds costs and needs to be minimised. If the flashing is thick, then physical removal can damage the forging as the mechanical strength of the thick flashing becomes significant. Prediction of defects, such as flashing, can be performed using SPH because of its ability to solve problems involving very high strains in material. The correct prediction of workpiece size using simulation can assist in avoiding such potential forging problems arising from the use of incorrect size for the workpiece. 3.3. Effect of material hardening on plastic deformation Next we investigate the effect of hardening on the material flow during forging. One of the advantages of forging over other manufacturing processes, such as casting, is to be able to control the grain direction in the metal microstructure of the manufactured part through control of the applied stresses. It is useful to colour the initial workpiece in horizontal strata. This allows us to monitor the deformation of regions and the changes that occur to the strata as the workpiece deforms under the applied force from the top die. Fig. 7 shows the same forging using two different work hardening moduli for the plasticity. The hardening modulus used on the left is 0.0167 MPa whilst on the right is 1.67 MPa. The material is coloured in strata to show the difference in the deformation patterns. The effects of deformation on the grain microstructure can be deduced from the orientations of the interfaces between colours. For the case of lower hardening modulus, the surface of the workpiece (blank) remains smooth during compression. For metal with a high hardening modulus, the surface of the workpiece strain hardens and becomes rippled with small scale rounded structures protruding from the surface. These lead to lap formation where surface material (including oxide) is folded back into the forging leading to regions that will be mechanically weaker. These are issues one might expect in forging harder metals such as titanium. The ability of SPH to automatically track free surface evolution allows complex surface and defect structures such as lap formation to be directly predicted as part of the solid flow solution. This is an attractive aspect of SPH for these types of manufacturing processes. There are other non-trivial differences in the filling produced by the strong hardening of the material, specifically relating to the distribution of the blue material (the first strata to be extruded). In the core of the forging, the strata boundaries are flatter and more horizontal indicating that the deformation is more localised near the initial contact points with the lower die and that the effect of shear in the direction of motion of the die is insignificant. Under the compression from the top die,
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3847
Fig. 7. Effect of hardening parameter on material deformation under compaction. The hardening modulus for the right column is 100 times that of the left column. Particles are coloured in initially horizontal strata to show the pattern of metal deformation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
the blue harder material near the contact points pushes straight down (refer to 750 and 950 ms frames of Fig. 7) whereas the softer material bulges down smoothly. This could be interpreted as there being greater resistance to bending and the existence of higher shear strain gradients for the harder material than the softer material.
3848
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
This example also demonstrates prediction of one of the issues relating to improper design of the die, that being the ability of the workpiece to move laterally within the die during loading, leading to asymmetry in forged components. Here the workpiece is able to move slightly to the right as it is being pushed into the die, leading to an asymmetric filling pattern with the left side moving faster than the right. In more extreme cases this can lead to under filling of the die in some regions and flashing from an overfilled die on the opposite side of the workpiece. 4. Industrial scale complex forging processes: flow pattern and plastic deformation In this section, simulation of realistic forging processes is considered to establish the viability of implementing SPH modelling for improving the efficiency of forging operations and obtaining the desired product quality. The forging processes simulated incorporate complex three-dimensional geometries of the forged products and equipment that are representative of real components. These components were made of an Aluminium alloy (A6061) with properties specified in Table 1. Fig. 8 shows the geometry of the first of the 3D components. It consists of 4 cylinders of different radii placed concentrically on top of one another. The radii are 40, 60, 80 and 120 mm, respectively. In the same order, their heights are 30, 30, 50 and 100 mm. Each of the 8 pins is 25 mm high and of 7.5 mm radius. In the simulation, 51,283 SPH particles of 4 mm spacing were used to represent the die, piston and cylindrical workpiece (blank). This is a difficult test case with significant geometric complexity. It is a representative of many forgings where the metal is plastically forced into a complex die shape by a single linear movement of one die (on the right) moving directly towards the stationary die on the left. This is a considerably more geometrically complex version of the demonstration example shown in Fig. 1. Fig. 9 shows the forging process as the piston compresses the blank into the die. The workpiece is coloured by equal width strata which are orthogonal to the direction of motion of the moving die in order to track the material deformation. These strata could as easily be different materials (i.e. a composite) and the mixing and deformation of each of these layers could be investigated. Initially, the relative positions of these colour blocks are fairly constant and there is little deformation as the blank is pushed into the die due to the blank having a smaller diameter than the diameter of the lower cylindrical sections of the die. At 0.25 s, the front of the workpiece enters the second cylindrical region from the top. As this cylinder diameter is smaller than the workpiece, material outside this diameter and near the outer perimeter is blocked and comes to near rest below the second cylinder. The material from the central region of the workpiece is pushed up into the space of the second cylinder. Note the small amount of cyan coloured material near the top at time 0.55 and 0.60 s. By 0.65 s, some material reaches the top cylinder. As the top cylinder has a smaller radius than that of the second cylinder, only material from the central region of the workpiece flows into this. These successive geometric constrictions create increasing resistance to the flow of material upwards. This causes the material to slow down and then, when further compressed by the piston die below, to spread outward laterally into the spaces around the outside of the workpiece. By 0.70 s, material near the piston has spread out sideways to reach the outer surface of the largest cylindrical region of the die (green and yellow material). At 0.72 s, the leading material begins to fill the eight pins that extend from the largest diameter cylindrical region at the base. Material can also be observed starting to flow into the cross shaped extensions from the cylinder at the top level of the die. By 0.77 s, the forging is completed. Fig. 10 shows the progression of this 3D forging with the workpiece coloured by plastic strain where blue indicates 0 and red indicates a plastic strain of 4.0. At 0.25 s when the front of the workpiece has only just started entering the second cylindrical region from the top, the metal is still not plastically loaded and there is very little deformation. At 0.50 s the central portion of the workpiece has entered into the second cylindrical region from top and very weak plastic deformation is experienced at the periphery of this central region. At 0.55 and 0.60 s the circular peripheral portion of the workpiece around the
Fig. 8. Die geometry for the first of the realistic 3D forgings.
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3849
Fig. 9. Progress of the elastoplastic deformation of the metal during forging of the first of the 3D components. The material is coloured by its initial position into strata orthogonal to the piston motion.
second cylinder is further plastically deformed with the plastic strain levels reaching values of around 1.0. By 0.65 s as some metal starts moving into the top cylinder, a second high plastic strain section of smaller circumference around the top cylinder develops. At 0.70 s as the material near the piston spreads out sideways, plastic deformation is observed close to the end of the lowest cylinder. At the same time plastic strain levels increase around the periphery of the upper two cylinders. By
3850
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
Fig. 10. Progress of the elastoplastic deformation of the metal during forging of the first of the 3D components. The material is coloured by the value of plastic strain with blue for 0.0, green for 2.0 and red for 4.0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
0.72 s the metal starts filling the eight pins that extend from the cylindrical base. Due to high deformation, the plastic strains reach a value of around 2.0 in this region. The cross shaped extensions from the top cylinder also start getting filled at this time with plastic strain levels in excess of 2.0. By the time the forging is completed at 0.77 s, maximum plastic strain of
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3851
Fig. 11. Final shape of the first 3D forging. The material is coloured by (left) the initial position of the material using strata orthogonal to the piston motion, and (right) plastic strain, with dark and the light blue for low value (0–1.0), green for medium strains (1.0–2.0) and yellow for the highest plastic strains (2.0–3.0). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 12. Die geometry for the second of the realistic 3D forgings.
around 3.0 is experienced by the metal that has entered into the pins at the bottom of the cylinder and the cross shaped section in the top cylinder. Fig. 11 shows the shape of the final forging. The upper regions of the forging are composed entirely of the dark blue coloured material, whilst the top surface of the largest diameter base cylindrical section is a blended mix of green and yellow materials. It is entirely yellow at the outer edge, green at the inner edge and grades progressively between these. The material mixing is not radially symmetric as there are significant small scale mottled variations. This shows that the materials are being mixed by the plastic deformation process in a way that is not completely uniform at the smaller scales. Fig. 11 also shows the forging coloured by the plastic strain. Highest plastic strains of around 3.0 are observed in the pins extending from the cylindrical base plate and in the cross extensions from the top cylindrical section of the die. These small scale features create the highest strains which then considerably affects the material microstructure. Since mechanical strength is of paramount consideration in these attachment structures, the ability to predict the state of plastic deformation in these locations is important. Fig. 12 shows two views of the geometry of the die for forging the second 3D component. The geometry does not include the four tubes of radius 20 mm extending from the four arms of the die that guide the motion of the pistons for compressing the workpiece into the shape of the die. The length of each of the four arms of the die is 75 mm. The radius of each arm is 20 mm and the radius of the enlarged section of the arm is 30 mm. A cylindrical workpiece of length 640 mm and radius 18 mm is placed in the die at the start of forging. This component has significant geometrical complexity and is a challenge to forge and to model since it has four interacting dies with two being fixed and with two being pushed from non-aligned directions. In the SPH simulation, 31,952 particles of 4 mm spacing were used to represent the workpiece, die and pistons.
3852
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
Fig. 13. Different stages of forging of the second 3D component. Material is coloured in vertical bands to show the deformation patterns. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 13 shows the different stages of forging. The particles in the workpiece are coloured in ten vertical strata to allow the deformation to be tracked. As the workpiece is compressed by the two horizontal pistons, more material is pushed into the die. At 1 s, three bulges are observed in the centre of the workpiece. At 2 s, the material begins to fill the two vertical arms of the die. The wider diameter sections of the horizontal arms are now almost filled. For times up to 2 s, the material flow is fairly symmetric with respect to a vertical plane through the middle of the die. At 2.25 s, the material motion upwards and downwards in the vertical arms is halted and the material is forced sideways at the ends. At 2.5 s, the metal in the vertical arms is being pushed sideways towards the side walls and the die is almost full. Forging is completed at 2.8 s with all the voids filled. Fig. 14 shows the workpiece coloured by plastic strain. Blue indicates a plastic strain of 0.0, green indicates 30.0 and red indicates a high plastic strain of 60.0. The material in this forging case is deformed intensively resulting in plastic strains
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3853
Fig. 14. Different stages of forging of the second 3D component. Material is coloured by plastic strain with blue indicating 0.0, green 30.0 and red indicating 60.0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
reaching very high values. At 1 s when only the central region is starting to deform, the plastic strain levels across the entire workpiece are relatively low. The sections of the metal at the far ends from where the metal is pushed into the central region start getting plastically loaded, so do the regions around the two big central bulges. As the metal deforms further at 2 s the plastic strain levels increase dramatically reaching values of around 30.0 in most parts of the workpiece. The sections of the workpiece close to the two horizontal ends and the ends of the metal that have started filling in the vertical regions experience strain levels in excess of 50.0. At 2.25 s when the metal starts filling sideways into the vertical sections, very high
3854
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
strain levels are experienced in these regions. Very high plastic strain is also still experienced at the two horizontal extremities of the workpiece in contact with the piston. By 2.8 s when the forging process is complete, very high strain levels in excess of around 60.0 are observed at all four extremities of the workpiece. For such forged parts, post heat treatment (e.g. annealing) is essential to prevent cracking in regions that have experienced very high plastic strains. 5. Conclusions In this paper we have demonstrated that SPH can be used to model forging processes. SPH is established to be a useful simulation tool that can provide insights into the nature of the solid metal flow, the overall deformation of different parts of the forged component and the local deformation of the metal. SPH has several important attributes that make it very well suited to solid material forming processes including: Very strong stability and robustness due to the absence of issues relating to mesh shape and quality under high deformations because of its mesh-less particle based nature. Ability to continue to simulate plastic deformation processes at extremely high strains (up to a plastic strain of 400% which is well beyond the limits typically observed in forging). Ability to follow the surface deformation of workpiece including self-collision automatically without any special treatment. Ability to interact with any number of die components of any geometric complexity and performing any type of motion with either kinematic or dynamic motion controls. The inherent ability to carry material history with each particle. This includes plastic strain, current yield stress due to strain hardening and material composition (such as when multiple materials or metal composites are used). SPH can also directly predict several types of defects, such as incomplete die filling (due to inadequate forging pressure or workpieces that are too small), flashing due to workpieces that are too large, lateral movement of unstably positioned workpieces leading to unexpected and/or sub-optimal material flow and lap formation (leading to internal seam lines and surface oxide entrapment). It also has the potential to provide new types of information for forging optimisation because of its inherent history carrying capability, specifically direct predictions of defects such as: Surface oxidation (scale) that can form if the metal is heated too much prior to forging. Coarse grain structure: A final forging may contain coarse grains if the finishing temperature is too high or strain is too low, which adversely affects mechanical properties. Local grain structure models can be easily implemented to predict these properties. Fracture: When the stresses due to external loads exceed the failure strength, material fracture may occur in the workpiece or dies. Explicit damage models can be easily included in SPH to allow direct prediction of crack initiation and growth. Flash cracks: In flash regions rapid cooling occurs, so excessive pressures can build up in the workpiece, resulting in cracks propagating into the bulk material. Cold shuts: Poor material flow patterns lead to metal folding over on itself forming weak seams that are pressed into the surface.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
A.N. Bramley, D.J. Mynors, The use of forging simulation tools, Mater. Des. 21 (2000) 279–286. S.L. Semiatin, Introduction to Bulk-Forming Processes. ASM Handbook, Metalworking: Bulk Forming (#06957G), vol. 14A, 2005. K. Serope, R. Steven Schmid, Manufacturing Engineering and Technology, Prentice-Hall, Upper Saddle River, NJ, 2006. Y.H. Moon, C.J. Van Tyne, W.A. Gordon, An upper bound analysis of a process-induced side-surface defect in forgings: part 2: characteristics and criteria curves, J. Mater. Process. Technol. 99 (1–3) (2000) 179–184. T. Ishikawa, N. Yukawa, Y. Yoshida, Y. Tozawa, Analytical approach to elimination of surface micro-defects in forging, CIRP Ann. Manuf. Technol. 54 (1) (2005) 249–252. J.R. Cho, D.Y. Yang, Three-dimensional finite element simulation of a spider hot forging process using a new remeshing scheme, J. Mater. Process. Technol. 99 (2000) 219–225. G. Li, T. Belytschko, Element free Galerkin method for contact problems in metal forming analysis, Eng. Comput. 18 (2000) 62–78. O. Maillard, P. Montmitonnet, P. Lasne, Simulation of cold forging by a 2D elastoplastic finite element method with automatic remeshing, J. Mater. Process. Technol. 34 (1–4) (1992) 93–100. C. PavanaChand, R. KrishnaKumar, Remeshing issues in the finite element analysis of metal forming problems, J. Mater. Process. Technol. 75 (1–3) (1998) 63–74. D. Peric, M. Dutko, D.R.J. Owen, Aspects of adaptive strategies for large deformation problems at finite inelastic strains, Stud. Appl. Mech. 47 (1998) 349–363. M. Zhan, L. Yuli, Y. He, Research on a new remeshing method for the 3D FEM simulation of blade forging, J. Mater. Process. Technol. 94 (2) (1999) 231– 234. G. Barton, X. Li, G. Hirt, Finite-element modeling of multi-pass forging of nickel-base alloys using a multi-mesh method, Trans Tech Publications Ltd., Vancouver, Canada, 2007, pp. 2503–2508. C.C. Wang, Finite element simulation for forging process using Euler’s fixed meshing method, Mater. Sci. Forum 575-578 (2008) 1139–1144.
P.W. Cleary et al. / Applied Mathematical Modelling 36 (2012) 3836–3855
3855
[14] J.R. Cho, N.K. Lee, D.Y. Yang, Three-dimensional simulation for non-isothermal forging of a steam turbine blade by the thermoviscoplastic finite element method, Proc. Inst. Mech. Eng. B. J. Eng. Manuf. 207 (B4) (1993) 265–273. [15] C. Guedes, J. César de Sá, A proposal to deal with contact and friction by blending meshfree and finite element methods in forming processes, Int. J. Mater. Form. 1 (4) (2008) 177–188. [16] L.-c. Liu, X.-h. Dong, C.-x. Li, Adaptive finite element-element-free Galerkin coupling method for bulk metal forming processes, J. Zhejiang Univ. Sci. A 10 (3) (2009) 353–360. [17] S. Yoon, C.T. Wu, Hui-Ping Wang, J.S. Chen, Efficient meshfree formulation for metal forming simulations, J. Eng. Mater. Technol. 123 (2001) 462–467. [18] J.-S. Chen, C. Pan, C.M.O.L. Roque, H.-P. Wang, A Lagrangian Reproducing Kernel Particle Method for metal forming analysis, Comput. Mech. 22 (3) (1998) 289–307. [19] G. Li, K. Sidibe, G. Liu, Meshfree method for 3D bulk forming analysis with lower order integration scheme, Eng. Anal. Bound. Elem. 28 (10 SPEC. ISS.) (2004) 1283–1292. [20] H. Wang, G. Li, X. Han, Z.H. Zhong, Development of parallel 3D RKPM meshless bulk forming simulation system, Adv. Eng. Softw. 38 (2) (2007) 87–101. [21] H. Liu, Y. Yang, C. Li, Meshfree method numerical simulation of sheet metal forming process, Wuhan Ligong Daxue Xuebao/J. Wuhan Univ. Technol. 28 (Suppl. 1) (2006) 815–819. [22] H.S. Liu, Z.W. Xing, Y.Y. Yang, Simulation of sheet metal forming process using Reproducing Kernel Particle Method, Commun. Numer. Methods Eng. 9999 (9999) (2009). n/a. [23] Y.H. Park, Rigid-plastic meshfree method for metal forming simulation, in: ASME 2003 Pressure Vessels and Piping Conference, Paper no. PVP20031907, 2003, pp. 231–236. [24] Y.H. Park, Numerical study of metal forming simulation using elasto-plastic and rigid-plastic meshfree analysis, in: ASME 2005 Pressure Vessels and Piping Conference, Paper no. PVP2005–71006, 2005, pp. 385–395. [25] H. Guang-ru, W. Hong-yu, A. Li-qiang, Thermal-mechanical analysis of metal forming based on meshfree method, in: Materials Science Forum, 2008, pp. 273–281. [26] S. Xiong, P.A.F. Martins, Numerical solution of bulk metal forming processes by the Reproducing Kernel Particle Method, J. Mater. Process. Technol. 177 (1-3) (2006) 49–52. [27] S. Xiong, J. Rodrigues, P. Martins, On background cells during the analysis of bulk forming processes by the Reproducing Kernel Particle Method, Comput. Mech. 40 (2) (2007) 233–247. [28] J.J. Monaghan, Smoothed particle hydrodynamics, Ann. Rev. Astron. Astrophys. 30 (1992) 543–574. [29] J.J. Monaghan, Smoothed particle hydrodynamics, Rep. Prog. Phys. 68 (2005) 1703–1759. [30] P.W. Cleary, M. Prakash, J. Ha, N. Stokes, C. Scott, Smooth particle hydrodynamics: status and future potential, Prog. Comput. Fluid Dyn. 7 (2–4) (2007) 70–90. [31] P. Cleary, J. Ha, V. Alguine, T. Nguyen, Flow modelling in casting processes, Appl. Math. Model. 26 (2) (2002) 171–190. [32] Y. Imaeda, S.-I. Inutsuka, Shear flows in smoothed particle hydrodynamics, Astrophys. J. 569 (1) (2002) 501–518. [33] L.D. Libersky, A.G. Petschek, Smooth particle hydrodynamics with strength of materials, Lect. Notes Phys. 395 (1991) 248–257. [34] C.A. Wingate, H.N. Fisher, Strength Modeling in SPHC, Los Alamos National Laboratory, Report No. LA-UR-93-3942, 1993. [35] J.P. Gray, J.J. Monaghan, R.P. Swift, SPH elastic dynamics, Comput. Methods Appl. Mech. Eng. 190 (49-50) (2001) 6641–6662. [36] J.P. Gray, J.J. Monaghan, Numerical modelling of stress fields and fracture around magma chambers, J. Volcanol. Geoth. Res. 135 (2004) 259–283. [37] Z.S. Liu, S. Swaddiwudhipong, C.G. Koh, High velocity impact dynamic response of structures using SPH method, Int. J. Comput. Eng. Sci. 5 (2) (2004) 315–326. [38] P.W. Cleary, M. Prakash, J. Ha, Novel applications of Smoothed Particle Hydrodynamics (SPH) in metal forming, J. Mater. Process. Technol. 177 (1–3) (2006) 41–48. [39] J. Bonet, S. Kulasegaram, Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations, Int. J. Numer. Methods Eng. 47 (6) (2000) 1189–1214. [40] J.L. Wilkins, Calculation of elastic–plastic flow, in: Methods of Computational Physics, Academic Press, New York, 1964, pp. 211–263.