# Numerical Analysis of Joule Heating Behavior and Residual Compressive Stress around Crack Tip under High Electric Load

Authors: Thomas Jin-Chee Liu
Journal: Modelling and Simulation in Engineering
 (2017)

---

## Abstract

This paper discusses the Joule heating effect and residual compressive stress near the crack tip under the electro-thermo-structural coupling state. For the crack tip field, the compressive condition is important for retarding or stopping the crack growth.


---

## Body

## 1. Introduction

The Joule heating behavior and electro-thermo-structural coupled-field around the crack tip were firstly investigated and discussed by Russian scholars in the 1980s [1, 2]. They used mathematical methods to analyze the crack-tip field under the electric current. Following their pioneer research, many typical studies were done for solving this problem in the last decades [3–12].Based on the past references, the Joule heating phenomenon near the crack tip is well known. Figure1 illustrates the concept of this problem. Due to the Joule heating effect, the thermoelectric phenomenon occurs in the conductive material under the electric load. In other words, the material gets hotter when it is subjected to the electric current. If the material has cracks or fractures, the electric current density concentration occurs around the crack tip. Then this electric concentration causes a local hot region at the crack tip due to the Joule heating. In addition, the electric current density has r-1/2 singularity at the crack tip [6]. This is similar to the elastic stress field.Figure 1
Electric current concentration and hot spot at crack tip.Under the Joule heating, the compressive stresses can be produced around the crack tip [1, 2]. This crack-tip compressive stress field is important to reduce or stop the potential crack growth. However, the past references did not show the time-history of the crack-tip stress. During the electric loading, unloading, and cooling process, the crack-tip stress may present the tensile or compressive state. If the plastic strain or deformation occurs at high temperature, the residual stress around the crack tip will be an important topic for the fracture problem.In this paper, the Joule heating behavior and residual stress around the crack tip will be investigated using the electro-thermo-structural coupling finite element analysis. The temperature and electric current density fields will be also obtained for estimating the crack tip behavior. In particular, this primary study will discuss the residual compressive stress and its importance for stopping the crack growth.


## 2. Problem Statement

Figure2 shows the geometric and loading conditions of the steel plate in this study. The plate with an edge crack is subjected to the constant direct current (DC) i0. The electric loading time is te. The plate is made of the mild steel with the dimension W×L×e. The crack length is a. This problem will be simulated by the three-dimensional finite element analysis with the solid model.Figure 2
Geometric and loading conditions.To consider practical conditions, the temperature-dependent material properties in Table1 [9, 13] are adopted in the finite element analysis. Also, the elastoplastic material properties as shown in Figure 3 are considered. Under the Joule heating effect, the electric-current-induced thermo-structural problem is transient. The initial temperature is 21°C. The convection coefficients on all solid/air interfaces are set as 10 W/m2-°C.Table 1
Temperature-dependent physical properties of mild steel [9, 13].
Temperature (°C)
Young’s modulusE (GPa)
Yielding strengthSY (MPa)
Coefficient of thermal expansionα(1/°C)
Thermal conductivityk (W/m-°C)
Specific heatCp (J/kg-°C)
Resistivityρ(Ω-m)
21
206.8
248
10.98 × 10−6
64.60
444
0.14224 × 10−6
93
196.5
238
11.52 × 10−6
63.15
452.38
0.18644 × 10−6
204
194.4
224
12.24 × 10−6
55.24
511.02
0.26670 × 10−6
315.5
186
200
12.96 × 10−6
49.87
561.29
0.37592 × 10−6
426.7
169
173
13.50 × 10−6
44.79
611.55
0.49530 × 10−6
537.8
117
145
14.04 × 10−6
39.71
661.81
0.64770 × 10−6
648.9
55
76
14.58 × 10−6
34.86
762.34
0.81788 × 10−6
760
6.9
14
14.05 × 10−6
30.46
1005.3
1.0109 × 10−6
871
—
—
13.05 × 10−6
28.37
1005.3
1.1151 × 10−6
982
—
—
—
27.62
1005.3
1.1582 × 10−6
1093
—
—
—
28.52
1189.6
1.1786 × 10−6
1204
—
—
—
—
1189.6
1.2090 × 10−6
Poisson’s ratioν = 0.3, density β = 7861.2 kg/m3, melting point = 1521°C.Figure 3
Bilinear elastoplastic stress-strain properties (SIG: stress, units: Pa; EPS: strain, dimensionless;T1~T6: different temperatures, units: °C).The contact condition between crack surfaces is considered as the coupled-field problem. The electric current and heat flow can pass through the crack surfaces when the crack contact occurs. The detailed information of the contact condition will be described in the next section.


