Numerical Study of the Effect of Inclusions on the Residual Stress Distribution in High-Strength Martensitic Steels During Cooling

In high-strength martensitic steels, the inclusions significantly affect the material performance especially in terms of fatigue properties. In this study, a numerical procedure to investigate the effect of the inclusions types and shapes on the residual stresses during the cooling process of the martensitic steels is applied systematically based on the scanning electronic microscopy (SEM) and energy dispersive spectrometer (EDS) results of different types of inclusions. The results show that the maximum residual stress around the interface between Mg-Al-O inclusion and the matrix is the largest, followed by TiN, Al-Ca-O-S, and MnS when the inclusions are assumed as perfect spheres for simplicity. However, these results are proved to be 28.0 to 48.0% inaccurate compared to the results considering actual shapes of inclusions. Furthermore, the convex shape of inclusion will lead to stress concentration in the matrix while the concave shape of inclusion will lead to stress concentration in the inclusion. The residual stress increases with the increase of inclusion edge angle. The increase rate is the largest for TiN inclusions on the concave angle, which leads to extreme stress concentration inside TiN inclusion.


Introduction
Driven by the increasing demand on the higher mechanical performance of engineering structures, high-strength steels are steadily developed and widely applied in recent decades in multiple areas, e.g., automotive, high-speed trains, and aerospace. The high strength is mainly obtained by the in-depth and complicated design of the microstructure of steels, which also leads to new damage mechanisms for multiphase steels [1][2][3][4][5]. In terms of fatigue properties, a fatigue limit is normally evident for the conventional steels developed some decades ago [6][7][8]. Murakami et al. [9] reported that the fatigue limit increases with the increase of steel hardness and decreases with the increase of the defect area. Therefore, increasing the strength of materials is leading to a higher hardness and eventually a higher fatigue strength limit assuming the same defect size. However, for the recently developed high-strength steels, it is found that even when stress amplitude is lower than fatigue limit, the fatigue crack still will occur when fatigue life reaches a certain level [9][10][11][12]. For higher stress and lower fatigue life, fatigue crack usually occurs from the component surface; for lower stress and longer fatigue life, fatigue crack usually occurs from inner microstructure [13]. The deformation heterogeneity in microstructure level, which is responsible for the damage and failure in general for steels, can be generally caused by two different mechanisms: the matrix inhomogeneity and the matrix-inclusion inhomogeneity [14,15]. Many researchers have focused on the effect of matrix inhomogeneity on mechanical properties, their damage mechanism with experiments and modelling. Lian et al. investigated the damage mechanism of a DP600 steel under different stress states experimentally [16] and via micromechanical simulation [17] and found that damage is initiated at the boundary of ferrite and martensite due to their strong deformation incompatibility. It was further revealed that the damage initiation is by nature stress state dependent due to the strain localization dependency of stress states. Krupp et al. [18] studied the fatigue crack initiation process of duplex stainless steel and quantified the interactions of cracks with the first microstructural barrier. Alharbi et al. [19] investigated the stress state of martensite and ferrite before cracking of dual phase DP1000 steel and estimated for crack initiation in the martensite. He et al. [20] reported that for a ferritic-pearlitic steel, the initiation of the cleavage fracture is located in the pearlitic phase where the ferrite and the cementite show a lamella structure with different properties.
In addition to steel matrix inhomogeneity, inclusions also play a significant role during the failure process. Fairchild et al. [21] reported the influence of inclusions on cleavage initiation and the relation between inclusion crack and steel matrix crack. Yang et al. [22] studied 42CrMo steel with different amounts of inclusions and found that the fatigue property was better with fewer inclusions. Sakai et al. [23] found that internal inclusions could lead to the formation of the penny-shape fine granular layer around inclusions during long-term cyclic loadings, which would cause small cracks inside the steel matrix. While investigating inclusions in steel, it is also found that there are multiple types of inclusions in steels [24,25] generated in different processes during the steel production [26,27]. For example, the Al 2 O 3 inclusion is produced during melting of steel while the TiN inclusions are produced during the casting of steels. Besides, the shape, size, and composition of inclusions, affected by the production method, also differ, leading to different effects on the mechanical properties of steels [28,29].
The production processes additionally induce another level of inhomogeneity, residual stresses, based on the microstructure-level inhomogeneity. For the steel matrix, it is found that the mechanical deformation, such as rolling [30] and indentation or surface treatment [31] could introduce the grain-level residual stresses and they are grain orientation dependent. For the interaction between the steel matrix and inclusions, the thermal treatment, e.g., cooling, could induce large local residual stress, which would affect the mechanical properties of the steel significantly [32]. These residual stresses are generally generated due to the different thermal expansion properties of inclusions and steel matrix during temperature variation. Brooksbank and Andrews [33,34] developed a model to calculate the residual stress around spherical and cylindrical inclusions and reported the residual stresses around the boundaries between different inclusions and steel matrix were different. This model was widely applied by researchers to calculate residual stresses around inclusions in different steels [35,36]. Murakami and Uchida [37] claimed that the residual stresses induced by different thermal expansion properties around Al 2 O 3 inclusions will lead to debonding of inclusions from steel matrix based on calculations, which seriously affects the mechanical properties of steel. In addition to common inclusions, carbides, which are close to the surface, will also cause irregular stress distribution [38]. In terms of finite element-based modelling, there are researchers simulating the residual stresses around inclusions and analyze their effects on mechanical properties. Pineau and Forest [39] performed an elastic-plastic finite element calculation to determine the residual stresses in the vicinity of the inclusions. The results show that the residual equivalent plastic flow occurs all around the inclusion during cyclic loading, which plays an essential role in crack initiation. Gu et al. [40,41] proposed a microstructure-based model considering residual stress induced by temperature variation. This model offers a more arcuate prediction on crack initiation site during fatigue process compared to the model neglecting residual stress. Although all these above-mentioned studies offer good illustrations of the numerical residual stress distributions around inclusions and their strong effects on the mechanical behaviour, the gap is also evident: Actual inclusion shape is not considered. For simplicity, these models and calculations only considered the inclusion as a simple geometrical shape, such as spherical and cylindrical. For actual inclusions in steel, the shape varies dramatically and contains multiple details, which will greatly affect the residual stress distribution around inclusions. In terms of experiments, quantifying residual stress around inclusions still faces challenges due to the redistribution of residual stress while machining specimens. Therefore, it is aimed in this study to apply a numerical procedure to investigate the effect of the inclusions types and shapes on the residual stresses systematically during the cooling process of the martensitic steels. It is emphasized that by using a numerical approach, the shape effect and the effect of the inclusions type in terms of the mechanical properties can be distinguished and individually studied. The ideal spherical shape is first assumed for all the inclusions, and the numerical results are compared with the analytical model. Furthermore, models were established with the actual shapes of inclusions in steel, and the residual stress distributions were analyzed numerically to clarify the relation between shape variation and residual stress distribution of different inclusions. A further quantitative comparison study of the residual stress result between the inclusions with ideal shape and the actual shape was conducted. This study provides the basis for further researches on numerical prediction of crack behaviour and cracks mechanism modelling for high-strength steels.

