Accepted: June 24, 2019 | Published Online: June 26, 2019

Citation: Prawoto Y, Mazein PG, Kosterin AV, Ahmad Z (2019) Effect of Shot Size in Surface Improvement via Shot Peening; Analytical, Modeling and Experimental Approaches. Int J Metall Met Phys 4:033.

Copyright: © 2019 Prawoto Y, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

This paper discusses about the surface engineering improvement by means of shot peening. Shot peening improves the surface engineering quality by eliminating the tool marks, such as machining, grinding, stamping and other surface defects. Most importantly, the improvements of shot peening are produced by combination of compressive stress and cold work. Compressive stresses are beneficial in increasing resistances to fatigue failures, while the cold work effects of shot peening treatments can increase the surface hardness. It employs a steel shot media running at high speed to smash the material's surfaces systematically, which results in a compressed and condensed surface. Several factors influence the effectiveness of shot peening. They are namely shot nature, shot size, shot hardness, Almen intensity, shot angle and coverage. The focus of this research is finding the influence of the shot size. The approaches used are analytical, computational, and experimental studies. The center of the study is the influence of shot size in the generation of the compressive residual stress. Through analytical, computational and experimental approaches, it is confirmed that the shot size has influence on the depth of the peak residual stress rather than the surface residual stress. This research reaffirms the qualitative understanding with quantitative and comprehensive approaches that can also be used for further in depth research in surface engineering improvement especially via shot peening.

Keywords

Shot peening, Shot size, Residual stress, Surface engineering, Cold work

Introduction

Shot peening has been used for decades as a measure to overcome surface engineering problems in various industries. Its improvements are produced mainly by combinations of compressive residual stress and cold work. Compressive residual stresses are known to be beneficial in increasing resistances to fatigue failures and corrosion fatigue, while the cold work effects of shot peening treatments can increase the surface hardness of many materials [1]. It is believed to be the most economical and effective method of producing and making surface residual compressive stresses to increase the product life of treated metal parts. The increased strength of treated parts allows for lighter-weight parts that exhibit high wear and fatigue resistance.

The process can be defined as work hardening to the surface of components by propelling streams of spherical shots to the surface. The surface layer of material yields plastically to generate residual compressive stress. Among the practitioners, it has been known well that many parameters influence the efficiency of shot peening process. These are the peening coverage, saturation, shot material, shot size, speed, and peening time [2].

In this research, the focus is on the effect of the shot size. This topic is chosen due to the reason that fundamental understanding on this influence is low. For very long time, this subject has been regarded as a "black engineering" which the perceived of benefit using larger size is real but the reasons of those benefits were not disseminated by those who understand. In a way this paper is promoting to the academics for the science of shot peening, instead of keeping the "black art of shot peening" in the industry. This attempt is also on the line with the policy among the shot peeners [3] from long time ago.

Analytical Solution

The starting point for the construction of the shot peening model is the solution of the problem for a single impact of a ball on a half-space [4]. The initial data for the models are:

V - Speed of flying shot balls;

R - Radius of shot balls;

ρ - Density of the shot ball material;

E - Elastic modulus;

µ - Poisson's ratio;

${\sigma }_u=G_u({\epsilon }_u)$ - Material hardening curve.

A ball flying at a speed V hits the surface, contact forces develop in the contact zone, and the kinetic energy of the ball begins to be transformed into the elastic energy of the ball/half space and plastic deformation of the material to be shot peened. After the velocity of the ball becomes zero, the following condition is achieved:

$$\frac{m{v}^2}{2}=K(A\Pi +A\text{Ш})(1)$$

Where

$A\Pi$ - Is half-space energy;

$A\text{Ш}$ - Elastic energy of a ball;

$K$ - Correction factor, taking into account the influence of the roughness of the hydrodynamic film formed at the moment of impact.

Equation (1) is the starting point for the formation of boundary conditions for the contact problem with an inhomogeneous half-space. In order to calculate the work $A\Pi$, it is necessary to find the stress distribution in the half-space. To do this, we solve a system of differential equations describing the elastoplastic behavior of the material. Taking into account the axial symmetry, the set of equations is most conveniently written in the cylindrical coordinate system (r,θ,Z), where r - radius, Z - axis, directed along the axis of symmetry, θ - polar angle.

It is assumed that the mass forces are absent after collision, so the equations of equilibrium are as follows:

$$\frac{d{\sigma }_r}{dr}+\frac{d{\tau }_{rz}}{dZ}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=0(2a)$$

$$\frac{d{\sigma }_Z}{dZ}+\frac{d{\tau }_{zr}}{dr}+\frac{{\tau }_{zr}}{r}=0(2b)$$

In what follows, the following notations are used:

${\sigma }_r,{\sigma }_Z$ - Radial and axial stresses;

${\text{ε}}_r,{\text{ε}}_Z$ - Radial and axial deformations;

u, $\omega$ - Radial and axial movements;

${\tau }_{zr}=\tau$ - Tangential stresses;

${\gamma }_{Zr}=\gamma$ - Shear strain;

$\theta$ - Bulk deformation.

$$\frac{d({\sigma }_{r}r)}{dr}+\frac{d(r{\tau }_{rz})}{dZ}+{\sigma }_{\theta }=0(3a)$$

$$\frac{d({\sigma }_{Z}r)}{dZ}+\frac{d({\tau }_{zr}r)}{dr}=0(3b)$$

This is done to make it easier to approximate the system being solved on the Z axis. Indeed, the system (3a,3b) can be correctly solved only under the condition that:

