Modeling and Numerical Simulation of Material Science
Vol.09 No.02(2019), Article ID:91759,12 pages

Calculation of the Stress Intensity Factor in an Inclusion-Containing Matrix

Claude Lincourt1, Jacques Lanteigne2, Madhavarao Krishnadev3, Carl Blais3

1UTC Pratt & Whitney Canada, Longueuil, Canada

2Hydro-Quebec Research Institute, Varennes, Canada

3Department of Mining, Metallurgical and Materials Engineering, Laval University, Quebec, Canada

Copyright © 2019 by author(s) and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

Received: February 13, 2019; Accepted: April 9, 2019; Published: April 12, 2019


The intent of this paper is to propose an engineering approach to estimate the stress intensity factor of a micro crack emerging from an inclusion in relation with the morphology of the inclusion and its relative stiffness with the matrix. A micromechanical model, based on the FEA (finite element analysis) of the behavior of cracks initiated at micro structural features such as inclusions, has been developed using LEFM (Linear Elastic Fracture Mechanics) to predict the stress intensity factor of a micro crack emerging from an inclusion. Morphology of inclusions has important connotations in the development of the analysis. Stress intensity factor has been estimated from the FEA model for different crack geometries. Metallographic analysis of inclusions has been carried out to evaluate the typical inclusion geometry. It also suggests that micro cracks less than 1 µm behave differently than larger cracks.


Inclusions, Finite Element Analysis, Micro-Crack

1. Introduction

Fracture toughness of high strength steels is controlled by the onset of crack initiation, which depends on the local stress distribution in the vicinity of a micro structural stress raiser. Crack growth then occurs in a limited region, whose extent depends on the local stress variation resulting from material heterogeneities. The knowledge of the stress condition for a micro crack emerging from a heterogeneity is then of prime importance to determine its sensitivity to microstructure.

2. Materials and Experimental Procedure

2.1. Materials

Different steels were used to measure typical inclusion size to be used to build the FEA model. The steels used in this metallographic study are presented in Table 1 for reference purpose. All these steels are commonly used in general industrial applications. Steel A has been subjected to rare earth inclusion shape control practice and presents a microstructure with fine and rounded inclusions as shown in Figure 1. Steels B to J are spring steels and present typical elongated manganese sulfide inclusions as seen on Figure 2.

2.2. Inclusion Analysis and Calculation of Stress Concentration Factor

Metallographic was used to characterize the inclusions. Polished and lightly etched specimens were used. A Clemex image analyzer linked to a Nikon Ephiphot optical microscope has been used to evaluate the morphology and distribution of inclusions. Inclusions in the transverse direction were considered due to their elongated shape. The stress concentration factors at inclusions have been calculated using the following formula [1] :

K t = 1 + 2 c ρ (1)

Kt = stress concentration factor

c = semi-major axis of the inclusion

ρ = radius of the notch

For simplicity, this equation may be written as:

K t = 1 + 2 L t (2)

Table 1. Chemical composition of the steels used in this study.

Figure 1. Round inclusion (rare earth shape control).

Figure 2. Elongated inclusion.

L = total length of the inclusion

t = thickness of the inclusion

The results are shown in Table 2.

3. Calculating the Stress Intensity Factor of Micro-Crack Using Finite Element Analysis

A 2D finite element mesh idealization of a 40 mm by 80 mm steel bar is shown in Figure 3(a). In Figure 3(b), the central portion is magnified. This mesh refinement is intended to represent a sequence of voids in the fracture plane (Y = 0), either circular or elongated, depending on the selected element material. For example, in Figure 4 and Figure 5, grey colored elements are high Modulus elements, as compared to low Modulus elements in white. Therefore, Figure 4 is

Table 2. Typical inclusion Kt for the steels used in this study.

(a) (b)

Figure 3. (a) General view of the model; (b) Portion of the model at a micro-scale level.

Figure 4. Detailed view of a round void.

Figure 5. Detailed view of an elongated void.

the representation of a circular void surrounded by the metal matrix while Figure 5 represents a void elongated along the X-axis of the model. This FEA model consists of 20,831 6-node triangles and 8-node quadrilateral elements, both types quadratic.