## 3. Principles and Finite Element Modelling

### 3.1. Basic Principles

In this paper, the analysis is the electro-thermo-structural coupled-field problem. First, for the electric current field, it obeys the following equations [14]: (1)E=-∇ϕ,J=1ρE,∇·J=0,where E, J, ϕ, and ρ are the electric field (V/m), electric current density (A/m2), electric potential (V), and resistivity (Ω-m), respectively.For the transient thermal analysis, the rules are as follows [14, 15]: (2)q′′=-k∇T,k∇2T+q˙=βCp∂T∂t,q˙=ρJ2,where q′′, k, T, q˙, β, Cp, and t are the heat flux (W/m2), thermal conductivity (W/m-°C), temperature (°C), heat generation (W/m3) of Joule heating, mass density (kg/m3), specific heat (J/kg-°C), and time (s), respectively.The thermoelastic analysis couples the thermal and elastic stress fields as follows [16]:(3)σji,j+Xi=βu¨i,i,j=x,y,zεij=12ui,j+uj,i,i,j=x,y,zεij=1E1+νσij-νI1-EαΔTδij,i,j=x,y,z,where σij, εij, Xi, ui, u¨i,E, ν, I1, α, ΔT, and δij are the stress (Pa = N/m2), strain (dimensionless), body force (N/m3), displacement (m), acceleration (m/s2), Young’s modulus (Pa), Poisson’s ratio (dimensionless), stress invariant (I1=σxx+σyy+σzz), coefficient of thermal expansion (1/°C), temperature difference (ΔT=T-T0 where T0 is the reference temperature) and Kronecker delta, respectively.In this study, the elastoplastic stress-strain behavior is considered. The von Mises yield criterion [17, 18] is used to analyze the plastic stress and strain. The yield surface is defined as(4a)σeqv=SY(4b)σeqv=12σ1-σ22+σ1-σ32+σ3-σ22,where SY, σeqv, and σi (i=1,2,3) are the yielding strength, von Mises equivalent stress, and principal stress, respectively. In this study, the stress-strain relationships in Figure 3 are used.The finite element equations of the electro-thermo-structural coupled-field analysis are as follows [18]: (5)M00000000U¨T¨V¨+C00CtuCt0000U˙T˙V˙+KKut00Kt000KvUTV=FQI,whereU,T,V,F,Q, andIare the vector forms of the displacement, temperature, electric potential, force, heat flow rate, and electric current, respectively. The material constant matrices M, C, Ct, Ctu, K, Kt, Kut, and Kv are the structural mass, structural damping, thermal specific heat, thermo-structural damping, structural stiffness, thermal conductivity, thermo-structural stiffness, and electric conductivity, respectively. The coupled heat flow matrix Q contains the effects of the thermal loading and Joule heating. Ctu and Kut are thermo-structural coupled terms. Equation (5) is a directly coupled nonlinear equation which is solved using the Newton-Raphson iterative method [18].


### 3.2. Boundary and Initial Conditions

Referring to Figure2, the boundary and initial conditions are listed as follows: (6)Jx,-L2,t=i0teW,ϕx,L2,t=0,0≤x≤WTx,y,0=T0=21°Cuix,y,0=u˙ix,y,0=u¨ix,y,0=0.The convective heat transfer on all structure/air interfaces is considered. The convection condition h=10 W/m2-°C and T∞=21°C is applied.