$$\underset{r\to 0}{\lim }r{\sigma }_r=0,\underset{r\to 0}{\lim }r{\sigma }_Z=0(4)$$

These limiting relations are also established when finite-difference formulas are realized. To solve the problem, the finite volume of the half-space in the form of a cylinder of radius R and height H are selected. On the upper base of the cylinder, the boundary conditions are given in the displacements (Figure 1), which are: a) Within the contact area u = u₀, $\omega ={\omega }_0$, where (u0, ${\omega }_0$) is the introduction of an indenter into a half-space at specified distances from the Z axis; b) Outside the contact area it is free, that is, at the corresponding nodes the system (Equation 4) is solved.

In addition, in connection with the lack of friction, the tangential stresses on the surface are zero $\tau =0$. On the lateral face of the cylinder r = R and on the lower base Z = H the boundary conditions are given in the stresses, which are calculated from the formulas of the theory of elasticity, namely, on the lateral surface:

$${\sigma }_r={\sigma }_r^{y\Pi p}(R,Z);\tau ={\tau }^{y\Pi p}(R,Z)(5a)$$

at the bottom:

$${\tau }_r={\tau }_r(r,H);{\sigma }_z={\sigma }_z^{y\Pi p}(r,H),(5b)$$

Where ${\sigma }_r^{y\Pi p},{\sigma }_z^{y\Pi p},{\tau }^{y\Pi p}$ are calculated from the formulas for solving the problem of indentation by a spherical indenter into an elastic half-space. In this case, the dimensions of the cylinder must be chosen so that the area of plastic deformation is inside the cylinder. On the axis of symmetry r = 0 as boundary values we take the following conditions: u_r = 0, since on the axis of symmetry the radial displacement is zero; ${\tau }_{zr}=0$ similarly, by virtue of symmetry.

To describe the plasticity processes, we use the equations of the theory of small elastoplastic deformations and the method of elastic solutions, which, after some transformations, can be written in the form of Hooke's law:

$${\text{ε}}_r=\frac{1}{E^{*}}\left[{\sigma }_r-V^{*}\left({\sigma }_{\theta }+{\sigma }_r\right)\right];(6a)$$

$${\text{ε}}_z=\frac{1}{E^{*}}\left[{\sigma }_{\theta }-V^{*}\left({\sigma }_r+{\sigma }_z\right)\right];(6b)$$

$${\text{ε}}_{\theta }=0,\gamma =\frac{1}{{\sigma }^{*}}\tau ,(6c)$$

Where

$$E^{*}=\frac{\frac{{\sigma }_u}{{\text{ε}}_u}}{1+\frac{3-2V}{3E}.\frac{{\sigma }_u}{{\text{ε}}_u}};{\sigma }^{*}=\frac{{\sigma }_u}{3{\text{ε}}_u};(7a)$$

$$V^{*}=\left(\frac{1}{2}-\frac{1-2V}{3E}.\frac{{\sigma }_u}{{\text{ε}}_u}\right)/\left(1+\frac{1-2V}{3E}.\frac{{\sigma }_u}{{\text{ε}}_u}\right),(7b)$$

Here ${\sigma }_u$ - Stress intensity;

${\text{ε}}_u$ - Strain intensity.

Then the solution of the problem in the theory of plasticity is reduced into solving a problem of the contact theory of elasticity with variable elasticity parameters, determined by the formula described in Equations (6) and (7), and the relationship between the elasticity parameters:

$${\sigma }^{*}=E^{*}/2(1+V^{*})(8)$$

The solution of contact tasks is organized as follows:

-The initial introduction of the ball into the half-space is specified and, with the given boundary conditions, the method of variable elasticity parameters determines the stresses and deformations in the material; the energy of the deformed half-space is then determined from the formulas:

$$A={\displaystyle {\int }_0^{H}2\pi }{\displaystyle {\int }_0^{R}rWdrdz,}(9)$$

Where W - Energy density.

$$W=W_o+W_{\Phi },(10)$$

Where $W_o$ - Bulk deformation energy; and

$W_{\Phi }$ - Forming energy;

$$W_o=\frac{1-2V^{*}}{3E^{*}}{\left({\sigma }_r+{\sigma }_Z+{\sigma }_{\theta }\right)}^2;(11)$$

$$W_{\Phi }={\displaystyle {\int }_0^{\text{ε}}{\sigma }_u}{\text{ε}}_{u}d{\text{ε}}_u(12)$$

The verification of the fulfillment of the condition is by taking into account the assumption that the ball is perfectly rigid and Aш = 0. If the condition does not hold, then we continue to insert the ball into the half-space, until the condition is met.

The method of variable elasticity parameters is as follows. In the first approximation, solve the usual problem of the theory of elasticity, when the variable elasticity parameters are constant E^* = E, V^* = V since the system becomes ordinary Hooke's law. At the first stage, we use the analytical formulas. From the obtained value, we find the deformations ${\text{ε}}_z^{(1)},{\text{ε}}_r^{(1)},{\text{ε}}_{\theta }^{(1)}$ and the stresses ${\sigma }_z^{(1)},{\sigma }_r^{(1)},{\sigma }_{\theta }^{(1)}$. According to the last values, at each point we determine the stress intensity ${\sigma }_u^{(1)}$ and the strain intensity ${\text{ε}}_u^{(1)}$. From the deformation curve we find the stress intensity ${\sigma }_u^{{(1)}^{*}}$, which corresponds to the calculated value ${\text{ε}}_u^{(1)}$, then we set ${\sigma }_1^{*}={\sigma }_u^{{(1)}^{*}}/3{\text{ε}}_u^{(1)}$, and find the modules $E_1^{*}\text{ И }{V}_1^{*}$, according to Equations (7a, 7b).

