Using the FEM Method in the Prediction of Stress and Deformation in the Processing Zone of an Elastic/Visco-Plastic Material during Diamond Sliding Burnishing

Abstract

This article concerns the application of the FEM method for the prediction of stress and deformation states in the workpiece during diamond sliding burnishing (DSB). An updated Lagrange (UL) description was used to describe the phenomena at a typical incremental step. The states of strain and strain rate are described by non-linear relationships without linearization. The material parameters were estimated during tensile tests to determine the characteristics of the 41Cr4 steel. Its hardness was also tested. Its aim was to prepare a table with the material properties of the above-mentioned steel and its implementation for numerical analyses. A Cowper–Symonds material model was used to model the displacement process of the wedge on an elastic/visco-plastic body reflecting the DSB process. The computer model was validated, and a good convergence of the results was obtained. Applications in the ANSYS/LS-Dyna program were developed to simulate the process of DSB. The results of numerical analyses were presented, among others, to explain the influence of the rake angle on the condition of the surface after machining, as well as the phenomenon of chip formation. The results of numerical simulations were verified experimentally on a test stand.

1. Introduction

1.1. Modeling of Selected Phenomena in the Diamond Sliding Burnishing Process

Diamond sliding burnishing is a quite well-known technological process. The main advantages include, among others, improved quality of the surface layer, improved fatigue strength or improved corrosion resistance. Numerous papers describe different advantages of using diamond burnishing as an alternative method in relation to the turning process. Below, several papers are shown as an example of the huge popularity of diamond sliding burnishing. The Authors generally research the influence of technological parameters or shape of the tool on the state of the surface layer after research. However, all the papers have the same common features. The Authors will mention these after a brief analysis of the literature.

In paper [1], the Authors present the geometric changes taking place on the surface of an object made of 16HG steel, hardened and burnished with a sliding diamond. They conducted real and simulation tests in the ABAQUS program on the basis of a two-dimensional section of the surface of a real object.

In works [2,3], it was found that the multiplicity of the load (the number of passes) during burnishing should not be too high, as it may cause excessive crushing with all its harmful consequences. Increasing the number of cycles above several is not advisable for economic reasons either, because as a result of the phenomenon of metal hardening, each cycle brings worse and worse effects [2,3].

The model of material deformation using a non-deformable rigid ball, known as the Hertz model, is quite often and with good results, adopted for an analysis of deformation processes in the contact zone of two elements, e.g., in an analysis of strengthened burnishing [4,5]. However, this model is not useful in the case of burnishing and hardening, or even in the case of burnishing itself. The formulas obtained on its basis in comparison with the values of forces used in experimental studies [6,7] indicate significant differences in the results.

The model of material deformation by means of a ball rolling on the machined surface is a better approximation of the burnishing processes, and due to its high similarity to sliding burnishing, it can be used to calculate the tool clamping force. However, attempts to calculate the tool clamping force during smooth burnishing provide inflated results; the values obtained differ significantly from the experimental data [8,9,10].

The model of material deformation by means of a ball sliding on the surface of the workpiece was considered by many authors, among which the studies described in [11] belong to the classics. The effects of burnishing, as practice shows, depend not only on the elastic-plastic properties of the workpiece and the geometry of the contact between the tool and the workpiece, but also on the initial surface roughness, and this factor is not taken into account by classical models, treating output surfaces as perfectly smooth. This is confirmed by the experimental formulas presented in the literature [12,13] that are used to select the burnishing parameters, but do not take into account the effect of the initial roughness of the surface machined. Therefore, it is necessary to take this factor into account in theoretical considerations.

For the reasons mentioned above, in an analysis of burnishing, it is necessary to adopt a different model than the one found in the literature. Such a model should take into account the stereometric state of the surface and, in particular, the height of the initial surface roughness, as the roughness of the surface after burnishing depends to a large extent on it [14,15,16,17,18,19,20].

Work hardening of the surface layer increases the resistance to the development of fatigue processes, especially in materials without very high hardness, such as aluminum alloys [21]. The hardening depth should be up to about 0.25 of the cross-sectional dimension of the part; at greater depths, no improvement in fatigue strength is observed [22].

In the work [23], it was found that anti-corrosion properties after burnishing are not the best, although in some cases sliding burnishing improves them. On the other hand, in [24,25,26,27,28,29,30] it was shown that the presence of compressive stresses may reduce the intensity of tribological wear. After slide burnishing, the depth of hardening is much smaller than in the case of the strengthening burnishing method (e.g., burnishing rolling), usually, hardening occurs within the height of unevenness—this makes it impossible to meet the anti-fatigue condition [31].

After sliding burnishing, compressive stresses are formed in the surface layer. Their maximum value depends on the burnishing parameters and the properties of the material processed. This was demonstrated in study [32], where the presence of compressive stresses in the surface layer is generally considered to have a positive effect on the operational properties, especially on the fatigue strength of machine elements as well as on tribological wear.

In works [33,34] it was shown that the structural changes after sliding burnishing are small and similar to changes after rolling pressure shot peening. In most cases, these changes are invisible even under an optical microscope.

In the work [35], it was stated that sliding burnishing increases the strength of power machine rotors. In the paper [36], a special device was designed and produced for burnishing the curvilinear contours of injection molds made of AISI P20 steel with burnishing forces reaching several hundred Newtons.

There are many studies related to testing the burnishing force of rough surfaces [25,37,38]. The Authors demonstrate how important it is to prepare the surface before diamond sliding burnishing. Papers [39,40] present the decorative qualities of burnished surface layers, which additionally gain increased fatigue strength and resistance to abrasion and corrosion.

An innovative approach was proposed in [41], where AISI 4140 steel was locally heated before using the burnishing tool. This resulted in the use of a lower burnishing force as well as a greater hardness of the surface machined, along with an improvement in the quality of the surface layer.

It has been shown in many publications that burnishing will result in favorable changes in the geometric structure of the surface through the plastic deformation of micro-unevenness and shaping the structure with flat tops with a high load-bearing capacity [32,42,43,44,45].

Diamond sliding burnishing is a very popular process. However, in every case we could observe very precise descriptions of technological parameters or tool geometry, but the beginnings of the burnishing process are discussed nowhere. Some tool geometry could cause different technological processes. It is important to explain in what time we can observe transition from one process to another. It is the main goal of this paper to describe phenomena during the sliding burnishing process transition to the turning process for two tool depths.

1.2. Modeling in Continuum Mechanics

Mathematical models describing a continuous medium usually constitute systems of integral or partial differential equations. This is due to the fact that in a continuous medium, variables are dependent on at least one spatial variable and time (if the object tested is variable in time). The mathematical model is supplemented with appropriate boundary and initial conditions. As a result, a boundary or boundary-initial problem is obtained. Only in simple problems is it possible to obtain a precise analytical solution to a given mathematical problem. Most often, the solution requires the use of approximate methods. In practice, approximate numerical methods are used most often.

Numerical methods of an approximate solution of continuous boundary-initial problems include the finite difference method MRS [46], the method of boundary elements MEB [47], the finite element method FEM [48,49,50,51,52,53], various meshless methods, such as the fuzzy particle method (smoothed particle hydrodynamics, SPH) [54], the free point method [55] and the material point method, also known as the material point method or the particle-in-cell-method [56,57].

The essence of the finite difference method is to replace derivatives with appropriate difference quotients determined on a discrete set of points. The finite difference method was initially applied to differential equations (local formulation), later the range of its applications was extended to problems formulated in variational forms (global formulation). The standard MRS formulation used regular (square, rectangular or cubic) node grids.

In advanced formulations, an arbitrary mesh of nodes is used. Difference quotients can be built on arbitrarily generated points, owing to which this MRS formulation can be classified as a gridless method [58]. It is also used in the mechanics of solids, e.g., a well-known program for solving problems in geomechanics FLAC [59] is based on the finite difference method.

In this paper, the finite element method (FEM) was used to discretize a continuous problem. The essence of FEM is a division of the continuous area under consideration into a finite number of sub-areas (finite elements) connected to each other at node points and an approximation of the solution in the area of elements using interpolation functions (shape functions) and values in nodes [60]. FEM equations are obtained by means of a formulation of an integral (global) problem using the principle of variation or the method of weighted residuals.

The method of boundary elements [61] consists in reducing the system of differential equations with given boundary conditions to the system of integral equations defined on the boundary of the area considered. Contrary to the finite element method and the finite difference method, its application does not require any discretization of the interior of the area, but only its edge. This is the main advantage of this method.

One of the disadvantages of using FEM is the need to generate an appropriate mesh of elements, which in the case of complex three-dimensional geometries is not easy for some types of elements. In an analysis of large deformations, the mesh of elements is subject to excessive distortions, which often makes it necessary to re-mesh the computational domain many times. Avoiding these problems is one of the reasons for the growing popularity of meshless methods. There are many meshless methods. The meshless finite difference method was mentioned above. Other meshless methods include the smoothed particle hydrodynamics method (SPH), the free point method, the finite point method [62] and the element-free Galerkin method [63]. The fuzzy particle method was initially developed for discrete problems in astrophysics [64]. Later, it was applied as a discretization method for continuous problems.