### 3.3. Finite Element Modelling

The software ANSYS is adopted to perform the finite element modelling and calculation. In Figure4, it shows the finite element model of the steel plate with the typical dimensions: L = 240 mm, W =10 mm, e = 1 mm, and a=5 mm. The mesh is constructed by the 20-node solid element SOLID226 which has the capability of the electro-thermo-structural coupled-field analysis. Due to the symmetry, only half thickness is analyzed. In addition, the quarter-point elements are used to simulate the r-1/2 singularity at the crack tip [9–12, 19].Figure 4
Finite element model. (a) Entire view of plate, (b) local region near crack, and (c) three-dimensional view.
(a)
(b)
(c)


### 3.4. Contact Condition on Crack Surfaces

The electro-thermo-structural contact condition on both crack surfaces is considered in the numerical analysis. In ANSYS, two contact element types, TARGE170 and CONTA174, are adopted. Aside from the contact stress, the electric and thermal contacts are considered as the following equations [18, 20]: (7)J=ηcelϕ1-ϕ2q′′=ηcthT1-T2,where ηcel and ηcth are, respectively, the electric conductance and thermal conductance of the contact surfaces. The terms (ϕ1-ϕ2) and (T1-T2) are, respectively, the electric potential difference and temperature difference between both contact surfaces. In this study, ηcth is assumed to be very large so that the thermal contact resistance on crack surfaces can be ignored. However, ηcel is defined as follows [20]:(8)ηcel=1ρclc,where ρc and lc are the electric contact resistivity and characteristic length, respectively. For mild steel, the typical values ηcel=6.29×1081/(Ωm2), ρc = 6.2586 × 10−5 Ωm, and lc = 2.54 × 10-5m (0.001 inch) [13] are used. In addition, the coefficient of friction on crack surfaces is neglected.


## 3.1. Basic Principles

In this paper, the analysis is the electro-thermo-structural coupled-field problem. First, for the electric current field, it obeys the following equations [14]: (1)E=-∇ϕ,J=1ρE,∇·J=0,where E, J, ϕ, and ρ are the electric field (V/m), electric current density (A/m2), electric potential (V), and resistivity (Ω-m), respectively.For the transient thermal analysis, the rules are as follows [14, 15]: (2)q′′=-k∇T,k∇2T+q˙=βCp∂T∂t,q˙=ρJ2,where q′′, k, T, q˙, β, Cp, and t are the heat flux (W/m2), thermal conductivity (W/m-°C), temperature (°C), heat generation (W/m3) of Joule heating, mass density (kg/m3), specific heat (J/kg-°C), and time (s), respectively.The thermoelastic analysis couples the thermal and elastic stress fields as follows [16]:(3)σji,j+Xi=βu¨i,i,j=x,y,zεij=12ui,j+uj,i,i,j=x,y,zεij=1E1+νσij-νI1-EαΔTδij,i,j=x,y,z,where σij, εij, Xi, ui, u¨i,E, ν, I1, α, ΔT, and δij are the stress (Pa = N/m2), strain (dimensionless), body force (N/m3), displacement (m), acceleration (m/s2), Young’s modulus (Pa), Poisson’s ratio (dimensionless), stress invariant (I1=σxx+σyy+σzz), coefficient of thermal expansion (1/°C), temperature difference (ΔT=T-T0 where T0 is the reference temperature) and Kronecker delta, respectively.In this study, the elastoplastic stress-strain behavior is considered. The von Mises yield criterion [17, 18] is used to analyze the plastic stress and strain. The yield surface is defined as(4a)σeqv=SY(4b)σeqv=12σ1-σ22+σ1-σ32+σ3-σ22,where SY, σeqv, and σi (i=1,2,3) are the yielding strength, von Mises equivalent stress, and principal stress, respectively. In this study, the stress-strain relationships in Figure 3 are used.The finite element equations of the electro-thermo-structural coupled-field analysis are as follows [18]: (5)M00000000U¨T¨V¨+C00CtuCt0000U˙T˙V˙+KKut00Kt000KvUTV=FQI,whereU,T,V,F,Q, andIare the vector forms of the displacement, temperature, electric potential, force, heat flow rate, and electric current, respectively. The material constant matrices M, C, Ct, Ctu, K, Kt, Kut, and Kv are the structural mass, structural damping, thermal specific heat, thermo-structural damping, structural stiffness, thermal conductivity, thermo-structural stiffness, and electric conductivity, respectively. The coupled heat flow matrix Q contains the effects of the thermal loading and Joule heating. Ctu and Kut are thermo-structural coupled terms. Equation (5) is a directly coupled nonlinear equation which is solved using the Newton-Raphson iterative method [18].