In the second step we solve the problem of the theory of elasticity with the obtained elasticity parameters, we determine in the second approximation the displacement u⁽²⁾ and ω⁽²⁾, then ${\text{ε}}_Z^{(2)},{\text{ε}}_r^{(2)},{\text{ε}}_{\theta }^{(2)}\text{ И }{\sigma }_Z^{(2)},{\sigma }_r^{(2)},{\sigma }_{\theta }^{(2)}$, from them we find the intensities ${\sigma }_u^{(2)},{\text{ε}}_u^{(2)}$ at each point of space, we calculate ${\sigma }_u^{{(2)}^{*}}$ by the deformation curve, we assume ${\sigma }_2^{*}={\sigma }_u^{{(2)}^{*}}/3{\text{ε}}_u^{(2)}$ etc.

The calculations are continued until the obtained results of the approximation calculations are different from the results (n - 1) of approximations by a given amount with the required accuracy. As the main criterion in the program, the condition, $\left|{\sigma }_n^{*}-{\sigma }_{n-1}^{*}\right|<\text{ ε}$, which means virtually invariance of the elastic parameters, thus process converges.

The impact of a ball on a surface is a complex process, for the description of which it is necessary to use the equations of thermoelasticity, plasticity, impact theory, hydrodynamics. In addition, the presence of a complex surface profile, which is formed due to roughness, makes this task difficult to resolve. At the same time, as the analysis has shown, it is possible to identify the determining equations on which the behavior of the model depends and discard the remaining nonessential bonds. Since the speed of flying ball in the process of hydrobasting is not high, then we neglect the dynamic effects, and assume that the shock is quasistatic. It is quasistatic if:

• The deformations are considered to be concentrated in the vicinity of the contact area and are determined by the static theory, the wave motion in the bodies is neglected;

• Each body moves at any time with the velocity of its center of mass. The quasistatic conditions remain valid also in the case of plastic deformations, since the presence of plastic flow reduces the intensity of the contact pressure and, consequently, the energy going to the elastic wave motion. In the shot peening, the impact speed is known to be up to 70 m/s, it is possible to use the relations for inelastic contact stresses under static conditions with the yield stress is replaced by a dynamic yield strength. For this reason, we neglect the influence of thermal stresses, since even assuming that the entire energy of the ball is spent on heating, it still does not suffice to exert a significant influence on the distribution of residual stresses. In practice, we assume that the surface is perfectly smooth, and the expenditure of energy expended on the deformation of the scallops will be taken into account in formula (3a, 3b) by introducing corresponding corrections in the coefficient K.

Meanwhile, it is known from experiments and numerical calculations that the stress intensity is maximal on the axis of symmetry of the imprint and gradually decreases, tending to zero with increasing distance from the axis of symmetry of the print. With this in mind, for stress intensity, we can write the expression of:

$${\sigma }_u=\beta (Z).e^{-\text{ε}{r}^2}(13)$$

Applying similar arguments for residual stresses and taking into account that they essentially depend on the yield strength of the material, we obtain the expression:

$${\sigma }_r^{\text{OCT}}=\theta {}_r(Z)(1-{\mu }_r{r}^2)/k_r{\sigma }_T;(14a)$$

$${\sigma }_{\theta }^{\text{OCT}}=\theta {}_{\theta }(Z)(1-{\mu }_{\theta }{r}^2)/k_{\theta }{\sigma }_T,(14b)$$

Where ${\sigma }_T$ - Yield strength;

$k_r$ and $k_{\theta }$ - Co-efficients of sensitivity of residual stresses to the yield strength of the material.

The functions $\beta (Z),\theta {}_r(Z),\theta {}_{\theta }(Z)$ are arbitrary, but based on physical meaning it is necessary to demand that:

$$\underset{r\to \infty }{\lim }\beta (Z)=0;(15a)$$

$$\underset{r\to \infty }{\lim }\theta {}_r(Z)=0;(15b)$$

$$\underset{r\to \infty }{\lim }\theta {}_{\theta }(Z)=0(15c)$$

Figure 2 schematically shows the distribution of intensity of load and residual stresses in the area of prints, which are calculated from formulas (13) and (14). These formulas reflect the qualitative picture of the distribution of residual stresses.

We expand the function (13) into the Taylor series.

$${\sigma }_u=\beta (Z)\left[1+{\displaystyle \sum {(-1)}^2\frac{1}{n!}}{\text{ε}}^n{r}^{2n}\right](16)$$

If the yield stress is subject to the condition

$${\sigma }_T\le {\sigma }_u^{\max },(17)$$

Then in Equation (16) we can drop all terms except the first with a small error, as a result we obtain:

$${\sigma }_u=\beta (Z)\left[1-\text{ ε}{r}^2\right](18)$$

Using the expression in Equation (18) we find the current radius of the hardening area:

$$\overline{{\sigma }_T}=\beta (Z)\left[1-\text{ ε}{l}^2\right],(19)$$

Hence

$$l^2=\frac{1}{\text{ε}}(1-\frac{\overline{{\sigma }_T}}{\beta (Z)})(20)$$

Now, taking into account Equation (18), we obtain:

