A novel multiscale material plasticity simulation model for high-performance cutting AISI 4140 steel

The achievable machined surface quality relies significantly on the material behavior during the high-performance cutting process. In this paper, a multiscale material plasticity simulation framework is developed to predict the deformation behaviors of AISI 4140 steel under various high-performance cutting conditions. The framework was built by coupling a three-dimensional discrete dislocation dynamic (3D-DDD) model with a finite element method (FEM) through the optimization of a dislocation density-based (DDB) constitutive equation (compiled as a user-defined subroutine in ABAQUS). The movement of edge and screw dislocations such as generation, propagation, siding, and their interactions, was performed by 3D-DDD, and the statistical features of dislocations were used to optimize the critical constants of the DDB constitutive equation. For validation, a classic FEM cutting model (Johnson-Cook constitutive equation) was employed as a reference. The simulation results indicated that the proposed multiscale model not only can precisely predict the stress, strain, cutting force, and temperature as those predicted by the classic FEM simulations, but also capture the microstructure characteristics such as grain size and dislocation density distributions under the tested cutting conditions. Severe dynamic recrystallization phenomena were found at the core shear zones. The recrystallization process reached a dynamic equilibrium at the machined surfaces when the cutting speed is larger than 280 m/min or the external-assisted temperature is between 200 and 350°, indicating an optimal range of machining parameters for improved surface integrity.


Introduction
High-performance manufacturing is a promising technology that comprises new developments to improve machining abilities through integrated in-process inspection and optimization of the processing chain [1,2]. Many mechanical machining methods, such as high-speed cutting [3], hybrid machining [4], and ultra-precision machining [5], have been continuously developed to produce functional surfaces. The complexity of material behaviors under different machining conditions brings extreme challenges in developing a "microstructural alteration control" method to address the link between surface quality and functional performance [6,7]. The dislocation evolution and grain refinement at the core shear zones and machined surface zone were reported to affect the surface functionality and reliability of machined components. For instance, Courbon et al. [8] studied the dynamic recrystallization phenomena of chips when cutting ferritic-pearlitic steels and observed the stain accumulation phenomenon at the shear deformation region. Ginting et al. [9] reported the dependence of the surface integrity properties, including roughness, defects, microhardness, and microstructure on machining parameters in dry milling titanium alloy. Bai et al. [10] studied the dislocation distribution at subsurface damage layer in micro-milling Ti-6Al-4V, and pointed out that the work hardening effect might attribute to the dislocations piling-up zone and persistent slip band observed at the machined subsurface. However, it is recognized that the fundamental research on the microstructural evolution processes in high-performance cutting requires many experimental matrices and cumbersome procedures of nano-characterization tests.
In contrast, modelling and simulation methods are continuously developed and constitute a computational way to reveal the underlying mechanisms of material microstructural evolution under various machining conditions. According to the differences in physical/mechanical hypotheses, those cutting models could be broadly classified as atomic-scale simulation, mesoscopic calculation, multiscale modelling, and finite element methods (FEM). The typical examples and highlights of their advantages and disadvantages are summarized in Table 1. For the atomistic and mesoscopic models, the detailed information of dislocations activities in a cutting process can be predicted from the process simulation. Tong et al. [11,12] demonstrated how dislocation activities (nucleation, propagation, and interaction) could alter the material removal mechanism in nanometric cutting when using a single tip and a multi-tip nanoscale diamond cutting tools. To enlarge the simulation size scale, Bai et al. [13,14] proposed a new mesoscopic modelling method according to DDD to study the microstructure characteristics of Ti-6Al-4V under various conditions. Based on quasicontinuum (QC), a hybrid MD-FEM model [15] has been proposed to address the nanoscale deformation phenomena in the meso-macroscale machining process. However, the computational load of those methods remains huge, which significantly hinders their applications in practice.
The FEM-based cutting model has been widely used to predict material deformation behaviors in various cutting conditions but limited by its continuum mechanics hypothesis to provide any detailed microstructural information. In the past decade, a set of phenomenological constitutive equations were adopted to consider the grain size evolution in the machining process, such as the Zener-Holloman model [22], Kocks-Mecking model [23], and Johnson-Mehl-Avrami-Kolmogorov model [24]. The crystal plasticity finite element method (CPFEM) also provides a tool for describing the mechanical response in low-speed cutting [20]. For phenomenological constitutive equations, modifications of constitutive relations are necessary to account for the response rate of material deformations and their non-linear behaviors. Most recently, many novel physical-based constitutive equations have been developed to consider the dislocation activities and grain deformations in machining. It provides a promising approach to reduce the computational load while improving the reliability of FEM-based cutting simulations. Typically, Melkote et al. [24] presented an enhanced physical-based material model to predict the chip formation and grain refinement process in the core shear band. A microstructural alteration-induced flow stress softening effect was observed. Based on the dislocation density-based (DDB) constitutive equation [28], Ding and Shin revealed the grain size change and grain misorientation mechanisms in the machined surface of commercially pure titanium [25], aluminum alloy [26], and AISI 52100 steel [27]. Akhondzadeh et al. [29] further modified the DDB strain hardening model through the systematic coarse-graining analysis of 3D-DDD simulation under uniaxial tension. Their modified model successfully obtained the strain hardening rate against loading orientations and compared with the experimental results.
In this paper, a multiscale material plasticity simulation framework, coupling of a 3D-DDD model with a DDB constitutive equation, is proposed to simulate the microstructure transformation in high-performance cutting AISI 4140 steel. The 3D-DDD model is adopted to calculate several critical values that used to optimize the DDB constitutive equation coefficients pertaining to controlling dislocation multiplication and annihilation activities. Subsequently, a DDD-FEM hybrid cutting model is built to simulate the machining processes, and the influences of processing parameters, e.g., cutting speeds and external-assisted temperatures, on the machined surface integrity are studied and discussed in details.