Materials and Methods
The material in the present study is high carbon martensitic bearing steel with 1.12% C and 1.41% Cr. The main chemical composition is shown in Table 1. This material was hot rolled into a steel bar with a diameter of 35 mm, vacuum oil quenched after holding for 20 min at 835 • C and tempered for 120 min at 180 • C. The mechanical properties were tested with the specimen represented in Figure 1. The specimen was perpendicular to the cross section of the steel bar. According to recent studies [42][43][44][45], the manufacturing methods of the specimen can affect the results of the tensile test. In this study, the cutting method of a wire electro discharge machining was adopted. With this method, the geometric-dimensional requirements can be easily met without any significant change in the mechanical properties. The extension rate of the tensile test was 1.0 mm/min. A metallographic specimen was also taken from the steel bar and polished by SiC paper and diamond suspensions in the cross-section of steel bar. Nonmetallic inclusions were observed with scanning electronic microscopy (SEM) and the compositions of inclusions were analyzed with energy dispersive spectrometer (EDS).

Modelling and Simulation
The geometrical profiles of typical inclusions in the steel observed with SEM were captured as coordinates with MATLAB (version 2016a, MathWorks, Natick, MA, USA, 1984). Due to the differences in contrast between inclusion and steel matrix in the SEM picture, the MATLAB codes can assist to recognize the details of inclusion shape in a desirable way. Thereafter, the coordinates of inclusion shape were processed by a Python script using the Abaqus (Dassault Systèmes Simulia Crop., Providence, RI, USA, 1994) extension of the Python programming language to create the inclusion component of the Abaqus model. The version of Abaqus is Abaqus 2017. The steel matrix was set as a square area (30 × 30 μm 2 ) and assumed to be a homogeneous material with a single phase and elastic only. The inclusions were positioned in the middle of the matrix area so that the residual stress around the inclusion was not affected by the matrix interface.
The simulated process was the cooling procedure during oil quenching. The temperature variation was from 835 °C to the room temperature, which was set at 20 °C. The temperature was defined with pre-defined field in Abaqus/Standard. The thermal expansion coefficients and the mechanical property parameters of different inclusions and steel matrix were obtained from references and experiments in the present study.
During the cooling process, the thermal stress of every node during the elastic deformation in the inclusion-matrix system was calculated with Equation (1).
where i is the material definition number, which represents steel matrix and different inclusions; αi is the thermal expansion coefficient; T0 is the initial temperature, and T is the current temperature.
To avoid over constraining the thermal expansion of the system and guarantee the accuracy of the simulation, it is important to define proper interface conditions. In the present study, the left side of the square steel area was fixed in the x-direction and still can move in the y-direction. The element type is set as "continuum-plane stress-3 nodes (CPS3)". The simulation type is static.