## 3.2. Boundary and Initial Conditions

Referring to Figure2, the boundary and initial conditions are listed as follows: (6)Jx,-L2,t=i0teW,ϕx,L2,t=0,0≤x≤WTx,y,0=T0=21°Cuix,y,0=u˙ix,y,0=u¨ix,y,0=0.The convective heat transfer on all structure/air interfaces is considered. The convection condition h=10 W/m2-°C and T∞=21°C is applied.


## 3.3. Finite Element Modelling

The software ANSYS is adopted to perform the finite element modelling and calculation. In Figure4, it shows the finite element model of the steel plate with the typical dimensions: L = 240 mm, W =10 mm, e = 1 mm, and a=5 mm. The mesh is constructed by the 20-node solid element SOLID226 which has the capability of the electro-thermo-structural coupled-field analysis. Due to the symmetry, only half thickness is analyzed. In addition, the quarter-point elements are used to simulate the r-1/2 singularity at the crack tip [9–12, 19].Figure 4
Finite element model. (a) Entire view of plate, (b) local region near crack, and (c) three-dimensional view.
(a)
(b)
(c)


## 3.4. Contact Condition on Crack Surfaces

The electro-thermo-structural contact condition on both crack surfaces is considered in the numerical analysis. In ANSYS, two contact element types, TARGE170 and CONTA174, are adopted. Aside from the contact stress, the electric and thermal contacts are considered as the following equations [18, 20]: (7)J=ηcelϕ1-ϕ2q′′=ηcthT1-T2,where ηcel and ηcth are, respectively, the electric conductance and thermal conductance of the contact surfaces. The terms (ϕ1-ϕ2) and (T1-T2) are, respectively, the electric potential difference and temperature difference between both contact surfaces. In this study, ηcth is assumed to be very large so that the thermal contact resistance on crack surfaces can be ignored. However, ηcel is defined as follows [20]:(8)ηcel=1ρclc,where ρc and lc are the electric contact resistivity and characteristic length, respectively. For mild steel, the typical values ηcel=6.29×1081/(Ωm2), ρc = 6.2586 × 10−5 Ωm, and lc = 2.54 × 10-5m (0.001 inch) [13] are used. In addition, the coefficient of friction on crack surfaces is neglected.


## 4. Results and Discussions

### 4.1. Validation of Finite Element Model

The accuracy of the finite element model must be valid. For the validation, the coupled-field analysis is simplified to the pure electric problem. Furthermore, the SOLID226 elements are degenerated to SOLID231 elements for the ANSYS model in Figure4. In this case, the resistivity is 7.2∗10-7 Ω-m. Also, the distributed current J0=6.25∗106 A/m2 is used to replace the concentrated current i0.Using the limited electric potential extrapolation technique (LEPET) [9], the electric current density factor KJ at the crack tip can be obtained from the finite element results. In this case, the numerical result of KJ is 8.857∗105 Am-3/2. According to the analogy and analytical methods [8, 9, 21], the analytical solution of the edge crack problem under electric load J0 is(9)KJ=BJ0πa,where B=1.128 for the case of a/W=0.5. Substituting all values to (9), KJ is 8.836∗105 Am-3/2. By comparing both values of KJ, the finite element result has good agreement with the analytical solution with the small numerical error 0.24%.


