A sequential DEM-FEM coupling method for shot peening simulation

    A sequential DEM-FEM coupling method for shot peening simulation Fubin Tu, Dorian Delbergue, Hongyan Miao, Thierry Klotz, Martin Levesque PII: DOI: Reference:

S0257-8972(17)30274-8 doi: 10.1016/j.surfcoat.2017.03.035 SCT 22204

To appear in:

Surface & Coatings Technology

Received date: Revised date: Accepted date:

7 October 2016 12 February 2017 15 March 2017

Please cite this article as: Fubin Tu, Dorian Delbergue, Hongyan Miao, Thierry Klotz, Martin Levesque, A sequential DEM-FEM coupling method for shot peening simulation, Surface & Coatings Technology (2017), doi: 10.1016/j.surfcoat.2017.03.035

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

RI P

T

A sequential DEM-FEM coupling method for shot peening simulation Fubin Tua , Dorian Delberguea,b , Hongyan Miaoa , Thierry Klotza , Martin Levesquea,∗ Laboratory for Multiscale Mechanics (LM2 ), CREPEC, Department of Mechanical Engineering Ecole Polytechnique de Montreal C.P. 6079, succ. Centre-ville, Montreal, Quebec, H3C 3A7, Canada b Laboratoire d’Optimisation des Procedes de Fabrication en Avances (LOPFA), Department of Mechanical Engineering Ecole de Technologie Superieure 1100 rue Notre-Dame Ouest, Montreal, Quebec, H3C 1K3, Canada

MA

NU

SC

a

ED

Abstract

AC CE

PT

Shot peening is a cold-working process widely used to form and enhance the fatigue life of metallic components. The process consists of projecting high-velocity particles onto a metallic surface. This study introduces a new, and experimentally validated, sequentially coupled Discrete Element Model (DEM) - Finite Element Model (FEM) to predict the process’ effects in terms of residual stresses and roughness. A shot stream was first simulated in DEM to obtain the velocity distribution of impacting shots. The target’s progressive hardening was accounted for by adjusting the Coefficients of Restitution (CoRs) for shot-target interactions as the number of impacts evolved through a meshless method. The extracted impacting shots were then impinged onto a representative cell in a dynamic FEM model to evaluate the shot peening effects. The simulations were compared against experimentally measured roughness and residual stresses at full coverage. The study shows that using a constant average CoR yields results that are quite similar to those with an evolving CoRs, for a fraction of the computational cost. Moreover, cases where shot-shot interactions were accounted for yielded lower ∗

Corresponding author Email address: [email protected] (Martin Levesque) URL: http://www.polymtl.ca/lm2 (Martin Levesque)

Preprint submitted to Surface and Coatings Technology

March 18, 2017

ACCEPTED MANUSCRIPT

RI P

T

residual stresses and roughness than cases where such interactions were not accounted for. Nevertheless, all simulated cases delivered simulations that were in good agreement with the experiments, which validates the proposed approach.

SC

Keywords: Shot peening, Residual stress, DEM-FEM coupling, Discrete Element Method, Finite Element Method, Inconel 718 1. Introduction

AC CE

PT

ED

MA

NU

Shot peening is a cold-working process widely used to form and enhance the fatigue life of metallic components. During the process, a large amount of small sphere-like particles, called shots, are bombarded onto a metallic surface, at high velocities (20∼150m/s) [1, 2]. These impacts induce tensile plastic strains and introduce beneficial compressive residual stresses in a thin surface layer. For a thick part, the residual stresses are approximately equal to the process-induced stresses since the bulk prevents any springback. For a thin shot peened part like a strip or a wing panel, the residual stresses are those at the end of the process, when all mechanical constraints have been removed. These stresses can be estimated by analytical means [3] if the process-induced stresses are known. The shot peening process simulation is a challenging topic, for a number of reasons. First, the shot stream is constituted of numerous shots of different sizes traveling at different velocities. Second, the stream interacts with the treated part and shots can also collide between each other before and/or after hitting the treated surface. The actual shot size-velocity distribution hitting the surface to be treated is usually unknown. Finally, the residual stresses estimation involves dynamic elasto-plastic computations that are sensitive to the constitutive theory and inputted material parameters. The Finite Element Method (FEM) has been used over the last decades to predict shot peening-induced residual stresses. For example, earlier attempts assumed deterministic impacting shots distributions (i.e., where shots positions, sizes and velocities are known a priori) on axisymmetric [4, 5], periodic [6, 7] and non-periodic [8, 9, 10] cells. A number of authors [11, 12, 13, 14, 15, 16, 17] adopted stochastic 3D FE models where the shots positions were randomly determined to capture this aspect of the process. Rayleigh damping [11, 12, 18] or infinite elements [9, 13, 19] were used in the simulations to prevent shock waves reflection and hence obtain a time stabilized residual stress 2

ACCEPTED MANUSCRIPT

AC CE

PT

ED

MA

NU

SC

RI P

T

field. To reduce the computational cost associated with such simulations, Mylonas et al. [14] optimized the peening time interval between two consecutive impacts while Kang et al. [10] implemented a 3D explicit dynamic analysis for shot-target impacts and an implicit static springback analysis for quickly stabilizing the impact waves. Moreover, shot-shot interactions are typically neglected to further reduce the simulation’s computational time. Finally, the surface over which the impacts are accounted for is typically of the order of 1∼3 times the shots diameters [11, 12, 16], again for computational reasons. More recent works focused on shot stream simulation through the discrete element method (DEM) [2, 20, 21, 22]. In the DEM [23], particles in motion obey Newton’s second law and their contact interactions are governed by specific force-displacements laws. Although several elasto-plastic contact models have been proposed [24, 25], most authors rely on the Coefficient of Restitution (CoR) e [26], defined as the ratio of relative velocities after and before an impact, to account for the energy dissipated during an impact. Up to now, the DEM has been widely used to describe the static or dynamic behavior of granular materials and an open-source DEM software Yade [27, 28] has been developed by a community of users. Several authors [2, 20, 21, 22] relied on the DEM to simulate shot peening streams. In that context, Bhuvaraghan et al. [29] relied on average and constant CoRs for shot-target interactions while Murugaratnam et al. [21] adjusted a facet-wise CoR for every shot-target contact as the number of impacts increased. Nguyen et al. [30] used computational fluid dynamics code ANSYS-FLUENT to simulate the air-shot-solid wall interactions and combined the results with a single impact FE model for predicting coverage. Badreddine et al. [31] proposed a CAD-based model to simulate shot dynamics in ultrasonic shot peening without considering shots rotations. The above-mentioned models are suitable for simulating a large number of shot-shot and shot-target interactions. However, they cannot accurately predict the target’s hardening behavior as well as the resulting residual stresses and hardness. A number of authors have coupled DE models, where the shot stream is realistically accounted for, to FE models for computing residual stresses [1, 32, 33, 29, 21, 22]. According to the concurrent coupling scheme [1, 32, 33], the DE and FE models communicate with each other at every simulation time step. By opposition, in the sequential coupling scheme [29, 21], a shot stream is first simulated with the DEM, and shots velocities or contact forces as well as their positions are sequentially transferred into a dynamic FE model. Han et al. [1, 32, 33] developed a concurrent DE-FE coupling method in 3

ACCEPTED MANUSCRIPT