Modelling and Simulation
The geometrical profiles of typical inclusions in the steel observed with SEM were captured as coordinates with MATLAB (version 2016a, MathWorks, Natick, MA, USA, 1984). Due to the differences in contrast between inclusion and steel matrix in the SEM picture, the MATLAB codes can assist to recognize the details of inclusion shape in a desirable way. Thereafter, the coordinates of inclusion shape were processed by a Python script using the Abaqus (Dassault Systèmes Simulia Crop., Providence, RI, USA, 1994) extension of the Python programming language to create the inclusion component of the Abaqus model. The version of Abaqus is Abaqus 2017. The steel matrix was set as a square area (30 × 30 µm 2 ) and assumed to be a homogeneous material with a single phase and elastic only. The inclusions were positioned in the middle of the matrix area so that the residual stress around the inclusion was not affected by the matrix interface.
The simulated process was the cooling procedure during oil quenching. The temperature variation was from 835 • C to the room temperature, which was set at 20 • C. The temperature was defined with pre-defined field in Abaqus/Standard. The thermal expansion coefficients and the mechanical property parameters of different inclusions and steel matrix were obtained from references and experiments in the present study.
During the cooling process, the thermal stress of every node during the elastic deformation in the inclusion-matrix system was calculated with Equation (1).
where i is the material definition number, which represents steel matrix and different inclusions; α i is the thermal expansion coefficient; T 0 is the initial temperature, and T is the current temperature.
To avoid over constraining the thermal expansion of the system and guarantee the accuracy of the simulation, it is important to define proper interface conditions. In the present study, the left side of the square steel area was fixed in the x-direction and still can move in the y-direction. The element type is set as "continuum-plane stress-3 nodes (CPS3)". The simulation type is static.

Properties of Inclusion and Steel Matrix
The stress vs. strain curve of the martensitic steel is shown in Figure 2 based on the tensile tests described in Section 2. According to the figure, the mechanical properties of the steel matrix are shown in Table 2.
The stress vs. strain curve of the martensitic steel is shown in Figure 2 based on the tensile tests described in Section 2. According to the figure, the mechanical properties of the steel matrix are shown in Table 2.  inclusions usually precipitate in the grain interface as spherical, cylindrical or irregular. All these different shape features will contributed to the variation of residual stress distribution around inclusions, which will be discussed in Section 4.3. The thermal expansion coefficients and the mechanical properties of different inclusions, as well as the steel matrix, are presented in Table 3, which are the required material parameters of the simulation in the following sections.    There are four main types of inclusions in the martensitic steel: spinel (Mg-Al-O), calcium aluminum (Al-Ca-O-S), titanium nitride (TiN), and manganese sulfide (MnS). The typical morphologies are shown in Figure 3. The Mg-Al-O inclusions are usually polygon. The Al-Ca-O-S inclusions are usually spherical. The TiN inclusions are usually irregular in shape. The MnS inclusions usually precipitate in the grain interface as spherical, cylindrical or irregular. All these different shape features will contributed to the variation of residual stress distribution around inclusions, which will be discussed in Section 4.3. The thermal expansion coefficients and the mechanical properties of different inclusions, as well as the steel matrix, are presented in Table 3, which are the required material parameters of the simulation in the following sections. The stress vs. strain curve of the martensitic steel is shown in Figure 2 based on the tensile tests described in Section 2. According to the figure, the mechanical properties of the steel matrix are shown in Table 2.  inclusions usually precipitate in the grain interface as spherical, cylindrical or irregular. All these different shape features will contributed to the variation of residual stress distribution around inclusions, which will be discussed in Section 4.3. The thermal expansion coefficients and the mechanical properties of different inclusions, as well as the steel matrix, are presented in Table 3, which are the required material parameters of the simulation in the following sections.