2. Materials and Methods

2.1. Modeling of Diamond Sliding Burnishing and Micro-Cutting Processes

This article concerns the use of the FEM method to predict stress and strain states and the geometric structure of the surface of the workpiece after the diamond sliding burnishing process. The updated Lagrange description was used to describe the phenomena at a typical incremental step. Non-linear incremental models without any linearization were used to describe physical phenomena. The results of the numerical simulation were experimentally verified on the example of selected parameters of the surface roughness profile.

Sliding burnishing, similarly to micro-cutting, can be treated as a process of moving a wedge (tool) over a workpiece. For a given tool geometry (vertical angle β and radius of cutting edge rounding r_d) and machining depth a_p, there is a limit rake angle ${\mathsf{\gamma }}_{gr}$ at which material flash is formed in front of the tool, but no chip is formed. Thus, for rake angles $\mathsf{\gamma }\le {\mathsf{\gamma }}_{gr}$, the sliding burnishing process takes place, while for $\mathsf{\gamma }>{\mathsf{\gamma }}_{gr}$, the flash turns into a chip and then the micro-cutting process takes place (Figure 1).

The processes of sliding burnishing and micro-cutting were considered as a geometrically and physically non-linear boundary-initial problem, in which there are non-linear, moving and variable in time and space boundary conditions. The boundary conditions in the tool–workpiece contact areas are unknown.

A mathematical description of non-linear phenomena requires the use of principles of formulating boundary-initial problems, which are different than those in linear problems, and require more complex methods of solving them.

Geometric non-linearity results from a non-linear relationship between strain and displacement. It is caused by a change in the geometry of the workpiece during the machining process (deformation, cracking).

Physical (material) non-linearity is caused by the non-linear mechanical properties of the workpiece material, leading to a non-linear stress–strain relationship.

2.1.1. The Concept of Incremental Description

To describe the phenomena in a typical incremental step, the updated Lagrange description was used, assuming three systems of orthogonal Cartesian coordinates ({x} movable associated with the object, {y} movable associated with the tool and {z} fixed reference system). The purpose of an incremental analysis is to determine the geometry of the body and the states of incremental displacements, displacement velocities, accelerations, strain, strain rates, stress, etc., existing in this body, in subsequent, discrete moments of time τ = 0, Δt, 2Δt, …, corresponding to a small increment of time. It is assumed that solutions are known for all the time steps in the interval [0, t] (including time t), and solutions will be searched for at time $\mathsf{\tau }=\mathrm{t}+\Delta \mathrm{t}$. By repeating this procedure for subsequent moments Δt, a solution is obtained in a given time interval [t_p, t_k].

Issues related to incremental description were presented in some works, e.g., in the paper by M. Kleiber, various concepts of incremental descriptions in non-linear mechanics were analyzed. An updated Lagrange description (Figure 2) was adopted to describe the wedge displacement process, which is characterized by high universality, based on the assumption of an additivity of individual deformation components on the original or current configuration. This concept of incremental description and algorithms of solutions developed for typical incremental step time has already been used by the Authors for numerical analyses of other technological processes of plastic forming, i.e., cutting [65,66,67,68,69,70,71,72], thread rolling [73], burnishing [74,75,76,77,78], drawing [79,80,81], grinding with a single abrasive grain [82] and for analyses of large plastic deformations and for predicting fatigue strength and product life [83,84,85], and it has been positively verified.

According to the updated Lagrange description, at times τ < t, the values of all the functions appearing in the problem are known. According to the assumption referred to above, at the current moment t, the initial (original) configuration of the body ${}_{}{}^0\mathrm{C}$ and the current configuration ${}_{}{}^{\mathrm{t}}\mathrm{C}$ are known. On the other hand, the next equilibrium configuration ${}_{}{}^{\mathrm{t}+\Delta \mathrm{t}}\mathrm{C}$ is sought, at the moment $\mathsf{\tau }=\mathrm{t}+\Delta \mathrm{t}$ (Figure 2).

The incremental description formulated in this way in a non-linear process analysis is associated with a number of significant problems. They concern the selection of appropriate coordinate systems, the definition of strain and stress measures and their increments, as well as the determination of the rules for their accumulation at each incremental step.

2.1.2. Measures of Increments of Physical Quantities

The non-linearity of the object material was described using an incremental model, taking into account the influence of the strain history and the strain rate. The object is treated as a body in which elastic deformations (in the range of reversible deformations) and viscous and plastic deformations (in the range of irreversible deformations) may occur with non-linear hardening. This body was abbreviated as E/VP. The non-linear Huber–Mises–Hencky (abbreviation: HMH) plasticity condition, associated flow law and mixed reinforcement (isotropic-kinematic) were used to build a material model. The state of the material after the preceding treatments was also taken into account by introducing the initial states: displacements, stresses and strains.

An increment ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta \mathrm{f}\left({}_{}{}^{\mathsf{\tau }}\mathrm{x},\mathsf{\tau }\right)$ of any function ${}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{f}\left({}_{}{}^{\mathrm{t}}\mathrm{x},\mathrm{t}\right)$ at a typical step t→τ = t+Δt is defined using the formalism differential, in the following form:

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta \mathrm{f}\left({}_{}{}^{\mathsf{\tau }}\mathrm{x},\mathsf{\tau }\right)\stackrel{\scriptscriptstyle\mathrm{def}}{=}{}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{f}\left({}_{}{}^{\mathsf{\tau }}\mathrm{x},\mathsf{\tau }=\mathrm{t}+\Delta \mathrm{t}\right)-{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{f}\left({}_{}{}^{\mathrm{t}}\mathrm{x},\mathrm{t}\right),$$

(1)

where ^τx is a vector (with components ^τx_i, i = 1, 2, 3) of the position of the body particle X at any time τ. The notation ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta \mathrm{f}$ indicates an increment of the function f from time t to time τ (the left superscript) expressed in the system {tx} and related to the configuration ${}_{}{}^{\mathrm{t}}\mathrm{C}$ at time t (the left subscript).

In the process analyzed of moving the wedge, the measures of increments of displacement, strain, strain rate and stress need to be defined. The movement of the body is described by smooth mapping:

$${}_{}{}^{\mathsf{\tau }}x={}_{}{}^{\mathsf{\tau }}x\left({}_{}{}^0\mathrm{x},\mathsf{\tau }\right).$$

(2)

The particle displacement vector ${}_{\mathrm{t}}{}^{\mathsf{\tau }}u$ and its components ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{i}}$, and the increment displacement ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta u$ vector and its components ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i}}$ take the following form:

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}u={}_{\mathrm{t}}{}^{\mathsf{\tau }}x-{}_{\mathrm{t}}{}^{\mathrm{t}}x, {}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{i}}={}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{x}{}_{\mathrm{i}}-{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{x}{}_{\mathrm{i}}, \mathrm{i}=1,2,3,$$

(3)

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta u={}_{\mathrm{t}}{}^{\mathsf{\tau }}u-{}_{\mathrm{t}}{}^{\mathrm{t}}u, {}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i}}={}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{i}}-{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i}}, \mathrm{i}=1,2,3,$$

(4)

where all the components of vectors are related to a common frame of reference ${}_{}{}^{\mathrm{t}}\mathrm{C}$. An adequate measure of the strain state of the workpiece, in the adopted updated Lagrange description, is the Green–Lagrange strain tensor ${}_{\mathrm{t}}{}^{\mathrm{t}}T{}_{\mathsf{\epsilon }}$, while the measure of the strain rate state is the Green–Lagrange strain rate tensor ${}_{\mathrm{t}}{}^{\mathrm{t}}T{}_{\stackrel{·}{\mathsf{\epsilon }}}$.

The tensors listed are defined and related to the current configuration ${}_{}{}^{\mathrm{t}}\mathrm{C}$. The measure of strain increment, strain rate and stress increment are the increments of strain tensors ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{T}}_{\mathsf{\epsilon }}$ and strain rate ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{T}}_{\stackrel{·}{\mathsf{\epsilon }}}$ using Green–Lagrange, respectively. The ${}_{\mathrm{t}}{}^{\mathrm{t}}\mathsf{\epsilon }{}_{\mathrm{ij}}$ components of the ${}_{\mathrm{t}}{}^{\mathrm{t}}T{}_{\mathsf{\epsilon }}$ tensor are related to ui displacements according to the formula in Cauchy notation [48]:

$${}_{\mathrm{t}}{}^{\mathrm{t}}\mathsf{\epsilon }{}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{j},\mathrm{i}}+{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{j}}\right),$$

(5)

where ${}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}$ is the partial derivative ${}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}=\partial \left({}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i}}\right)/\partial \left({}_{}{}^{\mathrm{t}}\mathrm{x}{}_{\mathrm{i}}\right)$. The components ${}_{\mathrm{t}}{}^{\mathrm{t}}\mathsf{\epsilon }{}_{\mathrm{ij}}$ can be decomposed into a linear part