Multiscale material plasticity simulation framework DDD-FEM
The framework of the proposed DDD-FEM is shown in Fig. 1. The initialization parameters including workpiece initial microstructure, mechanical properties, and machining conditions are numerical inputs. To ensure the simul a t i o n a c c u r a c y , t h e d e t e r m i n a t i o n o f t h e dislocation-related coefficients of DDB constitutive equation in FEM module is crucial. Normally, the values of α * , β * , and k 0 , are calibrated by comparing the DDB constitutive equation predictions of the flow stress data with the classic model predicted stress-strain relationships. This however will induce parametric uncertainties. It has been reported that dislocation features calculated from DDD can provide an effective approach to test the validity of microstructure-based constitutive relations [29]. Hence, the DDD simulation module is used to collect detailed physics-based dislocation density evolution under various boundary conditions in the proposed framework.
The multiscale coupling was realized by tracking the statistic features of dislocation nucleation, motion, collision, annihilation, junction zipping, and unzipping, etc., in DDD module, based on which the coefficients of DDB constitutive equation are optimized. The unified constitutive equation was then embedded into the VUHARD subroutine of ABAQUS software by compiling the FORTRAN script and calculated the dislocation fields. A multi-physics thermos-mechanical coupling model was established to validate the effectiveness of present coupling for the specific workpiece material. The proposed simulation framework can then be used to simulate the high-performance cutting process to gain the essential evolution features of material mechanical behaviors and microstructure characteristics of the chip and machined surface. The grain size distribution, dislocation density evolution, crystal orientation change, cutting force, etc., can be exported as outputs.

Dislocation density-based constitutive equation
The fundamental assumption of present model is to include dislocation cell structures in the DDB constitutive equation [30]. The cell consists of dislocation cell interior structure and dislocation wall structure. The increase rate of cell interior dislocation density follows the basic physics-based rule: where ρ c and ρ w denote the dislocation density for cell interior and cell wall, respectively.ρ c is the variation rate of ρ c ; α * governs the multiplication process of Frank-Read sources;γ r w denotes the shear strain rate of grain boundary, and b is the Burgers vector. Some discrete dislocations will leave the cell interior during the deformation and become part of the boundary. Therefore, the loss rate of dislocation density within the cell interior can be expressed as follows: where β * governs the interaction process between cell interior and wall, f represents the volume fraction of cell wall, andγ r c is considered as the shear strain rate for cell interior. Meanwhile, the dynamic recovery of grain includes the dislocation annihilation, junction, and cross-slip activities and are descript as below: where k 0 governs the loss rate of dislocations,γ 0 represents t h e r e f e r e n c e s t r a i n r a t e , a n d n d e n o t e s t h e temperature-sensitive coefficient. The dislocation density evolution for cell interior is governed by: Similarly to the Eq. (1), the dislocation nucleation rate in wall structure can be summarized as below: whereρ w is the variation rate of ρ w . As shown in Eq. (2), the discrete dislocations leaving the cell interior are integrated into cell boundary as follows: The difference between Eqs. (2) and (6) results from the physical difference between the volumes of cell interior and wall structure. Furthermore, the loss rate of dislocations in boundary can be taken into account by the following: The above Eqs. (4), (5), and (6) result in an integrated equation: The grain size d, total dislocation density ρ tol , and shear strain rateγ can be obtained as follows: where f 0 and f ∞ denote the initial volume fraction and the saturation volume fraction. M is Taylor factor andε is von Mises strain rate. Thus, the resolved stress τ r is as follows: Fig. 1 The framework of DDD-FEM multiscale simulation.
where m represents the strain rate sensitivity coefficient, receptivity. Based on the DDB constitutive equation, the twinning [31], phase transformation [32], and grain misorientation [26] can all be further calculated.