Materials are linear elastic and matrix elements are assumed to be 200,000 times as rigid as void elements. Plain strain situation is assumed and uniform normal (+Y) stress is applied at the remote nodes illustrated in Figure 3(a). The remote stress is normalized per unit thickness. Restriction of nodes against Y motion is applied at the specimen middle line Y = 0, which is the fracture plane, except at cracked element nodes. Owing to this plane of symmetry, the effective specimen size is actually 40 mm by 160 mm.

The no-crack, single circular void situation was first analyzed. Mesh was further refined at the tip of the void until convergence of the local stress was attained. The process finally ended with the smallest element size equal to 0.023 µm in Figure 4 and Kt = 3.03, which is close to the stress concentration factor Kt = 3 for a circular hole in an infinite plate [2] .

Cracked geometries were analyzed using the same mesh refinement. Restriction against X motion was applied at a single node of the model (generally at the crack tip) to prevent the X rigid body motion. Different crack lengths were analyzed by detaching nodes at the crack tip along the plane of symmetry. For each crack configuration, a complete solution was performed. Depending on the boundary conditions applied, along with element stiffness, the crack length becomes the independent variable of the analysis, and various asymmetrical cracked geometries can be analyzed.

3.1. Simulation of Voids and Inclusions

Based on results of micro structural analysis, it was decided to simulate two extreme inclusion morphologies, namely, the round and the elongated shapes. For round a void, the radius was set at 1 µm (refer to Figure 4). For elongated voids, two round voids were combined with their centers set 22 µm apart, thus creating a 24-µm stringer (refer to Figure 5). This void shape resembles the elongated manganese sulfide (MnS) stringers that are commonly observed in steels. This assessment is based on quantitative metallography.

3.2. Stress Concentration Factor Kt

The stress concentration factor is estimated as follows:

K t = σ max σ (3)

where σ max is the maximum stress calculated with FEA and σ is the remote stress.

Analytical expressions can be found for the stress concentration factor for the two types of void and, therefore, may be compared to FEA with no-crack configurations. As mentioned earlier, a round void behaves as a circular hole in an infinite medium, for which Kt = 3 [2] . This study found a convergent solution with Kt = 3.03 (a difference of only 1%). For the elongated void shown in Figure 5, the FEA result is Kt = 8.26. Rice [3] proposed the following equation for the stress concentration factor of a flat type void:

K t = 2.43 c ρ (4)

which gives Kt = 8.42 with c = 12 µm and ρ = 1 µm as in Figure 5. Here again the discrepancy between FEA and analytical results is very small (a difference of about only 2%, refer to Table 3).

Simulation of inclusions poses some questions. First, it is well known that typical inclusions in steel are considered weakly bonded to the matrix [4] while, in FEA models, the displacement at an interface is fully compatible between the two materials. Similar questions have been dealt with in FEA simulation of graphite nodules in cast iron [5] . In this particular case, the nodules have been considered to behave like voids because of the weak bond and the low stiffness of the graphite nodules. In this study, inclusions were considered as voids in the matrix. Attributing an extremely low stiffness value to an element of material data is a way to introduce a void in the model where otherwise it would be a steel element, thus allowing the mesh to represent various configurations.

Table 3. Stress concentration factor calculated with FEA model.

3.3. Stress Intensity Factor KI

Stress intensity factors were determined using the quadratic element mid-side node re-localized at the quarter position. Such crack elements were located along the path of the crack. The only effective one, in any simulation, is the one located at the crack tip. The stress intensity factor was calculated from [6] , assuming that the state of stress is plane strain:

K = E 4 ( 1 ν 2 ) ( 2 π L ) ( 4 u B u C ) (5)

where uB is the crack opening displacement at the quarter node, uC is the crack opening displacement at the corner of the element (position C), E is the Young’s Modulus, ν the Poisson’s ratio and L is the element size.

The stress intensity form factor, defined as:

F ( a ) = K I σ π a (6)

with KI the SIF in mode I, as calculated from Equation (5), σ , the remote stress and a, the crack length, was calculated from FEA results.

4. Results and Discussion

4.1. Stress Intensity form Factor