$$F(Z,l)={\displaystyle {\int }_0^l\alpha \beta }(Z)r(1-\text{ε}{r}^2)dr=\alpha \beta (Z)l^2(\frac{1}{2}-\frac{\text{ε}{l}^2}{4}),(21)$$

While taking into account Equation (19):

$$F(Z,l)=\alpha \beta (Z)\frac{1}{4\text{ε}}\left(1-\frac{\overline{{\sigma }_T}}{\beta (Z)}\right)\left(1+\frac{\overline{{\sigma }_T}}{\beta (Z)}\right)(22)$$

Further, we take the integrals for $F_r$ and $F_{\theta }$

$$F_r=2\pi {\displaystyle {\int }_0^{R}r}{\theta }_r(Z)(1-{\mu }_r{r}^2)/k_r\overline{{\sigma }_T}(Z)dr(23)$$

Obviously, in order to take this integral, it is necessary to specify the limits of integration. The radius of integration R will be found from the expression:

$$\theta {}_r(Z)(1-{\mu }_r{r}^2)/k_r\overline{{\sigma }_T}=0,(24)$$

Hence:

$$R=\sqrt{\frac{1}{{\mu }_r}}(25)$$

Integrating (23) while taking into account (25), we obtain:

$$F_r(Z,l)=\frac{\theta {}_r(Z)}{4{\mu }_r{k}_r\overline{{\sigma }_T}(Z)}(26)$$

Obviously, by carrying out similar actions for $F_{\theta }$ we obtain:

$$F_{\theta }(Z,l)=\frac{\theta {}_{\theta }(Z)}{4{\mu }_{\theta }{k}_{\theta }\overline{{\sigma }_T}(Z)}(27)$$

We now substitute the obtained expression in Equations (22), (26), (27) into Equation (3).

$$\frac{\overline{d{\sigma }_T^{\Sigma }}(z,t)}{dt}=2\pi q\left[\alpha \beta (Z)\frac{1}{4\text{ε}}(1-\frac{\overline{{\sigma }_T^2}}{{\beta }^2(Z)}-\frac{1}{2\text{ε}}(1-\frac{\overline{{\sigma }_T}}{\beta (Z)}))\right];(28a)$$

$$\frac{\overline{d{\sigma }_X^{\Sigma \text{ OCT}}}}{dt}=\frac{2q}{\pi }\left[\frac{\theta {}_r(Z)}{4{\mu }_r{k}_r\overline{{\sigma }_T}}+\frac{\theta {}_{\theta }(Z)}{4{\mu }_{\theta }{k}_{\theta }\overline{{\sigma }_T}}\right](28b)$$

We now consider in more detail the first differential equation, since it does not depend on ${\sigma }_{(Z,t)}^{\Sigma \text{ OCT}}$, then it can be solved separately. For the convenience of the solution, we introduce the following coefficients:

$$k=2\pi q;a=\frac{k(r-\alpha )}{4\text{ε}\beta (Z)};b=\frac{k}{2\text{ε}};c=\frac{\alpha k\beta (Z)}{4\text{ε}}(29)$$

Then the first Equation (28) is transformed as follows:

$$\frac{d{\overline{\sigma }}_T(z,l)}{dt}=a\overline{{\sigma }_T^2}(z,t)-b{\overline{\sigma }}_T(z,t)+c(30)$$

We shall seek a solution of (30) in the form of:

$${\overline{\sigma }}_T=y_1+\frac{1}{y^{\prime }}(31)$$

Where $y_1$ - Particular solution.

We assume that $y_1$ is a certain value that does not depend on the parameter t.

$$y_1=Y(32)$$

Then substituting Equation (32) into Equation (31) we obtain:

$$a{Y}^2+bY+c=0(33)$$

Solving this equation for Y we obtain two roots:

$$Y_{1,2}=\frac{\beta (Z)}{2-a}(1\pm \sqrt{1-a(2-a)}(34)$$

Now (30) it can be reduced to a linear differential equation:

$$\frac{d\overline{y}}{dt}+(2aY-b)\overline{y}=-a(35)$$

Substituting here (34) and taking into account (28), (29), we obtain the equation:

$$\frac{d\overline{y}}{dt}+\frac{k}{2\text{ε}}(A_0-1)\overline{y}+\frac{k(r-a)}{4\beta (Z)\text{ε}}=0,(36)$$

Where

$$A_0=1\pm \sqrt{1-a(2-a)}(37)$$

For simplicity, suppose that the material is ideally hardened, then $\alpha =1$ and (37) is transformed as follows:

$$\frac{d\overline{y}}{dt}+\frac{k}{4\beta (Z)\text{ε}}=0(38)$$

Integrating (38), we obtain the expression:

$$\overline{y}=-(\frac{kt}{4\beta (Z)\text{ε}}+c_1)(39)$$

Then (31) takes the following form:

$${\sigma }_T=\beta (Z)-\frac{1}{\frac{kt}{4\beta (Z)\text{ε}}+c_1}(40)$$

The constant $c_1$ is found from the boundary condition:

$${{\overline{\sigma }}_T|}^{t=0}={\sigma }_T^{\text{ИCX}};(41a)$$

$${\sigma }_T(Z,t)=\beta (Z)-\frac{1}{c_1};(41b)$$

Hence $c_1$ is found:

$$c_1=\frac{1}{-{\sigma }_T^{\text{ИCX}}+\beta (Z)}(42)$$

Substituting (42) into (41), we obtain:

$${\sigma }_T=\beta (Z)-\frac{1}{\frac{kt}{4\beta (Z)\text{ε}}+\frac{1}{\beta (Z)-{\sigma }_T^{\text{ИCX}}}}(43)$$

Study the solution obtained by us. If the hardening time increases, then: ${\sigma }_T(Z,t)\to \beta (Z)$ in the limit:

$$\underset{t\to \infty }{\lim }{\sigma }_T(Z,t)=\beta (Z)(44)$$

Now consider the relationship between residual stresses and deflection of control plates $\Delta f$.

For simplicity, we take in (30) that ${\mu }_r={\mu }_{\theta },k_r=k_{\theta },$ then we obtain:

$$\frac{d\overline{{\sigma }_X^{\Sigma \text{ OCT}}}(Z,t)}{dt}=\frac{2q}{\pi }\left(\frac{{\theta }_r(Z)-{\theta }_{\theta }(Z)}{4\mu k{\sigma }_T(Z,t)}\right)(45)$$

We substitute into the equation (45):

$$d=\frac{2q}{\pi }(\frac{{\theta }_r(Z)-{\theta }_{\theta }(Z)}{4\mu k}),(46)$$

Then the equation takes the form:

$$\frac{d\overline{{\sigma }_X^{\Sigma \text{ OCT}}}(Z,t)}{dt}=\left(\frac{d}{{\sigma }_T(Z,t)}\right)(47)$$

Equation (40) is reduced to the form:

$${\sigma }_T=\frac{k\beta (Z)(\beta (Z)-{\sigma }_T^{\text{ИCX}})t+4\beta (Z)\text{ε}{\sigma }_T^{\text{ИCX}}}{kt(\beta (Z)-{\sigma }_T^{\text{ИCX}})+4\beta (Z)\text{ε}}(48)$$

Taking into account (48), the equation (47) is becomes:

$$\frac{d\overline{{\sigma }_X^{\Sigma \text{OCT}}(Z,t)}}{dt}=\frac{d[kt(\beta (Z)-{\sigma }_T^{\text{ИCX}})+4\beta (Z)\text{ε}]}{kt(\beta (Z)-{\sigma }_T^{\text{ИCX}})+4\beta (Z){\sigma }_T^{\text{ИCX}}\text{ε}}(49)$$

We make the substitution in the Equation (49):

$$a=k(\beta (Z)-{\sigma }_T^{\text{ИCX}});b=4\beta (Z)\text{ε};(50a)$$

$$f=k(\beta (Z).(\beta (Z)-{\sigma }_T^{\text{ИCX}});g=4\beta (Z){\sigma }_T^{\text{ИCX}}\text{ε};(50b)$$

$$\frac{d\overline{{\sigma }_X^{\Sigma \text{OCT}}}}{dt}=d\frac{at+b}{ft+g}(50c)$$

We shall carry out a qualitative analysis of Equation (50). At the initial time, the increment of residual stresses should be maximal, since there is no hardening, in the future it should decrease to a certain limit, which can be easily found from the relation in Equation (50) (Figure 3).

$$\underset{t\to \infty }{\lim }\frac{d\sigma (Z,t)}{dt}=d\frac{a}{f}(51)$$

The presence of this limit is easily explained by the existence of a maximum hardening limit:

$${\sigma }_T<{\sigma }_T^{\max },0<t<\infty (52)$$

This limit exists in many materials and is explained by its physical properties. We integrate now to equate to (50).

$$\overline{{\sigma }_X^{\Sigma \text{OCT}}}=d\left[\frac{at}{f}+\frac{fb-ag}{f^2}(\ln (ft+g)-\ln g)\right](53)$$

The analysis shows that the calculated curve can be divided into 3 phases:

I. Fast initial growth phase;

II. Transitional period;

III. Saturation phase.

Since the solution of the adaptability equations at each i-th step of the integration requires the solution of the one-shot problem, it is most advantageous to apply the interactive solution methods, since the resolution of the equations at each new step begins with the previous solution, which increases the computational speed. In particular, here we have chosen a two-layer integration scheme for solving the problem. The main decision points listed below.

1) Based on the preset depth of ball penetration, determine the contact radius:

$$B=\sqrt{R_{\text{Ш}}.D_e,}(54)$$

Where $R_{\text{Ш}}$ - Ball radius;

$D_e$ - Depth of implementation.

2) Determine the number of nodes N1 in the contact area and check whether it falls into the required range, if not, the calculation is terminated. This is done to ensure that the number of grid nodes does not exceed the specified range, and, on the other hand, that the accuracy of the solution is sufficient.

3) Determine the coordinates of the grid nodes X (1), Y (1).

4) Determine the force of the ball pressing:

$$P=\frac{8{\sigma }_{0}B{D}_e}{3(1-v{)}^{\prime }}(55)$$

Where ${\sigma }_0$ - Initial shear modulus.

5) By the given force P and the radius of the print B determine the elastic stresses according to:

$${\sigma }_r^{\text{yΠp}}={\sigma }_r^{\text{yΠp}}(r,Z);{\sigma }_{\text{Z}}^{\text{yΠp}}={\sigma }_Z^{\text{yΠp}}(r,Z);{\sigma }_{\theta }^{\text{yΠp}}={\sigma }_{\theta }^{\text{yΠp}}(r,Z);(56)$$

6) On the basis of (55) find the boundary conditions:

On the axis $r=0$ we obtain:

$${\sigma }_r={\sigma }_{\theta }=K\left[(1+V)(1-Zarctg\frac{1}{Z}-\frac{1}{2(1+Z^2)})\right];(57a)$$