### 4.2. Electric Concentration and Hot Spot at Crack Tip

The Joule heating effect makes high temperature field in the steel plate. Due to the local electric concentration at the crack tip, a hot spot exists. In Figure5, it shows the electric current density and temperature fields around the crack under the following conditions: L=240 mm, W=10 mm, e=1 mm, a=5 mm, and i0=950 A. The electric loading time is te=0.3333 s. The numerical results obviously prove the phenomena of the electric concentration and hot spot at the crack tip.Figure 5
(a) Electric current density field (units: A/m2). (b) Temperature field (units: °C).
(a)
(b)


### 4.3. Residual Compressive Stress at Crack Tip

In this section, the geometric and loading conditions of the steel plate are the same as the above section. As shown in Figure6, the unloading and cooling processes are added in the analysis to obtain the residual stress around the crack tip. The loading process is during t = 0~0.3333 s. Then the electric current is removed and it remains the cooling process till t=600 s.Figure 6
Electric loading, unloading, and cooling history.In Figure7, the time-history of the crack tip temperature is shown. It can be seen that the crack tip temperature reaches the maximum value (969°C) at the end of the loading time (te=0.3333 s). After the unloading process, the temperature decays quickly and approaches the steady state. At t=600 s, the crack tip temperature decreases to 25°C.Figure 7
Time-history of crack tip temperature (t = 0~600 s).Figures8–10 show the time-history of the normal stress σy at the point B near the crack tip. The distance from B to the crack tip is 0.2 mm. In Figure 8, it shows the negative values of σy during the cooling process till t=600 s. Also, it approaches the steady state. At t=600 s, the value of σy is −301 MPa. The negative value implies the compressive normal stress. In other words, the residual compressive stress exists near the crack tip after the electric load is removed.Figure 8
Time-history of stress near crack tip (t = 0~600 s).Figure 9
Time-history of stress near crack tip (t = 0~30 s).Figure 10
Time-history of stress near crack tip (t = 0~2.5 s).In Figures9 and 10, the detailed stress values during shorter time intervals are shown. The stress causes complicated fluctuation during the electric loading and unloading processes. Then it decays and approaches the steady state with sufficient cooling time.Figure11 shows the contour of the stress σy near the crack tip at the end of the cooling process (t=600 s). There are negative stress values around the crack tip. It means that the residual compressive stresses occur in this area. According to the results of Figure 11, Figure 12 shows the stress distribution in front to the crack tip. The compressive stress field is also investigated near the crack tip.Figure 11
Contour of stressσy near crack tip (t=600 s).Figure 12
Stress distribution in front to crack tip (t=600 s).


### 4.4. Effects of Electric Load

Figure13 shows the time-history of the stress at point B under i0=950 A and i0=400 A. Comparing two curves in the figure, it is found that there is no residual stress under i0=400 A. It implies that the magnitude of the electric load is an important parameter for producing the residual compressive stress near the crack tip. If the electric load is not sufficient, the temperature and stress will be too small to make the yielding (plastic) strain and residual stress.Figure 13
Time-history of stress near crack tip under different electric loads.


## 4.1. Validation of Finite Element Model

The accuracy of the finite element model must be valid. For the validation, the coupled-field analysis is simplified to the pure electric problem. Furthermore, the SOLID226 elements are degenerated to SOLID231 elements for the ANSYS model in Figure4. In this case, the resistivity is 7.2∗10-7 Ω-m. Also, the distributed current J0=6.25∗106 A/m2 is used to replace the concentrated current i0.Using the limited electric potential extrapolation technique (LEPET) [9], the electric current density factor KJ at the crack tip can be obtained from the finite element results. In this case, the numerical result of KJ is 8.857∗105 Am-3/2. According to the analogy and analytical methods [8, 9, 21], the analytical solution of the edge crack problem under electric load J0 is(9)KJ=BJ0πa,where B=1.128 for the case of a/W=0.5. Substituting all values to (9), KJ is 8.836∗105 Am-3/2. By comparing both values of KJ, the finite element result has good agreement with the analytical solution with the small numerical error 0.24%.