The stress intensity form factor for round and elongated voids were calculated using Equation (6). The results are presented in Table 4.

With the result of Table 4, we may then plot the graph of the geometrical correction curve against crack length. The results are presented in Figure 6.

As it can be seen from this figure, the F curve exhibits a sudden exponential drop at the first three steps of propagation, after which it stabilizes as the crack behaviour becomes only dependent on the crack length. One particularity of this curve is the relatively high value of F for very short cracks.

4.1.1. Comparison with the Results of Breiznitskii [7] and Tweede-Rooke [8] for Round Voids

In this section, we want to compare the F function given by our FEA model with models from literature. The two references chosen are from Breitznitskii and Tweede-Rooke.

The model of Breitznitskii [7] is shown in Figure 7.

Table 4. Geometrical correction factor F against crack length “a” for round and elongated voids.

Figure 6. Geometrical correction curve against crack length for round void/inclusions.

Figure 7. Model from Brezhnitskii [7] .

The second model used here for purpose of comparison is proposed by Tweede-Rooke [8] and is defined by:

K I = σ π a F I (7)


FI = geometrical correction factor (in mode I)

The geometrical factor FI is defined as:

F I = F 1 / 2 F 2


F 1 / 2 = ( 2 + ε ) / ( 2 + 2 ε ) [ 1 + ( 0.2 ε ) / ( 1 + 3 ε ) 3 ] ; ε = a / R ,

where R is the void radius


F 2 = 1 + [ 1 / ( 2 ε 2 + 1.93 ε + 0.539 ) + 1 / 2 ( ε + 1 ) ]

Using the applicable formula for each of these models, we may calculate the geometrical correction factors for all five analyzed crack configurations. These are presented in Table 5 and compared to the FEA results. Figure 6 is a plot against the crack length a.

Though FEA results well agree with these models at large crack lengths, a significant discrepancy is observed when it is shorter than 0.2 micron.

At the shortest crack configuration analyzed (0.023 micron), F from FEA is at least eight times greater then that obtained from references while it reduces to twice at 0.098 micron and to almost nothing at 0.199 microns, where both F are almost equal. Hence, based on our FEA result, the F factor for a microscopic crack less than 0.2 micron emerging from a round void (with a radius of 1 micron) is significantly higher than the one predicted by the models.

4.1.2. Comparison with the Results of Breiznitskii [7] for Elongated Voids

We will use here the same approach and the same model proposed by Breitznitskii as presented in the previous section, the only difference being that the stress concentration factor Kt of the elongated void is higher. Now if we use the same formula given by the Breitznitskii model, we can calculate the geometrical correction factor for a micro-crack. These values are presented in Table 6 with the F factor calculated with the FEA model for an elongated void.

Figure 8 presents the curve of the geometrical correction factor F against crack length “a” plotted with the results from Table 6.

As for a round void, this figure reveals a significant difference for short cracks between the F factor calculated with FEA and the reference one. This difference is observed for the first 0.6 micron of growth after which the two curves merge together. The span over which the difference exists with the Breitznitskii model is longer than for the round void just because here the stress concentration factor is significantly higher. In comparison with the Breitznitskii model, the F factor calculated with the FEA model is at 0.023 micron of crack extent three times

Table 5. Geometrical correction factor F calculated from FEA and the models from Breitznitskii [7] and Tweede-Rooke [8] for a round void.

Table 6. Geometrical correction factor F calculated from FEA and the model from Breitznitskii for an elongated void.

Figure 8. Comparison between the geometrical correction factor F versus crack length calculated with FEA and reference models from Breitznitskii [7] and Tweede-Rooke [8] for a round void.

greater while it is 50% more at 0.2 micron and become finally equal around 0.6 micron of crack growth. The difference in the F functions between the round and the elongated void is definitely due to the higher stress concentration factor magnitude for the elongated void. Except for the first propagation step where the F values are almost equivalent in both cases, the stress intensity factor needed to propagate a crack emerging from an elongated void will require a lower remote stress than for the round void (Figure 9).

4.2. The Effect of Stress Raiser end Radius on the Geometrical Correction Factor F