$${\sigma }_Z=K\frac{1}{1+Z^2};\tau =0(57b)$$

On the surface $Z=0$ we have $\tau =0$

$${\sigma }_r|r\le b=K\left[\sqrt{1-r^2}+\frac{1-2v}{3r^2}\left({\left(1-r^2\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-1\right)\right];(58a)$$

$${\sigma }_{\theta }|r\le b=K\left[\sqrt{1-r^2}+\frac{1-2v}{3r^2}\left({\left(1-r^2\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}-1\right)\right];(58b)$$

$${\sigma }_Z|r\le b=K\sqrt{1-r^2};(58c)$$

$${\sigma }_{\theta }|r>b=-{\sigma }_r,{\sigma }_Z|r>b=0(58d)$$

7) Find the boundary conditions on the lateral face:

$$F_{Gb}(b,Z)={\sigma }_r^{\text{yΠp}}(b,Z);(59a)$$

$$F_{\tau b}(b,Z)={\tau }^{\text{yΠp}}(b,Z);\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}(59b)$$

$$F_{GN}(r,H)={\sigma }_r^{\text{yΠp}}(r,H);(59c)$$

$$F_{\tau N}(r,H)={\tau }^{y\text{Πp}}(r,H),(59d)$$

Where b, H - Height and width of the calculation area.

8) On the basis of (58), find the elastic displacements at the nodes.

9) To solve the interatomic problem, take the calculated elastic stresses as initial conditions, in addition, assume at the boundary that the boundary conditions correspond to the conditions found from the elastic solution (58) and (59).

10) Determine the intensity of stresses taking into account the residual stresses ${\sigma }_X^{\overline{\text{OCT}}}$ :

$${\sigma }_u=\frac{1}{\sqrt{2}}\sqrt{{({\sigma }_r+{\sigma }_X^{\overline{\text{OCT}}}-{\sigma }_Z)}^2+{({\sigma }_{\theta }+{\sigma }_X^{\overline{\text{OCT}}}-{\sigma }_Z)}^2+6{\tau }^2}\text{ (60)}$$

To calculate the stress state in the plastic region, we use the variable elasticity method. According to this, the shear modulus σ and Poisson's ratio V at each point of the area are calculated (Figure 4). By σ^* and V^* their values depending on the value σ_u are indicated. σ and V are the values from the previous iteration. Find the intensity of the deformations by the formula:

$${\text{ε}}_u=\frac{{\sigma }_u}{3}\left(\frac{1}{\sigma }+\frac{1-2V}{E}\right)\text{ (61)}$$

Further, we find the intensity along the deformation curve:

$${\sigma }_u{}^{*}=a_n+b_n{\text{ε}}_u,\text{ (62)}$$

Where n - Number of the segment of the broken line, which approximates the curve, the deformation. After this, we determine the new value of the shear modulus by the formula.

$${\sigma }^{*}=\frac{1}{\frac{3{\text{ε}}_u}{{\sigma }_u^{*}}+\frac{1-2v}{E}}\text{ (63)}$$

11) Find the deformations from the displacements, then calculate the stresses at the nodes from the obtained elastic deformations and the recalculated elastic modulus.

$${\sigma }_r=\lambda \theta +2\sigma {\text{ε}}_r;\text{ (64a)}$$

$${\sigma }_{\theta }=\lambda \theta +2\sigma {\text{ε}}_{\theta };\text{ (64b)}$$

$${\sigma }_Z=\lambda \theta +2\sigma {\text{ε}}_Z;\text{ (64c)}$$

12) To determine the displacements at the grid nodes, it is necessary to compile a system of linear equations with fixed elastic parameters:

$$AY=F,\text{ (65)}$$

Where Y - Array of radial and axial displacements;

F - Boundary conditions, which are determined from the equations (59)

a) In order to find the matrix A, it is necessary to express the equilibrium equations in the form of a displacement function. Substituting (64) into the equilibrium equations, we obtain a system of second-order differential equations:

$${\Delta }^{2}u+\frac{1}{1-2V}.\frac{d}{dx}\left(\frac{du}{dx}+\frac{dV}{dy}+\frac{dw}{dz}\right)+\frac{F_x}{\mu }=0;\text{ (66a)}$$

$${\Delta }^{2}V+\frac{1}{1-2V}.\frac{d}{dy}\left(\frac{du}{dx}+\frac{dV}{dy}+\frac{dw}{dz}\right)+\frac{F_y}{\mu }=0;\text{ (66b)}$$

$${\Delta }^{2}w+\frac{1}{1-2V}.\frac{d}{dz}\left(\frac{du}{dx}+\frac{dV}{dy}+\frac{dw}{dz}\right)+\frac{F_z}{\mu }=0\text{ (66c)}$$

The derivatives in (66) are found by three points:

$$y^{\prime }(x_1)=\frac{1}{H}\left(-y_3+4y_2-3y_1\right);\text{ (67a)}$$

$$y^{\prime }(x_i)=\frac{1}{H}\left(y_{i+1}-y_{i-1}\right),1<i<n;\text{ (67b)}$$

$$y^{\prime }(x_n)=\frac{1}{H}\left(3y_n-4y_{n-1}+y_{n-2}\right),\text{ (67c)}$$

Here $$y_i=y(x_i),$$ H - Double grid spacing.

Substituting the expression (67) into the equilibrium equation (66) and representing it as the product of matrices (65), we define the matrix A.

b) To solve the equation (65), we use the three-layer iterative method. According to this method, it is necessary to calculate the remainder by the formula:

$$r_k=A{Y}_k-F\text{ (68)}$$

c) Calculate the vector column. $A_{rk}$ for this we use the equations (65) and (66), but instead of displacement we substitute the remainder.

d) Calculate the coefficients A, Е, F and if $K>1$ , then calculate the coefficients B, C, D:

$$A=(A_{rk},r_k);B=(A_{rk},r_{k-1});\text{ (69a)}$$

$$C=(r_k,r_k-r_{k-1});D=(r_{k-1},r_k-r_{k-1});\text{ (69b)}$$

$$E=(A_{rk},A_{rk});F=(r_k,r_k);F_{st}=(r_{k-1},r_{k-1})\text{ (69c)}$$

e) Calculate the new vector of the solution Y at the -th integration by the formulas:

$$Y^{(1)}=Y^{(0)}-\tau .r_0,\text{ ΠpИ }k=1\text{ (70a)}$$

$$Y^{(k+1)}=Y^{(k-1)}+\alpha \left(Y^{(k)}-Y^{(k-1)}\right)-\alpha \tau r_k,\text{ ΠpИ }k\ge 2,\text{ (70b)}$$

Where

$${\alpha }_k=\frac{\left(A-B\right)B-DE}{(C-D)E-{\left(A-B\right)}^2};{\alpha }_1=1;\text{ (71a)}$$

$${\tau }_1=\frac{A}{E};{\tau }_k=\frac{B}{E{\alpha }_k}+\frac{\left(A-B\right)}{E},k\ge 2\text{ (71b)}$$

Since the contact task is being solved, it is necessary to fix the movement of the nodes in the contact area. In order for their values to remain unchanged in the iteration process, the discrepancy in these nodes is forcibly equated to zero, which leads to automatic fixation of displacements in these nodes.

The iterative process proceeds until the residual is sufficiently small, namely:

$$\left|\frac{F}{F_{st}}-1\right|<\text{ ε},\text{ (72a)}$$

Where $\text{ε}$ - Relative accuracy.

13) Calculate the new elastic modulus from formulas (61), (62), (63). Find the change in the elastic modulus at the (k + 1) iteration.

Further we check the condition: $\Delta E<\Delta E_g,$ where $\Delta E$ - Maximum increment of the shear modulus; $\Delta E_g$ - Maximum increment error.

If the condition is not fulfilled, go to (10) and repeat the calculation again, otherwise go to the next step. This cycle limits the solution of the elastoplastic problem for a given depth of ball penetration.

14) Calculate the energy $∍$ of the half-space;

$$∍={\displaystyle {\int }_V({∍}_{\text{o}}+{∍}_{\Phi })}dxdydz;\text{ (73a)}$$

$${∍}_{\text{o}}=\frac{1-2V^{*}}{6E^{*}}{({\sigma }_r+{\sigma }_{\theta }+{\sigma }_z)}^2;\text{ (73b)}$$

$${∍}_{\Phi }=\frac{{\sigma }_u{\text{ε}}_u}{2},\text{ (73c)}$$

Where ${∍}_{\text{o}}$ - Bulk strain energy;

${∍}_{\Phi }$ - Forming energy.

Check whether the kinetic energy of the ball $A_{\text{o}}$ is equal to the energy of the half-space. If not, go to step 17, otherwise go to the next step.

16) Define a new depth of the introduction of balls into a half-space, while:

$$h=h+\Delta h,if∍<A_o;\text{ (74a)}$$

$$h=h-\Delta h,if∍>A_o,\text{ (74b)}$$

Where h - Penetration depth of the ball;

$\Delta h$ - Increment of penetration depth.

To achieve the convergence, the step is calculated using the following the formula:

if $∍(h+\Delta h)>A_{\text{o}}$ and $∍(h)<A_o$ , then $\Delta h=\raisebox{1ex}{$\Delta h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$    (75)

After recalculating the depth of implementation, go to step 1 of the algorithm.

17) Calculate the residual stresses by the formula:

$${\sigma }_{ij}^{\text{OCT}}={\sigma }_{ij}^H-{\sigma }_{ij}^{y\prod \text{p}},\text{ (76)}$$

Where ${\sigma }_{ij}^H$ - Load stresses;

${\sigma }_{ij}^{y\prod \text{p}}$ - Stresses, discharges, obtained from the elastic solution.

The solutions obtained from the algorithm are based on the fact that at a sufficiently large distance from the contact zone the elastoplastic solution converges to the elastic solution. This assumption allows to limit the calculated area to a cylinder of the radius B and the height H, and at the nodes located on the outer boundary of the cylinder, the displacements obtained from the elastic solution.

Figure 5 shows the distribution of stresses with ${\text{K}}_{\text{ΠpoΠ}}\text{ = 1,}$ the ball radius $R_{\text{Ш}}=1,25mm.$ The figure shows the example of the residual stress distribution with a dummy mechanical property of the steel materials. For our experiment, the comparable stress would be the right bottom stress distribution.

Finite Element Modeling