Determination of the dislocation-related coefficients
In the DDB constitutive equation, a series of constants are used to describe the level of dislocation interactions. Excessive optimization objects and a limited set of experimental data will result in over-fitting issue. In the present work, the critical conditions associated with the onset of dynamic recrystallization are decided by dislocation-related parameters in the DDB constitutive equation. Indeed, the flow stress and microstructure transformation are governed by the multiplication and interactions of dislocations in response to the applied loadings [33]. The dislocation dynamics simulation can accurately calculate the plastic response from the collective movement of dislocation networks [34,35]. An open-source dislocation dynamics model-ParaDiS (developed at Lawrence Livermore National Laboratory) was employed in this work to simulate the dislocation activities.
In the 3D-DDD, the dislocation activities are controlled by a series of governing rules. The curved dislocation lines are discretized into a set of nodes connected by straight segments, and the velocities and forces on the segments are considered. The position of the segment is represented by vector r(s, t), in which s and t denote the location along the line and simulation time, respectively. Multiple categories forces are exerted into the dislocation segment, such as self-force, dislocations interaction force, externally applied force, and drag force. The integrated governing equation of dislocation interaction and movement is as follows: where g(·) represents the operational procedure and can be calculated by g ¼ M F tot drive r s ð Þ; σ applied ; … Â Ã À Á . To conduct the simulation of the microstructural responses of AISI 4140 steel under high-strain rate loading, a novel box geometry of 24000 b × 24000 b × 24000 b (b is the length of the Burgers vector) was developed. On the basis of previous experimental observations [36], random dislocations with the order of 1.7 × 10 12 m −2 were used as the initial computational condition. A relaxing process was performed to allow the simulation system reaching equilibrium configuration before imposing uniaxial compression ofε =10 6 s −1 . The 3D-DDD simulation results of dislocations evolution and spatial distribution within the representative cell are shown in Fig. 2. The morphologies of dislocations at different physical simulation times (from 3 to 6 ns) indicated that the dislocation density in the unit cell increased massively under the high strain rate. The high shear strain rate activated cross-slip on all possible slip systems of AISI 4140 steel. Figure 3 shows the dislocation density and engineering flow stress as a function of strain. The increase of dislocation density, in principle, results in the subsequent processes of plastic deformation and recrystallization. The yield point at the strain level of about 0.219% was 413 MPa, from which the flow stress-strain curve start deviating from the linear elastic response. The predictive value is consistent with the steady-state value of 415 MPa for AISI 4140 steel as reported in [37]. The dislocation density increases after the yield point, and a parabolic-like work hardening feature was observed. With the further increase of dislocation density, it reached a critical condition of dynamic recrystallization (DRX). In DRX process, the dislocation density increases dramatically to 0.241 × 10 14 m −2 when the stain increases to 0.6%, being one order of magnitude higher than that in activated point. The material is then in a stage of large plastic deformation. Therefore, the DDD results can be used to calculate the related dislocation coefficients, i.e., α * and k 0 , of DDB constitutive equation. Table 2 is a summary of the constants of DDB constitutive equation for AISI 4140 steel. Based on the DDD calculation, the nucleation rateρ nuc c is 1.431 × 10 22 m −2 s −1 , and the average annihilation rateρ ann c is 4.667 × 10 19 m −2 s −1 , and the average cell interior dislocation density ρ c is 1.427 × 10 13 m −2 . Theγ c andγ w are assumed to be equal to the resolved shear strain rate, and the reference strain rateγ 0 is 10 5 s −1 [38]. The dislocation-related coefficients α * and k 0 for AISI 4140 steel are 0.515 and 3.425 at 25°C according to Eqs. (1) and (3). The constant β * is set as 0.0017. Other parameters, not demonstrated here, are all derived from previous experiments and conclusions [39][40][41], as list in Tables 2 and 3. The flow  (15). Figure 4 compares the stress-strain curves predicted by DDB constitutive equations with those predicted by classic Johnson-Cook (JC) constitutive equation under various thermal conditions, and a good accordance has been achieved. By further calibrating the above curves at various temperatures, the k 0 was determined to be temperature-dependent as follows: where T is the cutting temperature (°C). When the cutting temperature increases from ambient 25 to 600°C, k 0 would increase from 3.425 to 5.536. Figure 5 is the data flow chart for the multiscale modelling process. Commercial FEM software with the explicit module was employed to build the cutting model, which can simulate the chip formation and surface generation during the high-performance manufacturing process. The 8-node thermally coupled displacement and temperature finite elements were used to divide the workpiece mesh. Mesh convergence study has been carried out to obtain the best possible values by referring to chip formation and computational efficiency.  Minimum element edge lengths for chip bulk and cutting edge were 0.0183 mm and 0.0102 mm, respectively. The cutting tool properties specified in the simulation were based on an uncoated tungsten carbide tool. The workpiece was assumed to be fixed on the back and bottom sides, and the tool was allowed to move horizontally from the right to the left while keeping vertically. To simulate the thermo soften effects observed in laser-assisted machining, in this cutting model, the initial temperatures of workpiece can be set to the respective temperatures. The detailed simulation parameters have been summarized in Table 4. .