For very short cracks, FEA yielded SIF significantly higher than that predicted from LEFM [7] [8] [9] . For example, a typical LEFM calculation [9] gave a form

Figure 9. Comparison between the geometrical correction factor F versus crack length calculated with FEA and the reference model from Breitznitskii [7] for an elongated void.

factor F(a, Kt) of 3.3 at a 0.023 µm crack emerging from a 2 µm diameter round void while the FEA result was 27. For crack lengths above 0.2 µm, however, both methods led to practically the same results. The same observation holds for elongated voids (with Kt = 8.34) for which a form factor of 28 was obtained. Therefore, it seems that, for a very short crack, the radius of the void appears to be the controlling SIF parameter, whereas, beyond 1 µm of propagation, this phenomenon disappears and the crack tip singularity combined with the stress concentration factor takes over.

In view of this unexpectedly steep F function seen at very short cracks, SIF were alternatively estimated from the Hellen-Parks [10] virtual crack extension method. With this method, the J-Integral is first estimated and hence the SIF is derived. These numerical SIF values, either based on the crack tip opening displacement formula of Leslie-Banks [6] or on the potential energy difference of Hellen-Parks, differed by only 8%, thus providing the same estimation whatever the calculation technique.

5. Summary and Conclusions

It has been determined that FEA can actually be applied at a micro structural level to describe the behavior of a micro-crack emerging from micro structural features that act as stress raisers. The major findings are as follows:

1) Stress intensity factors for micro-cracks emerging from metallurgical stress concentrations calculated from plane strain finite element analysis (FEA) shows higher values in comparison with models from literature.

2) The analysis suggests that micro-cracks smaller than 1 µm behave differently from that of larger cracks. For a very short crack, the radius of the void appear to be the controlling SIF parameter, whereas, beyond 1 µm of propagation, this phenomenon disappears and the crack tip singularity combined with the stress concentration factor takes over.


The authors are grateful to NSERC (Natural Sciences and Engineering Research Council of Canada) for the financial support provided. They also express their sincere thanks to Pratt and Whitney Canada, IREQ (Research Institute of Quebec Hydro) and Timken (USA), for the provision of materials and services. They are also grateful to Dr Real Bouchard, CANMET for mechanical testing and Mme. Maude Larouche, Laval University for metallographic.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Lincourt, C., Lanteigne, J., Krishnadev, M. and Blais, C. (2019) Calculation of the Stress Intensity Factor in an Inclusion-Containing Matrix. Modeling and Numerical Simulation of Material Science, 9, 17-28.


  1. 1. Knott, J.F. (1973) Fundamentals of Fracture Mechanics. Butterworth, London.

  2. 2. Peterson, R.E. (1974) Stress Concentration Factors. John Wiley & Sons.

  3. 3. Rice, J.R. (1968) Fracture, an Advance Treatise. Academic Press.

  4. 4. Venkatasubramanian, T.V. and Baker, T.J. (1982) Role of Elongated MnS Inclusions in Hydrogen Embrittlement of High-Strength Steel. Metal Science, 16.

  5. 5. Dahlberg, M. (1997) Micromechanical Modeling of Nodular Cast Iron, a Composite Material. International Journal of Cast Metals Research, 9, 319-330.

  6. 6. Bank-Sills, L. and Bortman, Y. (1984) Reappraisal of the Quarter-Point Quadrilateral Element in Linear Elastic Fracture Mechanics. International Journal of Fracture, 25, 169-180.

  7. 7. Brezhnitskii, L.T. (1966) Propagation of a Crack Terminating at the Edge of a Curvilinear Hole in a Plate. Soviet Materials Science, 2, 16-23.

  8. 8. Tweed, J. and Rooke, D.P. (1973) The Distribution of Stress near the Tip of a Radial Crack at the Edge of a Circular Hole. International Journal of Engineering Sciences, 11, 1185-1195.

  9. 9. Tada, I., Paris, P.C. and Irwin, G.R. (1973) The Stress Analysis of Cracks Handbook. Del Research Corporation, Hellertown.

  10. 10. Parks, D.M. (1974) A Stiffness Derivative Finite Element Technique for Determination of Crack Tip Stress Intensity Factors. International Journal of Fracture, 10, 487-502.