AC CE

PT

ED

MA

NU

SC

RI P

T

which the shot peened part was simulated with FEM, shots were modeled with DEM and their interactions were predicted by a penalty-based approach, to simulate a shot peening process. Only a limited number of shots were simulated and shot-shot interactions were not accounted for. Wang et al. [34] concurrently combined smoothed particle hydrodynamics (SPH) where the interactions among adjacent shots were accounted for to a FE model where shot peening effects were simulated. Hong et al. [2] independently employed the DEM to study a shot stream and the FEM to simulate the elastic-plastic dynamic shot-target impacts. Bhuvaraghan et al. [29] converted the shot-target contact forces obtained from a DE model to contact pressures, which were sequentially input into a dynamic FEM model. Although applying contact pressures in a FE model is more straightforward than running contact analyses, this method must deal with time dependent pressure application areas, an approximate method for converting contact forces into pressures as well as a way to transfer circular contact areas onto a triangular or quadrilateral mesh. All these challenges have a detrimental effect on the method’s accuracy. By opposition, shots impacting velocities and positions can be straightforwardly prescribed as initial conditions into a dynamic FE model. Murugaratnam et al. [21] predicted shots velocities and positions from a DE model and sequentially transmitted them into a dynamic FEM model. Rousseau et al. [22] adopted the same sequential DEM-FEM coupling scheme to analyze the effect of bead quantity on a peened surface in ultrasonic shot peening. The concurrent DEM-FEM coupling scheme is physically motivated but has been shown to require lengthy computations. Moreover, the different length scales in both DE and FE models represent a challenge because the dimension of a DE model must be significantly reduced to achieve a seamless coupling with a FE model [35]. In the sequential DEM-FEM coupling scheme, the DE and FE simulations are running independently of each other, which results in a complete separation of both length and time scales [35]. This feature, although less physically representative than the concurrent scheme, is computationally efficient. This survey reveals that a number of limitations are still existing in the shot peening treatment prediction: 1. In all the FE-based methods [11, 12, 13, 14, 15, 16, 17], the shots impacting velocity was typically assumed as a constant and the distributed shots impacting velocities induced by shot-shot and shot-target interac4

ACCEPTED MANUSCRIPT

NU

SC

RI P

T

tions were not accounted for. It was consequently deemed to be of interest to study the effects of accounting for, or not, the effects of the impacting shots distribution on the residual stresses and roughness. 2. Murugaratnam et al. [21] adjusted the CoRs as the number of impacts increased, for each discrete facet. Only the facet containing the contact point affected a shot-target interaction. This approach required a very fine facet mesh before reaching convergence. Moreover, the authors did not study the CoR’s influence on the impacting velocities distribution. More importantly, the CoR’s influence on the residual stress profile was not studied in detail. 3. No experimental comparisons were used to assess the relevance of all the presented DE-FE coupling works published so far.

AC CE

PT

ED

MA

Coverage is a shot peening process control variable that is defined as the ratio of the area covered by peening dimples to the total treated surface area. Coverage is directly related to the peening time and is either assessed in industrial application through visual or automated inspection. A part where 98% of the surface is covered by dimples is known to be fully covered. Coverage evaluation from FE simulations has typically been evaluated by computing the surface area on the treated part having an accumulated plastic strain greater than an arbitrarily set threshold [11, 17]. This technique, in addition to be arbitrarily set, does not account for material that would visually seem as covered but has no plastic strain. Surface roughness is defined as the discrepancy between the peak and valley along a sampling line, which may provoke stress concentration. Finally, the peening induced stresses can bend a thin strip up to a certain height so as to equilibrate the residual stress profile. Peening intensity is defined as the point, on the curve of arc height against peening time beyond which the arc height increases by 10% when the peening time doubles. Intensity is measured on standardized steel strips called ‘Almen strips’. Peening intensity has been evaluated from FE simulations by combining simple beam theories [3, 11] or by more advanced peen forming models [12]. The main objective of this paper is to propose a sequential DE-FE coupling scheme for predicting the shot stream and the responses of a shot peened part in terms of induced residual stresses, coverage, surface roughness and peening intensity. A new meshless method for adjusting the CoR as the number of impacts at the same area increases has been developed in the DE portion of the proposed scheme. In addition, the FE model’s optimal 5

ACCEPTED MANUSCRIPT

NU

SC

RI P

T

dimensions and numerical efficiency has been investigated. The influence of the CoR on the shots impacting velocity distribution, induced stress profiles, surface roughness and peening intensity for IN718 has been studied in detail. The predicted residual stress profiles at full coverage were compared to an experimentally measured profile. This manuscript is organized as follows. Section 2 presents the equipment and procedures used for generating the experimental results. Section 3 deals with the DE and FE models, as well as their sequential coupling scheme. Section 4 presents numerical results while Section 5 concludes this work. 2. Experiment setup and procedures

AC CE

PT

ED

MA

2.1. Shot peening experiment Shot peening was carried out with an air-pressurized shot peening machine manufactured by Canablast and Genik. Figure 1(a) shows the peening cabinet while Figure 1(b) is a schematic representation of the peening equipment. The nozzle was moved by a FAUNC Robot M-20iA manipulator arm at a velocity of 22mm/s along x direction (see Figure 1). The nozzle’s radius rn was 6.35mm, while the stand-off distance d was 300mm. Spherical conditioned cut wire (CW14, density ρs =7800kg/m3 , nominal diameter ds =0.3556mm) shots were propelled at an average velocity v0 of 75m/s (measured with a ShotMeterG3 provided from Progressive Technology, USA) and with a mass flow rate of 6.80kg/min (controlled by MagnaValves provided from Electronics Inc). These parameters led to an Almen intensity of 8A, which is a shot peening intensity typically used in aerospace. Nickel SuperAlloy Inconel 718 (IN718) samples having dimensions of 76.2mm×50.8mm×10.2mm were submitted to this process until full coverage, as is typically done in aerospace. EBSD imaging confirmed that the material was untextured. 2.2. Residual stresses measurements The sin2 ψ X-ray diffraction (XRD) technique was adopted for measuring the macroscopic residual stress profile induced in the shot peened part. This technique uses the crystal lattice as a strain gauge and provides a local measurement of the current residual stresses state. XRD measurements were made by assuming a non textured material having an equi-bi-axial stress

6

AC CE

PT

ED

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Fig. 1. Shot peening equipment: (a) Picture of the actual setup; (b) Schematic diagram.

state. The measurements were carried out using a Proto iXRD diffractometer (Proto Manufacturing Ltd.) equipped with two linear detectors. A manganese (Mn) tube was used and the {311} family of diffraction planes was chosen for its high angle of measurements and its better accuracy on Δ θ measurements. The X-ray Elastic Constant (XEC) 1/2S2 was determined by a micro-tensile test following the same procedure as [36] for the diffraction planes. A specimen having the same chemical composition and microstructure as the shot peened samples was used. The parameters adopted for all the residual stresses measurements are summarized in Table 1. 7

ACCEPTED MANUSCRIPT

2mm 1/2S2 = 4.84×10−6 (MPa−1 )

Bragg angle

˚ Mn (λ=2.103 A)

7 (−28◦ ≤ ψ ≤ 28◦ ) {311} 2θ = 151.8◦