Residual Stress Distribution Around Spherical Inclusions
The residual stress distribution around inclusions was already calculated by Brooksbank and Andrews [33]. In their model, inclusions were simplified as perfect spherical. The radial stress in the surface between inclusion and steel matrix can be calculated with Equations (2) and (3).
where α M and α are the coefficients of linear expansion of the steel matrix and the inclusion, respectively; v M and v are the Poisson's ratios of the steel matrix and the inclusion, respectively; E M and E are the Young's modulus of the steel matrix and inclusion, respectively; ∆T is the difference between the holding temperature before vacuum quenching and the room temperature; R is the radius of the inclusion, and R M is the radius of the steel matrix around a single inclusion, which is half of the distance between two inclusions with the same type.
To clarify the influence of inclusion type on residual stress distribution and make a comparison of the calculation results in the present study with the model of Brooksbank and Andrews, perfect spherical Mg-Al-O, Al-Ca-O-S, TiN, and MnS inclusions were inserted in the square steel matrix. The geometry model and the calculation results of residual stress distribution are shown in Figure 4. It is obvious that the residual stress in the steel matrix around the inclusion decreases when the distance to inclusion surface increases. The position of the maximum stress is on the surface between the inclusion and steel matrix. When the inclusion is spherical, the stress on the surface between inclusion and steel matrix keeps constant around the inclusion. This constant maximum stress varies with inclusion type. The value of Mg-Al-O inclusion is the largest, followed by TiN and Al-Ca-O-S inclusions while the value of MnS inclusion is the smallest. These simulation results are also compared with the calculation results with Equations (2)

Residual Stress Distribution Around Actual Inclusions
When inclusions are assumed as spherical, the shape details are actually neglected. Nevertheless, as shown in Figure 3, the inclusions of different types usually possess different geometry shapes, which will lead to important stress concentration around interfaces between the inclusions and the steel matrix. In this Section, four types of inclusions in Figure 3 were inserted in the square steel area, and the residual stress distributions were simulated. The simulation results are shown in Figure 6. Generally, the residual stress also decreases when the distance to inclusion surface increases, which is consistent with the case of the perfect spherical inclusion, while the maximum residual stresses for these four types of inclusions all increased compared to the maximum stress of spherical inclusions. Besides, the residual stress on the surface between the inclusion and the steel matrix changes with the shape of inclusions.

