Thermographic evaluation of metal crack propagation during cyclic loading

Journal Pre-proofs Thermographic evaluation of metal crack propagation during cyclic loading B. Hajshirmohammadi, M.M. Khonsari PII: DOI: Reference:

S0167-8442(19)30306-4 https://doi.org/10.1016/j.tafmec.2019.102385 TAFMEC 102385

To appear in:

Theoretical and Applied Fracture Mechanics

Received Date: Revised Date: Accepted Date:

6 June 2019 7 October 2019 9 October 2019

Please cite this article as: B. Hajshirmohammadi, M.M. Khonsari, Thermographic evaluation of metal crack propagation during cyclic loading, Theoretical and Applied Fracture Mechanics (2019), doi: https://doi.org/ 10.1016/j.tafmec.2019.102385

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

Thermographic evaluation of metal crack propagation during cyclic loading B. Hajshirmohammadi and M. M. Khonsari1 Department of Mechanical and Industrial Engineering Louisiana State University 3283 Patrick Taylor Hall Baton Rouge, LA 70803, USA Abstract

An analytical relationship between the stress intensity factor and temperature rise for a moving crack tip is proposed. Comparison of the predictions with experimental results shows good agreement in both trend and magnitude. The results reveal the potential of applying theromographical technique for assessing the behavior of crack growth.

Nomenclature 𝜔(𝜃) 𝑊𝑝 𝑇 𝑇0 𝜎𝑖𝑗(𝑛,𝜃) 𝜖𝑖𝑗(𝑛,𝜃) ∆𝜎

Absolute ambient temperature (K) HRR non-dimensional functions for stress distribution HRR non-dimensional functions for strain distribution Stress range (MPa) Elastic modulus (GPa) Strain range

𝐼𝑛′ 𝑟 ∆𝐾𝐼 𝜔0

Cyclic plastic strain hardening exponent HRR dimensionless coefficient Distance to crack tip(cylindrical coordinate) (m) First mode stress intensity factor (MPa.m1/2) Plastic region radius at angle 0 (m)

𝑎

Crack length (m)

𝜈

Poisson’s ratio

𝑡 𝑘 𝜆

𝜎′0

Cyclic strength coefficient (MPa) Cyclic yield stress (MPa)

𝜉

Moving coordinate

𝜎𝑒𝑞

Equivalent stress (MPa)

𝑣

Crack moving speed (m/s)

𝜀𝑝𝑒𝑞

Equivalent strain

𝛽

Frequency

𝛼

Ratio of Heat dissipated to plastic work Dimensionless constant

𝑘′

𝑓

Corresponding author [email protected]

1|Page

𝑛′

Specimen thickness (m) Thermal conductivity (W/mK) Thermal diffusivity (m2/s)

𝐸 ∆𝜖

1

Plastic region radius at angle θ (m) Cyclic plastic strain energy density (J/m3K) Absolute temperature (K)

1. Introduction Prediction of fatigue failure has long been an area of interest in the design of mechanical components that experience cyclic loading with rich volumes of research papers from various fronts. In particular, fatigue crack growth (FCG), as an inextricable part of fatigue failure, has captured the attention of many investigators since the end of the 20th century in an attempt to determine a relationship between the applied load and the propagation rate. While models are available that relate these parameters, the problem remains untouched when the loading of components under fatigue changes unpredictably during the operation. To this end, thermographic technology, as a non-destructive approach, can provide insight into the physical processes that take place during FCG. In the present study, a method is presented to predict a relationship between the crack propagation and the temperature rise near a moving crack tip by utilizing an infrared (IR) camera. Many research papers have provided models for quantifying crack propagation rate during the cyclic loading. Paris and Erdogan [1] were the first researchers who used the stress intensity factor to determine a relationship between the loading condition, crack length and the propagation speed. Their method was later modified by Walker [2] who considered the effect of mean stress. When a cracks start to propagate, a cyclic plastic zone is formed ahead of its tip. Due to the presence of the cyclic plastic strain in this area, heat generation is concentrated on the crack tip with a noticeable temperature rise in its vicinity. A part of the plastic dissipation, which is the source of this heat generation, is converted to heat. According to Zehnder et al [3], the ratio of the dissipated plastic energy which transforms to heat depends on a number of factors such as the yield strain, plastic strain, and the strain hardening of the material. According to Kallivayalil [4] and Wong and Kirby [5], the ratio of plastic dissipation to that converted to heat is in the range of 0.8 to 0.95. More recently, thermal imaging techniques have been used to attain more clear insight into the propagation mechanism. For example, Meneghetti et al. [6-9] have studied energy dissipation around the crack and found a relationship between energy dissipation and propagation rate in a specimen made of stainless steel. Studies pertaining to other materials such as Titanium alloys are reported by Vshivkov et al. [10, 11]. Another noteworthy example is the work of Zhang et al. [12] who investigated the cyclic plastic zone and found four different stages of propagation in 4003 ferritic stainless steel. Plastic dissipation mechanism was first proposed by Rice [13] as a method to formulate FCG and was solved using FEM by Klingbeil [14]. Later, Smith [15] studied the tensile overload in stainless steel 304 by applying the dissipated energy criterion to predict the crack propagation rate. The application of thermography for fatigue study is becoming widespread [16-23]. Recently, Huang et al. [24] proposed new methods for rapid determination of fatigue limit using thermal analysis. Manquin et al.[25] studied the initial stage of fatigue loading thermally and showed that in small stress levels, the heat generated was almost constant with cycle number. Yang et al. [26] used thermography for rapid evaluation of the reliability of high-cycle fatigue. Haung et al. [27] showed the usefulness of infrared measurements in life prediction of CFRP composites. Jürgen et al.[28] found thermoelastic effect to be useful in determination of dissipated energy of fatigue. 2|Page

In the present paper, the method introduced by Pandy et al. [29] is used to derive a formula between FCG and the temperature rise at the crack tip during cyclic loading. This approach is based on the classic work of HRR (Hutchinson-Rice-Rosengren) [30, 31] on the singularity solution for stressstrain angular distribution around the crack tip. The predictions of the present model are experimentally verified using stainless steel SS304 and SS316.

2. Theory and formulation According to the Ramberg-Osgood [32], the stress and strain behavior of metals can be described by the following relationship (1)

1/𝑛′

( )

∆𝜖 ∆𝜎 ∆𝜎 = + 2 2𝐸 2𝑘′

where

∆𝜖

and

2

∆𝜎 2

are the range of strain and stress amplitudes, 𝑛′ represents the cyclic plastic strain

hardening exponent, and 𝑘′ is cyclic strength coefficient. For materials that obey the power law hardening [32] the crack tip stress-strain field can be evaluated using the theory put forward by H.R.R. [30, 31]. (2)

𝑛′

(

𝛥𝜎𝑖𝑗 = ∆𝜎0

∆𝐾2𝐼

)

𝑛′ + 1

𝜎𝑖𝑗(𝑛′,𝜃)

𝛼𝛥𝜎20𝐼𝑛′𝑟

1

𝛥𝜖𝑖𝑗 =

(

)

𝛼∆𝜎0

∆𝐾2𝐼

𝐸

𝛼𝛥𝜎20𝐼𝑛′𝑟

𝑛′ + 1

(3) (

𝜖𝑖𝑗 𝑛′,𝜃

)

where 𝜎𝑖𝑗 and 𝜖𝑖𝑗 are non-dimensional functions of 𝑛′ and 𝜃. In this equation, 𝐼𝑛′ is an integration constant, 𝜎0 is the cyclic yield stress (∆𝜎0≅2𝜎0) and 𝐸 is the modulus. The parameter 𝛼 is a dimensionless constant given by Eq. 4. 𝛼=

2𝐸 1

(4) 𝑛′ ― 1

(2𝑘′)𝑛′𝛥𝜎0 𝑛′ where 𝑟 and 𝜃 are polar coordinates shown in Fig. 1 with the crack tip located on the origin (𝑟 = 0).