Numerical validation of DDD-FEM simulation framework
To demonstrate the effectiveness and stability of the proposed DDB-FEM model, a FEM cutting model using the classic JC constitutive equation was performed with the same cutting parameters, and the simulation results from the two models were compared in terms of stress distribution, cutting force, and cutting heat. Figure 6 a and b show the stress distribution predicted by the DDB-FEM model at cutting speed of 40 m/min and workpiece temperature 25°C. The continuous chip forming indicated a localized, highly deformed stress concentration zone. Significant stress gradient was found not only in the tool-workpiece contact zone but also in the steady-state machined surface. Those results are in excellent agreement with the simulation results predicted by the reference JC constitutive equation (Fig. 6c, d). A slight difference was found on the peak value of von Mises Stress, which were both less than 1% and 2.7% in the initial stage and stabilization stage, respectively. The calculated cutting temperatures predicted by different constitutive equations are shown in Fig. 7. As expected, the temperatures predicted by the two models are in a similar distribution and with only a slight difference of 0.2% and 2.8% in the peak values. The cutting force is another important indicator for the high-performance machining process. The cutting force can be adopted to demonstrate the integral simulation state relative to the transient stress or temperature distribution features. As shown in Fig. 8, the curves of tangential force and feed force approached steady stages after the sharp increases at the initial cutting stage, and then fluctuating around constant values. The cutting forces from different constitutive equations are in a high consistency at all times. The agreement between the results predicted by two models indicates that the proposed multiscale numerical framework is stable and does not change the material property in macro scale.

Effect of cutting speed
To identify the particulars of microstructure transformation for chip and machined surface of AISI 4140 steel, the multiscale modelling framework was used to demonstrate the distribution and evolution of dislocation densities together with grain sizes in cutting. For a more detailed elucidating of dislocations accumulation and its contribution occurring in the localized Fig. 6 Stress distribution predicted by DDB and JC constitutive equations (v = 40 m/min, T w = 25°C). a DDB prediction at the cutting distance of 0.06 mm. b DDB prediction at the cutting distance of 0.7 mm. c JC prediction at the cutting distance of 0.06 mm. d JC prediction at the cutting distance of 0.7 mm Fig. 7 Comparison of cutting temperatures predicted by DDB and JC constitutive equations (v = 40 m/min, T w = 25°C). a DDB prediction at the cutting distance of 0.06 mm. b DDB prediction at the cutting distance of 0.7 mm. c JC prediction at the cutting distance of 0.06 mm. d JC prediction at the cutting distance of 0.7 mm deformation regimes, the microstructure transformations in the primary shear zone (PSZ), second shear zone (SSZ), and machined surface zone (MSZ) were analyzed. The entire chip deformation and dislocation densities distribution are shown in Fig. 9a, and the average dislocation densities at PSZ, SSZ, and MSZ are shown in Fig. 9b-d, respectively. The prediction of dislocation density reached the order of 10 9 mm −2 , which is in good agreement with the experimental observations by Gonza lez et al. in [42]. It was found that, except for the average dislocation density at MSZ, the average densities at PSZ and SSZ were both increased with the increase of cutting speed. The activities of discrete dislocations were mainly accompanied by energy dissipation, while the increase of cutting speed will result in the large local strain rate and deform Average dislocation density at PSZ. c Average dislocation density at SSZ. d Average dislocation density at MSZ energy at the core shear zones. Therefore, a mass of dislocation dipoles at PSZ and SSZ were generated and released to compromise the severe plastic deformations in those areas. In contrast, increased machining speed restricted the dislocation accumulation at MSZ, and further alleviating the machined surface hardening phenomenon.
The influence of cutting speed on the dynamic recrystallization process is shown in Fig. 10. The grains were severely deformed near the main shear zones. The minimum grain size was found to be 0.338 μm which is similar to the grain size characterized in selected shear localized regions by Courbon et al. [8]. A long narrow and highly recrystallization band was identified at PSZ as shown in Fig. 10a. The recrystallization phenomena at SSZ were severer than those at PSZ. This is due to the relatively severe impact contributed by the cutting tool.
The accumulation of the cutting energy contributed to the generation of non-equilibrium microstructure within the coarser grains so that the continuous recrystallization phenomena occurred in the tool-workpiece contact zone. On the basis of the Hall-Petch strengthening rule, excessive refined grains would aggravate strain hardening and result in early tool wear/ damage as well as degraded machined surface quality.