MA

XEC

Number of inclinations Diffraction planes

Cr

Aperture size

Target

Description

RI P

Filter

Item

SC

Detectors

Description PSSD (Position sensitive scintillation detector) 20◦ 2θ range

NU

Item

T

Table 1 X-ray diffraction parameters for residual stress measurements.

PT

ED

A layer removal technique relying on electropolishing was combined with XRD measurements to extract a residual stress profile through the thickness since the X-Ray penetration is relatively small in metals. The Moore and Evans correction [37] was used for the stress relaxation and redistribution caused by the layer removal process. This correction reads Z H Z H σm (y1 ) σm (y1 ) σ(y) = σm (y) + 2 dy1 − 6y dy1 (1) y1 y12 y y

AC CE

where σ is the corrected residual stress, y is the distance from the sample’s bottom surface to the measured point, σm is the measured stress and H the part’s initial thickness. 2.3. Surface roughness measurement Profile roughness measurement setup provided from Mitutoyo was adopted to measure the surface roughness of a shot peened part. The touch arm translated 15, 000μm along a sampling line and the length between each two sampling points was 1μm for every measurement. The sampling lines were randomly chosen in both longitudinal and transverse directions of the shot peened surface. The initial surface roughness of the IN718 part was 2.97 ± 1.46μm. 2.4. Mechanical properties measurement To obtain the mechanical properties of IN718, symmetric strain-controlled cyclic tests (strain amplitude Δ = 2%, strain ratio R = −1) were per8

ACCEPTED MANUSCRIPT

SC

3. Sequential DE-FE coupling scheme

RI P

T

formed on a MTS 318.25 uni-axial test machine at a constant strain rate of 1.0×10−3 s−1 . The cross-section area of a sample was 11.4mm×11.4mm and the gauge length was 9.0mm. Epsilon 3542-025M-025-ST extensometers were installed in both the longitudinal and transverse directions for measuring the deformations to determine the strains.

MA

NU

3.1. Overall strategy The proposed DE-FE sequential coupling method is a variant of that introduced by Murugaratnam et al. [21]. First, the shot stream was simulated in Yade to obtain the shots impacting velocities and spatial coordinates. The extracted shots impacting velocities and positions were then imported into ABAQUS/Explicit for computing the target’s response. Figure 2 schematizes the various simulation steps.

AC CE

PT

ED

3.2. DEM model 3.2.1. Governing equations The shots were described by spherical discrete elements and the metal part was represented by a rigid facet. The shot’s motion was governed by Newton’s second law as per [23, 27, 28]: m¨ x=

NC X

Fi

(2a)

i=1

N

C X 2 Mi mrs 2 ω˙ = 5 i=1

(2b)

where m is a shot’s mass, x is its current mass center, an overdot indicates time derivative, NC is the number of contacts around a shot, Fi is the ith contact force, ω is the shot’s angular velocity, and Mi is the momentum computed from Fi . The contact force was related to the relative displacement Ui (subscript i is omitted in sequel for simplifying the notation) through a contact law. The elastic parts of the shot-shot and shot-target interactions obeyed HertzMindlin nonlinear model [38] in which the normal and tangential contact forces were described by Feln = K n U n (3a) 9

T

ACCEPTED MANUSCRIPT

RI P

Input shot peening parameters

MA

Record the shots impacting velocities and locations just before they impact the target

NU

Build a dynamic FE model featuring a zone of interest on the target

Simulate the shot stream

No

SC

Build a model for shots and target (Shots are described by spherical particles and target is discretized by rigid facets)

Prescribe the extracted impacting velocities and locations as initial condition, execute dynamic impacting analyses

ED

Is the peening time reached? Yes

Extract the compressive residual stresses, plastic strain and displacement End

PT

End

AC CE

Fig. 2. Proposed sequential DE-FE coupling scheme.

Fels =