3|Page

𝑦

𝑦

𝑃

𝑠 𝛾

𝐵 𝑅

𝑟

𝜔(𝜃)

𝜃

𝑥

𝛹

𝐴

𝑟𝑝

𝑥

𝜃 𝜔(𝜃)

(a)

(b)

Figure 1. Cyclic plastic zone in front of the crack using (a) von Mises criteria and (b) simple circle as the plastic region. As it is shown in Fig.1, The plastic zone size, which can be found according to the von Mises criterion (Fig 1 (a)), can be replaced by a circle as illustrated in Fig.1 (b). Using the equivalent stress and plastic strain, one can arrive at the following relationship to describe the cyclic plastic strain energy per unit volume per cycle. 𝑊𝑝 =

( ) 𝑛′ ― 1 𝑛′ + 1

(5)

𝜎𝑒𝑞𝜀𝑝𝑒𝑞

Substituting Eqs. 2 and 3 in Eq. 5, yields the following equation for heat generation ahead of the crack tip in terms of the frequency f. 𝑛′ ― 1 ∆𝐾2𝐼 𝜎𝑒𝑞(𝑛′,𝜃)𝜀𝑒𝑞(𝑛′,𝜃) 𝑊𝑝 = 𝑓𝛽 𝐸𝐼𝑛′𝑟 𝑛′ + 1

( )

(6)

where 𝛽 represents the portion of plastic work converted to heat and ∆𝑘 is the range of stress intensity factor defined as 𝑘𝑚𝑎𝑥 ― 𝑘𝑚𝑖𝑛. Parameters 𝜎𝑒𝑞(𝑛′,𝜃) and 𝜀𝑒𝑞(𝑛′,𝜃) represent the equivalent stress and plastic strains defined as follows. 1

𝜎𝑒𝑞 = (

𝜎2𝑟

𝜀𝑝𝑒𝑞 =

+

𝜎2𝜃

2 2 (𝜀 + 3 𝑟

)

+ 𝜎𝜃𝜎𝑟 + 3𝜎2𝑟𝜃 2 1 2 2 2 𝜀𝜃 ― 2𝜀𝑟𝜀𝜃 + 3𝜀𝑟𝜃

)

(7) (8)

where 𝜎𝑟,𝜎𝜃 and 𝜎𝑟𝜃 are angular distribution functions of stress defined by Hutchinson et al. [30] and 𝜀𝑟,𝜀𝜃 and 𝜀𝑟𝜃 are angular distribution function of strain around the crack tip in plain strain condition. If the plastic region is estimated by a circle of radius 𝑟𝑝, then boundary of the plastic zone can be defined by Eq. 9. 𝜔(𝜃) = 2𝑟𝑝cos (𝜃)

where the plastic zone radius 𝑟𝑝 is given by the following expression [33]. 4|Page

(9)

𝑟𝑝 =

1 ∆𝐾2𝐼 8𝜋 𝜎20

2.1.

(10)

Thermal field around the crack

The equation governing the heat conduction with provision for heat source 𝑊𝑝 is: ∂2𝑇 ∂𝑥

+ 2

∂2𝑇

(11)

𝑊𝑝 1∂𝑇 + = 𝑘 𝜆 ∂𝑡 ∂𝑦2

where 𝑇 is temperature, 𝜆 is the thermal diffusivity, 𝑘 is the conductivity, and 𝑡 represents time. For a moving heat source with the velocity 𝑣 the heat transfer equation is represented by Eq. 12 ∂2𝑇

𝑣∂𝑇 ∂2𝑇 𝑊𝑝 + + =― 2 2 𝑘 𝜆 ∂𝑡 ∂𝜉 ∂𝑦

(12)

In this relationship, the 𝜉 coordinate is attached to the tip of crack and moves with the same speed. That is, (𝜉 = 𝑥 ― 𝑣𝑡). Referring to Fig. 1, the differential temperature rise at point 𝑃 due to the point heat source 𝐵 is found by Eq. 13. 𝑑𝑇 =

𝑊𝑝 2𝜋𝑘

(

𝑒𝑥𝑝 ―

(13)

) ( )

𝑣𝑠 𝑣𝑠 cos 𝛾 𝐾0 2𝜆 2𝜆

where 𝐾0 is the modified Bessel function of the second kind and zero order. The total temperature rise due to all the point sources in plastic zone can be found by integrating eq. 13. The result is given by Eq. 14 [29]. 𝑇 ― 𝑇0 =

𝜋 2

∫ ∫ ―

𝜋 2

𝜔(𝛹) 𝑊𝑝 0

2𝜋𝑘

(

𝑒𝑥𝑝 ―

(14)

) ( )

𝑣𝑠 𝑣𝑠 cos 𝛾 𝐾0 𝑑𝑟𝑑𝛹 2𝜆 2𝜆

where 𝑇0 is the ambient temperature. By considering 𝑠 = 𝑅2 + 𝑟2 ― 2𝑅𝑟𝑐𝑜𝑠(𝛹 ― 𝐴) in Eq. 14, the temperature can be predicted by the following equation. 𝑇(𝑅,𝐴) = 𝑇0 +

∫

1 ― 𝑛′ ∆𝐾2𝐼 𝑓𝛽 1 + 𝑛′2𝜋𝑘𝐸𝐼𝑛′

𝜔(𝛹) = 2𝑟𝑝cos (𝜃)

∫

(15)

𝜋 2 ―

𝜋 2

(

𝜎𝑒𝑞(𝑛′,𝛹)𝜖𝑒𝑞(𝑛′,𝛹)𝑒𝑥𝑝 ―

0

) (

)

𝑣 𝑣 12 (𝑅𝑐𝑜𝑠𝐴 ― 𝑟𝑐𝑜𝑠𝛹) 𝐾0 (𝑅2 + 𝑟2 ― 2𝑅𝑟𝑐𝑜𝑠(𝛹 ― 𝐴)) 2𝜆 2𝜆

𝑑𝑟𝑑𝛹

To determine the temperature rise around the crack, a double integration on the plastic zone area needs to be performed. This is a time-consuming process in addition to having a singularity issue when the value of 𝑠 in Eq. 14 approaches zero. Typically, the integration is divided into two regions with the singular point on the border of the two. This approach has been long practiced, yet it still

5|Page

has numerical difficulties. In order to solve these problems, an attempt is made to simplify the integration to a form that is more convenient to apply. The temperature rise on the crack tip can be determined by considering 𝐴 = 0 and 𝑅 = 0 in Eq. 15. Since the zero order Bessel function of the second kind can be replaced with natural logarithm function ( 𝐾0 (𝑥)≅ ― 𝑙𝑛(x)) for 𝑥 values close to zero and exp (𝑥)≅1 + x for small values of argument 𝑥≪1

𝑥, Eq. 15 is simplified to Eq. 16. 1 ― 𝑛′ ∆𝐾2𝐼 𝑓𝛽 ∆𝑇𝑡𝑖𝑝 = ― 1 + 𝑛′2𝜋𝑘𝐸𝐼𝑛′

𝜋 2

∫ ∫ ―

𝜋 2

𝜔(𝛹)

(

𝜎𝑒𝑞(𝑛′,𝛹)𝜖𝑒𝑞(𝑛′,𝛹) 1 +

0

(16)

)( )

𝑣 𝑣𝑟 𝑟𝑐𝑜𝑠𝛹 𝑙𝑛 𝑑𝑟𝑑𝛹 2𝜆 2𝜆

Setting the limit of integration into the plastic zone and using Eq. 10 yields 1 ― 𝑛′ ∆𝐾2𝐼 𝑓𝛽

𝜋

𝜔(𝛹) = 2𝑟𝑝cos (𝛹)

∆𝑇𝑡𝑖𝑝 = ― 1 + 𝑛′2𝜋𝑘𝐸𝐼 ∫2 𝜋∫0 𝑛′