Effect of external-assisted temperature
Laser-assisted machining has been reported to improve the machinability of hard metals effectively. The effect of external-assisted temperature on material deformation was studied with the proposed model. Figure 11 a is the representative picture showing the dislocation densities distribution and microstructure characteristics at cutting speed of 120 m/ min and workpiece temperature of 350°C. The increased workpiece temperature inhibited the dislocation multiplication, and a decrease of the average dislocation density was observed at PSZ and MSZ (Fig. 11b). The reason is that the increased workpiece temperature boosted enormously dislocation activities so that enough plastic strain could be produced by a small number of dislocations to meet the significant plastic deformation demands. The high workpiece temperature facilitated the cell dynamic recovery and decelerated the continuous recrystallization as well (Fig. 12). The average grain size went up from 0.389 to 0.589 μm at PSZ when the workpiece temperature increased from 25 to 400°C. An important notice is that the averaged grain sizes were stable in chip and machined surface when the workpiece temperatures were kept between 200 and 350°C. The results indicate that a dynamic equivalent status has been achieved for the applied strain, strain rate, and cutting temperature.
To find the optimal range of processing parameters, the combined effects of various processing conditions on the machined surface integrity were further studied. Figure 13 shows the grain size evolution as a function of cutting speeds and external-assisted temperatures. It is found that a lower cutting speed and a higher workpiece temperature result in a coarser grain size after machining. The maximum mean value of grain size was about 0.636 μm for the tested conditions. Under the three applied cutting speeds, the grain sizes increased by approx. 60% when the workpiece temperature increased from 25 to 400°C. This indicates that the external-assisted temperature is the major factor for the growth of grains. Indeed, the large grain size indicated the thermal soften effect of AISI 4140 steel in the main shear zones, thereby contributing to an approximate 20% lower cutting force as indicated in Fig.  13d. An important notice is that with the increase of cutting speed, the refined grains can facilitate the strain hardening effect. However, there is a sign of slowing for grain refinement phenomenon when the cutting speed is higher than 280 m/min. The mean values of grain size tend to stable with the further increase of cutting speed, indicating an optimal range for the selection of cutting speed.

Conclusions
In crystals, slip induced by discrete dislocation evolution is identified as the dominant mechanism of plastic deformation in most circumstances. To establish a physical link between microstructural alteration and the strain hardening phenomenon, a DDD-FEM multiscale simulation framework for high-performance manufacturing was presented in this work. The real-time tracking of dislocations multiplication, interaction, propagation, and annihilation, was conducted in the DDD calculation platform. These statistic microstructural evolution features were adopted to optimize the dislocation-related coefficients in the DDB constitutive equation. A flow stress model incorporating the physics of dislocation accumulation, dynamic recrystallization, and dynamic recovery was developed to describe the material deformation and capture the microstructure characteristics in the FEM module. The present framework has been validated through comparisons with classic simulation results (chip morphology, stress-strain distribution, temperature distribution, and cutting force) with a good consistency. Significant dynamic recrystallization phenomena have been found in the primary shear zones, and the grain size at the machined surface is affected by the external-assisted temperature and cutting speed. A stable microstructural organization was identified while the machining speed is larger than 280 m/min, or the external-assisted temperatures lie between 200 and 350°, indicating an optimal range of processing conditions for improving machined surface integrity. The present multiscale cutting model will be further extended to a wide range of metals and verified by the electron backscatter diffraction (EBSD) experimental tests for the evolution of grain size and misorientation angle of chip and machined surface in the future work.

Ethics approval and consent to participate
Research participants were not subjected to harm in any ways whatsoever. The respect for the dignity of all the research participants had been prioritized, and full consent was obtained from the participants prior to the research study.

Consent for publication
The authors consent to Springer the non-exclusive publication rights and warrant that their contribution is original and that they have full power to make this grant.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.