(

Fel0s + k s 4 U s μFeln

if sticking if sliding

(3b)

where subscript “el” indicates the elastic contact force, U n is the relative normal displacement, K n is the normal stiffness, Fel0s is the tangential force at the previous time step, k s is the tangent tangential stiffness and μ is the friction coefficient. K n and k s are expressed as [28] √ ! 2 hGi 2˜ r √ n (4a) Kn = U 3 (1 − hνi)

10

ACCEPTED MANUSCRIPT




13  2 hGi 3 (1 − hνi) r˜  |Feln | ks =  1 − hνi

(4b)

RI P

T

2

NU

SC

where hGi is the average shear modulus of two adjacent shots, hνi is their average Poisson’s ratio, 1r˜ = r11 + r12 with r1 and r2 the radii of two adjacent shots. The normal and tangential viscous contact forces were calculated by [39] r 2 ln e 5 n ˙n n n ˙n ˜U (5a) Fvis = C U = − p K m 2 ln e + π 2 4 r

5 n ˙s s K m Fvis = C s U˙ s = − p ˜U (5b) ln2 e + π 2 4 where C n and C s are respectively the normal and tangential viscous coefficients, m1˜ = m11 + m12 with m1 and m2 the masses of two adjacent shots and e is the CoR.

ED

MA

2 ln e

AC CE

PT

3.2.2. General model setup According to standard AMS2431-3D, CW14 shots diameters are comprised within 0.36 ± 0.025mm. This tolerance is quite tight and shots having identical diameters of 0.3556mm were simulated. The shot peening nozzle was translated along the x axis (see Figure 1(b)). Assuming that the nozzle’s plane center was (x0 , d, 0), a shot’s initial position was randomly generated according to   x = x0 + 2 (rn − rs ) rand (0, 1) − (rn − rs ) (6) y=d   z = 2 (rn − rs ) rand (0, 1) − (rn − rs )

where rs = ds /2 is the shot’s radius, rand (0, 1) is a uniform distribution random number generator delivering numbers in the range [0,1]. The generated shots lied inside a circular nozzle through the constraint (x − x0 )2 + z 2 ≤ (rn − rs )2

(7)

The shots initial velocities were assumed to be identical and oriented along y (see Figure 1(b)). This assumption had to be made since it was not possible 11

ACCEPTED MANUSCRIPT

MA

NU

SC

RI P

T

to measure the actual shots velocities and trajectories distributions. Note that the proposed algorithm could account for these distributions. CW14 shots Young’s modulus and Poisson’s ratio were respectively set to Es =200GPa and νs =0.3. The friction coefficient between the shots and the target was set to μ = 0.2, as was done in the FE simulations. Shots hardening was neglected, similar to [21] and a constant CoR es s = 0.4 was used for shot-shot interactions. The whole simulation area was 40mm×40mm in Yade. The CoR for a shot-target interaction, es t , is typically a function of parameters like shot velocity, shot material, angle of impact, target material hardness, etc [29]. As in [21], the hardening effect of a shot peened part was simulated in terms of the number of times an area was hit by shots. The meshless method for adjusting the shot-target CoRs is presented in the next subsection. Shots initial positions and velocities, as well as impacting positions and velocities when no shot-shot interactions were considered, were generated in Matlab.

PT

ED

3.2.3. Variation of CoR for shot-target interactions The target’s hardening effect on the CoR was accounted for through a new meshless approach (see Figure 3). The procedure was as follows:

AC CE

1. A set of points {x0i } corresponding to the FE nodes, were used to discretize the target’s surface. The number of times a point i was submitted to an impact was denoted by Ii while the CoR at point i after a number of impacts Ii was denoted by es ti (Ii ). 2. The CoR for an impact position x0 was calculated according to P |x0 −x0i |≤h es ti (8) es t (x0 ) = P 0 −x0 ≤h 1 x | | i

where h is the influence domain radius for one impact, as shown in Figure 3. Equation (8) can be interpreted as the average of CoRs within a circular area of radius h. 3. The number of impacted times was incremented for each point i located at x0i for which |x0 − x0i | ≤ h was met.

Parameters h and es ti (Ii ) were computed from single impact FE models and were tabulated. Section 4.2 provides more details. 12

ACCEPTED MANUSCRIPT

x'

h

RI P

x

SC

z

T

y

NU

Fig. 3. Schematic diagram of a meshless method for adjusting the CoR for impact position x0 .

PT

ED

MA

3.2.4. Simulation time optimization A simple method was devised so as to delete rebounding shots that no longer interacted with any incoming shot. A shot was deleted if any of the following conditions were met   x < −20mm      x > 20mm (9) z < −20mm    z > 20mm    y > 40mm, |v| ≤ 0.8v 0

AC CE

where |v| is the norm of v.

3.3. FEM model 3.3.1. Meshing and element choices The dynamic FE simulations were conducted with ABAQUS v6.14-1 commercial software. The shots were modeled as rigid spherical surfaces or deformable balls. The target was meshed by 8-node reduced integration 3D elements (C3D8R). Figure 4 shows an example of the mesh. The upper surface was divided into 5 regions shown in Figure 5. Region 1 was constituted of infinite elements for absorbing reflection elastic waves in the far field and served as displacement boundary conditions. Region 2 was the area occupied by coarse finite elements to transition from finite elements to infinite elements. Region 3 was the zone allowing impact indentations to be separated from region 2. text and tint have been optimized and more details are found in section 4.1. Region 4 was the zone confining all the impact centers 13

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

Fig. 4. Mesh.

AC CE

PT

ED

and its length was set to I = 3rs . The elements in region 4 were cubes whose sides had a length of rs /16. Region 5 was the representative zone that was used for extracting the average induced residual stresses. In this manuscript, R = I was used. x and z displacements were constrained on the bottom surface. The simulations were conducted until full coverage was reached in order to make direct comparisons with the experiments. 3.3.2. Material properties IN718 properties were obtained from cyclic tests and led to a Young’s modulus E = 205GPa and Poisson’s ratio ν = 0.32. The mass density was extracted from the literature as ρ = 8100kg/m3 [29]. A combined isotropickinematic hardening model implemented into ABAQUS was adopted to capture the material behavior under the action of cyclic loading. Stabilized cycle data was used to define the kinematic component. A single element simulation was performed to validate the tabular data extracted from cyclic test as instructed in [40]. Figure 6 shows the comparison between the cyclic test data and the simulation curve. A kinematic contact model [40] was used to capture the shot-target interactions. A number of authors [41, 42] have shown that the resulting residual stresses are relatively insensitive to friction coefficient for values ranging from 0.1 to 0.5. The friction coefficient was arbitrarily set to 0.2. 14

T

ACCEPTED MANUSCRIPT

RI P

text tint

R

SC

z

Region 2

Region 4 Region 5

MA

Region 3

NU

Region 1

I

x

1500 1200

PT

900

ED

Fig. 5. Schemetic view of the regions defined in the z-x plane at y = 0.

300

AC CE

Stress (MPa)

600

0

−300 −600 −900

−1200

−1500 −0.02

Experiment Simulation −0.012

−0.004 0.004 Strain

0.012

0.02

Fig. 6. Comparison between the cyclic test data and a single element simulation for IN718 (strain amplitude Δ = 2%, strain ratio R = −1).

15

ACCEPTED MANUSCRIPT

4. Numerical results

MA

NU

SC

RI P

T

3.3.3. Numerical strategies for reducing the computation time Different strategies were investigated to reduce the computational time. First, the analysis physical time duration was taken to be as short as possible to limit the number of explicit computation steps. This was done by assuming that no shots could hit the surface simultaneously and by reducing the time interval between two successive impacts. Second, rather than meshing all the shots at once, it was decided to exploit the *IMPORT function of ABAQUS [40] through which only a number of meshed shots were impacted and the final results were used as initial conditions in a subsequent analysis where new shots would hit the target. Such a procedure decreased the number of elements and contact pairs involved in a given analysis but required information to be transferred from the random access memory (RAM) to the hard disk, which decreased the computations performance. Section 4.1 presents the results of these analyses.

AC CE

PT

ED

4.1. FE model optimization The residual stress profiles for one and two impacts with deformable (elastic) shots are compared with those obtained from rigid shots. The mechanical properties of the wire steel were extracted from the literature [43]. It is worth mentioning that these parameters are those before conditioning. The conditioned shots were simulated as elastic balls since the plastic deformation of a shot leads to a non-spherical shot shape that is eliminated through the sieving system of the shot peening machine (in other words, only spherical shots hit the target). Figure 7 shows the residual stress profiles for elastic and rigid shots after one and two impacts. The residual stress profile for elastic shots is close to that for rigid shots. However, it should be noted that the computational time for two impacts with elastic shots is 4.6 times that for rigid shots. Figure 8 shows P EEQs (equivalent plastic strains) and residual stress σz profiles after a single shot impact at 75m/s, which is the highest velocity simulated. The figure shows that both the impact dimple radius and the yield zone radius were smaller than 0.5rs . It should be noted that the target hardens with the number of impacts, which should reduce the dimple size for successive impacts. tint was roughly set to 0.5rs for that reason. Figure 8 shows that using rs /16 as the elements sizes led to dimples meshed by more than 10 elements, which was judged sufficient [12]. Single impact analysis 16

ACCEPTED MANUSCRIPT

T

300

RI P

-300

SC

-600 -900

Rigid ball, One impact Elastic ball, One impact Rigid ball, Two impacts Elastic ball, Two impacts

-1200 -1500 0.1

0.2

MA

0.0

NU

Residual stress (MPa)

0

0.3

0.4

Depth (mm)

ED

Fig. 7. Residual stress profiles for elastic and rigid shots after one and two impacts.

AC CE

PT

revealed that a single shot-target impact duration was about 0.5 μs for a shot velocity of 75m/s, which should be the shortest interaction time in this study since the maximum shot velocity was 75m/s. The peening time interval between two successive impacts was arbitrarily set to Δtint = 1.5μs in order to avoid any interactions between two successive impacts. With tint fixed, 40 randomly positioned shots were impacted for text = {0.1778, 0.3556, 0.5334}mm. The induced stress profiles determined by averaging the σx and σz at every element layer in y direction for different text are shown in Figure 9. The figure shows that the residual stress profiles converged as text increased. text = 2rs =0.3556mm was chosen so as to limit the computational time. The *IMPORT option in ABAQUS was used to simulate 20 randomly distributed shots impacts in {1, 2, 4, 5, 10 and 20} analyses (for example, the *IMPORT option was used 19 times on single shot analysis). Each run was performed with 2 CPUs on a computer having an Intel(R) Core(TM) i7 4790 CPU and 32GB RAM. The computational time versus the number of ABAQUS/Explicit analyses is shown in Figure 10. It can be seen that using two separate runs led to the minimal computational time, which represents a decrease of 10% with respect to the single run. This result stems, most

17

ACCEPTED MANUSCRIPT

Z

MA

Residual stress (MPa)

NU

(a)

SC

RI P

T

0.5rs

ED

(b)

PT

Fig. 8. Contours of P EEQs and σz for a single impact: (a) P EEQs; (b) σz .

likely, from a trade-off between the reduced computational time of smaller models and the time used to access the hard disk.

AC CE

4.2. Coefficients of Restitution determination Multiple impacts at the same location were run to obtain the CoRs as a function of the number of impacts. Table 2 lists the CoRs es t for normally impacting shots at 75m/s as a function of the number of impacts. The number of impacts was set as Ii ≤ 4 based on the fact that each point is likely to be hit 3–4 times, at most [29]. Moreover, Table 2 suggests that the CoR increased with the number of impacts and stabilized after 4 impacts. This behavior is mostly due to the material’s hardening. Analyses with constant CoRs of {0.33, 0.4, 0.44} were also carried out, for comparison purposes. Table 2 CoRs for normally impacting with different numbers of impacts.

Number of impacts CoR es t

1 0.33

2 0.40

3 0.43

18

4 0.44

T

ACCEPTED MANUSCRIPT

RI P

400

SC

0 -200 -400

NU

Induced stress (MPa)

200

-600

text=0.3556mm text=0.5334mm

MA

-800

text=0.1778mm

-1000 0.0

0.1

0.2

0.3

0.4

Depth (mm)

ED

Fig. 9. Induced stress profiles for different text after 40 impacts at 75m/s when neglecting shot-shot interactions.

AC CE

PT

4.3. DE simulations results 32,000 shots were simulated with evolving shot-target CoRs as well as with CoRs of {0.33, 0.4, 0.44} on a computer having an Intel(R) Core(TM) i7 4790 CPU and 32GB RAM, and according to the methodology described at Section 3.2. Table 3 lists the computational time for each CoR. It can be seen that the computational time for evolving CoRs was about 18 times of that for the constant CoR=0.4. The results also shows that the simulation time decreased with increasing CoR, which might be explained by the fact that less shots remained near the surface, and hence were deleted, with higher CoRs than lower CoRs. Table 3 Computational time for different CoRs.

CoR Computational time (h)

Evolving 72

0.33 4.4

0.4 4.0

0.44 3.8

Figure 11 shows a snapshot of a shot stream near the shot peened part for evolving CoRs. The shots rebounding from the target are colored in blue 19

4

3.6x10

T

ACCEPTED MANUSCRIPT

RI P

nCPUs=2

SC

4

NU

2.8x10

4

2.4x10

2.0x10

4

MA

Computational time (s)

4

3.2x10

1

2

4

8

16

Number of ABAQUS/Explicit analyses

ED

Fig. 10. Computational time versus number of multiple ABAQUS/Explicit analyses for 20 randomly distributed impacts.

AC CE

PT

while those about to hit the target are colored in red. The figure suggests that some shots never hit the target due to their interactions with either the rebounding or incoming shots. Figure 12 provides the impacting shots velocities distribution where an impacting velocity in the range of 0 ∼7.5m/s implies that some of these shots never reached the target. The figure shows that, as the CoR increases, the number of shots never hitting the surface or hitting the surface with a velocity greater than 0.9v0 increases. The result for evolving CoRs lies between those of constant CoRs of 0.33 and 0.44. Figure 13 compares the impacting velocity distributions for simulations involving a fixed CoR of 0.4 with a fixed nozzle, a moving nozzle traveling at 22mm/s and 200mm/s along x. The figure shows the velocities distribution are insensitive to the robot’s velocity in this range, mostly due to the relatively small robot’s velocity, when compared to that of the shots. 4.4. Coverage assessment from the FE simulations A new coverage evaluation procedure was devised. The procedure tried to reproduce what the human eye would see from the deformed surface predicted from the FE simulations. The main idea was to compute the surface 20

MA

NU

SC

RI P

T

ACCEPTED MANUSCRIPT

ED

Fig. 11. A snapshot of a shot stream near a shot peened part: Shots rebounding from the target are in blue; Shots about to hit metal part or rebounding from other shots are in red.

AC CE

PT

normal at every position on the deformed surface and then compute the angle between that normal and the vertical axis y. The normals were computed from analytical surfaces generated on each element thanks to their shape functions. The surface normal n was calculated as per n=

∂x ∂x × ∂ξ ∂η

(10)

where x is a vector containing the (x, y, z) positions of the peened surface, (ξ, η) are the parent coordinates and × indicates the cross product. x was obtained from the classical relationship as x=

4 X

Ni (ξ, η)xi

(11)

i=1

i η) where Ni (ξ, η) = (1+ξi ξ)(1+η are the shape function, (ξi , ηi ) are the coordi4 th nates of the i node in the parent element (i.e., 1 or −1) and xi is a vector containing the (xi , yi , zi ) coordinates of the ith node in physical element. It

21

ACCEPTED MANUSCRIPT

40

RI P SC

25 20 15

NU

Distribution (%)

30

T

CoR=0.33 CoR=0.4 CoR=0.44 Evolving CoRs

35

10

MA

5 0 0~7.5

7.5~15 15~22.5 22.5~30 30~37.5 37.5~45 45~52.5 52.5~60 60~67.5 67.5~75 Impacting velocity (m/s)

ED

Fig. 12. Distributions of impacting velocities for CoR es t = 0.33, CoR es t = 0.4, CoR es t = 0.44 and Evolving CoRs.

AC CE

PT

was assumed in the sequel that the deformed surface was defined by the triplet (x, uy , z), where x and z are the current positions and uy is the y displacement. The ξ0 and η0 coordinates corresponding to specific x0 and z0 coordinates were generated by  4 P   Ni (ξ, η)z0i = 0 z0 − i=1 (12) 4 P   Ni (ξ, η)x0i = 0  x0 − i=1

through a Newton-Raphson method, where (z0i , x0i ) is the initial position of the ith node. The angle between the current normal n and (0, 1, 0) at (z0 , x0 ) was calculated by   n ∙ (0, 1, 0)T θ = arccos (13) |n|

Figure 14 (a) depicts the grey scale map for surface normals in Region 4 (see Figure 5) after 12 impacts for evolving CoRs. The figure was generated 22

ACCEPTED MANUSCRIPT

40

RI P

T

CoR=0.4, fixed nozzle CoR=0.4, moving nozzle 22mm/s CoR=0.4, moving nozzle 200mm/s

35

SC

25 20 15

NU

Distribution (%)

30

10

0

0~7.5

MA

5

7.5~15 15~22.5 22.5~30 30~37.5 37.5~45 45~52.5 52.5~60 60~67.5 67.5~75 Impacting velocity (m/s)

ED

Fig. 13. Comparison of the impacting velocity distribution: CoR es t = 0.4; fixed nozzle and moving nozzle.

AC CE

PT

by extracting the normals on an array where each point was separated by rs /64. Figure 14 (b) shows a binary image where areas for which θ < 5◦ were set to white. To obtain the coverage, the white zones under the impacting positions were manually filled in black by Paint as shown in Figure 14 (c). The corrected map shows that around 42.5% of the pixels are in black, which means that around 42.5% of the surface is covered by dimples. Although not perfect, and somewhat arbitrary since the binarizing angle was set to 5 ◦ , this procedure has the potential to better capture physical coverage. Future works will address its refinement. Figure 15 presents the coverage against number of impacts for evolving CoRs and without shot-shot interaction. Equation of the form Cf (nI ) = 100(1 − e(−mnI ) )

(14)

was used to fit the data, where Cf is the fitting coverage, nI is the the number of impacts and m is a fitting parameter. For evolving CoRs, m = 0.047, it became m = 0.074 when neglecting shot-shot interactions. For all the cases accounting for shot-shot interactions, full coverage was obtained after 84 23

SC

RI P

T

ACCEPTED MANUSCRIPT

(b)

(c)

NU

(a)

MA

Fig. 14. Grey scale map for surface normal after 12 impacts for evolving CoRs: (a) Contour plots for the θ angle; (b) Binarized contour plot with a threshold of θ = 5◦ ; (c) Corrected binarized contour.

100

ED

90 80

50 40

PT

60

AC CE

Coverage (%)

70

30 20

Evolving CoRs Without shot-shot interaction

10

0 0

20

40 60 Number of impacts

80

100

Fig. 15. Coverage versus number of impacts for cases: (a) Evolving CoRs; (b) Without shot-shot interaction. * and + denotes the coverage of one of six simulations.

impacts while it was reached after 60 impacts for case neglecting shot-shot interactions.

24

ACCEPTED MANUSCRIPT

RI P

T

400

0

SC

-200

CoR=0.33 -400

CoR=0.4

CoR=0.44

NU

Induced residual stress (MPa)

200

-600

Evolving CoRs

-800

MA

Without shot-shot interaction

-1000 0.0

0.1

Experiment

0.2

0.3

0.4

Depth (mm)

PT

ED

Fig. 16. Average induced residual stresses (σx and σz ) profiles at full coverage for different cases: CoR=0.33, CoR=0.4, CoR=0.44, Evolving CoRs, Without shot-shot interaction and Experiment.

AC CE

4.5. Induced stress profiles and comparison with XRD measurements Figure 16 depicts the average induced residual stresses (σx and σz ) profiles at full coverage for different CoRs. Table 4 lists the maximum compressive residual stresses and thicknesses of compressive residual stresses. Comparison with experimentally measured residual stresses reveals the case accounting for shot-shot interaction led to a better estimation of the maximal compressive residual stress, while neglecting shot-shot interaction captured the thickness of compressive residual stresses better. The simulated tensile residual stresses are larger than XRD measurements because infinite elements around the representative cell acted as constraints, while that was not the case for the experiments. It should be noted that Yang et al. [16] found that the tensile residual stresses approached 0 when imposing periodic boundary conditions, which should be more representative of the process. Figure 16 and Table 4 also suggest that the residual stresses field is fairly insensitive to the range of CoRs simulated.

25

ACCEPTED MANUSCRIPT

SC

Thickness of compressive residual stresses (mm) 0.15 0.15 0.15 0.15

-906

-877±44

0.16 0.19

MA

CoR=0.33 CoR=0.4 CoR=0.44 Evolving CoRs Without shot-shot interaction Experiment

Maximal compressive stress (MPa) -889 -888 -890 -885

NU

Case

RI P

T

Table 4 Maximal compressive stresses and thicknesses of compressive residual stresses for different cases: CoR=0.33, CoR=0.4, CoR=0.44, Evolving CoRs, Without shot-shot interaction and Experiment.

4.6. Surface roughness and comparison with experimental measurements The peak-low roughness parameter RPV is defined as

ED

RPV = max(uy ) − min(uy ).

(15)

AC CE

PT

where max(uy ) and min(uy ) are the peak and valley displacements along the measured line. 25 sampling lines were used for extracting roughness on the samples and 98 sampling lines were used to extract RPV from the simulations. The mean values and the standard deviations for the studied cases are listed in Table 5. Except for the simulation neglecting shot-shot interactions, the simulation results are in good agreement with the experimental measurements. 4.7. Peening intensity investigation Shot peening saturation curves of IN718 strips were simulated according to a procedure that is quite similar to the standardized Almen intensity measurement. The arc heights of IN718 strips were calculated from computed residual stress fields and simple beam theory for increasing number of simulated impacts. The strips arc heights were calculated from the bending moment equilibrating the residual stress field [11] Z t t Mx = − (16) σx ( − y)bdy 2 0 26

ACCEPTED MANUSCRIPT

and the arc height given by

NU

26.79±5.13

SC

RI P

Surface roughness RPV (μm) 20.03±1.74 19.31±3.40 19.98±3.21 20.71±4.69 22.78±4.25

MA

Case Experimental CoR=0.33 CoR=0.4 CoR=0.44 Evolving CoRs Without shot-shot interaction

T

Table 5 Surface roughness for different cases: Experiment, CoR=0.33, CoR=0.4, CoR=0.44, Evolving CoRs and Without shot-shot interaction.

Ah =

3Mx l2 2Ebt3

(17)

AC CE

PT

ED

where E is the strip’s Young’s modulus; l, b and t are strip length, width and thickness. Figure 17 presents the arc heights calculated from the induced residual stress profiles for evolving CoRs as a function of number of simulated impacting shots nI . The arc height was fitted according to Ah(nI ) =

p1 n I 2 + p 2 n I nI + p 3

(18)

where p1 , p2 and p3 are fitting parameters. The arc height at saturation was 0.1019mm for evolving CoRs and was reached for 46 impacts. Figure 18 shows the average arc heights versus number of impacting shots for the studied CoRs and the case without shot-shot interaction. The figure shows that, when compared to the simulations with evolving CoRs, the arc height was higher for the case without shot-shot interactions and slightly lower for CoR = 0.33, which is consistent with the impacting velocities distributions extracted from the DE simulations. Table 6 lists the arc heights at saturation and the corresponding impacting shots numbers. Table 6 reveals that the peening intensity for the case without shot-shot interaction was 5.2% higher than that resulting from evolving CoRs while that for a CoR=0.33 was 1.8% lower. It is worth noting that accounting for shot-shot interactions, 27

ACCEPTED MANUSCRIPT

T

0.12

RI P SC

0.08

0.06

0.04

NU

Arc height (mm)

0.1

0 0

20

MA

0.02

Saturation Curve Saturation time Double Saturation time

40 60 80 Number of impacting shots

100

120

ED

Fig. 17. IN718 strip arc height versus impacting shots number for evolving CoRs: 6 simulations. * denotes the arc height of one of the six simulations.

AC CE

PT

irrespectively of the CoRs, led to similar arc heights at saturation. Figure 19 presents the induced residual stress profiles at saturation, for all cases. The figure demonstrates that induced residual stress profile for peening inensity relies on shot-shot interaction. 5. Discussion

Aluminum alloy AA7050 (yield stress 464MPa) and copper alloy C11000 (yield stress 255MPa) were simulated as target materials to further investigate the influence of CoR on residual stress profile. These materials were selected as they are much more ductile that the simulated IN718. The mechanical properties of AA7050 were extracted from monotonic tensile tests while those of C11000 were extracted from [44]. Through the FE simulation for multiple impacts on the same position, the CoR for AA7050 was 0.30 after the first impact, while it was 0.42 after the fourth impact. For C11000, the CoR was 0.13 after the first impact and 0.2 after the sixth impact. Figure 20 and Figure 21 respectively show the residual stress profiles for AA7050 and for C11000 with the lowest and highest CoRs at full coverage. 28

ACCEPTED MANUSCRIPT

RI P

T

0.14 0.12

SC

0.08 0.06

Without shot-shot interaction CoR=0.33 CoR=0.4 CoR=0.44 Evolving CoRs

NU

Arc height (mm)

0.10

0.04

0.00 0

10

20

MA

0.02

30

40

50

60

70

80

90

100

Number of impacting shots

ED

Fig. 18. IN718 strip arc heights versus impacting shots numbers for different cases: CoR=0.33, CoR=0.4, CoR=0.44, Evolving CoRs and Without shot-shot interaction.

AC CE

PT

For these two softer material, the residual stress profiles are also insensitive to the CoR that is in the range of the lowest and highest CoRs. It could be possible to devise a fictitious material where the significance of using an evolving CoR could be evidenced. However, this study dealt with a wide range of real materials that did not exhibited a significant effect. 6. Conclusion

A velocity based sequential DE-FE coupling scheme was proposed to simulate the shot peening process. The shot stream was modeled in Yade to obtain impacting shots velocities and positions. The extracted shots were impinged upon an IN718 target in ABAQUS/Explicit to study peening effects like induced residual stress profile, surface roughness, arc height at saturation, etc. This work brought significant improvements to existing simulation techniques, namely: (1). A meshless method was proposed to adjust Coefficients of Restitution for shot-target interactions to account for the target’s hardening. 29

ACCEPTED MANUSCRIPT

T

Table 6 Peening intensities and corresponding numbers of impacts for different cases: CoR=0.33, CoR=0.4, CoR=0.44, Evolving CoRs and Without shot-shot interaction.

0.1001 0.1023 0.1025 0.1019

SC

CoR=0.33 CoR=0.4 CoR=0.44 Evolving CoRs Without shot-shot interaction

Number of impacts 45 50 50 46

RI P

Peening intensity (mm)

Case

18

MA

NU

0.1072

400 200

CoR=0.33

PT

-400

ED

-200

CoR=0.4

-600

AC CE

Induced stress (MPa)

0

CoR=0.44 Evolving CoRs

-800

-1000 0.0

Without shot-shot interaction 0.1

0.2

0.3

0.4

Depth (mm)

Fig. 19. Induced stress profiles related to peening intensities for all the cases: CoR=0.33, CoR=0.4, CoR=0.44, Evolving CoRs and Without shot-shot interaction.

(2). The shots that had no influence on others were deleted to reduce computational time in DEM. (3). Multiple ABAQUS/Explicit analyses are employed to reduce FEM simulation cost. (4). A new physically motivated method for extracting coverage from FE 30

ACCEPTED MANUSCRIPT

T

200

RI P

0

SC

-100 -200 -300 -400 0.1

0.2

0.3

MA

0.0

NU

Residual stress (MPa)

100

0.4

CoR=0.3 CoR=0.42

0.5

0.6

Depth (mm)

PT 0

AC CE

Residual stress (MPa)

100

ED

Fig. 20. Residual stress profiles for AA7050 with the lowest and highest CoRs.

-100

CoR=0.13 CoR=0.20

-200

0.0

0.1

0.2

0.3

0.4

0.5

0.6

Depth (mm)

Fig. 21. Residual stress profiles for C11000 with the lowest and highest CoRs.

simulations was proposed. The conclusions of this work are as follows: 31

ACCEPTED MANUSCRIPT

ED

MA

NU

SC

RI P

T

(1). The CoR had an effect on the shots impacting velocities distribution: the higher the CoR, the more the shots maintaining their initial velocity. Some shots never reached the target. (2). A constant CoR close to the average CoR delivered results that were close to the simulations with evolving CoRs. Simulations with constant CoRs were 18 times faster than those with evolving CoRs, which suggests that using an average CoR might be an efficient means for reducing computation time. (3). At full coverage, the induced residual stress profile was insensitive to CoRs when accounting for shot-shot interactions, which also suggests that using a constant CoR could provide means to reduce the computation costs. (4). At full coverage, the case neglecting shot-shot interactions overestimated the maximal compressive stress and the surface roughness. The arc height at saturation and the corresponding number of impacts were also overestimated. (5). The residual stress profiles and surface roughness for CoRs of {0.4, 0.44} and evolving CoRs matched experiment measurments.

AC CE

PT

These works suggest that the proposed new DEM-FEM coupling methodology provides reasonable results. It should be noted, however, that all simulations, accounting for shot-shot interactions or not, are very close to each other and to the experimental data. For this very specific case of a flat geometry, accounting for, or not, shot-shot interactions yields results that are in close agreement with experimental data. Different conclusions could be drawn for peening more complex geometries where shot-shot interactions would play a more significant role. Acknowledgements The work was supported by the research funding from the project CRIAQ MANU508 supported by Natural Sciences and Engineering Research Council of Canada (NSERC), Pratt & Whitney Canada, Bell Helicopter Textron Canada Ltd., L-3 Communications MAS and Heroux-Devtek Inc.

32

ACCEPTED MANUSCRIPT

References

RI P

T

[1] K. Han, D. Peric, A. Crook, D. Owen, A combined finite/discrete element simulation of shot peening processes-Part I: studies on 2D interaction laws, Engineering Computations 17 (2000) 593–620.

SC

[2] T. Hong, J. Ooi, B. Shaw, A numerical simulation to relate the shot peening parameters to the induced residual stresses, Engineering Failure Analysis 15 (2008) 1097–1110.

NU

[3] M. Guagliano, Relating almen intensity to residual stresses induced by shot peening: a numerical approach, Journal of Materials Processing Technology 110 (2001) 277–286.

MA

[4] K. Mori, K. Osakada, N. Matsuoka, Finite element analysis of peening process with plasticity deforming shot, Journal of Materials Processing Technology 45 (1994) 607–612.

ED

[5] K. Schiffner, C. Droste gen. Helling, Simulation of residual stresses by shot peening, Computers & Structures 72 (1999) 329–340.

AC CE

PT

[6] T. Kim, J. Lee, H. Lee, S. Cheong, An area-average approach to peening residual stress under multi-impacts using a three-dimensional symmetrycell finite element model with plastic shots, Materials & Design 31 (2010) 50–59. [7] F. Yang, Z. Chen, S. Meguid, 3D FE modeling of oblique shot peening using a new periodic cell, International Journal of Mechanics and Materials in Design 10 (2014) 133–144. [8] Y. Al-Obaid, Three-dimensional dynamic finite element analysis for shot-peening mechanics, Computers & Structures 36 (1990) 681–689. [9] M. Klemenz, V. Schulze, I. Rohr, D. Lohe, Application of the FEM for the prediction of the surface layer characteristics after shot peening, Journal of Materials Processing Technology 209 (2009) 4093–4102. [10] X. Kang, T. Wang, J. Platts, Multiple impact modelling for shot peening and peen forming, Proceedings of the Institution of Mechanical Engineers, Part B: Journal of Engineering Manufacture 224 (2010) 689–697.

33

ACCEPTED MANUSCRIPT

RI P

T

[11] H. Miao, S. Larose, C. Perron, M. Levesque, On the potential applications of a 3D random finite element model for the simulation of shot peening, Advances in Engineering Software 40 (2009) 1023–1038.

SC

[12] A. Gariepy, S. Larose, C. Perron, M. Levesque, Shot peening and peen forming finite element modellingtowards a quantitative method, International Journal of Solids and Structures 48 (2011) 2859–2877.

NU

[13] S. Bagherifard, R. Ghelichi, M. Guagliano, A numerical model of severe shot peening (SSP) to predict the generation of a nanostructured surface layer of material, Surface and Coatings Technology 204 (2010) 4081– 4090.

MA

[14] G. Mylonas, G. Labeas, Numerical modelling of shot peening process and corresponding products: residual stress, surface roughness and cold work prediction, Surface and Coatings Technology 205 (2011) 4480– 4494.

PT

ED

[15] S. Gangaraj, M. Guagliano, G. Farrahi, An approach to relate shot peening finite element simulation to the actual coverage, Surface and Coatings Technology 243 (2014) 39–45.

AC CE

[16] F. Yang, Z. Chen, S. Meguid, Effect of initial surface finish on effectiveness of shot peening treatment using enhanced periodic cell model, International Journal of Mechanics and Materials in Design 11 (2015) 463–478. [17] M. Jebahi, A. Gakwaya, J. Levesque, O. Mechri, K. Ba, Robust methodology to simulate real shot peening process using discrete-continuum coupling method, International Journal of Mechanical Sciences 107 (2016) 21–33. [18] T. Kim, H. Lee, H. Hyun, S. Jung, Effects of rayleigh damping, friction and rate-dependency on 3D residual stress simulation of angled shot peening, Materials & Design 46 (2013) 26–37. [19] P. Sanjurjo, C. Rodriguez, I. Penuelas, T. Garcia, F. Belzunce, Influence of the target material constitutive model on the numerical simulation of a shot peening process, Surface and Coatings Technology 258 (2014) 822831. 34

ACCEPTED MANUSCRIPT

T

[20] N.-L. C., M. Zarwel, C. Diviani, D. Hertz, H. Zahouani, T. Hoc, Surface impact analysis in shot peening process, Wear 302 (2013) 10581063.

RI P

[21] K. Murugaratnam, S. Utili, N. Petrinic, A combined DEMFEM numerical method for shot peening parameter optimisation, Advances in Engineering Software 79 (2015) 13–26.

NU

SC

[22] T. Rousseau, T. Hoca, P. Gillesb, C. Nouguier-Lehona, Effect of bead quantity in ultrasonic shot peening: Surface analysis and numerical simulations, Journal of Materials Processing Technology 225 (2015) 413– 420.

MA

[23] P. Cundall, O. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65.

ED

[24] L. Vu-Quoc, X. Zhang, L. Lesburg, Normal and tangential forcedisplacement relations for frictional elasto-plastic contact of spheres, International Journal of Solids and Structures 38 (2001) 6455–6489.

PT

[25] M. Brake, An analytical elastic plastic contact model with strain hardening and frictional effects for normal and oblique impacts, International Journal of Solids and Structures 62 (2015) 104–123.

AC CE

[26] C. Wu, L. Li, C. Thornton, Rebound behaviour of spheres for plastic impacts, International Journal of Impact Engineering 28 (2003) 929– 946. [27] J. Kozicki, F. Donze, A new open-source software developed for numerical simulations using discrete modeling methods, Computer Methods in Applied Mechanics and Engineering 197 (2008) 4429–4443. [28] Smilauer V. et al., Yade Documentation 2nd ed, The Yade Project, 2015. [29] B. Bhuvaraghan, S. Srinivasan, B. Maffeo, Shot peening simulation using discrete and finite element methods, Advances in Engineering Software 41 (2010) 1266–1276. [30] V. Nguyen, H. Poh, Y. Zhang, Predicting shot peening coverage using multiphase computational fluid dynamics simulations, Powder Technology 256 (2014) 100–112.

35

ACCEPTED MANUSCRIPT

RI P

T

[31] J. Badreddine, S. Remy, M. Micoulaut, E. Rouhaud, V. Desfontaine, P. Renaud, Cad based model of ultrasonic shot peening for complex industrial parts, Advances in Engineering Software 76 (2014) 3142.

SC

[32] K. Han, D. Peric, D. Owen, J. Yu, A combined finite/discrete element simulation of shot peening processes-Part II: 3D interaction laws, Engineering Computations 17 (2000) 680–702.

NU

[33] K. Han, D. Owen, D. Peric, Combined finite/discrete element and explicit/implicit simulations of peen forming process, Engineering Computations 19 (2002) 92–118.

MA

[34] F. Wang, J.and Liu, F. Yu, G. Zhang, Shot peening simulation based on sph method, The International Journal of Advanced Manufacturing Technology 56 (2011) 571578.

ED

[35] J. Wernik, S. Meguid, Coupling atomistics and continuum in solids: status, prospects, and challenges, International Journal of Mechanics and Materials in Design 5 (2009) 79110.

PT

[36] D. Delbergue, D. Texier, M. Levesque, P. Bocher, Comparison of two x-ray residual stress measurement methods: Sin 2 ψ and Cosα, through the determination of a martensitic steel x-ray elastic constant, Materials Research Proceedings 2 (2017) 55–60.

AC CE

[37] M. Moore, W. Evans, Mathematical correction for stress in removed layers in x-ray diffraction residual stress analysis, SAE Technical Paper 66 (1958). [38] C. Thornton, S. Cummins, P. Cleary, An investigation of the comparative behaviour of alternative contact force models during elastic collisions, Powder Technology 210 (2011) 189–197. [39] D. Antypov, J. Elliott, B. Hancock, Effect of particle size on energy dissipation in viscoelastic granular collisions, Physical Review E 84 (2011) 021303. [40] ABAQUS Analysis User’s Manual, ABAQUS Analysis User’s Manual, version 6.14, Dassault Systemes, 2015.

36

ACCEPTED MANUSCRIPT

RI P

T

[41] S. Meguid, l. G. Shaga, J. Stranart, 3D FE analysis of peening of strainrate sensitive materials using multiple impingement model, International Journal of Impact Engineering 27 (2002) 119–134.

SC

[42] M. Meo, R. Vignjevic, Finite element analysis of residual stress induced by shot peening process, Advances in Engineering Software 34 (2003) 569–575.

NU

[43] S. J441, (R) Cut wire shot, Technical Report, surface vehicle recommended practice, 1993.

AC CE

PT

ED

MA

[44] A. copper, A guide to working with copper and copper alloys, Technical Report, Copper development association, 2008.

37

ACCEPTED MANUSCRIPT

T

Highlights

RI P

• A new sequentially coupled DEM-FEM was introduced to predict the shot peening effects.

SC

• The influence of shot-shot interactions on the shot peening effects were studied. • The DEM-FEM simulation results were compared against experimental measurement.

NU

• The outward normal of the shot peened surface was used to predict the coverage.

AC CE

PT

ED

MA

• Techniques for reducing the computational cost were used in DE and FE simulations.

38