Influence of hard and soft inclusions inside a ferritic matrix

. In the case of an interior fracture mode in fatigue, material impurities play a fundamental role in the generation and the propagation of cracks. In general, the presence of hard or soft inclusions, voids, materials imperfections, etc., alter the local stress state in the matrix generating the accumulation of plastic deformations which lead to the crack formation. In terms of inclusions, high strength steels are often characterized by the presence of aluminum oxides Al 2 O 3 or Manganese sulfide MnS. The experimental works reported by several authors pointed out that the critical location for the crack nucleation is, often, at the edge of the inclusion and it subsequently propagates on a plane orthogonal to the loading direction. In a previous work, the authors investigated the stress distribution around a spherical inclusion inside a ferritic matrix, pointing out the role of the material anisotropy and the different crystallographic orientation of the matrix. However, the investigations dealt with simple loading conditions and linear elasticity. The present paper aims to extend the field of investigations to the elasto-plastic domain, focusing the attention on the role of the crystallographic orientation and comparing the results with the previous study.


Introduction
The fatigue life of high strength steel components is highly influenced by the presence of materials defects (voids, inclusions, heterogeneity in the material matrix, etc.), which often represent the locus for the micro-crack formation. Several studies [1][2][3] investigated the role of the material defects in the crack formation, focusing the attention on the individuation of the main factors responsible for the limitation of the fatigue life. An exhaustive overview by Murakami [4] individuated the defect characteristics that influence the most the fatigue life of metallic components in: type of defect (inclusion, pore, shrinkage, etc.), shape (spherical, ellipsoidal, complex, etc.), position (internal, surface, etc.), size (related or not to the loading direction).
Alfredsson and Olson [3] presented a work where the fatigue initiation is predicted by means of the Finley multi-axial critical plane criterion [5]. Their investigations supported 1 the conclusions carried out by other authors that inclusions, with a higher Young's modulus (i.e. Al2O3) than the surrounding matrix, are generally more harmful than soft inclusions (i.e. CaS, MnS). In their paper, Alfredsson and Olson carried out both experimental and numerical observations, focusing the attention on inclusions of spherical shape, motivating the choice by the experimental evidence. Another interesting work was carried out by Vincent et al. [1] extending the work presented in a previous paper [6]. In [1], the authors presented an upgrade of the Defect Stress Gradient (i.e. DSG) approach, using the Equivalent Inclusion Method (EIM) [7], and considering cavities of different shapes inside a metallic matrix. Their work showed how the EIM can be successfully used to predict the stress around the defect, saving the computation time required by the FEM. Two loading conditions were analysed: uniaxial tension loading and torsional loading, with the defect perfectly centred inside the specimen. Cerullo and Tvergaard [2] analysed the effect of small inclusions on the fatigue life of wind turbine bearings by means of micromechanical analyses. A modified version Dang Van fatigue criterion [8,9] was adopted to individuate the load level below which the material failure is avoided. In detail, the analyses were carried out by means of 2D FE simulations under the assumption of plane strain. Two types of inclusion were considered: a circular Al2O3 particle and a squared TiN particle, positioned at the centre of the sample. Cerullo and Tvergaard correlated the volume fraction occupied by the inclusion with the damage-like parameter derived by the modified Dang Van criterion, finding that a higher damage factor was always reached at the interface between the inclusion and the matrix, regardless the volume fraction of the inclusion. The possible explanation for this result was due to the difference in the elastic modulus between the particles and the surrounding matrix. In [10], the authors carried out a series of numerical investigations by means of FE, at a mesoscale level, on a cylindrical sample with a spherical inclusion at the centre. Two different types of inclusion were considered: a harder Al2O3 particle and a softer MnS inclusion, both of spherical shape. In addition, anisotropic elasticity was considered for the ferritic matrix and isotropic elasticity was assumed for the inclusions. The field of investigations was limited to the elastic regime and to uniaxial loading conditions. The present paper aims to extend the material response to the elasto-plastic regime considering a crystal plasticity framework at finite strain. Moreover, the loading conditions are here enriched considering the superimposition of a hydrostatic component of the stress in order to evaluate the effect of the mean stress.
The paper is organized as follows. Section 2 presents the theoretical framework and the set of equations used for the numerical analyses. The constitutive equations were implemented via user subroutine for the commercial code ABAQUS. Section 3 deals with the numerical results. An initial overview of the sample in terms of geometry, mesh, material parameter and boundary condition is shown. Subsequently, the results are organized by loading condition, starting from the uniaxial extension, and then considering the new loading conditions. At last, the conclusions are reported.