All computer modeling in the present study was conducted using the nonlinear finite element (FE) code LS-DYNA 971. A model of the workpiece surfaces was developed using a constant stress solid element formulation which suitably model the impacted surface response of the workpiece. The impacting balls of two different sizes namely 0.8 mm and 1.3 mm were also modelled using solid elements with reduced integration techniques in combination with hourglass control. In a way it is similar to that of Purohit’s [5] with multiple shot balls. The workpiece is modelled with four geometrical parts with same nodes between the internal surfaces. The model was discretised using solid element with two mesh sizes to capture better simulation results. Stiffness-based hourglass control was employed to avoid spurious zero energy deformation modes and the use of reduced integration was to avoid volumetric locking. The elastic and elasto-plastic material models were used to model the impacting ball and the impacted surfaces, respectively. The material properties included for the impacted surface are the yield stress of 1433 MPa with a plastic hardening behavior. For the impacting ball, a very high value of stiffness has been assigned to represent high rigidity of the structure when contacting with other surfaces. The finite element discretization of the shot peening simulation is depicted in Figure 6. Shot peening loading was simulated by defining an impact velocity of 90 m/s to the impacting balls over a three repeated impact loading. Figure 7 shows the schematic diagram of the impacting points and sequence on the workpiece model. The evaluation of residual stress and impacting depth was carried out underneath of Point 5. The simulation was done in two separate phases namely explicit and implicit simulations. Explicate simulation was employed to simulate the impact loading on the workpiece. From this simulation, impact analysis may be carried out resulting in obtaining the initial stress field in the workpiece. Subsequently, the output obtained from explicit simulation was employed to simulate implicit simulation. Thereafter, the residual stress in material can be obtained from the implicit analysis on work piece stress after impact. The impacting ball was constrained to translate vertically only along the z axis, and the bottom and edges nodes of the workpiece surfaces were fully constrained as shown in Figure 7. A reliable contact modeling for the interaction between the impacting ball and workpiece is paramount and needs to be accurately established. Automatic surface to surface contact was defined for the contact interfaces of the impacting ball with the workpiece in order to avoid any penetration occurred. The friction coefficients associated with contact between the balls and workpiece were defined as 0.3 and 0.2 for the static and dynamic respectively.

Figure 8 and Figure 9 show the result of the FE models. One can see that as the number of the pass increases, the maximum compressive residual stress also increases. At the same time, the depth of the peak also increases as the size of the ball increases. These were in agreement with both our analysis and experiment.

Experiment

The experiment was done by shot peening on two different plates, namely plate A and plate B. On the plate A, shot ball size was 0.8 mm and on plate B, the shot ball size 1.3 mm. The experiment was performed with the same machine with the parameters to match the analysis. The coverage of the shot peening was approximately 98%. This is to mimic what is usually called saturated shot peening. Beyond this coverage, more shot peening amount that is put into the system will not affect the residual stress significantly. The samples were than profiled with an x-ray diffraction device. Figure 10 and Figure 11 show the results. The 0.8 mm ball produces slightly less than 300 MPa at the surface, while the 1.3 mm ball produces more than 500 MPa. The peak for the shot ball size of 0.8 mm produces 0.1 mm maximum depth of compressive residual stress, while the 1.3 mm ball produces peak at 0.19 mm. This trend is similar with both the analytical and finite element results as well as the finding of other researchers [6].

Discussions

Figure 12 and Figure 13 are the direct comparison between the simulation and the experiment. One can see that while the direct comparison do not match in detail, but they are at the right tendency. The disagreement between the analysis, finite element model and the experiment was expected from the beginning. The computation limitation was also prevented us to do until saturation, in our case, 3 passes. However, since this is our first attempt to academically materialize the concept of the shot peening in calculation such discrepancy is highly anticipated. Other researchers also pursue similarly [2,7-9]. Among the practitioners, the saturated residual stress at the surface is usually predicted as:

$${\text{Surface σ}}_{RS}\text{[MPa] }=-276\frac{Ah}{R}+7.1{\gamma }_R-\text{ 0}.59{\sigma }_{pre}-451\text{ (76)}$$

At the surface, where Ah is the arc height and R is the radius of shot, while ${\gamma }_R$ is the retain austenite. For the peak, usually it becomes

$$\text{Max}{\text{. σ}}_{RS}\text{[MPa] }=-172\frac{Ah}{R}+7.1{\gamma }_R-\text{ 0}.54{\sigma }_{pre}-882\text{ (77)}$$

The unit is in MPa. The above empirical equations are quite famous in the industrial world although it is lack of scientific basis. Based on Equations (76) and (77), and our measurement that arc-height are 0.54 and 0.28 respectively, our prediction of surface residual stress would be -530 MPa and -512 MPa respectively. The peak residual stress would be 917 MPa and 906 MPa. Again, direct comparison of the prediction using these empirical equations with experiment and with computational results are still difficult. We can only say that the trend is similar.

Figure 14 shows the coverages at different passes. It shows clearly that our first attempt in this computation is still an early stage and can be modified into broader area and more variation of the coverage.

In the fatigue based industrial application, it is very critical that a product has both sufficient surface residual stress and depth of peak to prolong the fatigue life. This trend explain why some products use double shot peenings, which are the first using large shot ball and followed by smaller shot ball. This is thought to guaranty delaying both the fatigue crack initiation and early stage of fatigue crack propagation.

Conclusions

Through analysis, finite element modeling and experiment, it is proved that shot ball size influences mainly on the depth of the peak residual stress. This conclusion is nothing new among the practitioners. However, this research proved that academically this field has so much promise and potential to develop further. It would be an interesting topics to incorporate the chaos theory, since while the number of individual impacts is linearly proportional to shot flow, exposure area, and exposure time, the shot peening coverage is not linearly proportional because of the random nature of the process.

Acknowledgment

Part of the work was performed during the first author’s stay as a visiting professor at the South Ural State University, Russia. The financial support and contribution from the International office of South Ural State University is highly appreciated.

---