Residual Stress Distribution Around Actual Inclusions
When inclusions are assumed as spherical, the shape details are actually neglected. Nevertheless, as shown in Figure 3, the inclusions of different types usually possess different geometry shapes, which will lead to important stress concentration around interfaces between the inclusions and the steel matrix. In this Section, four types of inclusions in Figure 3 were inserted in the square steel area, and the residual stress distributions were simulated. The simulation results are shown in Figure 6. Generally, the residual stress also decreases when the distance to inclusion surface increases, which is consistent with the case of the perfect spherical inclusion, while the maximum residual stresses for these four types of inclusions all increased compared to the maximum stress of spherical inclusions. Besides, the residual stress on the surface between the inclusion and the steel matrix changes with the shape of inclusions. To quantitatively analyze the effect of the inclusion shape on residual stress changes, two parameters were introduced: the edge angle θsv,i and the local residual stress σsv,i. These two parameters represent the smooth of the inclusion edge and the corresponding local residual stress for each θsv,i. The calculation method of θsv,i is shown in Equation (4) and the schematic is shown in Figure 7.
where i is the identification number of digital point on the inclusion edge. All the digital point were sequenced continually in an anticlockwise direction. θi is the angle between the x-axis and the straight line through the (i-1) th and the i th digital point on the inclusion edge, which can be calculated with Equation (5).
where xi and yi are the coordinates of the digital i th point on the inclusion edge. When θi is a positive value, the local part of inclusion is convex; when θi is a negative value, the local part of inclusion is concave. The smaller the absolute value of θsv,i is, the smoother the inclusion edge is. When θsv,i is To quantitatively analyze the effect of the inclusion shape on residual stress changes, two parameters were introduced: the edge angle θ sv , i and the local residual stress σ sv,i . These two parameters represent the smooth of the inclusion edge and the corresponding local residual stress for each θ sv , i . The calculation method of θ sv , i is shown in Equation (4) and the schematic is shown in Figure 7.
where i is the identification number of digital point on the inclusion edge. All the digital point were sequenced continually in an anticlockwise direction. θ i is the angle between the x-axis and the straight line through the (i-1) th and the i th digital point on the inclusion edge, which can be calculated with Equation (5).
where x i and y i are the coordinates of the digital i th point on the inclusion edge. When θ i is a positive value, the local part of inclusion is convex; when θ i is a negative value, the local part of inclusion is concave. The smaller the absolute value of θ sv , i is, the smoother the inclusion edge is. When θ sv , i is zero, the local shape of the inclusion is a straight line. For a right angle, θ sv , i is 90 • . It is to be noted that, the distance between each neighboring point will affect the result of θ sv , i . Therefore, the distance between each neighboring point is kept consistent in this study which offers a value of 3.6 • for θ sv , i in a perfect circle.
Appl. Sci. 2018, 8, x FOR PEER REVIEW 9 of 13 zero, the local shape of the inclusion is a straight line. For a right angle, θsv,i is 90°. It is to be noted that, the distance between each neighboring point will affect the result of θsv,i. Therefore, the distance between each neighboring point is kept consistent in this study which offers a value of 3.6° for θsv,i in a perfect circle.  Figure 6 were extracted and calculated. The relation between these two parameters is shown in Figure 8. To avoid the interaction of neighboring points, the data with θsv,i lower than 5° were neglected. As shown in the figure, the local residual stress σsv,i increases with the absolute value of edge angle θsv,i for all four types of inclusions, which indicates that no matter the local shape of the inclusion is convex or concave, the sharp edge will lead to stress concentration. When the local shape is convex, the stress concentration is in the steel matrix side; when the local shape is concave, the stress concentration is in the inclusion side (see Figure 6).   Figure 6 were extracted and calculated. The relation between these two parameters is shown in Figure 8. To avoid the interaction of neighboring points, the data with θ sv , i lower than 5 • were neglected. As shown in the figure, the local residual stress σ sv,i increases with the absolute value of edge angle θ sv , i for all four types of inclusions, which indicates that no matter the local shape of the inclusion is convex or concave, the sharp edge will lead to stress concentration. When the local shape is convex, the stress concentration is in the steel matrix side; when the local shape is concave, the stress concentration is in the inclusion side (see Figure 6). zero, the local shape of the inclusion is a straight line. For a right angle, θsv,i is 90°. It is to be noted that, the distance between each neighboring point will affect the result of θsv,i. Therefore, the distance between each neighboring point is kept consistent in this study which offers a value of 3.6° for θsv,i in a perfect circle.  Figure 6 were extracted and calculated. The relation between these two parameters is shown in Figure 8. To avoid the interaction of neighboring points, the data with θsv,i lower than 5° were neglected. As shown in the figure, the local residual stress σsv,i increases with the absolute value of edge angle θsv,i for all four types of inclusions, which indicates that no matter the local shape of the inclusion is convex or concave, the sharp edge will lead to stress concentration. When the local shape is convex, the stress concentration is in the steel matrix side; when the local shape is concave, the stress concentration is in the inclusion side (see Figure 6).  θ sv , i . However, for the concave part of inclusion, large θ sv , i on the surface of TiN inclusion will lead to extremely large stress concentration on the inclusion side. This is also the reason that TiN inclusion itself often cracks in steel. These cracks inside inclusions can extend into steel matrix during the service of steel, causing failure of component [4]. Therefore, to obtain better performance it is important to avoid sharp edges of TiN inclusion.
To clarify the effect of the shape on the maximum residual stress around different types of inclusions, the maximum residual stresses of perfect spherical inclusions and inclusions with actual shape are shown in Figure 9. The maximum residual stresses of perfect spherical inclusions of all types are obviously smaller than the cases considering inclusion shapes. When we assume the inclusion shape as a perfect sphere for simplicity, there will be a non-negligible error of around 28.0 to 48.0%.
In addition, θsv,i, inclusion type also plays a significant role in the effect on local residual stress variation with θsv,i. For the convex part of inclusion, the local residual stress of MnS inclusion increases the fastest with θsv,i among these four types of inclusions. When there is a convex angle of 25° on the inclusion surface, the local residual stress will increase about 1.64 times for MnS, 1 time for Al-Ca-O-S, 0.75 times for Mg-Al-O, and 0.53 times for TiN. For concave part of inclusion, the local residual stresses of Mg-Al-O and TiN inclusions increase the fastest, followed by MnS and Al-Ca-O-S. When there is the concave angle of 25° on the inclusion surface, the local residual stress will increase about 0.48 times for Al-Ca-O-S, 0.51 times for MnS, 1.32 times for Mg-Al-O, and 1.25 times for TiN. Although the local residual stress increases the fastest for convex part of MnS inclusion, the concentrated residual stress on the steel matrix side still does not exceed the other three types of inclusion under the same value of θsv,i. However, for the concave part of inclusion, large θsv,i on the surface of TiN inclusion will lead to extremely large stress concentration on the inclusion side. This is also the reason that TiN inclusion itself often cracks in steel. These cracks inside inclusions can extend into steel matrix during the service of steel, causing failure of component [4]. Therefore, to obtain better performance it is important to avoid sharp edges of TiN inclusion.
To clarify the effect of the shape on the maximum residual stress around different types of inclusions, the maximum residual stresses of perfect spherical inclusions and inclusions with actual shape are shown in Figure 9. The maximum residual stresses of perfect spherical inclusions of all types are obviously smaller than the cases considering inclusion shapes. When we assume the inclusion shape as a perfect sphere for simplicity, there will be a non-negligible error of around 28.0 to 48.0%.