The theoretical framework
Following the formulation proposed by of Asaro and Lubarda [11], the description of the material deformation in a elasto-plastic crystallographic framework is introduced by defining the deformation gradient    F φX . Where ( , ) t φX is a map associating the spatial position in the current configuration at the time t with the material position X in the reference configuration. Then, it is assumed that the total deformation gradient can be divided into two parts: the elastic one e F (responsible for the distortion of the crystal lattices) and the plastic deformation due to crystallographic slip p F , which can be 2 combined multiplicatively as ep  F F F . The initial configuration is a configuration before the motion, the current one is considered after the total motion whereas the intermediate configuration refers to a configuration excluding the elastic deformation.
By the multiplicative decomposition, the velocity gradient can be expressed as: where d and w are the rate of deformation tensor and the spin tensor respectively. The definition of the previous quantities can be done focusing on the single contributes.
The crystallographic slip is activated when the resolved shear stress on each slip surface of the slip system reaches the critical value (i.e. CRSS). Since the constitutive theory for the crystal plasticity has been well established, it will be briefly recalled in this section, accordingly to the theory presented by Asaro and Lubarda [11] and by Pierce, Asaro and Needleman [12]. It is assumed () 0  s the unit vector defining the slip direction and () 0  m the unit vector that defines the normal direction of the slip surface at the initial or reference configuration for the  slip system. Considering the crystal lattice as nondeformable during the plastic slip it is possible to write the previous vectors and the resolved shear stress as: ; : Where J is the determinant of the deformation gradient and σ is the Cauchy stress. Considering the total plastic power as the sum of the plastic power of the single slip systems, the last expression in Eq. (4) can be used to write the definition of the plastic rateof-deformation and the plastic spin tensors as follows: Here ()  μ and ()  ω are the symmetric and skew-symmetric components of   On the other hand, the definition of the crystallographic slip has been assumed in accordance with Asaro [13]: Where a is the reference slip system, () g  is the slip resistance of the  system, where n controls the strain rate dependency, h  is the latent hardening moduli and h  is the selfhardening, function of the initial hardening modulus 0 h and of 0  and s  , respectively the initial CRSS and the saturation stress.
Where represents the elasticity constants tensor, * σ is the corotational stress rate evaluated in the intermediate configuration and whose definition can be given using the Jaumann rate of the Cauchy σ stress as follows: (8). Combining Eqs. (7) and (8) (9).
In the previous expression, the definition for the elastic transformation and plastic rates of deformation and spins can be obtained inverting Eqs.  dd μ w w ω (10).

Numerical analyses
The set of constitutive equations, briefly introduced in section 2, was implemented via user subroutine for the commercial code Abaqus (ver. 6.14-4). The FE analyses carried out in this work deal with a series of loading tests on a metallic sample containing a spherical material inclusion. Two different non-metallic materials were considered for the inclusions Al2O3 and MnS, respectively. The effect of the size, as well as, different shapes or positions of the inclusion inside the matrix were not considered in this study. The attention is focused on understanding the stress localization around the inclusion that leads to the formation of plastic deformations as a function of the difference of Young's moduli between the matrix and the inclusion. Moreover, the influence of the elastic anisotropy for the ferritic matrix is also taken into consideration. It has to be mentioned that a perfect bond between the inclusion and the matrix is assumed. This aspect requires future investigations since it represents a strong assumption on the influence of the inclusion. Partial debond or a weaker bond could, in fact, significantly alter the stress concentrations around the inclusions.

Mesh, geometry and material parameters
The sample consists of a cylindrical ferric body with a spherical inclusion at the centre. The ratio between the sphere radius and the cylinder radius is 0.20, giving a volume fraction of 0.005. According to Cerullo and Tvergaard [2], increasing sizes of the volume fraction of the inclusion reduced the damage factor estimated with the modified Dang Van criterion, probably due to the fact that a bigger inclusion reduces the stress peak for the matrix in the reference cell. Future works will investigate this aspect.
The focus of the paper is to evaluate the effect of the Young's modulus ratio between the ferritic matrix and the inclusion, therefore, two non-metallic inclusions are considered: a hard inclusion (i.e. Al2O3) and a soft inclusion (i.e. MnS). The cylinder is constrained at the base along the direction y and rotations are prevented. A prescribed displacement along y is imposed on top of the cylinder, giving a maximum nominal axial strain of 0.25% (see Fig. 1a). Different values of uniform pressure are applied on the lateral and top surface of the sample, in order to study the effect of a hydrostatic component of the stress (see Fig. 4). In detail the pressure values are: -100 MPa, -50 MPa, 0 MPa, +50 MPa and +100 MPa. The mesh consists of 126144 hexahedral elements with reduced integration (i.e. Abaqus element C3D8R), for a total amount of 135501 nodes (see Fig. 1b). For sake of comparison, the ferrite has been considered alternatively isotropic and anisotropic to evaluate the influence of an anisotropic elastic tensor. The elastic constants are assumed accordingly to Tjahjanto et al. [14] for the anisotropic elastic tensor and for the crystal plasticity constants (Table 1). In detail, perfect plasticity was assumed. The constants for the isotropic behaviour of the ferrite and the Al2O3/MnS inclusions were taken from [15] (see Table 2).