$${}_{\mathrm{t}}{}^{\mathrm{t}}\overline{\mathsf{\epsilon }}{}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{j},\mathrm{i}}\right),$$

(6)

and a non-linear part

$${}_{\mathrm{t}}{}^{\mathrm{t}}\tilde{\mathsf{\epsilon }}{}_{\mathrm{ij}}=\frac{1}{2}{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{j}}.$$

(7)

The increment of the components ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathsf{\epsilon }}_{\mathrm{ij}}$ of the tensor ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{T}}_{\mathsf{\epsilon }}$ will be derived based on the definition of the increment of any function. Thus, writing the components according to Equation (5) at times τ and t and substituting them into Equation (1), after applying the incremental decomposition $({}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}={}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i},\mathrm{j}}$) and performing the appropriate algebraic operations, we obtain the following:

$$\begin{gathered} {}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathsf{\epsilon }}_{\mathrm{ij}}={}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathsf{\epsilon }{}_{\mathrm{ij}}-{}_{\mathrm{t}}{}^{\mathrm{t}}\mathsf{\epsilon }{}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{j},\mathrm{i}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{k},\mathrm{j}}\right)-\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{j},\mathrm{i}}+{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{j}}\right) \\ =\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{j},\mathrm{i}}+{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{j}}\right), \end{gathered}$$

(8)

where ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i},\mathrm{j}}=\partial \left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i}}\right)/\partial \left({}_{}{}^{\mathrm{t}}\mathrm{x}{}_{\mathrm{i}}\right)$ denotes the partial derivative.

Equation (8) is valid for various incremental descriptions. For a specific description, an appropriate simplification should be used. For the updated Lagrange description adopted here, the expression ${}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathrm{t}}\mathrm{u}{}_{\mathrm{k},\mathrm{j}}=0$. The increment of the ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathsf{\epsilon }}_{\mathrm{ij}}$ components can be decomposed into a linear part

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}\overline{\mathsf{\epsilon }}{}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{j},\mathrm{i}}\right)$$

(9)

and a non-linear part

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\tilde{\mathsf{\epsilon }}}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{u}}_{\mathrm{k},\mathrm{j}}\right),$$

(10)

In a similar manner, the following formulas for the increments of the components ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\mathsf{\epsilon }}}_{\mathrm{ij}}$ of the strain rate increase tensor ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\mathrm{T}}_{\dot{\mathsf{\epsilon }}}$ for the linear part ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\overline{\mathsf{\epsilon }}}}_{\mathrm{ij}}$ are

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\overline{\mathsf{\epsilon }}}}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\mathrm{u}}}_{\mathrm{i},\mathrm{j}}+{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\mathrm{u}}}_{\mathrm{j},\mathrm{i}}\right)$$

(11)

as is a non-linear part

$${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\tilde{\mathsf{\epsilon }}}}_{\mathrm{ij}}=\frac{1}{2}\left({}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\mathrm{u}}}_{\mathrm{k},\mathrm{i}}·{}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta {\dot{\mathrm{u}}}_{\mathrm{k},\mathrm{j}}\right).$$

(12)

The introduction of non-linear terms means that exact formulas will be used, without any linearization.

The stress state and stress increment were described by the second (symmetric) Pioli–Kirchhoff stress tensor ${}_{\mathrm{t}}{}^{\mathrm{t}}T{}_{\mathsf{\sigma }}$ and its increment ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta T_{\mathsf{\sigma }}$, respectively.

The components ${\Delta \mathsf{\sigma }}_{\mathrm{ij}}$ of the stress increase tensor ${}_{\mathrm{t}}{}^{\mathsf{\tau }}\Delta T_{\mathsf{\sigma }}$ for the body E/VP were calculated from the following formula:

$${\Delta \mathsf{\sigma }}_{\mathrm{ij}}={\mathrm{C}}_{\mathrm{ijkl}}^{\left(\mathrm{E}\right)}·\left({\Delta \mathsf{\epsilon }}_{\mathrm{kl}}-{\mathrm{A}}_{\mathrm{kl}}/\mathrm{B}\right),$$

(13)

where ${\mathrm{A}}_{\mathrm{kl}}$ are tensor components and B is a scalar, according to the following formulas:

$${\mathrm{A}}_{\mathrm{kl}}={\tilde{\mathrm{S}}}_{\mathrm{kl}}\left\{{\tilde{\mathrm{S}}}_{\mathrm{ij}}·{\mathrm{C}}_{\mathrm{ijkl}}^{\left(\mathrm{E}\right)}·{\Delta \mathsf{\epsilon }}_{\mathrm{kl}}-\frac{2}{3}·{\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)·{\dot{\mathrm{E}}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}·{\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right\},$$

(14)

$$\mathrm{B}={\tilde{\mathrm{S}}}_{\mathrm{ij}}·{\mathrm{C}}_{\mathrm{ijkl}}^{\left(\mathrm{E}\right)}·{\tilde{\mathrm{S}}}_{\mathrm{kl}}+\frac{2}{3}·{\mathsf{\sigma }}_{\mathrm{Y}}^2\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)·⟨{\dot{\tilde{\mathrm{C}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)+\frac{2}{3}·\mathrm{E}}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}\rangle .$$

(15)

