## 1 Introduction
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. [PERSON] and [PERSON] [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]. [PERSON] and [PERSON] [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. [PERSON] and [PERSON] [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. [PERSON] 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.