Uniaxial tensile test
The first set of analyses does not consider the effect of a superimposed hydrostatic pressure and the only loading condition is a prescribed displacement along the y-direction applied on top of the cylinder. The response of the sample is checked in terms of nominal stress vs. nominal strain curves, varying the crystallographic orientation of the ferritic matrix. The sets of Euler angles are the same as in [10], due to the highly symmetric configuration of the BCC structure. The same approach is adopted for the analyses in the subsequent sub-section 3.3. Perfect elasto-plastic behaviour is considered since the interest of the work is to check whenever the first plastic slip takes place. In fact, the generation of local plastic deformations in the surrounding of the inclusion gives useful information on the fatigue life in terms of nominal stress and crack location site. As a comparison, the graphs in Fig. 2a and b report an additional dashed curve obtained considering the sample without inclusion and assigning to each element a different crystallographic orientation, in order to simulate a polycrystal stress-strain curve. The first plastic slip in the polycrystal analysis is indicated with a yellow marker on the curve.
The first step is to compare the material response without the inclusion considering isotropic elasticity or anisotropic elasticity (i.e. black lines in Fig. 2a and Fig. 2b, respectively). Comparing Fig. 2a with Fig. 2b, it is possible to notice that the elastic anisotropy varies the inclination of the elastic parts of the stress-strain curves, depending on the crystallographic orientation, which consequently modify of the size of the elastic domain in terms of strain. This last aspect can be better understood by observing Fig. 2b. Here, the black markers indicate the nominal yield stresses which, compared with Fig. 2a, take place at different nominal strains.
The effect of the harder and softer inclusions can be observed in Fig. 2b. The differences in Young's modulus between the ferritic matrix and the Al2O3 or MnS inclusions alter the stress field in the surrounding of the inclusion, generating some stress localization. This triggers the activation of some slip systems, giving a different response compared with the analyses without the inclusion. Fig. 3 summarises the results obtained in the analyses. In detail, the graph reports the normalized nominal stress at the first local plastic slip σy against the nominal yield stress of the isotropic elastic analyses without the inclusion (i.e. Fig. 2a). It can be concluded that the Al2O3 inclusions are in general more harmful since they trigger the generation of plastic deformation at lower stress regime, confirming the conclusion reported in [2,3]. However, depending on the Euler angles set, the differences in nominal stress between the Al2O3 and the MnS inclusions seem to become smaller when the BCC cubic structure is rotated to give an alignment between cubic diagonal and the pulling direction. In case of perfect alignment, the MnS inclusion seems to be more harmful (see Fig. 3).

Uniaxial tensile test with hydrostatic pressure
The last set of analyses modify the boundary conditions by means of uniform pressure on the surfaces of the sample. Four different values were investigated: -100MPa, -50 MPa, +50MPa and -100MPa. The magnitude of the confining pressures was chosen to be maximum as the 0  in Table 1. The following Fig. 4 reports a sketch of the new loading conditions. A prescribe displacement condition was applied on top of the cylindrical sample as in section 3.2. Moreover, the same crystallographic orientations were again considered. The following Fig. 5 summarizes the results obtained in the numerical analyses. In detail, the two contour fields report the axial stress, at the first plastic slip around the inclusions, normalized against the yield stress of the homogeneous cylinder without inclusion. In the y-axis are displayed the different crystallographic orientations whereas the x-axis reports the different values of the confining pressures. Observing the two graphs it is possible to notice the following main aspects. The hard inclusion can be considered the most harmful for the material with the exception of the (35.26˚,45˚,0˚) orientation, confirming again the results of section 3.2. Moreover, the MnS inclusion seems to be more sensitive against the confining pressure for the orientations (0˚,0˚,0˚), (0˚,15˚,0˚) and (0˚,30˚,0˚), showing a wider range of normalized nominal stress along the x-axis. In general, a tensile hydrostatic pressure (i.e. -50 and -100 MPa) results in early plastic slips for both the inclusions.

Conclusions
The present work aimed to investigate the effect of material inclusions inside a metallic matrix in order to understand the mechanism that leads to the formation of localized plastic deformations and subsequently to the material failure during the service life of structures or components. The main results can be summarized in the following points. The crystallographic orientation affects the formation of the plastic deformation (see Fig. 2a).
The anisotropy of the elastic tensor influences the size of the elastic domain in terms of strain (see Fig. 2a). Moreover, the presence of an inclusion alters the stress state in the surrounding matrix causing the formation of plastic slips. In general, the Al2O3 inclusion induces plastic slips at lower nominal stress compared with the MnS inclusion. This tendency varies with the Euler angles. In particular, when the cubic diagonal of the BCC structure is aligned with the pulling direction, the MnS particle triggers the formation of plastic deformations at lower nominal stress compared with the Al2O3 inclusion (see Fig. 3). The MnS inclusion seems to be more sensitive to the variation of confining pressures. In general, the addition of a tensile hydrostatic stress component facilitates the formation of plastic slips. Vice versa, the addition of a compressive hydrostatic stress component delays the formation of plastic slips (see Fig. 5).