―2

( ( )+

𝜎𝑒𝑞(𝑛′,𝛹)𝜖𝑒𝑞(𝑛′,𝛹) 𝑙𝑛

𝑣𝑟 2𝜆

𝑣 2𝜆𝑟𝑐𝑜𝑠𝛹

(17)

( ))𝑑𝑟𝑑𝛹

𝑙𝑛

𝑣𝑟 2𝜆

Operating the first integration, yields the following equation for temperature rise as a single integration shown by Eq. 18. 1 ― 𝑛′ ∆𝐾2𝐼 𝑓𝛽 ∆𝑇𝑡𝑖𝑝 = ― 1 + 𝑛′2𝜋𝑘𝐸𝐼𝑛′

∫

(

―

(

(

′ ′ 𝜋𝜎𝑒𝑞(𝑛 ,𝛹)𝜖𝑒𝑞(𝑛 ,𝛹) 2𝑟𝑝cos (𝛹)𝑙𝑛 2

𝑣𝑟𝑝𝑐𝑜𝑠𝛹 𝑣 1(18 + (𝑟𝑝𝑐𝑜𝑠𝛹)2 𝑙𝑛 ― 𝑑𝛹 𝜆exp (1) 𝜆 2𝜆 2 )

𝑣𝑟𝑝cos 𝛹

((

)

) ― ) can be neglected in comparison with 2𝑟 cos (𝛹)𝑙𝑛(

𝑣𝑟𝑝𝑐𝑜𝑠𝛹

The term 𝜆(𝑟𝑝𝑐𝑜𝑠𝛹)2 𝑙𝑛( 𝑣

𝜋 2

2𝜆

1 2

𝑝

) ))

𝑣𝑟𝑝cos 𝛹 𝜆exp (1)

) in the

last equation since it is of second order of 𝑟𝑝. Simplifying Eq.18 yields: 1 ― 𝑛′∆𝐾2𝐼 𝑓𝛽𝑟𝑝 ∆𝑇𝑡𝑖𝑝 = ― 1 + 𝑛′ 𝜋𝑘𝐸𝐼𝑛′

∫

𝜋 2

𝑣𝑟𝑝cos 𝛹

(

(𝑛′,𝛹)𝜖𝑒𝑞(𝑛′,𝛹)cos (𝛹)𝑙𝑛

𝜋𝜎𝑒𝑞 ― 2 𝑣𝑟𝑝cos 𝛹

𝜆exp (1)

)

(19) 𝑑𝛹

𝑣𝑟𝑝

In the last relation 𝑙𝑛( 𝜆exp (1) ) can be replaced with 𝑙𝑛(𝜆exp (1)) ― 𝑙𝑛(cos 𝛹). Further, it can be shown that 𝑙𝑛(cos 𝛹) has a small contribution to the total integration in the range of values of engineering material properties normally used. The validity of these simplifications is investigated in the next section. Therefore, the final relation is shown by Eq. 20. ∆𝑇𝑡𝑖𝑝 = ―

𝑣𝑟𝑝 1 ― 𝑛′𝑓𝛽∆𝐾2𝐼 𝑟𝑝𝐼′𝑛′ 𝑙𝑛 𝜆exp (1) 1 + 𝑛′ 𝜋𝑘𝐸𝐼𝑛′

(

)

(20)

where 𝐼′𝑛′ is given by Eq. 21. 𝐼′𝑛′ =

∫

𝜋 2 ―

′ ′ 𝜋𝜎𝑒𝑞(𝑛 ,𝛹)𝜖𝑒𝑞(𝑛 ,𝛹)cos 𝛹𝑑𝛹

(21)

2

It is seen that 𝐼′𝑛′ is only a function of 𝑛′ which is a material property that can be found knowing the two angular distribution functions 𝜎𝑒𝑞 and 𝜖𝑒𝑞. The approach is proposed by Hutchinson et al. [30] . Values for 𝐼′𝑛′ are given in Table 1. Table 1. Value of 𝐼′𝑛′ given by Eq.21 for different cyclic work hardening values.

6|Page

0.1 1.18

𝑛′ 𝐼′𝑛′

2.2.

0.15 1.17

0.2 1.14

0.25 1.12

0.3 1.10

0.35 1.06

0.4 1.03

0.5 0.98

Numerical procedure

To evaluate the accuracy of Eq. 20, its predictions are compared with the numerical solutions of the original complete integration of Eq. 15. The double integration of Eq. 15 is performed using the two-dimensional Simpson’s integration rule. The values for 𝜎𝑒𝑞(𝑛′,𝛹) and 𝜖𝑒𝑞(𝑛′,𝛹) are determined using the approach suggested by Nikishkov [34] (see Appendix A).The propagation speed 𝑣 in Eq. 20 is a function of stress the intensity factor range (∆𝐾). Therefore, an iterative approach is needed to solve Eq. 20. The Paris law is used to determine 𝑣. Convergence is achieved normally after several trials. The temperature rise given by Eq. 20 and Eq. 15 are compared for different values of stress intensity factor in Fig. 2. The material properties used for the simulation is found in Table. 2 for SS 304 when frequency is 20 Hz. 80 Eq. 15 Eq. 20

70

∆𝑇 (𝐾)

60 50 40 30 20 10 0

0

20

40

60

80

100