In Equations (13)–(15) ${\Delta \mathsf{\epsilon }}_{\mathrm{kl}}$ are the components of the increment of the tensor of total strain (elastic, plastic and viscous), ${\mathrm{C}}_{\mathrm{ijkl}}^{\left(\mathrm{E}\right)}$ is the fourth-order tensor characterizing the elastic properties of the object material, ${\tilde{\mathrm{S}}}_{\mathrm{ij}}$ is a cumulative component of the stress deviator at the time t considered, ${\mathrm{E}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}$ is the instantaneous strain-hardening modulus, ${\dot{\mathrm{E}}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}$ is the instantaneous strain rate-hardening modulus, $\tilde{\mathrm{C}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)$ is the instantaneous modulus of translational hardening and ${\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)$ is the instantaneous plasticizing stress of the object material, depending on the instantaneous values of the effective deformation ${\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$ and visco-plastic strain rate ${\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$.

The yield stress model ${\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)$ takes the following form:

$${\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)=\mathrm{K}{\left({\mathsf{\epsilon }}_0+{\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)}^{\mathrm{n}}·{[\mathrm{a}·{\left({\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)}^3+\mathrm{b}·{\left({\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)}^2+\mathrm{c}·\left({\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)+\mathrm{d}]}^{{\mathrm{n}}_1}$$

(16)

The constants a, b, c, K, n, n₁, ${\mathsf{\epsilon }}_0$ in the yield stress model (16) were determined in experimental studies and for 41Cr4 steel are listed in

Table 1. Material parameters used in the simulation.

ρ [kg·m−3]	E [GPa]	Re [MPa]	ET [MPa]	ν [-]	K [MPa]	ε0 [-]	C [s−1]	P [-]
7830	207	310	763	0.27	938	0.05	40	5
n [-]	n1 [-]	a [-]	b [-]	c [-]	d [-]	β [-]	m [-]	${\mathsf{\epsilon }}_f$[-]
0.1694	1	3·10−10	−10−6	7·10−4	1.0	1	0.2	2.0

where ρ—mass density of the object material, E—Young’s modulus, ν—Poisson’s ratio, K—hardening parameter, Re—initial yield stress, E_T—strain-dependent hardening modulus, ε₀—initial strains, n—hardening factor, n₁—hardening coefficient depending on the strain rate, m, P—material constants determining the sensitivity of the material to the rate of plastic deformation, C—parameter dependent on the strain rate, a, b, c, d—coefficients depending on the strain rate, ${\mathsf{\epsilon }}_{\mathrm{f}}$—failure strain.

.

The increment of the yield stress at a typical step t → τ = t + Δt is defined by the following equation [74]:

$${\Delta \mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)=\frac{\partial {\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)}{\partial {\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}}·{\Delta \mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}+\frac{\partial {\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)}{\partial {\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}}·{\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$$

(17)

where ${\Delta \mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$ and ${\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$ are the increments of the visco-plastic effective strain and the effective strain rate, respectively, and are calculated from the following formulas:

$${\Delta \mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}=\sqrt{\frac{2}{3}·\Delta D_{\mathsf{\epsilon }}⨷\Delta D_{\mathsf{\epsilon }}} ,$$

(18)

$${\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}=\sqrt{\frac{2}{3}·\Delta D_{\dot{\mathsf{\epsilon }}}⨷\Delta D_{\dot{\mathsf{\epsilon }}}}$$

(19)

where the sign “⨷” means the multiplication of tensors, $\Delta D_{\mathsf{\epsilon }}$ and $\Delta D_{\dot{\mathsf{\epsilon }}}$ are increments of deviators of strain and strain rates, respectively.

Partial derivatives $\partial {\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)/\partial {\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$ and $\partial {\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)/\partial {\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}$ are calculated using a geometric interpretation of the derivative. Then the partial derivative at a given time t is equal to the tangents of the angle of inclination of the tangent to the corresponding hardening curve at time t.

2.1.3. The Equation of the Motion of a Discrete Object on Typical Step Time

An incremental functional of the total energy of the system, which consists of a tool and an object, has been developed. From the condition of the stationarity of this functional, a variational, non-linear equation of motion and deformation of the object was derived for a typical step time. This equation was solved using spatial discretization employing the finite element method, obtaining discrete systems of the equations of motion and the deformation of the object in the process of displacement of the wedge, for a typical finite element and for a typical time step $\mathrm{t}\to \mathsf{\tau }=\mathrm{t}+\Delta \mathrm{t}$. By writing specific equations of motion for all the finite elements, separated from the tool and the workpiece, after their summation, the equation of motion of a discrete object is obtained, in the updated Lagrangian description, on a typical incremental step, of the following form [74]:

$$M·\Delta \ddot{r}+C·\Delta \dot{r}+\left(K+\Delta K\right)·\Delta r=\Delta R+\Delta F+F+R$$

(20)

where $M$ is the instantaneous mass matrix, $C$ represents theinstantaneous damping matrix, $K$, $\Delta K$, respectively, indicate the instantaneous stiffness matrix and its increment, $F$, $\Delta F$, respectively, indicate the instantaneous column vector of the internal loads of nodes and its increment, $R$, $\Delta R$, respectively, represent the instantaneous column vector surface loads and its increment, $\Delta r$ indicates the column vector of node displacement increments, $\Delta \dot{r}$ indicates the node velocity increment column vector and $\Delta \ddot{r}$ indicates the node acceleration increment column vector.

The Matrix Equation (20) being a system of N differential equations of the second order with constant (incremental step) coefficients and appropriate initial conditions $\left\{r\left(\mathrm{t}=0\right)\right\}=\left\{r_0\right\}$, $\left\{\dot{r}\left(\mathrm{t}=0\right)\right\}=\left\{{\dot{r}}_0\right\}$ and boundary values constitutes an incremental formulation of the dynamic equilibrium of the deformed solids in contact, for the case of sliding the wedge on the surface of the rough rollers.

The System of Equation (20) contains N equations with 2N known elements of the vector of internal forces $F$ and external forces $R$ and $3{\mathrm{N}}^2$ elements of the matrix $M$, $C$ and $K$. However, it also contains unknown elements, i.e., the components of the vectors: an increment of displacements of nodes $\Delta r$, an increment of the velocity of nodes $\Delta \dot{r}$, an increment of node accelerations $\Delta \ddot{r}$, an increase of the internal loads of the object $\Delta F$ and ${\mathrm{N}}^2$ unknown elements of the object stiffness increment matrix $\Delta K$. In this equation, there is also a part of the components of the vector of the increase in external loads $\Delta R$ (in contact areas) that is unknown.

Equation (20) was solved using partial linearization and the method of central differences (explicit), in which the differential approximation of the derivatives of partial displacements is assumed according to the following equations:

$${\dot{r}}^{\mathrm{t}}=\frac{1}{2·\Delta \mathrm{t}}·\left(r^{t+\Delta t}-r^{t-\Delta t}\right)$$

(21)

$${\ddot{r}}^{\mathrm{t}}=\frac{1}{\Delta {\mathrm{t}}^2}·\left(r^{t+\Delta t}-2·r^t+r^{t-\Delta t}\right)$$

(22)

then Equation (20) takes the following form

$$\tilde{M}·r^{t+\Delta t}=\tilde{Q}$$

(23)

where $\tilde{M}$ is the effective mass matrix and $\tilde{Q}$ the effective load vector.

$$\tilde{M}=\frac{1}{\Delta {\mathrm{t}}^2}·M+\frac{1}{2·\Delta \mathrm{t}}·C$$

(24)

$$\tilde{Q}=R+F+\frac{2}{\Delta {\mathrm{t}}^2}·R·r^t-\left(\frac{1}{\Delta {\mathrm{t}}^2}·M+\frac{1}{2·\Delta \mathrm{t}}·C\right)·r^{t-\Delta t}.$$

(25)

From Equation (23), the unknown column vector of the displacements of the nodes of the discrete object at the end of the considered time step $r^{t+\Delta t}$ is calculated.

In the already formulated constitutive relationships (material models) Equations (8) and (11), the Cauchy index notation was used to facilitate mathematical operations. In order to use these relationships for the Object Motion Equation (23) derived and to develop computer applications for a numerical simulation of the wedge displacement process, it is convenient to write them in a matrix form.

Relationship (13) now takes a matrix form. The increment of the stress tensor $\Delta T_{\mathsf{\sigma }}$ is calculated from the following formula:

$$\Delta T_{\mathsf{\sigma }}=C^{\left(E\right)}·\left(\Delta \mathsf{\epsilon }-A/\mathrm{B}\right)=C^{\left(E\right)}·⟨\left(1-{\tilde{\mathrm{S}}}^{**}\right)·\Delta \mathsf{\epsilon }-\Delta {\mathsf{\epsilon }}^{**}\rangle ,$$

(26)

where $C^{\left(\mathrm{E}\right)}$ is the Hook matrix, $\Delta \mathsf{\epsilon }$ and $\Delta {\mathsf{\epsilon }}^{**}$ are the column vectors of the Green–Lagrange strain increment, ${\tilde{\mathrm{S}}}^{**}$ is the scalar, A is the matrix and B is the scalar, according to the following formulas:

$$A=\tilde{S}·\left\{\tilde{S}·C^{\left(\mathrm{E}\right)}·\Delta \mathsf{\epsilon }-\frac{2}{3}·{\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)·{\dot{\mathrm{E}}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}·{\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right\},$$

(27)

$$\mathrm{B}=\tilde{S}·C^{\left(\mathrm{E}\right)}·\tilde{S}+\frac{2}{3}·{\mathsf{\sigma }}_{\mathrm{Y}}^2\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)·\langle {\dot{\tilde{\mathrm{C}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)+\frac{2}{3}·\mathrm{E}}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}\rangle .$$

(28)

The increment of the total strain tensor $\Delta T_{\mathsf{\epsilon }}$ is calculated from the following formula:

$$\Delta T_{\mathsf{\epsilon }}=\frac{1}{1-{\tilde{\mathrm{S}}}^{**}}·\left(D^{\left(\mathrm{E}\right)}·\Delta \sigma -G/\mathrm{B}\right)=\frac{1}{1-{\tilde{\mathrm{S}}}^{**}}·\left(\Delta {\mathsf{\epsilon }}^{\left(\mathrm{E}\right)}-\Delta {\mathsf{\epsilon }}^{**}\right),$$

(29)

where $D^{\left(\mathrm{E}\right)}={[C^{\left(\mathrm{E}\right)}]}^{-1}$ is the inverse of Hook’s matrix, $\Delta \mathsf{\epsilon }$ is the column vector of the Green–Lagrange increment of total strain and G is a matrix according to the following formula:

$$G=\frac{2}{3}·{\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{VP}\right)},{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}\right)·{\dot{\mathrm{E}}}_{\mathrm{T}}^{\left(\mathrm{VP}\right)}·{\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{VP}\right)}·\tilde{S}.$$

(30)

In numerical calculations, the matrix notation of Formula (13) obtained after its discretization by the finite element method is also used:

$$\Delta \sigma =C^{\left(\mathrm{E}\right)}·\overline{B}+\left[\Delta q\right]·\tilde{B}·\left\{\Delta q\right\}+\Delta {\sigma }^{**}=\left(1-{\tilde{\mathrm{S}}}^{**}\right)·\left(\overline{S}+C^{\left(\mathrm{E}\right)}\right)·\left[\Delta q\right]·\tilde{B}·\left\{\Delta q\right\}+\Delta {\sigma }^{**},$$

(31)

where $\overline{B}$ and $\tilde{B}$ are the instantaneous matrices of the linear and non-linear dependence of the strain increment $\Delta \mathsf{\epsilon }$ on the displacement increment of the nodal points $\Delta r$, respectively, $\overline{S}=C^{\left(\mathrm{E}\right)}·\overline{B}$ is the stress matrix of the discrete system, $\left[\Delta q\right]$ and $\left\{\Delta q\right\}$ are, respectively, a matrix and a column vector of the displacements increment of the nodal points of the system on a step time.

In the calculation process, it is necessary to accumulate incremental quantities, while the problem of an accumulation of the components of the vectors of the displacement increase ${\Delta \mathrm{u}}_{\mathrm{i}}$, and the components of the strain increment tensors ${\Delta \mathsf{\epsilon }}_{\mathrm{ij}}$ and the components of the strain rate increment tensors ${\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{ij}}$ at any time τ is trivial. It consists of adding their respective components derived in the previous steps:

$${}_{}{}^{\mathsf{\tau }}\mathrm{u}{}_{\mathrm{i}}={{{\displaystyle \sum }}^{ }}_{\mathrm{t}=0}^{\mathsf{\tau }}{\Delta \mathrm{u}}_{\mathrm{i}}, \mathrm{i}=1 ÷ 3 (6),$$

(32)

$${}_{}{}^{\mathsf{\tau }}\mathsf{\epsilon }{}_{\mathrm{ij}}={{{\displaystyle \sum }}^{ }}_{\mathrm{t}=0}^{\mathsf{\tau }}{\Delta \mathsf{\epsilon }}_{\mathrm{ij}}, \mathrm{i}=1, 2, 3,$$

(33)

$${}_{}{}^{\mathsf{\tau }}\dot{\mathsf{\epsilon }}{}_{\mathrm{ij}}={{{\displaystyle \sum }}^{ }}_{\mathrm{t}=0}^{\mathsf{\tau }}{\Delta \dot{\mathsf{\epsilon }}}_{\mathrm{ij}}, \mathrm{i}=1, 2, 3.$$

(34)

The accumulation of yield stress increments ${\Delta \mathsf{\sigma }}_{\mathrm{Y}}$ and the increments of stress tensors $\Delta T_{\mathsf{\sigma }}$ and their components ${\Delta \mathsf{\sigma }}_{\mathrm{ij}}$ is not trivial; it requires the use of appropriate algorithms named “Yield stress accumulation” and “Stress accumulation”.

2.2. Experimental Research

2.2.1. Research Position

Experimental tests of the diamond sliding burnishing process were carried out on the stand shown in Figure 3. The NEF 400 CNC machining center with the CNC 3D Fanuc 210 is the control, manufactured by the German GILDEMEISTER Drehmaschinen GmbH from Bielefeld, DMG, was used for the tests. The device is a multifunctional machining center. It ensures a high quality and dimensionality of products.

The machining process was performed with a special tool for sliding burnishing with glued diamond tips, made of a diamond composite with ceramic 512 bonding phase, namely Ti3SiC2 [86], commissioned by the Authors to be performed by the Institute of Advanced Manufacturing Technologies in Krakow (IAMT). The construction of the tool is shown in Figure 4, while Figure 5 shows a general view of the tool and three spherical tips with different values of the radius r_d. In proper tests, a burnishing tip with a radius of ${\mathrm{r}}_{\mathrm{d}}=3.5 \mathrm{mm}$ was used. The required radius r_d was calculated from the following formula: ${\mathrm{r}}_{\mathrm{d}}=4.9-0.0055\mathrm{HV}=3.5 \mathrm{mm}$, where HV was the microhardness of the workpiece (the hardness of the samples made of 41Cr4 steel was ${\mathrm{HV}}_{30}=255$). During the burnishing process, the rollers were lubricated with Prowadol VG-68 universal lubricating oil.

The actual radius r_d of the diamond burnishing tip was measured on the Werth Video-Check-IP 250 coordinate measuring machine (Figure 5b). The measurements were made at five planes every 60°, obtaining: ${\mathrm{r}}_{\mathrm{d}}=3.335;3.347;3.340;3.334 \mathrm{and} 3.333 \mathrm{mm}$. Thus, the average rounding radius of the active surface of the diamond tip (with an accuracy of 0.1) is ${\mathrm{r}}_{\mathrm{d}}=3.3 \mathrm{mm}$. This value of the radius was assumed in numerical computer simulations.

2.2.2. Test Samples

For the experimental tests of the burnishing process, cylindrical samples (Figure 6) made of rolled bars of 41Cr4 steel in a normalized state were prepared, in accordance with the PN-EN 10002-1+AC1 standard [87]. The chemical composition of the material is shown in

Table 2. Chemical composition of 41Cr4 steel.

C [%]	Mn [%]	Si [%]	P [%]	Cr [%]	S [%]	Ni [%]	Mo [%]	W [%]	V [%]	Cu [%]
0.42	0.83	0.32	0.03	1.05	0.028	0.24	0.08	0.16	0.04	0.23

.

Before sliding burnishing, the outer surfaces of the samples were ground with a Korund grinding wheel type 05/0/8 38A80-LVBE 500 × 50 × 203 with a tool speed of 1500 [rev∙min⁻¹] with hydraulic feed. After grinding, each sample had 3 working zones with an average roughness of $\mathrm{Ra}=0.63 \mathsf{\mu }\mathrm{m}$.

The material parameters of the 41Cr4 steel used in the numerical calculations were determined on the Zwick Roell Z400E in the uniaxial tensile test of one- and three-stage samples (Figure 7 and Figure 8). The tests were carried out at an ambient temperature of 18 °C. The tests were carried out for 5 standardized cylindrical samples and for 5 three-stage cylindrical samples. The strain rate was the default value in accordance with the B standard for conducting static tensile tests. The stress–strain relationship diagram for 41Cr4 steel is shown in Figure 9.

The material constants K and n in the yield stress model (16) were determined using the R. H. Heyer method [88].

$$\mathrm{n}=\frac{\ln \left(\frac{{\mathrm{b}}_{\mathrm{C}0}}{{\mathrm{b}}_{\mathrm{B}0}}·\frac{{\mathrm{l}}_{\mathrm{C}0}}{{\mathrm{l}}_{\mathrm{B}0}}·\frac{{\mathrm{l}}_{\mathrm{B}}}{{\mathrm{l}}_{\mathrm{C}}}\right)}{\ln \left[\frac{\ln \left(\frac{{\mathrm{l}}_{\mathrm{B}}}{{\mathrm{l}}_{\mathrm{B}0}}\right)}{\ln \left(\frac{{\mathrm{l}}_{\mathrm{C}}}{{\mathrm{l}}_{\mathrm{C}0}}\right)}\right]}=\frac{\ln \left(\frac{19.8}{18.18}·\frac{20}{20}·\frac{22.46}{20.952}\right)}{\ln \left[\frac{\ln \left(\frac{22.46}{20}\right)}{\ln \left(\frac{20.952}{20}\right)}\right]}=0.1694,$$

$$\mathrm{K}=\frac{{\mathrm{F}}_{\max }}{{\mathrm{b}}_{\mathrm{C}0}·{\mathrm{g}}_0·\frac{{\mathrm{l}}_{\mathrm{C}0}}{{\mathrm{l}}_{\mathrm{C}}}·\ln \frac{{\mathrm{l}}_{\mathrm{C}}}{{\mathrm{l}}_{\mathrm{C}0}}}=\frac{164000}{19.8·307,9·\frac{20}{20.952}·\ln \frac{20.952}{20}}=938 \left[\mathrm{MPa}\right]$$

In the computer simulations, the elastic/visco-plastic Cowper–Symonds model of plastic stresses was also used. That model does not include any temperature parameters (e.g., the Johnson–Cook model requires temperature parameters) and it is appropriate for calculating the plastic strain rate which occurs during turning or diamond sliding burnishing, which was why it was used in the simulations. This model takes into account mixed isotropic-kinematic or kinematic plastic hardening and the effect of equivalent plastic strain and an equivalent rate of plastic strain, according to the Power Relation [89]:

$${\mathsf{\sigma }}_{\mathrm{Y}}={[1+\left({\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{P}\right)}/\mathrm{C}\right)]}^{\mathrm{m}}·\left({\mathrm{R}}_{\mathrm{e}}+\mathsf{\beta }·{\mathrm{E}}_{\mathrm{p}}·{\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{P}\right)}\right)$$

(35)

where:

-

B—the hardening parameter, isotropic-kinematic $\left(\mathsf{\beta }=1\right)$ or kinematic $\left(\mathsf{\beta }=0\right)$ plastic hardening,

-

${\mathrm{R}}_{\mathrm{e}}$—initial, static yield stress [MPa],

-

${\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{P}\right)}$—effective plastic strain rate [s⁻¹],

-

C—material parameter determining the influence of the strain rate of plastic deformation [s⁻¹],

-

$\mathrm{m}=1/\mathrm{P}$—material constant determining the sensitivity of the material to the plastic strain rate [-],

-

${\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{P}\right)}$—effective plastic strain [-],

-

${\mathrm{E}}_{\mathrm{p}}={\mathrm{E}}_{\mathrm{T}}·\mathrm{E}/\left(\mathrm{E}-{\mathrm{E}}_{\mathrm{T}}\right)$—material parameter dependent on the modulus of plastic hardening, where ${\mathrm{E}}_{\mathrm{T}}=\partial {\mathsf{\sigma }}_{\mathrm{Y}}\left({\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{p}\right)};{\dot{\mathsf{\epsilon }}}_{\mathrm{e}}^{\left(\mathrm{p}\right)}\right)/\partial {\mathsf{\epsilon }}_{\mathrm{e}}^{\left(\mathrm{p}\right)}$ is the strain-hardening modulus and E is the Young’s modulus of elasticity.

The material parameters in the Plastic Stress Models (16) and (35) were determined on the basis of the results of experimental tests, and it was proven that they have a significant impact on the accuracy of the simulation; these are summarized in Table 1.

—failure strain. The verification of the material parameters determined was carried out by means of a numerical analysis of the tensile strength of a standardized cylindrical sample and a comparison of the sample geometry obtained numerically with those measured experimentally. For this purpose, an original application known as “Stretching a cylindrical sample” in the APDL language in the ANSYS/LS-Dyna program was developed. As input data for the simulation, the material parameters previously determined from the experiment for the power model were introduced. Figure 10 shows selected stages of a numerical analysis of the sample stretching process. Figure 10a already shows the stage of the formation of the so-called “necks”. It is the next stage after the stage of elastic deformation. The material flowed and the neck formed. This is followed by plastic hardening of the material. During the narrowing of the so-called “neck” there is a stress concentration in this area. The strengthening capacity of the material is no longer high enough to overcome the loss in the width and thickness of the material. After exceeding the permissible deformations, the sample breaks in its smallest cross-section (Figure 10b). The greatest stresses occur under the surface of the sample at the fracture site. It is in this area that cracking occurs first, and then it moves outward from the sample.

A numerical measurement of the sample geometry after fracture was also carried out. For this purpose, a grid with a pitch of 1 mm was applied. The diameter of the sample in the neck measured in the experiment was $\mathrm{d}=14.98 \mathrm{mm}$, while the numerically calculated one was $\mathrm{d}=14.50 \mathrm{mm}$. Statistical tests of the significance of dimension differences showed that they are insignificant at the significance level of α = 0.05. This proves the correctness of the models, algorithms and numerical application developed.

2.2.3. Experiment Plan and Methodology of Statistical Elaboration of Research Results

Experimental tests were carried out in accordance with the five-level rotatable experiment plan. The experimental plan along with the actual and coded values of the parameters are listed in

Table 3. The five-level rotatable experiment plan.

Level Plan	Investigated Factors
	Real Values		Coded Values
	Feed	Force	Feed	Force
	f [mm·rev−1]	F [N]	f [-]	F [-]
1	0.05	8	-1	-1
2	0.09	8	+1	-1
3	0.05	22	-1	+1
4	0.09	22	+1	+1
5	0.10	15	α = +1.414	0
6	0.04	15	α = −1.414	0
7	0.07	25	0	α = +1.414
8	0.07	5	0	α = −1.414
9	0.07	15	0	0
10	0.07	15	0	0
11	0.07	15	0	0
12	0.07	15	0	0
13	0.07	15	0	0

.

Ranges of the variability of technological parameters in the sliding burnishing process were determined. The assumed area of research in the considered case is as follows:

(1)

feed: f = 0.04 ÷ 0.1 [mm∙rev⁻¹],

(2)

force: F = 5 ÷ 25 [N].

Three-fold repeatability of the tests was used. This task required the following steps to be conducted [90]:

• Determination of the variability range of the parameters studied,
• Choice of the class of the mathematical model,
• Coding the parameters analyzed,
• Collection of the experiment results,
• Elimination of results with gross mistakes,
• Calculation of inter-row variance and standard deviation,
• Checking the homogeneity of variance,
• Calculating the coefficients of regression function,
• Statistical analysis of the regression function,
• Examination of the significance level of the correlation coefficient,
• Checking the adequacy of the mathematical model,
• Decoding the regression function.

The average values of the outputs of the object $\widehat{\mathrm{Ȳ}}$ were approximated with a polynomial of the second degree with double interaction, obtaining the regression equation as two-parameter functions:

$$\widehat{\mathrm{Ȳ}}={\mathrm{b}}_{\mathrm{o}}+{\mathrm{b}}_1·{\overline{\mathrm{x}}}_1+{\mathrm{b}}_2·{\overline{\mathrm{x}}}_2+{\mathrm{b}}_{11}·{\overline{\mathrm{x}}}_1{}^2+{\mathrm{b}}_{22}·{\overline{\mathrm{x}}}_2{}^2+{\mathrm{b}}_{12}{\overline{·\mathrm{x}}}_1·{\overline{\mathrm{x}}}_2, \mathsf{\mu }\mathrm{m}$$

(36)

where

${\mathrm{b}}_{\mathrm{o}},{ \mathrm{b}}_1,{\mathrm{b}}_2,{ \mathrm{b}}_{11},{ \mathrm{b}}_{22}{ \mathrm{and} \mathrm{b}}_{12}$—unknown coefficients of the Regression Equation,

${\overline{\mathrm{x}}}_{\mathrm{i}}$—input variables: ${\overline{\mathrm{x}}}_1=\overline{\mathrm{f}}$ [mm·s⁻¹] and ${\overline{\mathrm{x}}}_2=\overline{\mathrm{F}}$ [N].

Using matrix calculus, the column vector {b} of the unknown coefficients in Equation (36) was calculated from the following matrix formula:

$$\left\{\mathrm{b}\right\}={({\left[\overline{\mathrm{X}}\right]}^{\mathrm{T}}\left[\overline{\mathrm{X}}\right])}^{-1}{\left[\overline{\mathrm{X}}\right]}^{\mathrm{T}}\left\{\overline{\mathrm{Y}}\right\},$$

(37)

where:

-

$\left[\overline{\mathrm{X}}\right]$—input variable matrix of dimension N×L, for data N = 13 and L = 6,

-

${\left[\overline{\mathrm{X}}\right]}^{\mathrm{T}}$—transposed matrix $\left[\overline{\mathrm{X}}\right]$,

-

${({\left[\overline{\mathrm{X}}\right]}^{\mathrm{T}}\left[\overline{\mathrm{X}}\right])}^{-1}$—covariance matrix,

-

$\left\{\overline{\mathrm{Y}}\right\}$—column vector of the average values of the experimental results.

The boundaries of the confidence region for Regression Function (36) were determined from the following formula:

$$\widehat{\overline{\mathrm{Y}}}\pm {\mathrm{t}}_{\mathrm{kr}\left(\mathsf{\alpha };\mathrm{f}=\mathrm{N}-\mathrm{L}\right)}·\frac{{\mathrm{S}}_{\mathrm{R}}}{\sqrt{\mathrm{N}-\mathrm{L}-1}}·\sqrt{{\left\{\overline{\mathrm{x}}\right\}}^{\mathrm{T}}{({\left[\overline{\mathrm{X}}\right]}^{\mathrm{T}}\left[\overline{\mathrm{X}}\right])}^{-1}\left\{\overline{\mathrm{x}}\right\}},$$

(38)

where:

$\widehat{\overline{\mathrm{Y}}}$—regression equation according to Equation (36),

${\mathrm{t}}_{\mathrm{kr}\left(\mathsf{\alpha };\mathrm{f}\right)}$—critical value of Student t test for significance level α = 0.05 and the number of the degrees of freedom f = N − L,

N—number of measurement points in the experimental design: N = 13,

L—number of unknown coefficients in Regression Equation (36); here, L = 6,

$\left\{\overline{\mathrm{x}}\right\}$ and ${\left\{\overline{\mathrm{x}}\right\}}^{\mathrm{T}}$—column vector of the functions of input variables (test factors in real form) and its transposition: ${\left\{\overline{\mathrm{x}}\right\}}^{\mathrm{T}}=\left[1 {\overline{\mathrm{x}}}_1 {\overline{\mathrm{x}}}_2 {\overline{\mathrm{x}}}_1{}^2 {\overline{\mathrm{x}}}_2{}^2 {\overline{\mathrm{x}}}_1·{\overline{\mathrm{x}}}_2\right]$,

${\mathrm{S}}_{\mathrm{R}}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{i}=\mathrm{N}}{\left({\widehat{\overline{\mathrm{y}}}}_{\mathrm{i}}-{\overline{\mathrm{y}}}_{\mathrm{i}}\right)}^2$—residual variance,

${\widehat{\overline{\mathrm{y}}}}_{\mathrm{i}}$—average values of model outputs for plan points calculated from Equation (36),

${\overline{\mathrm{y}}}_{\mathrm{i}}$—average values of experimental results.

The test results after statistical processing according to the algorithm presented in the study [90] were used to develop Regression Equations (38) and (41).

3. Results

3.1. Results of Numerical Analyses

The results of the displacement of a non-deformable wedge (E → ∞) on an elastic/visco-plastic substrate, which has the material characteristics of 41Cr4 steel after normalization annealing, are presented. For this steel, a Cowper–Symonds non-linear material was selected, according to Equation (35), obtaining a model dependent on the equivalent plastic strain and the equivalent rate of plastic strain according to the power law.

The value of the static friction coefficient was ${\mathsf{\mu }}_{\mathrm{s}}=0.1$, and the dynamic coefficient ${\mathsf{\mu }}_{\mathrm{d}}=0.05$. Wedge rounding radius ${\mathrm{r}}_{\mathrm{d}}=25 \left[\mathsf{\mu }\mathrm{m}\right]$, wedge internal angle β = 60 [°], wedge depth ${\mathrm{a}}_{\mathrm{p}}=25 \left[\mathsf{\mu }\mathrm{m}\right]$ and ${\mathrm{a}}_{\mathrm{p}}=250 \left[\mathsf{\mu }\mathrm{m}\right]$, wedge travel distance $\mathrm{s}=0.002 \left[\mathrm{m}\right]$, simulation time t= 0.0002 [s]. The displacement speed of the wedge is ${\mathrm{v}}_{\mathrm{c}}=10 \left[\mathrm{m}·{\mathrm{s}}^{-1}\right]$. The facility was divided into 60.000 PLANE162 finite elements, which were compacted in the contact zone. This number is based on multiple sensitivity analyses for different workpiece materials. For 41Cr4 steel and given values of ${\mathrm{a}}_{\mathrm{p}}$ and s, this number is the optimal value. The shape factor $\mathrm{WK}=1$ was assumed, resulting from the Authors’ own simulation studies and those performed by other researchers [65,68,79]. The Cowper–Symonds model contains a critical strain value above which the element is removed (the so-called element removal method). The value of this parameter was determined and it equals ɛ_f = 1.7.

The wedge in the contact area was discretized using 18.000 finite elements. A regular hexagonal grid with regular fields was used in the simulations. In the material machined, the mesh was refined in the contact zone and in the surface layer (the element’s size being 0.1 × 0.1 μm). The influence of the inclination of the wedge relative to the ground on the type of machining (micro-cutting or sliding burnishing) as well as on the occurring strains and stresses was examined. The range of the rake angle was $\mathsf{\gamma }=-7010 [°]$ with a step of 5 [°]. Selected steps from the simulation are shown in the screenshots. Some drawings have a 5 × 5 [μm] grid in order to better illustrate the size of the flash or chip created. The scale in the drawings shows the values of Huber–Mises–Hencky (HMH) stresses and is expressed in [Pa] or HMH strains expressed in a dimensionless number [-]. For each wedge inclination angle, an exemplary displacement of the node along the Y axis located just below the layer processed, located exactly 1 [mm] from the beginning of the sample, is also provided.

Sixteen numerical simulations were performed for wedge depressions ${\mathrm{a}}_{\mathrm{p}}=25 \left[\mathsf{\mu }\mathrm{m}\right]$ and ${\mathrm{a}}_{\mathrm{p}}=250 \left[\mathsf{\mu }\mathrm{m}\right]$, for different rake angles γ. In each simulation, the displacement along the Y axis of the same node, the number and location of which is shown in Figure 11, was analyzed.

A summary of the simulation results is shown in Figure 12. From the graph, one may read the limit rake angles of the wedge for a given recess at which the sliding burnishing process takes place. For depth ${\mathrm{a}}_{\mathrm{p}}=25 \left[\mathsf{\mu }\mathrm{m}\right]$ and rake angles ${\mathsf{\gamma }}_{\mathrm{gr}}<-43 [°]$, ${\mathrm{a}}_{\mathrm{p}}=250 \left[\mathsf{\mu }\mathrm{m}\right]$ and ${\mathsf{\gamma }}_{\mathrm{gr}}<-50 [°]$, the surface is smoothed, while rake angles are greater for micro-cutting. The influence of cutting conditions on the form of the chip created is also visible.

For example, the case for ${\mathrm{a}}_{\mathrm{p}}=25 \left[\mathsf{\mu }\mathrm{m}\right]$ and $\mathsf{\gamma }=-70 [°]$ is described in detail. Figure 13 shows stress and strain maps for $\mathsf{\gamma }=-70 [°]$. For this wedge rake angle, we observe the phenomenon of flash formation but without chip formation. This process may be compared to smooth burnishing. The maximum values of HMH stresses at the selected time step prevail in the contact area of the wedge with the workpiece and are ${\mathsf{\sigma }}_{\mathrm{e},\max }=3670 [\mathrm{MPa}]$ (Figure 13a). The largest deformations were at the point of contact of the rake face of the wedge with the workpiece, and they amounted to ${\mathsf{\epsilon }}_{\mathrm{e},\max }=1.99$ [-] (Figure 13b).

From the mesh applied (Figure 14a) it can be read that the height of the flash created was ca. 5 [μm]. The resulting burr due to the large rake angle had a size of approx. 20 [μm]. For this configuration of the wedge, the stresses propagate into the material. Figure 14b shows a diagram of the displacement of the node lying under the contact surface at a distance of 5 μm along the Y axis. It shows that in the initial phase of the process, up to ca. 8 [ms], the node was slowly displaced due to the approach of the tool. At 10 [ms], the tool was located exactly above the selected node, hence its rapid displacement to the value of −2.5 μm. After the tool had passed through, the elastic return of the material node to the value of −1.0 μm took place. Compared to the initial value, this result shows a WW reduction of 1.5 μm.

The exemplary results presented of numerical analyses allow, for a given tool geometry and the required burnishing depth, to determine the boundary rake angles for which the process of smoothing surface irregularities (burnishing) takes place.

3.2. Results of Experimental Studies

After the diamond sliding burnishing process, surface roughness parameters were measured: the root mean square height Sq [μm]; skewness Ssk [-]; kurtosis Sku [-]; the maximum peak heigh of the surface Sp [μm]; the maximum pit height of the surface Sv [μm]; the maximum heigh of the surface Sz [μm]; the arithmetical mean height of the surface Sa [μm]. Other factors measured were the height parameters of the roughness profile: the average maximum profile height of the ten largest peak-to-valley separations in the evaluation area Rz [μm], the arithmetical mean deviation of the profile Ra assessed [μm] and the total height of the profile Rt [μm] (

Table 4. Average values of surface and profile roughness parameters.

Plan Level	Surface Roughness Parameters							Profile Roughness Parameters
	Sq [μm]	Ssk [-]	Sku [-]	Sp [μm]	Sv [μm]	Sz [μm]	Sa [μm]	Rz [μm]	Rt [μm]	Ra [μm]
1	0.465	−0.028	3.14	2.39	2.54	4.93	0.369	1.98	2.30	0.326
2	0.769	−0.189	3.51	4.30	4.01	8.31	0.614	5.11	6.30	0.549
3	0.583	−0.255	3.70	5.20	2.94	8.14	0.444	2.74	3.22	0.384
4	0.721	0.098	3.31	3.82	4.38	8.20	0.568	4.65	5.17	0.557
5	0.753	0.003	3.18	3.74	3.96	7.70	0.599	4.23	4.52	0.556
6	0.567	−0.149	3.51	3.75	3.46	7.21	0.441	3.37	3.48	0.396
7	0.673	−0.187	3.36	3.66	3.82	7.47	0.528	3.87	4.40	0.479
8	0.957	−0.543	3.94	4.38	5.66	10.0	0.740	5.11	5.90	0.605
9	0.722	−0.212	3.30	3.84	4.48	8.32	0.571	3.59	4.06	0.529
10	0.810	−0.378	3.82	5.50	4.66	10.2	0.625	3.42	4.84	0.549
11	0.754	−0.245	3.44	3.98	4.51	9.20	0.613	3.61	4.13	0.473
12	0.796	−0.276	3.37	4.45	4.47	9.75	0.587	3.49	4.66	0.508
13	0.783	−0.306	3.67	4.38	4.58	8.64	0.596	3.54	4.28	0.504

). The measurements were carried out in the Laboratory of micro- and nanoengineering located in the Department of Precision Mechanics, the Koszalin University of Technology, using a Taylor Hobson Precision device with the TalySurf CCI6000 system. A Nikon lens with a magnification of 50× was used. The results were processed in the TalyMap Platinum program. Examples of the surface roughness results according to the Ra parameter after the sliding burnishing process are listed in

Table 5. Results of measurements of the Ra [μm] parameter of profile roughness.

Plan Level	Rresults of Measurements of the Ra [μm] Parameter of Profile Roughness
	Ra1 [μm]	Ra2 [μm]	Ra [μm]	Average Ra [μm]
1	0.295	0.350	0.335	0.326
2	0.58	0.504	0.565	0.549
3	0.39	0.364	0.398	0.384
4	0.572	0.557	0.544	0.557
5	0.535	0.556	0.577	0.556
6	0.395	0.37	0.424	0.396
7	0.477	0.481	0.479	0.479
8	0.641	0.623	0.551	0.605
9	0.528	0.527	0.534	0.529
10	0.453	0.540	0.656	0.549
11	0.458	0.473	0.49	0.473
12	0.479	0.513	0.533	0.508
13	0.490	0.521	0.502	0.504

. The view of the surface of the sample after grinding and after diamond sliding burnishing is shown in Figure 15.

Exemplary results of the study of the stereometry of the surface machined, which had the lowest value of the Ra parameter, are presented in Figure 16. All the surfaces analyzed have a degressive-progressive material share curve. The traces of diamond sliding burnishing are clearly visible in the autocorrelation function. The distribution of ordinates is close to normal with slight positive skewness. The distribution of the tops of the local hills is close to normal, consisting of two overlapping independent processes, as evidenced by the presence of two peaks.

The results of the research, after statistical elaboration, allowed the Authors to obtain the required mathematical model of the process. For example, the model describing the dependence of the arithmetic mean deviation of the profile from the mean line Ra of the surface after the process on technological parameters describes the equation for the coded variables ${\check{\overline{\mathrm{x}}}}_1$ and ${\check{\overline{\mathrm{x}}}}_2$:

$$\widehat{\overline{\mathrm{Y}}}={\widehat{\overline{\mathrm{R}}}}_a=0.5134+0.055·{\check{\overline{\mathrm{x}}}}_1-9.97·{10}^{-3}·{\check{\overline{\mathrm{x}}}}_2-3.61·{10}^{-6}·{\left({\check{\overline{\mathrm{x}}}}_1\right)}^2+3.9·{10}^{-4}·{({\check{\overline{\mathrm{x}}}}_2)}^2-6.17·{10}^{-3}·{\check{\overline{\mathrm{x}}}}_1·{\check{\overline{\mathrm{x}}}}_2,$$

(39)

where ${\check{\overline{\mathrm{x}}}}_i\in \left[-1;+1\right]; \mathrm{i}=1, 2$.

After decoding Equation (39) according to the code formulas:

$${\check{\overline{\mathrm{x}}}}_1=\left(\overline{\mathrm{f}}-0.07\right)/0.021,$$

(40)

$${\check{\overline{\mathrm{x}}}}_2=\left(\overline{\mathrm{F}}-15\right)/7.072,$$

(41)

the equation for the real values of the technological parameters f and F is obtained, along with the 95% confidence area of the following form:

$$\begin{gathered} \widehat{\overline{\mathrm{Y}}}={\widehat{\overline{\mathrm{R}}}}_a=-0.1464+14.9·\overline{\mathrm{f}}+0.003296·\overline{\mathrm{F}}-0.08227·\overline{\mathrm{f}}·\overline{\mathrm{F}}-71.38·{\overline{\mathrm{f}}}^2+1.56·{10}^{-5}·{\overline{\mathrm{F}}}^2 \\ \pm 0.0103·\sqrt{3.744·[{\left({\check{\overline{\mathrm{x}}}}_1\right)}^4+{\left({\check{\overline{\mathrm{x}}}}_2\right)}^4]+7.488·{\left({\check{\overline{\mathrm{x}}}}_1\right)}^2·{\left({\check{\overline{\mathrm{x}}}}_2\right)}^2-0.989·\left[{\left({\check{\overline{\mathrm{x}}}}_1\right)}^2+{\left({\check{\overline{\mathrm{x}}}}_2\right)}^2\right]+2.6}, \end{gathered}$$

(42)

where $\overline{\mathrm{f}}\in \left[0.04;0.1\right] \left[\mathrm{mm}·{\mathrm{rev}}^{-1}\right]$, $\overline{\mathrm{F}}\in \left[5;25\right] \left[\mathrm{N}\right],$ giving ${\check{\overline{\mathrm{x}}}}_1$ and ${\check{\overline{\mathrm{x}}}}_2$ according to formulas (40) and (41).

The graphs presented in Figure 17 show that the surface roughness Ra after burnishing treatment is affected by both the feed rate f and the burnishing force F. Figure 17b shows the special effect of the feed rate on the quality of the surface obtained. For small feed values, the surface obtained was characterized by the lowest value of the Ra coefficient. As the feed rate increased, the roughness also increased. Increasing the clamping force of the diamond burnishing tool in the case of 41Cr4 steel had a negative effect on the surface quality. In general, an increase in the clamping force resulted in an increase in the roughness coefficient Ra. Its lowest value was recorded for the lowest burnishing force.

Experimental research shows that the diamond sliding burnishing process is a very complex process. The technological quality of the product depends on many parameters related both to the condition of the semi-finished product and the conditions of the burnishing process. An incorrect selection of these leads to obtaining a product with defects, e.g., visible grooves in the material due to an incorrectly selected value of the burnishing force, etc. The most important parameters determining the quality of the surface, characterized by the Ra parameter include the tool geometry, the feed rate and the burnishing force. The oil mist lubrication used is also an important factor.

The regression equations developed can be used to determine the burnishing process conditions in terms of the required product quality, or vice versa, for given burnishing conditions, the quality of the surface machined can be predicted. Determining the optimal conditions for the process implementation requires solving the problem of a multi-criteria optimization with constraints. Solving the optimization task will allow one to obtain a product of the required quality at minimum production costs.

3.3. Validation of Models and Numerical Algorithms

The results of experimental research were also used to verify the developed models and algorithms of numerical calculations. For example, in Figure 18, Figure 19 and Figure 20, the results of numerical analyses on the surface roughness profile were compared with the results of experimental measurements.

In order to better compare the surface roughness from the experimental tests with the simulation tests, the points from the simulation in the ANSYS program were read and implemented in the Matlab program. The roughness profile generated from the ANSYS program is shown in Figure 19 and Figure 20. The comparison of the results of numerical analyses with the results of measurements shows that the differences at the adopted significance level of α = 0.05 are not significant. Therefore, both the process model and computer applications in the ANSYS/LS-Dyna system are developed correctly.

The comparison of the results of numerical analyses with the results of measurements shows that the differences at the adopted significance level of α = 0.05 are not significant. Therefore, it can be concluded that it is possible to predict the quality of diamond sliding burnished products with an accuracy sufficient for engineering practice on the basis of computer simulation results. The results of experimental research showed the usefulness of the developed mathematical models and algorithms for numerical calculations of the burnishing process. Numerical algorithms implemented in the ANSYS system can be used to assess the impact of burnishing technological conditions on product quality. They allow for a better understanding of the phenomena occurring in the contact and deformation zones and, therefore, they may constitute the basis for developing guidelines for the selection of burnishing conditions, due to the required technological quality of the product.

4. Summary and Conclusions from Numerical Analyses

(1)

The workpiece in the process of moving the wedge is subjected to a complex deformation process, characterized by large displacements and large deformations. Rendering complex changes in the outline of an object requires a sufficiently fine mesh of finite elements. The use of a very fine mesh of finite elements for the entire workpiece and throughout the process leads to a numerical model with a very large number of unknowns, requiring long calculation times. The developed effective discrete models and numerical algorithms allow for automatic mesh refinement only in areas with high curvature and large deformation gradients. In the numerical examples presented in this article, discretization with 4-node PLANE 162 elements was used. The finite elements used to discretize the object well modeled the complex deformation state, and at the same time, they are characterized by high computational efficiency.

(2)

The shape factor of the finite element (SF) significantly affects both the accuracy of mapping the outline of the tool and the value of stresses and strains as well as the course of cracking. The most favorable results of the stresses and strains calculated are obtained for $\mathrm{SF}=1$. Further increase of SF (up to nine times) has no significant impact on the accuracy of the stresses and strains calculated, defined as the difference of their values in elements and nodes.

(3)

The results of the computer simulation of the wedge displacement process on the effective model and real material confirmed the results of the experimental research. This proves the correctness of the developed mathematical model and the method of integrating the discrete equation of motion.

(4)

The results of numerical calculations confirm the possibility of a correct analysis of the burnishing or micro-cutting processes and the determination of the states of displacements, strains and stresses, using the discrete effective models developed.

(5)

The application developed in the ANSYS system enables a time analysis of deformation and stress states at any point of the material.

(6)

For example, Figure 14b shows the displacement of selected nodes in the contact area in the direction of the Y axis. It is possible to determine the separation point of the material moving with the chip or remaining in the surface layer of the product.

(7)

The basic problem in the design of the wedge displacement process is the development of appropriate tool constructions and the selection of process conditions in order to ensure the technical requirements and properties of the product obtained, while increasing tool life and process efficiency. The results of computer simulations may form the basis for the selection of technological parameters in the machining or burnishing process, and they can be used to develop a new technology with the necessary technological equipment.

(8)

As a result of the analyses carried out, it was found that the angle of inclination of the shear plane is not the same for each rake angle γ and it varies in the range of $\mathsf{\Phi }=37 [°]-55 [°]$, which is consistent with the values given in the professional literature. The values of maximum stresses, for almost every simulation, prevailed in the contact area of the wedge tip with the base material, as well as in the contact areas of the wedge rake surface with the surface of the forming chip, and in the area of the shear plane.

(9)

For given machining conditions and tool depth ${\mathrm{a}}_{\mathrm{p}}>{\mathrm{r}}_{\mathrm{d}}\left(1-\sin \mathsf{\gamma }\right)$ there is a limited rake angle ${\mathsf{\gamma }}_{\mathrm{gr}}$ below which the burnishing process takes place (no chip) and above it, the micro-cutting process (chip formation) occurs. For example, Figure 12 shows two cases of wedge displacement, for which the limit rake angle was, respectively, the depth of ${\mathsf{\gamma }}_{\mathrm{gr}}=-43 [°]$ for the depth of ${\mathrm{a}}_{\mathrm{p}}=25 \mathsf{\mu }\mathrm{m}$ and ${\mathsf{\gamma }}_{\mathrm{gr}}=-50 [°]$ for the depth of ${\mathrm{a}}_{\mathrm{p}}=250 \mathsf{\mu }\mathrm{m}$.