Conclusions
The cooling process during heat treatment will lead to residual stress around the interface between inclusion and steel matrix due to different thermal expansion coefficients and mechanical properties. In this study, the effect of the inclusion types and shapes on the residual stresses during the cooling process of the martensitic steels was investigated with a numerical procedure. The following results were obtained: (1) This residual stress varies with inclusion type and shape. For perfect spherical inclusion, the residual stress around the interface between Al-Mg-O inclusion and steel matrix is the largest, followed by TiN and Al-Ca-O-S inclusion while the residual stress around MnS inclusion is the smallest.

Conclusions
The cooling process during heat treatment will lead to residual stress around the interface between inclusion and steel matrix due to different thermal expansion coefficients and mechanical properties. In this study, the effect of the inclusion types and shapes on the residual stresses during the cooling process of the martensitic steels was investigated with a numerical procedure. The following results were obtained: (1) This residual stress varies with inclusion type and shape. For perfect spherical inclusion, the residual stress around the interface between Al-Mg-O inclusion and steel matrix is the largest, followed by TiN and Al-Ca-O-S inclusion while the residual stress around MnS inclusion is the smallest.
(2) When the inclusion shape is assumed as a perfect sphere for simplicity, there will be a non-negligible error of the maximum residual stress of around 28.0 to 48.0% compared with the residual stresses considering the actual shapes.
(3) The convex local shape of inclusion will lead to stress concentration in the steel matrix side, and the concave local shape of inclusion will lead to stress concentration in the inclusion side.
(4) The local residual stress of MnS inclusion increases the fastest with the increase of inclusion edge angle θ sv , i for the convex part of inclusion. For the concave part, the local residual stress of TiN inclusions increases the fastest, which leads to extreme stress concentration inside TiN inclusion and contributes greatly to material failure.

Conflicts of Interest:
The authors declare no conflict of interest.