∆𝐾 (𝑀𝑃𝑎 𝑚

Figure 2. Comparison of temperature rise on the moving crack tip given by Eqs. 15 and 20. for different stress intensity factors

It can be observed that the error is small in low-stress intensities and increases as it reaches the high values. Note that according to [35], ∆𝐾𝐼 = 100 𝑀𝑃𝑎 𝑚 is above the region where crack propagation can be found by applying Paris law for SS304. Thus, the maximum error happens using Eq. 20 instead of double integration of Eq. 15; within the region where the Paris law is applicable, the error is less than 2%.

7|Page

3. Experimental 3.1.

Material and specimen

To evaluate the efficacy of the model, the temperature rise is calculated around the crack tip for the specimens made of SS304 and SS316. The property data and chemical compositions of these materials are given in Tables 1 and 2, respectively. Table 2. Material properties of specimen Material

𝑬 (𝑮𝑷𝒂)

SS 304 SS 316

193 193

𝒌 𝑾 ) ( 𝒎𝑲

𝒏′

16.3 14.4

𝒌′ (𝑴𝑷𝒂)

0.26 0.27

1200 1150

(

𝝀 𝒎𝟐 𝒔

𝒕 (𝒎𝒎)

)

3.56 × 10-6 3.58 × 10-6

1.85 1.58

𝝈𝟎 (𝑴𝑷𝒂)

270 265

Table 3. Composition of material in Table 1. Composition SS 304 SS 316

3.2.

Fe

C

53.4874.5% 58.2373.61%

00.08% 00.08%

Cr

Cu

Mn

Mo

17.524% 1618.5%

0-1%

0-2%

0-2.5%

0-1%

0-2%

0-3%

Ni 02.5% 0-3%

Ni 00.1% 00.1%

P 0-0.2% 00.045%

S 0.35 % 0.35 %

Specimen preparation

Flat specimens are used with a pre-crack notch of length 5 mm as depicted in Fig. 3. All the specimens are cut from steel strips and coated with a thin layer of black paint sprayed on the surface to increase thermal emissivity to the known value of 0.98 for the particular color. The tests are done in ambient temperature of 298 𝐾. 5mm 5mm

𝑎

50mm

150 mm Figure 3. Flat specimen schematic used for the crack propagation test.

3.3 Experimental procedure

8|Page

Experimental data are collected using a servo-hydraulic test machine with the maximum capacity of 50 kN in tension-compression loading. The machine can apply cyclic load with frequency as large as 75 Hz. The surface temperature is detected by means of an infrared (IR) camera (FLAIR A615) capable of capturing 640 × 480 pixel microbolometer (Fig. 4). The sensitivity is 50 mK and 100 frames per second can be recorded in the 640 × 480 pixel mode and as high as 200 frames per second can be captured for 640 × 120 resolution. Crack tip temperature is measured by mounting a Close-Up IR Lens, (16mm × 12mm) 1.5X (25μm). A digital 40-1000X microscope is used to see the crack tip. The IR camera records the temperature from the front of the specimen and a microscope is placed behind the specimen where a set of LED lights make the surface clearer for the detection of the crack tip. As the crack moves during the test, the IR camera and microscope are moved using two separate sliding stages with electrical motors that makes it possible for both cameras to travel with variable speeds to track the movement of the crack. In order to accurately relate the crack length and propagation rate to stress intensity factor at any time during the test, each IR image has to be attributed to a specific microscope picture. This means that pictures have to be recorded and saved at the same rate. To perform this, a LabVIEW code is prepared that captures and saves images simultaneously from the two sources.

Top grip

Bottom grip

Moving microscope

Specimen Microscope stand Moving stage

Grip control IR camera Axial 50 KN Torque 2kN

Electric motor Moving Stage Adjustment screw

9|Page

Figure. 4. Hydraulic fatigue test machine and infrared camera with a close up lens mounted

4. Results and discussion The value of the stress intensity factor is found by a 2D finite element simulation using the ABAQUS software. The results of FE analysis are validated by a process of mesh refinement to ensure that results are independent of the mesh size. The J integral method is applied in this regards with second-order elements. For each crack length, the stress intensity factor is calculated in 10 contours for J integral computation. MPC Beam constraint is used on the top surface of the specimen to account for rotation restriction there. To compensate for the motion of specimen, for each thermal image, the maximum temperature location is flagged and the images are made stable by removing the vertical motion of maximum temperature location. Furthermore, each image is processed to reduce any unexpected temperature jumps coming from noises. For this purpose, a smoothening algorithm is used which is based on local regression using linear least squares and a second-degree polynomial. The captured image of the crack tip in cyclic loading using a close-up lens shows the temperature rise around the tip in Fig. 5. Also shown the temperature evolution in three stages. At the beginning of the crack propagation, a temperature rise is seen until the crack’s heat dissipation reaches a thermal equilibrium in the second region where the temperature is nearly steady. The third stage begins when the specimen reaches the time that around 90 percent of the total life (𝑁𝑓) is passed.

10 | P a g e

400 375

𝑲

365 355

380

345

𝑇 (𝐾)

335

1 mm

1 mm

325

360

310

340

𝑲

316

309

314

308

312

307 306

1 mm

𝑲

310 308

1 mm

305

306

320 300

0

0.6

0.2 0.3 0.4

0.6 5

𝑁 𝑁𝑓

0.8

0.98

1

Figure 5. Maximum temperature evolution for fatigue test of SS304 with f = 35 Hz and fully reversed 10 (𝐾𝑁) load. 𝑁

𝑁

𝑁

Temperature distribution around the crack is shown at 𝑁𝑓 = 0.3 , 𝑁𝑓 = 0.65 and 𝑁𝑓 = 0.98

The numerical solution of Eq. 15 is carried out using the MATLAB software. It is consisting of a shooting method algorithm to determine the stress-strain solution of HRR singularity using the Rung-Kutta approach. 𝑁

The temperature distribution solution contours around the crack tip is shown in Fig. 6 for 𝑁𝑓 = 0.3 shown in Fig. 5.

11 | P a g e

𝑌/𝑟𝑝

𝑇 = 308.0 K

Crack faces

308.5 309.0 309.5

𝑋/𝑟𝑝

Crack tip

𝑁

Figure 6. Solution of temperature distribution (Eq.15) around the crack tip for 𝑁𝑓 = 0.3 (∆𝐾𝐼 = 38 𝑀𝑃𝑎 𝑚) shown in Fig. 5

Fig. 7 shows the predicted crack propagation speed as a function of a range of stress intensity factor for SS304 and SS316. It is seen that the region where Paris law is applicable for these two materials falls within 40 to 70 (𝑀𝑃𝑎 𝑚) range. There is a resemblance between different stages of FCG in Fig. 7 and temperature evolution stages of Fig. 5. It is seen that the temperature rise at the end of the test corresponds to the 3rd stage of propagation in Figure 7. Physically, the final temperature rise is due to the fast accumulation of dissipated energy which makes the temperature unstable in this region. It is inferred from this figure that in a short period of time, the crack speed grows by 10-fold and this growth continues until the final fracture. After the fracture, the temperature drops rapidly (not shown in Fig. 5). The first stage of rapid growth after passing the threshold with a small amount of microplastic deformation is the reason for temperature rise. In the second stage, the inelastic effect and heat conduction reach an equilibrium and most of the life is spent there [36].

12 | P a g e

10-5

SS304 SS316

𝑉(𝑚/𝑐𝑦𝑐𝑙𝑒)

10-6

1.97

𝑉 = 1.22 × 10 ―10(∆𝐾)

10-7 𝑉 = 1.20 × 10 ―10(∆𝐾)1.98

30

40

50

60

70

80

90

∆𝐾(𝑀𝑃𝑎 𝑚)

Figure. 7. Crack propagation speed in different stress intensity factors for SS304 and SS316

The Paris law constants of the two steel specimens used in the study are found according to Fig. 7. The values of these constants are close to each other in the ambient temperature. In Fig. 8 the temperature distribution ahead of the crack tip as predicted by Eq. 15 is compared with the test results. The experimental crack propagation speed is used to solve the temperature rise relation. It is seen that maximum temperature occurs in a small distance ahead of the tip and its value is very close to the crack tip temperature. This indicates that without knowing the location of the crack tip on the thermal image, the maximum temperature value can be used instead in Eq. 20. These plots show good agreement between the test results and the numerical solution of the model. It can be observed that the temperature drop from the crack tip is higher in elevated ∆𝐾𝐼 values. For example, for SS304 temperature drop at 𝑟=8 mm when is around 1 degree when ∆𝐾𝐼 =40 while this number is 25 degrees for ∆𝐾𝐼=100. This shows the concentration of heat generation is concentrated in a small area ahead of the crack tip, which is the cyclic plastic zone.

13 | P a g e

400

Experiment Eq. 15 solution

380

𝑇 (K)

∆𝐾𝐼 = 100 (𝑀𝑃𝑎 𝑚)

90

360

80 340

70 60

320

50 40

300

0

1

2

3

4

5

6

7

8 10-3

𝑟 (m)

(a)

380 Experiment Eq. 15 solution

370 360 𝑇 (K)

∆𝐾𝐼 = 100 (𝑀𝑃𝑎 𝑚)

350 340

90

330 80

320 310 300

50

60

70

40

0

1

2

3

4 𝑟 (m)

(b)

5

6

7

8 -3

10

Figure. 8. Temperature evolution around the crack tip for different stress intensity factors and comparison with the experiment results for (a) SS304 and (b) SS316

14 | P a g e

Fig. 9 compares the experimental measurement of stress intensity factor (SIF), ∆𝐾, obtained thermographically with the value given by Eq. 20 for the two materials used in this study. In this figure, ∆𝐾 is plotted against the temperature rise at the crack tip. This is seen that temperature rise is a powerful parameter for effective prediction of the stress intensity. The temperature rise is small in the zone where the stress intensity factor is small. In Fig. 9 the simulation results are shown for two cases. For the first case, Eq. 20 is solved by using the crack speed (𝑣) predicted from the Paris law and in the second case the experimentally measured speed is applied in Eq. 20. It is seen that in when range is low, the Paris law gives reasonable results but in high speeds and that the thermography predicted SIF deviates from the experimental results. By using the experimentally measured speed in Eq. 20, the values predicted for stress intensity factor becomes close to the calculated SIF.

∆𝐾 (𝑀𝑃𝑎

160

160

140

140

120

∆𝐾

𝑚) 100

120

(𝑀𝑃𝑎 𝑚)

100

80

Experiment

60 40 20 0

Experiment

80

First case (Eq. 20 solution using Paris Law for crack speed

60

First case (Eq. 20 using Paris Law for crack speed)

Second case (Eq. 20 solution using experiment speed)

40

Second case (Eq. 20 solution using Experiment speed)

50

100

𝑇𝑒𝑚𝑝𝑒𝑟𝑎𝑡𝑢𝑟𝑒 𝑅𝑖𝑠𝑒 (K)

150

0

50

100

150

𝑇𝑒𝑚𝑝𝑒𝑟𝑎𝑡𝑢𝑟𝑒 𝑅𝑖𝑠𝑒 (K)

Figure. 9. Comparison between Stress intensity factor given by Eq. 20 and experimental results for (a) SS304 and (b) SS316

15 | P a g e

5. Conclusions A thermal imaging method applied in conjunction with the analysis of Pandy et al. [29] is proposed to determine the stress intensity factor of a moving crack in the first mode cyclic fatigue crack loading. The relation provides the stress intensity factor based on the temperature rise on tip of the crack. Experimental tests results show good agreement between the model and the empirical results. It is also shown that a simplified analytical model presented is capable of determining the stress intensity factor by measuring the temperature rise at the tip of the crack. In the proposed equation relating temperature rise and SIF, the crack speed could be plugged in from experimental results; however, implementing the Paris law for propagation speed yields reasonable results even for the range of SIF that are out of region where Paris law is applicable.

Appendix A Procedure for determination of 𝜎𝑒𝑞(𝑛′,𝛹) and 𝜖𝑒𝑞(𝑛′,𝛹) Ramberg-Osgood Equation for uniaxial case. 𝑛

(A.1)

()

𝜖 𝜎 𝜎 = +𝛼 𝜖 0 𝜎0 𝜎0

In which; 𝑛 = 1/𝑛′ 𝑛′: Cyclic plastic strain hardening exponent

𝜖: uniaxial strain 𝜎: uniaxial stress 𝛼: hardening coefficient 𝜎0: Yield stress 𝑛: Hardening exponent 𝜖0 = 𝜎0/𝐸

Ramberg-Osgood in Multi-axial case. 𝑛―1

()

𝜎𝑒 1 ― 2𝜈 3 1+𝜈 𝑆𝑖𝑗 + 𝜎𝑘𝑘𝛿𝑖𝑗 + 𝛼𝜖0 𝜖𝑖𝑗 = 3𝐸 2 𝜎0 𝐸

where 𝜎𝑒 is the equivalent stress. 16 | P a g e

𝑆𝑖𝑗 𝜎0

(A.2)

Equilibrium Equations in 3 directions; 𝒓 Direction: 𝜎𝑟 ― 𝜎𝜃 1 𝜎𝑟,𝑟 + 𝜎𝑟𝜃,𝜃 + + 𝜎𝑟𝑧,𝑧 = 0 𝑟 𝑟

(A.3)

𝜽 Direction: 1 2 1 𝜎 + 𝜎𝑟𝜃,𝑟 + 𝜎𝑟𝜃 + 𝜎𝑧𝜃 + 𝜎𝑧𝜃,𝑧 = 0 𝑟 𝜃,𝜃 𝑟 𝑟

(A.4)

𝒛 Direction 1

𝜎𝑧,𝑧 + 𝑟𝜎𝑧𝜃,𝜃 + 𝜎𝑧𝑟,𝑟 +

𝜎𝑟𝑧 𝑟

(A.5)

=0

Strain definition in respect with displacement field: 𝜖𝑟 = 𝑢𝑟,𝑟 1 1 𝜖𝜃 = 𝑢𝑟 + 𝑢𝜃,𝜃 𝑟 𝑟 𝜖𝑧 = 𝑢𝑧,𝑧 1 𝜖𝑟𝑧 = (𝑢𝑟,𝑧 + 𝑢𝑧,𝑟) 2 1 1 𝜖𝜃𝑧 = 𝑢 + 𝑢𝜃,𝑧 2 𝑟 𝑧,𝜃 1 1 1 1 𝜖𝑟𝜃 = 𝑢𝑟,𝜃 + 𝑢𝜃,𝑟 ― 𝑢𝜃 2 𝑟 𝑟 𝑟

( (

(A.6) (A.7) (A.8) (A.9)

)

(A.10)

)

(A.11)

𝑯𝑹𝑹 singularity denotes: 𝜎𝑖𝑗(𝑟,𝜃) 𝜖𝑖𝑗

𝜎0

(A.12)

= 𝑟𝑠𝜎𝑖𝑗(𝑟,𝜃) 𝑢𝑖

= 𝑟𝑠𝑛𝜖𝑖𝑗(𝑟,𝜃) which means 𝛼𝜖0 = 𝑟𝑠𝑛 + 1𝑢𝑖(𝑟,𝜃) where 𝑠 is a constant to be found as a part of solution. 𝛼𝜖0

(A.13)

Plastic strain from multiaxial Ramberg-Osgood Eq. A.2 reduces to: 𝜖𝑖𝑗 𝛼𝜖0

=

𝑛―1

()

3 𝜎𝑒 2 𝜎0

𝑆𝑖𝑗

(A.14)

𝜎0

The deviatoric stress is defined by Eq. A.15. 1 𝑆𝑖𝑗 = 𝜎𝑖𝑗 ― 𝜎𝑘𝑘𝛿𝑖𝑗 3

(A.15)

Plain stress condition (𝜎𝑧𝑧 = 𝜎𝑧𝜃 = 𝜎𝑧𝑟 = 0) Eq. A.15 reads: 1 1 1 𝑆𝑟𝑟 = 𝜎𝑟𝑟 ― (𝜎𝑟𝑟 + 𝜎𝜃𝜃 + 𝜎𝑧𝑧) = (2𝜎𝑟𝑟 ― 𝜎𝜃𝜃) = (2𝜎𝑟 ― 𝜎𝜃) 3 3 3 1 1 1 𝑆𝜃𝜃 = 𝜎𝑧𝑧 ― (𝜎𝑟𝑟 + 𝜎𝜃𝜃 + 𝜎𝑧𝑧) = (2𝜎𝜃𝜃 ― 𝜎𝑟𝑟) = (2𝜎𝜃 ― 𝜎𝑟) 3 3 3 1 1 1 𝑆𝑧𝑧 = 𝜎𝑧𝑧 ― (𝜎𝑟𝑟 + 𝜎𝜃𝜃 + 𝜎𝑧𝑧) = ― (𝜎𝜃𝜃 + 𝜎𝑟𝑟) = ― (𝜎𝑟 + 𝜎𝜃) 3 3 3

17 | P a g e

(A.16) (A.17) (A.18)

Equivalent stress The value of equivalent stress distribution function, 𝜎𝑒, for plain stress is: 1

(

)

2 3 𝜎𝑒 = 𝑆𝑖𝑗𝑆𝑖𝑗 2 𝑆𝑖𝑗𝑆𝑖𝑗 = 𝑆11𝑆11 + 𝑆12𝑆12 + 𝑆13𝑆13 + 𝑆21𝑆21 + 𝑆22𝑆22 + 𝑆23𝑆23 + 𝑆31𝑆31 + 𝑆32𝑆32 + 𝑆33𝑆33

(A.19)

𝑆𝑖𝑗𝑆𝑖𝑗 = 𝑆211 + 𝑆222 + 𝑆233 + 2𝑆12𝑆12 + 2𝑆13𝑆13 + 2𝑆23𝑆23

(A.20)

𝑆𝑖𝑗𝑆𝑖𝑗 = 𝑆2𝑟 + 𝑆2𝜃 + 𝑆2𝑧 + 2𝑆2𝑟𝜃 + 2𝑆2𝑟𝑧 + 2𝑆2𝜃𝑧

Therefore, for plain stress condition Equivalent stress is:

(

3 𝜎𝑒 = (𝑆2𝑟 + 𝑆2𝜃 + 2𝑆2𝑟𝜃) 2

[ {(

2

1 2

(A.21)

)

) (

3 1 (2𝜎𝑟 ― 𝜎𝜃) 𝜎𝑒 = 2 3

2

) }]

) (

1 + (2𝜎𝜃 ― 𝜎𝑟) 3

1 + ― (𝜎𝑟 + 𝜎𝜃) 3

2

1 2

(A.22)

Eq. A.22 in plain stress condition around the crack tip is given by Eq. A.23 𝜎𝑒 = [

𝜎2𝑟

+

𝜎2𝜃

1 2 2 3𝜎𝑟𝜃

(A.23)

]

― 𝜎𝑟𝜎𝜃 +

To find the shear strain distribution function, 𝜖𝑟𝜃, Eq. A.13 , Eq. A.14 and Eq. A.12 for 𝜎𝑒 are used as, 𝜖𝑟𝜃

𝑛―1

()

3 𝜎𝑒 = 𝑟 𝜖𝑟𝜃 = 𝛼𝜖0 2 𝜎0 𝑠𝑛

3 𝜎𝑒𝜎0𝑟𝑠 = 𝜎0 2 𝜎0

𝑆𝑟𝜃

( )

𝑛―1

𝜎𝑟𝜃𝜎0𝑟𝑠

(A.24)

𝜎0

Hence, 3 𝜖𝑟𝜃 = (𝜎𝑒)𝑛 ― 1𝜎𝑟𝜃 2

(A.25)

The same process is needed to find 𝜖𝑟 : 𝜖𝑟

𝑛―1

()

3 𝜎𝑒 = 𝑟 𝜖𝑟 = 2 𝜎0 𝛼𝜀0 𝑠𝑛

3 𝜎𝑒𝜎0𝑟𝑠 𝑠𝑛 𝑟 𝜖𝑟 = 2 𝜎0

(

(A.26)

𝜎0 1 (2𝜎𝑟 ― 𝜎𝜃)𝑟𝑠𝜎0 3 𝜎0

𝑛―1

( )

𝑆𝑟

1 𝜖𝑟 = (𝜎𝑒)𝑛 ― 1 𝜎𝑟 ― 𝜎𝜃 2

)

(A.27)

(A.28)

This is: Eq.A.13 is used. 𝑢𝑟 𝛼𝜀0 𝑢𝜃 𝛼𝜀0

= 𝑟𝑠𝑛 + 1𝑢𝑟 = 𝑟𝑠𝑛 + 1𝑢𝜃

18 | P a g e

(A.29) (A.30)

The derivative of 𝑢𝑟 and 𝑢𝜃 are found according to the last two relations. The result is shown in Eq. A.31 and Eq. A.32. 𝑢𝜃,𝑟 𝛼𝜀0 𝑢𝑟,𝜃 𝛼𝜀0

(A.31)

= (𝑠𝑛 + 1)𝑟𝑠𝑛𝑢𝜃

(A.32)

= 𝑟𝑠𝑛 + 1𝑢𝑟,𝜃

Using Eq. A.11 and substituting from Eqs. A.25, A.30, A.31 and A.32 in Eq. A.11 yields Eq. A.33. 𝜖𝑟𝜃𝑟𝑠𝑛 =

(

1 1 𝑠𝑛 + 1 1 𝑟 𝑢𝑟,𝜃 + (𝑠𝑛 + 1)𝑟𝑠𝑛𝑢𝜃 ― 𝑟𝑠𝑛 + 1𝑢𝜃 2 𝑟 𝑟

)

(A.33)

Simplifying yields: 1 𝜖𝑟𝜃𝑟𝑠𝑛 = 𝑟𝑠𝑛(𝑢𝑟,𝜃 + (𝑠𝑛 + 1)𝑢𝜃 ― 𝑢𝜃) 2

(A.34)

Thus, 𝑢𝑟,𝜃 + 𝑠𝑛𝑢𝜃 ― 2𝜖𝑟𝜃 = 0

(A.35)

The values for 𝜖𝑟 and 𝑢𝑟 are given by Eq. A.36 and Eq. A.37 as: 𝑟𝑠𝜖𝑟 =

𝜖𝑟

(A.36)

𝛼𝜖0

𝑟𝑠𝑛 + 1𝑢𝑟 =

𝑢𝑟

(A.37)

𝛼𝜖0

Substitution of 𝜖𝑟 and 𝑢𝑟 from Eqs. A.36 and A.37 into Eq.A.6 yields: (𝑠𝑛 + 1)𝑟𝑠𝑛𝑢𝑟,𝑟 ― 𝑟𝑠𝑛𝜖𝑟 = 0

(A.38)

which simplified to Eq. A.39. (𝑠𝑛 + 1)𝑢𝑟 ― 𝜖𝑟 = 0

(A.39)

Rewriting Eq. A.4 for plain strain condition (𝜎𝑧𝜃 = 0 and 𝜎𝑧𝜃,𝑧 = 0) regarding Eq. A.12 yields: 1 𝑠 2 𝑟 𝜎 + 𝑠𝑟𝑠 ― 1𝜎𝑟𝜃 + 𝑟𝑠𝜎𝑟𝜃 = 0 𝑟 𝜃,𝜃 𝑟

(A.40)

Therefore, 𝜎𝜃,𝜃 + 𝑠𝜎𝑟𝜃 + 2𝜎𝑟𝜃 = 0

(A.41)

or, 𝜎𝜃,𝜃 + (𝑠 + 2)𝜎𝑟𝜃 = 0

(A.42)

The same procedure is followed for Eq.A.3. 𝜎𝑟 ― 𝜎𝜃 1 𝜎𝑟,𝑟 + 𝜎𝑟𝜃,𝜃 + =0 𝑟 𝑟

(A.43)

Putting back values using Eq.A.12 𝑠𝑟𝑠 ― 1𝜎𝑟,𝑟 + 𝑟𝑠 ― 1𝜎𝑟𝜃,𝜃 + 𝑟𝑠 ― 1(𝜎𝑟 ― 𝜎𝜃) = 0

(A.44)

Thus, (𝑠 + 1)𝜎𝑟 + 𝜎𝑟𝜃,𝜃 ― 𝜎𝜃 = 0 𝜖𝜃

(A.45)

𝑠𝑛

𝜖𝜃 from 𝛼𝜀0 = 𝑟 𝜖𝜃 can be used in Eq. A.7. 𝑢𝑟 + 𝑢𝜃,𝜃 ― 𝜖𝜃 = 0

19 | P a g e

(A.46)

Finally, we arrive at the following set of differential equations using Eqs. A.35, A.39, A.42, A.45 and A.46.

{

𝑢𝑟,𝜃 + 𝑠𝑛𝑢𝜃 ― 2𝜖𝑟𝜃 = 0 (𝑠𝑛 + 1)𝑢𝑟 ― 𝜖𝑟 = 0 𝜎𝜃,𝜃 + (𝑠 + 2)𝜎𝑟𝜃 = 0 (𝑠 + 1)𝜎𝑟 + 𝜎𝑟𝜃,𝜃 ― 𝜎𝜃 = 0 𝑢𝑟 + 𝑢𝜃,𝜃 ― 𝜖𝜃 = 0

(A.47)

The second relation in Eq. A.47 is differentiated in respect with 𝜃: (A.48)

(𝑠𝑛 + 1)𝑢𝑟,𝜃 ― 𝜖𝑟,𝜃 = 0

Combining last equation with the first relation of the system of equations and applying Eq. A.50 yields Eq. A.51.

(

)

(

)

(

1 1 𝜖𝑟,𝜃 = (𝑛 ― 1)𝜎𝑒,𝜃𝜎𝑛𝑒 ― 2 𝜎𝑟 ― 𝜎𝜃 + 𝜎𝑛𝑒 ― 2 𝜎𝑟,𝜃 ― 𝜎𝜃,𝜃 2 2

(

(

1 1 1 (𝑛 ― 1)𝜎𝑒,𝜃𝜎𝑛𝑒 ― 2 𝜎𝑟 ― 𝜎𝜃 + 𝜎𝑛𝑒 ― 2 𝜎𝑟,𝜃 ― 𝜎𝜃,𝜃 2 2 (𝑠𝑛 + 1)

))

)

(A.49)

+ 𝑠𝑛𝑢𝜃 ― 3𝜎𝑛𝑒 ― 1𝜎𝑟𝜃 = 0

(A.50)

Solving Eq. A.51 for 𝜎𝜃,𝜃 results in Eq. A.51 ― 𝜎𝑟,𝜃 =

{

{

}

𝜎𝜃,𝜃 1 𝑛―1 1 𝑠𝑛 (2𝜎𝜃,𝜃𝜎𝜃 + 6𝜎𝑟𝜃,𝜃𝜎𝑟𝜃 ― 𝜎𝑟𝜎𝜃,𝜃)𝜎𝑛𝑒 ― 3 𝜎𝑟 ― 𝜎𝜃 ― ― 𝑛 ― 1(𝑢𝜃) + 3𝜎𝑟𝜃 𝑠𝑛 + 1 2𝜎2 2 2 𝜎𝑒 𝑒

{(

1 (𝑛 ― 1) (2𝜎𝑟 ― 𝜎𝜃 𝑠𝑛 + 1 2𝜎2𝑒

― 𝜎𝑟,𝜃 =

( ))(

) )

(A.51)

}

1 𝜎𝑟 ― 𝜎𝜃 + 1 2

{

}

𝜎𝜃,𝜃 1 𝑛―1 1 𝑠𝑛 (2𝜎𝜃,𝜃𝜎𝜃 + 6𝜎𝑟𝜃,𝜃𝜎𝑟𝜃 ― 𝜎𝑟𝜎𝜃,𝜃)𝜎𝑛𝑒 ― 3 𝜎𝑟 ― 𝜎𝜃 ― ― 𝑛 ― 1(𝑢𝜃) + 3𝜎𝑟𝜃 𝑠𝑛 + 1 2𝜎2 2 2 𝜎 𝑒 𝑒

{(

1 (𝑛 ― 1) (2𝜎𝑟 ― 𝜎𝜃 𝑠𝑛 + 1 2𝜎2𝑒

𝜎𝜃,𝜃 = ― (𝑠 + 2)𝜎𝑟𝜃 𝜎𝑟𝜃,𝜃 = 𝜎𝜃 ― (𝑠 + 1)𝜎𝑟

{(

𝑢𝜃,𝜃 = 𝜎𝑛𝑒 ― 1 𝜎𝜃 1 ―

) (

( ))(

) )

}

1 𝜎𝑟 ― 𝜎𝜃 + 1 2

(A.52)

)}

1 1 1 + 𝜎𝑟 ― + 2(𝑠𝑛 + 1) 2 𝑠(𝑠𝑛 + 1)

where 𝜎𝑒 = (𝜎2𝑟 + 𝜎2𝜃 ― 𝜎𝑟𝜎𝜃 + 3𝜎2𝑟𝜃)

1/2

𝑠 = ―1/(𝑛 + 1)

The system of differential equations given by Eq. A.52 is solved for the following unknowns. 𝜎𝑟,𝜎𝜃,𝜎𝑟𝜃,𝑢𝜃 (A.53) This is a system of nonlinear differential equations which is solved considering the following boundary conditions for displacement and stress around the crack.

20 | P a g e

𝑢𝜃(0) = 0 𝜎𝑟𝜃(0) = 0 𝜎𝜃(𝜋) = 0 𝜎𝑟𝜃(𝜋) = 0

(A.56)

This means: 𝑢𝜃(0) = 0 𝜎𝑟𝜃(0) = 0 𝜎𝜃(𝜋) = 0 𝜎𝑟𝜃(𝜋) = 0

(A.57)

The shooting method is applied for solution of Eq. A.69. To solve the differential equations using RungeKutta approach. To find initial values of 𝜎𝑟 and 𝜎𝑟𝜃 at 𝜃 = 0 there needs to be an initial guess for 𝜎𝑟|𝜃 = 0 and 𝜎𝑟𝜃|𝜃 = 0 to start with. The following values are helpful in this regard for different 𝑛 values given in (A.73) 𝑛 = [2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 15, 20]

(A.58)

𝜎𝑟|𝜃 = 0 = [0.842, 0.769, 0.721, 0.690, 0.669, 0.654, 0.643, 0.634, 0.628, 0.618, 0.608, 0.6]

(A.59)

𝜎𝑟𝜃|𝜃 = 0 = [1.053, 1.106, 1.125, 1.134, 1.139, 1.142, 1.144, 1.146, 1.147, 1.148, 1.1, 1.151]

(A.60)

HRR integration constant 𝐼𝑛′ and 𝐼′𝑛′ is defined by the following equations. 𝐼𝑛′ =

𝐼′𝑛′

∫

=

𝜋

{

―𝜋

∫

𝑛′ + 1

[

((

) )

]}

𝑛′ + 2 ′ 𝜎 𝑛 𝑐𝑜𝑠(𝜃) ― 𝑠𝑖𝑛(𝜃)(𝜎𝑟(𝑢𝜃 ― 𝑢𝑟,𝜃) ― 𝜎𝑟𝜃(𝑢𝑟 ― 𝑢𝜃,𝜃)) + 𝑛′ ― 2 + 1 (𝜎𝑟𝑢𝑟 + 𝜎𝑟𝜃𝑢𝜃)𝑐𝑜𝑠(𝜃) 𝑑𝜃 𝑛′ + 1 𝑛′ + 1 𝑒 1

𝜋 2 ―

(A.61) (A.62)

′ ′ 𝜋𝜎𝑒𝑞(𝑛 ,𝛹)𝜖𝑒𝑞(𝑛 ,𝛹)cos (𝛹)𝑑𝛹 2

The numerical solution of last equation is shown in the following plot for different values of 𝑛 which is 1

inverse of strain hardening coefficients (𝑛 = 𝑛′). The present simulation result is compared with Anderson, [37].

21 | P a g e

4.4 Present study Anderson [19]

4.2 4 3.8 3.6

𝐼𝑛′

3.4 3.2 3 2.8 2.6

2

4

6

8

10

12

14

16

18

20

𝑛 = 1/𝑛′ Fig.A.1. HRR integration constant 𝐼𝑛′ for different 𝑛 values and its comparison with Anderson [37].

1.15

1.1

1.05

𝐼′𝑛′ 1

0.95

0.9

2

4

6

8

10

12

14

16

𝑛 = 1/𝑛′ Fig.A.2. 𝐼′𝑛′ constant for different 𝑛 values.

22 | P a g e

18

20

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.

12. 13. 14. 15. 16. 17. 18. 19. 20.

Paris, P. and F. Erdogan, A critical analysis of crack propagation laws. Journal of basic engineering, 1963. 85(4): p. 528-533. Walker, K., Rethinking fatigue crack growth: Crack growth at constant R and strain energy. Fatigue and Fracture of Engineering Materials, Structures, 2018. 41(3): p. 700-707. Zehnder, A.T., Kallivayalil, Jacob A. Temperature rise at the tip of dynamically propagating cracks. in 1991 SEM Spring Conference on Experimental Mechanics. 1991. Kallivayalil, J.A. and A.T. Zehnder, Measurement of the temperature field induced by dynamic crack growth in Beta-C titanium. International journal of fracture, 1994. 66(2): p. 99-120. Wong, A. and G. Kirby III, A hybrid numerical/experimental technique for determining the heat dissipated during low cycle fatigue. J Engineering fracture mechanics, 1990. 37(3): p. 493-504. Meneghetti, G., M. Ricotta, and B.J.P.S.I. Atzori, The heat energy dissipated in a control volume to correlate the fatigue strength of bluntly and severely notched stainless steel specimens. 2016. 2: p. 2076-2083. Meneghetti, G. and M. Ricotta, The use of the specific heat loss to analyse the low-and high-cycle fatigue behaviour of plain and notched specimens made of a stainless steel. J Engineering Fracture Mechanics, 2012. 81: p. 2-16. Meneghetti, G. and M.J.I.J.o.F. Ricotta, Evaluating the heat energy dissipated in a small volume surrounding the tip of a fatigue crack. 2016. 92: p. 605-615. Meneghetti, G. and M. Ricotta, Evaluating the heat energy dissipated in a small volume surrounding the tip of a fatigue crack. International Journal of Fatigue, 2016. 92: p. 605-615. Vshivkov, A., Iziumova, A, Plekhov, O, Bär, J Experimental study of heat dissipation at the crack tip during fatigue crack propagation. J Frattura ed Integrità Strutturale, 2016. 10(35): p. 5763. Vshivkov, A., Iziumova, A Yum Panteleev, IA, Prokhorov, AE, Ilinykh, AV, Wildemann, VE and O. Plekhov. The study of the dissipation heat flow and the acoustic emission during the fatigue crack propagation in the metal. in IOP Conference Series: Materials Science and Engineering. 2017. IOP Publishing. Zhang, H., Wei, CY., Yan, ZF., Wang, WX., Guo, SF., Zhou, YG Research on fatigue crack propagation behaviour of 4003 ferritic stainless steel based on infrared thermography. J Fatigue Fracture of Engineering Materials Structures, 2016. 39(2): p. 206-216. Rice, J., Mechanics of crack tip deformation and extension by fatigue, in Fatigue crack propagation. 1967, ASTM International. Klingbeil, N.W., A total dissipated energy theory of fatigue crack growth in ductile solids. J International Journal of Fatigue, 2003. 25(2): p. 117-128. Smith, K., Application of the dissipated energy criterion to predict fatigue crack growth of Type 304 stainless steel following a tensile overload. J International Journal of Fatigue, 2011. 78(18): p. 3183-3195. Ricotta, M., Meneghetti, Giovanni., Atzori, Bruno., Risitano, Giacomo., Risitano, Antonino, Comparison of Experimental Thermal Methods for the Fatigue Limit Evaluation of a Stainless Steel. 2019. 9(6): p. 677. Corigliano, P., m Epasto, G., Guglielmino, E., Risitano, G Fatigue analysis of marine welded joints by means of DIC and IR images during static and fatigue tests. J Engineering Fracture Mechanics, 2017. 183: p. 26-38. Khonsari, M.M. and M. Amiri, Introduction to thermodynamics of mechanical fatigue. 2012: CRC press. Haghshenas, A. and M. Khonsari, Non-destructive testing and fatigue life prediction at different environmental temperatures. Infrared Physics Technology, 2019. 96: p. 291-297. Mehdizadeh, M. and M.J.I.P. Khonsari, On the application of fracture fatigue entropy to variable frequency and loading amplitude. J Theoretical Applied Fracture Mechanics, 2018. 98: p. 30-37.

23 | P a g e

21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37.

Jang, J., Khonsari, MM J Theoretical, On the evaluation of fracture fatigue entropy. Applied Fracture Mechanics, 2018. 96: p. 351-361. Corigliano, P., Cucinotta, F., Guglielmino, E., Risitano, G., Santonocito, D, Thermographic analysis during tensile tests and fatigue assessment of S355 steel. J Procedia Structural Integrity, 2019. 18: p. 280-286. Mehdizadeh, M. and M. Khonsari, On the role of internal friction in low-and high-cycle fatigue. International Journal of Fatigue, 2018. 114: p. 159-166. Huang, J., Pastor, Marie-Laetitia., Garnier, Christian., Gong, Xiaojing, Rapid evaluation of fatigue limit on thermographic data analysis. International Journal of Fatigue, 2017. 104: p. 293301. Maquin, F. and F. Pierron, Heat dissipation measurements in low stress cyclic loading of metallic materials: From internal friction to micro-plasticity. J Mechanics of Materials, 2009. 41(8): p. 928-942. Yang, W., Guo, Xinglin., Guo, Qiang., Fan, Junling, Rapid evaluation for high-cycle fatigue reliability of metallic materials through quantitative thermography methodology. International Journal of Fatigue, 2019. 124: p. 461-472. Huang, J., Pastor, Marie-Laetitia., Garnier, C., Gong, XJ, A new model for fatigue life prediction based on infrared thermography and degradation process for CFRP composite laminates. International Journal of Fatigue, 2019. 120: p. 87-95. Bär, J., Seilnacht, Luca., Urbanek, Ralf Determination of dissipated energies during fatigue tests on Copper and AA7475 with Infrared Thermography. J Procedia Structural Integrity, 2019. 17: p. 308-315. PANDEY, K.N., Chand, S, Analysis of temperature distribution near the crack tip under constant amplitude loading. J Fatigue Fracture of Engineering Materials Structures, 2008. 31(5): p. 316326. Hutchinson, J., Singular behaviour at the end of a tensile crack in a hardening material. J Journal of the Mechanics Physics of Solids, 1968. 16(1): p. 13-31. Rice, J., Rosengren, Gl F, Plane strain deformation near a crack tip in a power-law hardening material. Journal of the Mechanics Physics of Solids, 1968. 16(1): p. 1-12. Ramberg, W. and W.R. Osgood, Description of stress-strain curves by three parameters. 1943. Ranc, N., Palin-Luc, Thierry., Paris, Paul C, Thermal effect of plastic dissipation at the crack tip on the stress intensity factor under cyclic loading. J Engineering Fracture Mechanics, 2011. 78(6): p. 961-972. Nikishkov, G., An algorithm and a computer program for the three-term asymptotic expansion of elastic-plastic crack tip stress and displacement fields. J Engineering Fracture Mechanics, 1995. 50(1): p. 65-83. Meneghetti, G., Ricotta, Mauro., Atzori, Bruno The heat energy dissipated in a control volume to correlate the fatigue strength of bluntly and severely notched stainless steel specimens. Procedia Structural Integrity, 2016. 2: p. 2076-2083. Cui, Z., Yang, HW., Wang, WX., Yan, ZF., Ma, ZZ., Xu, BS., Xu, HY Research on fatigue crack growth behavior of AZ31B magnesium alloy electron beam welded joints based on temperature distribution around the crack tip. J Engineering Fracture Mechanics, 2015. 133: p. 14-23. Anderson, T.L., Fracture mechanics: fundamentals and applications. 2017: CRC press.

24 | P a g e

25 | P a g e

Highlights   

An analytical relationship between the stress intensity factor and temperature rise for a moving crack tip is proposed. The model is capable of determining the stress intensity factor by measuring the temperature rise at the tip of the crack. Experimental tests for crack speed for SS304 and SS316 are in good agreement with the model predictions.

26 | P a g e

LOUISIANA STATE UNIVERSITY ______________________________________________________________________________________

AND AGRICULTURAL AND MECHANICAL COLLEGE Prof. Michael M. Khonsari, Dow Chemical Endowed Chair in Rotating Machinery Department of Mechanical and Industrial Engineering, Baton Rouge LA 70803-6413 Telephone 225/578-9192 Fax: 225/578-5924 E-mail: khonsari@ lsu.edu

October 7, 2019

We confirm that this is an original research paper and has not been submitted elsewhere for publication. Further, we declare no conflict of interest.

Sincerely,

Michael Khonsari Dow Chemical Endowed Chair and Professor of Mechanical Engineering Fellow ASME, STLE, AAAS; Senior Member NAI Louisiana State University

27 | P a g e