## 4.2. Electric Concentration and Hot Spot at Crack Tip

The Joule heating effect makes high temperature field in the steel plate. Due to the local electric concentration at the crack tip, a hot spot exists. In Figure5, it shows the electric current density and temperature fields around the crack under the following conditions: L=240 mm, W=10 mm, e=1 mm, a=5 mm, and i0=950 A. The electric loading time is te=0.3333 s. The numerical results obviously prove the phenomena of the electric concentration and hot spot at the crack tip.Figure 5
(a) Electric current density field (units: A/m2). (b) Temperature field (units: °C).
(a)
(b)


## 4.3. Residual Compressive Stress at Crack Tip

In this section, the geometric and loading conditions of the steel plate are the same as the above section. As shown in Figure6, the unloading and cooling processes are added in the analysis to obtain the residual stress around the crack tip. The loading process is during t = 0~0.3333 s. Then the electric current is removed and it remains the cooling process till t=600 s.Figure 6
Electric loading, unloading, and cooling history.In Figure7, the time-history of the crack tip temperature is shown. It can be seen that the crack tip temperature reaches the maximum value (969°C) at the end of the loading time (te=0.3333 s). After the unloading process, the temperature decays quickly and approaches the steady state. At t=600 s, the crack tip temperature decreases to 25°C.Figure 7
Time-history of crack tip temperature (t = 0~600 s).Figures8–10 show the time-history of the normal stress σy at the point B near the crack tip. The distance from B to the crack tip is 0.2 mm. In Figure 8, it shows the negative values of σy during the cooling process till t=600 s. Also, it approaches the steady state. At t=600 s, the value of σy is −301 MPa. The negative value implies the compressive normal stress. In other words, the residual compressive stress exists near the crack tip after the electric load is removed.Figure 8
Time-history of stress near crack tip (t = 0~600 s).Figure 9
Time-history of stress near crack tip (t = 0~30 s).Figure 10
Time-history of stress near crack tip (t = 0~2.5 s).In Figures9 and 10, the detailed stress values during shorter time intervals are shown. The stress causes complicated fluctuation during the electric loading and unloading processes. Then it decays and approaches the steady state with sufficient cooling time.Figure11 shows the contour of the stress σy near the crack tip at the end of the cooling process (t=600 s). There are negative stress values around the crack tip. It means that the residual compressive stresses occur in this area. According to the results of Figure 11, Figure 12 shows the stress distribution in front to the crack tip. The compressive stress field is also investigated near the crack tip.Figure 11
Contour of stressσy near crack tip (t=600 s).Figure 12
Stress distribution in front to crack tip (t=600 s).


## 4.4. Effects of Electric Load

Figure13 shows the time-history of the stress at point B under i0=950 A and i0=400 A. Comparing two curves in the figure, it is found that there is no residual stress under i0=400 A. It implies that the magnitude of the electric load is an important parameter for producing the residual compressive stress near the crack tip. If the electric load is not sufficient, the temperature and stress will be too small to make the yielding (plastic) strain and residual stress.Figure 13
Time-history of stress near crack tip under different electric loads.


## 5. Conclusions

From the finite element results, the electric current density concentrates at the crack tip. Due to the Joule heating, it causes a hot spot at the crack tip. The residual compressive stress appears near the crack tip due to the high temperature and plastic deformation. Furthermore, the compressive condition can retard or stop the crack growth. The concept for stopping crack growth is shown in Figure14.Figure 14
Concept for stopping crack growth.This paper provides a primary study and conclusion of the residual compressive stress near the crack tip under the electro-thermo-structural coupling state. The compressive condition is practically important to the fracture mechanics problem.It is important to study the fatigue crack growth due to the Joule heating caused stresses. The extending work will be considered in the future research.



---

*Source: 3012949-2017-10-19.xml*