An investigation of the conditions leading to strip adhesion in industrial steel coils

IHIEI'iI!JlIiiE 1E);,I§IIllVilNhlRllWll.<G

ONllllE

R

1

By

Lulama Velile Wakaba

An Investigation of the Conditions Leading to Strip

Adhesion in Industrial Steel Coils

This thesis is submitted in accordance with the requirements for the degree

Magister Scientiae

In the Faculty of Natural and Agricultural Sciences

Department of Physics

at the

University of the Free State

Bloemfontein

Study leader: Dr M.F. Maritz

Co-study leader: Prof. H. C. Swart

Summary

During the production of steel strip, a significant amount of work hardening takes place when the steel is rolled into thin strips, which are stored in a coiled form. These steel coils are batch annealed in order to reduce the hardness and restore formability, before further processing takes place. The development of diffusion welds between spirals of steel coils, during batch annealing, is of particular interest because it prevents the coils from being unwound for further use. This problem is often referred to as strip adhesion or stickering.

In order for strip adhesion to develop, it is essential for some coil spirals to be in contact, while the inter-diffusion between spirals takes place. Furthermore, high temperatures also aid in the diffusion process. Itis therefore useful to study the temperature and resulting thermal stress distributions in the coil, during batch annealing. The temperature distribution allows for the calculation of thermal stress, which is the driving force for establishing contact between spirals, and in addition to this, the temperature distribution also provides some clues regarding the likelihood of inter-diffusion.

In this study, models of temperature and stress are presented. A two-dimensional finite difference model for temperature is presented and confirmed by an analytical solution of the same problem. Analysis of a three-dimensional temperature model in the third chapter shows that, as far as heat transfer is concerned, a cylindrical coil can be well approximated by a solid cylinder with a concentric hole. All the temperature modeling was done for the interior of a coil. Further discussion in later chapters shows that a cylindrical coil can also be treated as a solid cylinder for thermal stress modeling.

A long cylinder stress model in the fifth chapter provides some useful insight, as far as strip adhesion is concerned, even though it is a one-dimensional model that does not consider the effect of axial heat transfer. The radial compressive stress during the cooling

The thermal stress calculations are later extended to include a linear cooling temperature ramp and these results are arguably the most important findings of this study. According to these results, when the a cylinder is cooled in such a way that the outer edge lags behind the cooling of the inner edge, by a few hours, the compressive radial stress is greatly reduced. Consequently, the contact pressure between spirals at the most critical stage of batch annealing (where strip adhesion occurs) is decreased and the chance of strip adhesion developing is reduced.

• Dr M.F. Maritz for his invaluable expertise and assistance.

• Prof. H.C. Swart for his insight and guidance.

• Dr c.J. Greyling for his wealth of information and support.

• The University of the Free State and its staff for providing an academically nurturing

environment.

Acknowledgements

1.1BACKGROUND 4

1.2BATCH ANNEALING PROCESS 5

1.3DEFOX PLUS PROCESS 7

1.4OBJECTIVES OF THIS STUDY 8

1.5 SCOPE OF THESIS 9

TABLE OF CONTENTS

CHAPTER

1

4

INTRODUCTION

4

CHAPTER

2

10

TEMPERATURE

MODEL FOR A SOLID CYLINDER

I0

2.1 INTRODUCTION ; 10

2.2 SOLID CYLINDER MODEL Il

2.3 INITIAL AND BOUNDARY CONDITIONS 13

2.4 FINITE DIFFERENCES APPLIED TO THE SOLID CYLINDER 14

2.5 STABILITY CONDITIONS 16

2.6 TEMPERATURE DEPENDENCE OF THE RADIAL AND AxIAL THERMAL DIFFUSIVITY 20

2.7 NUMERICAL RESULTS FORA SOLID CYLINDER 22

2.7.1 Temperature ramp 26

2.8 ANALYTICAL SOLUTION FOR SOLID CYLINDER 26

2.9 ANALYTICAL RESULTS 37

2.10 COMPARISON OF ANALYTICAL AND NUMERICAL SOLUTIONS 38

CHAPTER

3

42

TEMPERATURE

MODEL FOR A CYLINDRICAL

COIL

42

3.1 INTRODUCTION 42

3.2 FINITE DIFFERENCE MODEL 42

3.2.1 Numerical Grid 43

3.2.2 Discretized Heat Equation 45

3.3 INITIAL CONDITIONS 47 3.4 NUMERICAL RESULTS 48 3.5 THETA DEPENDENCE 52

CHAPTER

4

60

ELASTICITY

THEORY

60

4.1 INTRODUCTION 60

4.5 DERIVATION OF EQUILIBRIUM CONDITIONS 67

4.6 LAGRANGlAN AND EULERIAN DESCRIPTIONS OF DEFORMATION 69

4.7 DEFORMATION AND STRAIN TENSORS 71

4.8 COMPATIBILITY EQUATION 74

4.9 HOOKE'S LAW FOR ISOTROPIC MEDIA 75

4.10 THERMAL STRESS 78

4.11 EQUILIBRIUM CONDITIONS IN POLAR COORDINATES 80

4.12 EQUATIONS OF STRAIN IN POLAR COORDINATES 84

CHAPTER 5

87

LONG CYLINDER STRESS MODEL

87

5.1 INTRODUCTION 87

5.2 FORMULATION OF MODEL 88

5 .3 INCORPORATION OF AN ANAL YTICAL TEMPERATURE DISTRIBUTION 93

5.4 CONSTANT HEATING RESULTS 96

5.5 CONSTANT COOLING RESULTS 99

5.6 STEADY STATE SOLUTION 103

5.7 LINEAR COOLING 108

CHAPTER 6

114

CONCLUSION

114

APPENDIX A

117

A.1 FINITE DIFFERENCE METHOD 117

APPENDIX B

119

B.l STARTER. M 119 B.2 BESSOLRO.M 120 B.3 FUNCTION BESSOLF.M 120 B.5 FINE2D.M 121 B.6 FUNCTION RAMPl.M 122 B.7 FINRAMP.M 122 B.8 TEMPRAMP.M 124 B.9 BESSOL321.M 124 B.I0 BESSOLCOMP.M 127

APPENDIX C

131

C.1 THERMD.M 131

APPENDIX D

135

D.1 PARTIAL DERIVATIVES 13 5

E.1 FUNCTION RAMPCOOL.M 137

Chapter 1

Introduction

1.1 Background

During the production of steel strip, a significant amount of work hardening takes place when the steel is rolled into thin strips, which are stored in a coiled form. These steel coils are batch annealed in order to reduce the hardness and restore formability, before further processing takes place. The development of diffusion welds between spirals of steel coils, during batch annealing, is of particular interest because it prevents the coils from being unwound for further use. This problem is often referred to as strip adhesion or stickering.

Most of the work done in this field deals with aspects related to heat transfer. Perrin and Johnson [1] discuss how the overall design of annealing systems impacts some aspects of aerodynamics and heat transfer of the circulating gases during, batch annealing. The thermal and physical properties of stacks of cold-rolled steel are investigated by Lisogor and Mitkalinnyi [2]. Some computer modeling and experimental results regarding the prediction of heating, soaking and cooling times in batch annealing, are presented by Rao, Barth and Miller [3]. Rovito [4] also provides a computer-based model for predicting the annealing time.

Leroy and Bouquegneau [5] discuss the possibilities and limits of batch annealing with 100% hydrogen. This reference also contains sections dealing with the gas-metal reactions and the development of radial thermal stress that take place throughout the annealing process. A more in-depth study of the gas-metal reactions can be obtained from a United States Patent by Zylla [6], which introduces a process for preventing the formation of adhesives when annealing steel coils (this is also known as the Defox Plus process).

1.2 Batch Annealing Process

In a typical batch annealing process, several coils are annealed in a bell-shaped furnace

(as shown in figure 1.1) and a reducing gas, i.e. hydrogen or a nitrogenlhydrogen mixture, is passed through the coils, in a circular fashion, to prevent oxidation [6]. The heat is supplied from outside the inner cover by means of a heater that fits over the system. The reducing gas circulates in the manner shown by the arrows in figure 1.1 (the dotted lines on each coil represent the inner ring). When the gas leaves the circulating fan, it moves along the heated walls of the inner cover and thereby receives heat. In the ensuing period, the gas then moves down through the coil center and interior, thus heating the coil from interior as well. The convector plates enhance gas circulation though the coils.

During the batch annealing process, heating occurs in the form of a temperature ramp, which increases to a maximum temperature of about 670°C before decreasing to room temperature. According to experimental findings, strip-adhesion usually takes place at the critical time interval shown on the temperature ramp in figure 1.2 [5]. The critical time interval is characterized. by a steep decreasing thermal gradient and high thermal stresses.

~~

~---~

~ I I

-t

I ~~ I I I

t

I I I I I I I I I I .... I I "" I I

t

i~

_~ II

t

I ~ I I I I I I I I I I

!~

~ I

t

t

I I I I I I

-

-

I lEE ~

I

I

i

Circulation fan

Figure 1.1: An illustration of the batch annealing process.

Coil

Convector

Inner

Critical time interval 700r--'-;~==~==~==~~~~~--'--'--~

_,

I I

,

---p---

I I \

,

,

_{,-_}

.... 500 100 600 200

O~--~--~--~--~--~----~--~--~--~--~

o 5 10 15 20 25 30 35 40 45 50 lime (hours)

Figure 1.2: A typical temperature ramp with an indicated critical time interval. The dotted line at 600°C has relevance to the Defox Plus process that is discussed in section

1.3.

1.3 Defox Plus Process

The Defox Plus process, which was developed by Peter Zylla [6], claims to prevent strip adhesion by circulating a gas mixture, during batch annealing, that oxidizes the coil spirals at temperatures above 600°C in the heating phase and reduces them below 600°C, in the cooling phase. This process was designed in such a way that an oxide layer is present on the spirals during the critical time interval, in order to discourage adhesion. After undergoing reduction, the coils emerge as desired. This process can be understood largely by studying the gas-metal reactions that take place during the oxidation and reduction stages.

1.4 Objectives of this study

-times in the temperature cycle. Compressive radial thermal stresses, resulting temperature gradients within the coils, can cause successive spirals to establish contact, thereby constricting the flow of gas and hindering particularly the reduction processes, in the critical time interval. If the coils were annealed in an oxidation-preventative gas atmosphere (i.e. an inert gas), radial compressive stresses would increase the chance of adhesion between spirals. It is therefore necessary to know the relevant thermal stress distributions for the annealing process. Furthermore, the thermal stress distributions can only be calculated once the temperature distribution is known.

This study was undertaken with three goals in mind. The first objective was to provide models of temperature and thermal stress that would shed more light on the conditions

• <

experienced by industrial steel coils during batch annealing. The second aim was to use the knowledge gained from the temperature and thermal stress models to formulate a suggestion that would reduce the likelihood of strip adhesion developing. The third, fmal and most important objective has been to expose the author to numerical and analytical techniques in a research environment. For this reason, no automatic software packages were used to solve equations.

With regards to the modeling of temperature, it should be noted that various models of batch annealing already exist. Some provide predictions of annealing times [4], whilst others consider such factors as: the gas flow rate, the thermal conductivity of the circulating gas, the effect of radiation (from the inner cover) and convection on the eventual temperature distributions in each coil [3]. The entire temperature topic, in this study, will focus on heat transfer modeling within a coil, based on supplied external boundary conditions.

1.5 Scope of thesis

A numerical solution to a solid cylinder (a solid cylinder with a concentric hole) model for heat transfer is developed in Chapter 2, as a first approximation to a cylindrical coil. This is subsequently followed by an analytical solution and comparison that confirms the numerical result. In Chapter 3, a numerical model for a three-dimensional coil is formulated and a theta-dependence test is performed, for the purpose of establishing the difference between heat transfer in a cylindrical coil and a solid cylinder. This has bearing on the computational speed at which temperature can be modeled in a cylindrical coil. Stability conditions for the numerical solutions are also derived.

Certain aspects of the general theory of elasticity are introduced in Chapter 4. Some fundamental concepts are treated in this chapter, and it forms a basis for the one-dimensional thermal stress model for a solid cylinder in Chapter 5. This is followed by a short discussion regarding temperature and stress environments that encourage and discourage the adhesion of spirals. Appearing in this chapter, as well, is a proposed reason for strip adhesion occuring during cooling and a suggestion that reduces the compressive radial stresses in the critical time interval. A final conclusion is made in

Chapter 2

Temperature Model for a Solid Cylinder

2.1 Introduction

In this chapter, the heat equation for axial symmetry is solved both analytically and numerically, with the intention of using the analytical solution to verify the numerical solution. The first section of this chapter begins with a description of how heat transfer in a solid cylinder (with a concentric hole) is treated. This is followed by a segment that deals with the relevant initial and boundary conditions.

Section 2.4 shows how the finite difference method can be applied to the heat equation. A derivation of the numerical stability conditions follows thereafter, just before the results are displayed. Included in the results, is a short discussion of the cold region that develops during heating. The analytical solution is derived in section 2.8 and the corresponding results are displayed in section 2.9. A comparison between numerical and analytical solutions is subsequently performed.

2.2 Solid Cylinder Model

The study of heat transfer in a solid cylinder can be considered as an approximation to heat transfer in a cylindrical coil. It is important to note that such an approximation depends on the assumption that there is no heat flow along the spirals of the coil (see figure 2.1). This assumption can also be re-stated by saying that the heat flow in the coil has no B-dependence.

(a) (b)

Figure 2.1: (a) Cylindrical coil. (b) Solid cylinder with a concentric hole and different thermal diffusivity values in the radial and axial directions.

A solid cylinder model has to account for the different thermal coefficients

k.

(in the z-direction of the material) and kr (in the radial direction). Since a steel coil, that undergoes

batch annealing, has a film of gas and oil between each spiral, it can be expected that k; and k= will be different [3].

In order to model the evolution of temperature in the steel coil during the batch annealing process, use will be made of the well-known heat equation [7]. In cartesian coordinates, the heat equation is,

where Tis the temperature, tis time and _ k; «, ---, pCp kv Kv =-"-, " pCp and k_ K =---z pCp.

Note that Kx, Kyand K=are referred to as the thermal diffusivity values for the x,y andz

directions. The density of the material is represented by p, Cp is the specific heat under constant pressure and kis the thermal conductivity.

For the purpose of the solid cylinder problem, it is more useful to transform (2.1) into cylindrical coordinates. In cylindrical coordinates, the following is obtained:

(2.2) where

kB

KB =--,

=.

(2.3) and k_ K_

=---.

- pCp

(2.4)

Since the solid cylinder is axially-symmetric, the temperature distribution of the whole cylinder will be known if the temperature distribution over cross section A is known (see figure 2.2). This reduces the 3-D problem to a 2-D problem

....

r-~

f'...

I ~ I I I I I ,

...

,--«.:_:"-~-~

A A

...

r

...

Figure 2.2: A diagram of an axially symmetric solid cylinder that has a temperature distribution which isfully described by cross section A.

2.3 Initial and Boundary Conditions

As a first approximation, assume that the batch annealing process occurs at a constant temperature of 670°C. Now, in order to test whether this model can predict the existance of a cold spot, the temperature 670°C will be applied on all external surfaces of the cylinder and the time evolution of the temperature will be calculated. Note that the initial temperature of the interior is at room temperature, which will be taken here as 23°C. The boundary condition will be taken as constant and equal to 670°C for any time greater than zero.

Boundary Temperature: 670°C

!

670'C

670°C

..

23°C (inside)

...

670°C

A

t

670'C

Figure 2.3: Initial and boundary conditions of solid cylinder.

The finite difference method will be used to solve equation (2.4) for the above mentioned initial conditions. An introduction of the finite difference method can be seen in appendix A.

2.4 Finite Differences Applied to the Solid Cylinder

The heat equation (2.4) needs to be solved for surface A (see figure 2.3). In order to use the finite difference technique, surface A needs to be discretized into a matrix of points. The step-length and index for the r-direction ish; and i,respectively. For the z-direction,

h= is the step-length andj is the index. The index for time is n, while trepresents the corresponding step-length. A temperature T\,j is therefore associated with each point (ij) at some time ID. Since ih;

=

r.jh,

=

zand nr= t, it follows that T\,j is an approximation of

Figure 2.4:A diagram of cross sectionA,after discretization.

When dealing with a function of two or more independent variables, the partial derivatives are approximated by partial differences. Therefore, by using partial differences,.the following approximations can be made:

a2T ::::::(Tni+l,J -2Tni.J +Tni_I,J]

ar

2 h 2 '

r

(2.5)

aZT ::::::(Tni'J+1 - 2Tni,J

+

Tni,j_l]

a

2 2 '

'Z h.

(2.6)

(2.7)

Tn+1"-Tn" [I_[Tni+l,J-Tni-l,J: [Tni+I,J-2Tni,J+Tni-l,JJ

I,} - I,}

+

T Kl'

+

Kl' 2

ri 2hl' hl'

(2.9)

[Tni,J+l -2Tni,J +Tni,J_1

JJ

+K_ 2

_ h:

According to the initial conditions, all the points on the boundary are at 670°C, while all the interior points are at room temperature (23°C). As discussed previously, the initial temperature distribution over surfaceA has a matrix of the form:

670 670 670 670 670

670 23 23 23 670

T

=

670 23 23 23 670

670 23 23 23 670

670 670 670 670 670

Equation (2.9) is applied to the interior matrix. elements, while the boundaries remain at 670°C. The calculation begins atn = 1,i =2 and j =2. All the interior matrix elements

are updated with each iteration. In other words, T\J is updated for each value of n, as time progresses. The stability conditions for this numerical scheme is investigated below.

2.5 Stability Conditions

The numerical scheme (2.9) is consistent with (2.4) meaning that as h-, h- and t:tend to zero, each of the difference replacements tend to the corresponding derivatives. However, this is not enough to ensure that a solution of (2.9) is a good approximation to a solution of (2.4), for the same initial and boundary conditions. An extra condition, the so-called "von Neumann stability condition" needs to be satisfied [8]. This is a condition that imposes a relationship between the magnitude of the time and space step sizes.

Kr

which is positioned in front of the term

[ Tili+l,} - Tili-I,}

J

2hr

maintains its largest value. With this assumption, (2.9) is a linear difference equation and will admit solutions of the type:

Til ;:" i(i(1),+}OJ_)

i,}

=

':> e . , . _(2.10)

A solution of the type (2.10) is called a mode. Von Neumann stability implies that T,h.,

and h- have such a relationship that no single mode will grow exponentially. Otherwise stated,

or

(2.11)

K j:11ei(i(U,+j{u:): :

+

"'='? (elllJ,. -2+e-IIlJ:) h,.-(2.12) K j:l1ei(iOJ,+j(U:): :]

+

z'=' (ei"': - 2+e-/(V:) h 2 z

After some simplification, (2.12) can be written as,

(2.13) If identities - -eieu _ e-iOJ sinew)

=

_

2i (2.14) and - -e'"

+

e-iOJ cosu»)

=

---2 (2.15)

are used, equation (2.13) can be expressed as:

C;=l+T[ K)

sin(w,.)+ 2~,. (cos(w,.)-l)+

2~z

(cOS(WJ-1)].

h~~ h k (2.16)

where (KI' K _) A

=

-11-

k2

(2.18) and B-

_---

KI' hrmin (2.19) According to (2.17), and

This implies that scheme (2.9) is stable if

(2.20)

Equation (2.20) shows that the stability condition also depends on Kr and K=. These

values, in general, depend on temperature. Some temperature dependence curves follow hereafter.

2.6 Temperature Dependence of the Radial and Axial Thermal

Diffusivity

A publication by Rao, Barth and Miller [3] shows some tabulated values of the thermal coefficients k.,

k,

and the specific heat Cp as a function of temperature. Graphical plots of these values can be seen in figures 2.5a, 2.5b and 2.5e. The calculated values for Kr and K= are plotted in figures 2.5d and 2.5e. These values were calculated using equation (2.3)

in conjunction with the above-mentioned data. A constant density of p = 7860 kg.m"

was used. x lOs x lO~ 2,1 1. 20 ~E ';'~

...

_•

~2,0 ~1.15 i- ₊ ~1.9 ~_1.10 .2

...

_.2

...

2 1.8 + 2 _{1. 05}

...

",.

...

.>t d 1.7 -~

...

'i

•

1. 00

..

...

:a 1.6 _{-~ 0,95}

...

..

"... d _1.5

•

-s 0,90

•

'" £7

_...

_£7

•

:~1.4 ~0,85

•

..

il1.3

..

1

8 0,80

•

+

i

1.2

i

0,75

i-...

~_{Ll ~} i-0,70 0 100 200 300 400 500 600 700 800 0 100 200 300 400 500 600 700 800 Temperature rC) Temperature fc) (a) (b)

0.055 0 0.05 0 0 0.045 0 1) 0 ~ 0.04 ₀ 2' 'l' 0 '" 0.035 ₀

I

0 0.03 0 0 0.025 0 0.02 0 0 0.015 0 100 200 300 400 500 600 700 800 Temperature ("C) X10.3 1000 ₊ 2.2 2.1 900 0 0 0 ₀ + 0 0 0 1) 1)r 2 0 0 .9 800 ₊ ~ 0 r 2' ~ 'r"1.9 0 ~ +

r

700 +

I

0 1.8 +

i

600 + + 1.7 + 0 + + 500 + + + 1.6 0 400 _1.5 0 100 200 300 400 500 600 700 800 _{0 100 200 300 400 500 600 700 800} Temperature ('C) _{Temperature (0C)} (C) (d) (e)

Figure 2.5: (a) Axial thermal conductivity as a function of temperature. (b) Radial thermal conductivity versus temperature. (c) Specific heat against temperature. (d) Radial thermal diffusivity as a function of temperature. (e) Axial thermal diffusivity versus temperature.

2.7 Numerical Results for a Solid Cylinder

A program-calledfine2d.m, that implements equation (2.9) for the initial and boundary conditions of section 2.3, was written in MATLAB (version 4.2c.l). This program IS

listed in appendix B.S. Mesh diagrams of the results are shown in figure 2.6.

lime= 0hours ). 600 "I' 400

;

e 200 ~ 0 lime=1hour (a) (b) lime =3 hours Time = 5 hours 700 600 700 500 600 ~ 400 _~ 500 400

i

300

i

300 200 ~ 200 ~ ₁₀₀ 100 0 1 0 1 0 0 _{0 0} z (c) (d)

lime=3 hours _{llme =3 hours}

~

J :

~-+-____;IN

560 1---4--I--'I.."\.

(e) (f)

Figure 2.6: (a) to (d) Temperature profiles for various time intervals. (e) A cross section in the z-direction after 3 hours. (f) A cross section in the r-direction after 3 hours.

It can be seen from figure 2.6 that the temperature distribution in the z-direction is symmetrical, while the r-direction exhibits a slightly anti-symmetrical temperature distribution. The symmetry in the z-direction can be expected, since the cylinder experiences the exact same conditions at the top and bottom. The reason for the lack of symmetry in the r-direction can be explained by the following physical arguments.

Top view of cylinder

Outer wall Inner wall

Region 1

Figure 2.7: Top view of a solid cylinder with imaginary sectors.

If one considers, say 8, equidistant regions on the outside wall of the cylinder and 8 equidistant regions on the inner wall (see figure 2.7): The amount of heat entering from

the inner regions is not the same as the heat entering from the outer regions.

The heat entering per unit area is the same on both inner and outer regions (both are at 670°C). But, the area of material comprising the inner region is less than that of the outer region, therefore, more heat enters from the outer region. This implies that a cold region, which is closer to the inner wall than the outer wall, will develop. So a temperature curve like the one in figure 2.6f will develop in the radial direction and the cold region will be closer to the inner ring (see figure 2.8).

Top view -of cylinder

Cold region

The extent of anti-symmetry is influenced by the dimensions of the coil and the magnitudes of Kr and K=. If K= is much larger than Kr, heat flow in the z-direction

dominates and the anti-symmetry of the r-direction becomes less significarit. On the other hand, if K,. is comparable or even larger than K=, the effect of anti-symmetry becomes

more pronounced. Note also that the cold region is closer to the inner ring for a thick coil than for a thin one.

Since the temperature distribution is symmetrical in the z-direction, with the minimum value at the mid-point of the cylinder, it is safe to say that the cold region develops as a ring that is situated in the center of the coil (see figures 2.8 and 2.9).

,--

...

~~~

. .

.

. 1 1 1 1 Cold nng-c; 1 1 1 1

r--:...~"·~···t·1

I

:~~."'··I

1 1 1 1 1 ... -

--

-.I, 1

...

)

,~::;

Figure 2.9: An illustration of the position of the cold ring.

The cold ring is at a position of less than a half of the distance from the inner to the outer wall, as shown in figure 2.8. The position of the cold ring also depends on how the cylinder is heated. In a typical batch annealing furnace, the outer wall of a cylindrical coil tends to be hotter than the- inner wall during heating and the converse is true during cooling. This would affect the position of the cold ring, particularly in the radial direction.

2.7.1 Temperature ramp

During an annealing cycle, the edges of the coil are directly exposed to a temperature ramp. Therefore, the edges remain at the same temperature as the ramp throughout the heating process. A temperature ramp was applied to the edges of a solid cylinder (this was implemented by the MATLAB program finramp.m, in appendix B.7). Figure 2.9 shows the time evolution of the hottest region (the edges) together with the cold region. Note, once again, that the temperature profile of the hottest region is identical to that of the applied temperature ramp.

Hottest region Coldes t region

5 10 15 20 25 30 35 40 45 50

"Time (hours)

Figure: 2.9: A comparison of the outer region and coldest region during a temperature cycle.

In the next section, an analytical solution will be obtained for a solid cylinder.

2.8 Analytical Solution for Solid Cylinder

This section seeks to obtain an analytical solution to the heat equation for axial symmetry by means of separation of variables. The heat equation for axial symmetry is given by (2.4),

This equation will be solved on the rectangular domain B={(r,z)lr E [rin,rouJ, z E [O,Hj}.

Let OB be the boundary of this domain, with rinandroutbeing the inner and outer radii of the coil, respectively. The coil height His in the z-direction.

Equation (2.4) will be solved subject to the boundary conditions

Tït.r.z) =0, t ~0, (r,z) E OB (2.22)

and the initial conditions

{o,(r,z) E

es

T(O,r,z) =

-To,{r,z) EB-aB (2.23)

This condition assumes a shifting in the temperature so that room temperature is at -To

and the externally applied temperature (670°C in this case) is 0.

In order to find a solution by means of separation of variables, T will be expressed as the following product: T(O,r,z)

=

ToR(r)Z(z) where

{a,

R(r)

=

-1, 0, r

=

r.In (2.24)

{O, Z(z)

=

-1, 0,

z=o

zE(O,H). z==H (2.25)

See figure 2.10, for a graphical illustration of the functions R(r) and Z(z).

...-c;..b __ ~••.•" 0.2 0

t

-0.2

l

-0.4

i

-0.6 -0.8 ~ -1 0 0.2 0.4 0.6 0.8 1.2 (a) 0.2 0 ~ -0.2

!

-0.4 "'-_.~~

l·~~~

.-:...,.._,. I-~ -0.8 -1 -0.2 0 0.2 0.4 0.6 0.8 1.2 z (b)

0.2 ... . ... . . . . .,,- ,;0-0.8 0.6 0.4 z (c)

Figure 2.10: (a) Initial condition (2.24). (b) Initial condition (2.25). (c) Initial condition for T(r,z,t).

Since a solution by separation of variables is desired, the solution can expressed as a product of three functions of single variables,

T(t,r,z)=U(t)R(r)Z(z) . (2.26)

~ow~ by. substituting (2.26) into (2.4), the_follo.wi~g expression is obtained:

1

Kr -U(t)R'(r)Z(z)

+

KrU(t)R"(r)Z(Z)

+

K=U(t)R(r)Z"(z)

=

U'(t)R(r)Z(z).

r

When the above expression is divided by U(t)R(r)Z(z), the following is obtained:

1R'(r) R"(r) Z"(z) U'(t)

K ---+K --+K --=--.

r r R(r) r R(r) = Z(z) U(t)

U'(t)

U(t) =-Jl (2.28)

Let-Jlbe the separation constant. Itthus follows that

and

K .!._ R'(r)

+

K R"(r) = _ _K Z"(z) =-A?

,. r R(r) ,. R(r) Jl. z Z(z) , (2.29)

The left hand side is function of r, while the right hand depends on z. This implies that both sides equal a constant value, say _,1,2, which is known as the separation constant. Consequently, equation (2.29) yields the following ordinary differential equations,

(2.30) where (2.31 ) K_ and K,.R'(r)

+

K,.r2 R"(r)

=

_r2 R(r)A2 . (2.32) Multiplication by r2R(r) reduces (2.29) to r2A2 r2R"(r)

+

rR'(r) +--R(r)

=

O. Kr

(2.33)

which is known as Bessel 's equation with n

=

0 [9].Itshould be noted from (2.31), that

(2.24)

Time Solution

The solution of (2.28) is

(2.34)

Z-Solution

The general solution of (2.30) is

Z(z)

=

Acos(wz)

+

Bsin(wz), (2.35)

where A and B are constants. The boundary condition Z(O) 0 requires that A=O, therefore,

Z(z)

=

Bsin(wz). (2.36)

The boundary condition Z(H)

=

0 requires that sin(wz)

=

O. This implies that wH =jn;

wherej is an integer. Therefore,

jtr

ar, =-,

J

H

is a set of valid values for to, and any linear combination of these modes will still be a solution of (2.30). As a result, the most general solution is:

CJ)

2(z)

=

IBj sin(())jz),

j=1

(2.37)

At this point, the coefficient Bj needs to be determined from the initial conditions (2.25). The initial conditions are enforced by expressing the initial condition in2(z) in the form of a Fourier Series expansion [10], yielding

4 C .

Bj

= ----:-'

lOr]

=

1,3,5, ... lf}

and

Bj = 0, for j=2, 4,6, ...

A superposition of the initial conditions and the approximating Fourier senes, for

FOURIER series 1.5 o (\" ,,-"Li V vV

-+--

I---

--

_-

--

,....--

-

--_

I

.r.

+-

--

---

I---L1\ ft

-

-

_Vvv

v~ 0.5 -0.5 -1 = -1.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 z

Figure 2.12: A Fourier series approximation of initial condition (2.25).

R-Solution

The most general solution of (2.33) is

(2.38)

where A and b are constants. According to boundary condition (2.24),

(2.39)

for all time. For convenience sake, the range of R(r) will be transformed from [Yin,rout] to

[0;1], where (J"= rin/rout, At this point, a set of a's and b's that satisfy (2.39) needs to be

found. The MATLAB routine bessolro.m performs this function (see appendix B.2).

Since it is desirable to have the quantity R(r)Z(z) equal to -1 in the interior, the initial condition (2.25) will be modified as follows:

{O,

R(r)

=

1, 0,

r =r.In

(2.40)

The functions that satisfy (2.39) are called eigenfuctions. The eigenfunctions of (2.39) are

(2.41)

From Sturm-Liouville theory [11], the eigenfunctions must be orthogonal with respect to a weight function w(r),

where} :;ék.

In a similar way as before (with Z(z), the function R(r) can be expressed in the form of a Bessel series. The Bessel series is composed of a linear combination of the eigenfunctions of (2.39). As a first step in building a Bessel series, the orthogonality of these eigenfunctions need to be tested.

Test for Orthogonality

Assume that for CXj and ea; }:;ék. In addition to this, let w(r)=r. Therefore, using an integration property of bessel functions [9],

Lrpj(ajr)Pk(akr)dr

=

2 r 2 [ajPk(akr)p'j (ajr)-akP;(ajr)p'k (akr)t

a, -aj

=

2 r 2 [ajPk(akr)p'j (ajr)-akPj(ajr)p'k (akr)

Bessel Series

A function R(r) defined on an interval [0'",I] can be approximated by the bessel series

Cf) R(r)

=

'LAkPk(akr) k=l (2.42) where (2.43) and

such thatpk(a/(Y) = pk(aJJ = 0, for k=l, 2, 3, ...

In this particular case, R(r) is given by [9], therefore,

is an expression for the bessel series coefficients. Figure 2.12 shows the initial condition for R(r) with a bessel series approximation superimposed.

BESSEL series 0.5 --- ---.--+--I---+---t--- --- ---o --- ----+--+--I----!--- --- ---0.5 --- --+----+---11---+--- -- t----1

---~ ""

r-; V V -1.5 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.1

Figure 2.13: A Bessel series approximation for condition (2.24).

Specific Solution

The final solution satisfying the initial and boundary conditions is therefore

2.9 Analytical Results

The program bessoI321.m, in appendix B.9, implements equation (2.44). The analytical results for a 21 x 21 grid are shown below.

Time=0 hours Time=2 hours

100 o 1 600 [)500 '<;;;'400 s ~ 300 ~ 200 t--(a) (b) Time = 4 hours 100 o 1 o 0 600 [)500 '<;;;'400 . s ~ 300 ~ 200 t--(c)

Note that at time =0 hours, there are some oscillations in the base of the diagram. This is due to the fact that only a finite number of terms in both Fourier and Bessel series were used. As time increased, the temperature profile behaved in a similar manner to the numerical solution obtained in the previous section.

2.10 Comparison of Analytical and Numerical Solutions

In section 2.4, the finite difference method was used to solve the heat equation (2.4), for the case of axial symmetry, over a discretized solid cylinder with a concentric hole. The same equation was also solved analytically in the previous section. If the numerical solution is correct, it must compare favorably with the analytical solution.

As the number of finite difference grid points for the cylinder increase (for a fixed dimension) and the time step dt decreases, equation (2.9) approximates equation (2.4) more accurately. Therefore, if the solution obtained in section 2.7 (for constant heating on the boundaries) is correct then it should approach the analytical solution, for a given time, as the number of grid points increase and as the time-step dt tends to zero. An increase in the number of grid points can also be viewed as a decrease in the step lengths

hy and h= (as the number of grid points increase hy and h- must tend to zero).

Tables 1, 2 and 3 show the temperature and percentage difference between the analytical and numerical solutions at the following times: 0.4 hours, 1.2 hours and 2 hours respectively. All three tables show the change between the two methods as the grid size increases and as the time-step decreases. The stability condition causes the time-step to decrease as the number of grid points increase.

The percentage difference was obtained by comparing the temperature difference to the temperature of the hottest region in the cylinder (in this case the boundaries experienced the highest temperature at 670°C).

Table 2.1: time

=

0.4 hours

Grid Size dt(m) Maximum Temperature Maximum Percentage

Difference (0C) Difference (%) 5x5 0.4 188.26 28.10 10x10 0.08 82.048 12.25 15x15 0.0363636 35.68 5.33 20x20 0.02 22.87 3.414 25x25 0.0125 13.75 2.052 30x30 0.0088888 9.669 1.443 40x40 0.005 5.0332 0.7512 50x50 0.0030769 2.907 0.4339

Table 2.2: time

=

1.2 hours

Grid Size dt(m) Maximum Temperature Maximum Percentage

Difference (0C) Difference (%) 5x5 0.4 162.8 24.30 10x10 0.08 39.39 5.88 15x15 0.0363636 18.35 2.739 20x20 0.02 10.07 1.503 25x25 0.0125 6.270 0.9358 30x30 0.0088888 4.345 0.6485 40x40 0.005 2.332 0.3480 50x50 0.0030769 1.361 0.2031

Table 2.3: time

=

2 hours

Grid Size dt (m) Maximum Temperature Maximum Percentage

.. Difference (DC) Difference (%) 5x5 0.4 105.1 15.69 . 10x10 0.08 24.12 3.6 15x15 0.0363636 11.23 1.675 20x20 0.02 6.212 0.9271 25x25 0.0125 3.868 0.5772 30x30 0.0088888 2.704 0.4035 40x40 0.005 1.479 0.2208 50x50 0.0030769 0.8833 0.1318 200 180

R

160 140

I

120 ~ ₁₀₀ Ijl

i

80 m 60 u rp ~ 40 M 20

"

Time

=

0.4 hours

•

Time

=

1.2 hours 0 Time

=

2 hours

"

•

"

0 "

•

0 _t;

_"

'é ~ j, o 5 10 15 20 25 30 35 40 45 50 Grid size

Figure 2.15: A graph of maximum change in temperature versus grid size.

A graph of the maximum change in temperature versus grid size, for tables 1,2 and 3 (for 0.4 hours, 1.2 hours and 2 hours), is shown in figure 2.15. At 0.4 hours, the maximum difference in temperature decreases rapidly until it becomes negligible at grid size of

a result of the fact that the temperature distribution moves closer to steady state as the time increases. At steady-state, the temperature of the entire cylinder reaches 670°C.

Chapter 3

Temperature

model for a Cylindrical Coil

3.1 Introduction

In the previous chapter, heat transfer in a solid cylinder was treated as an approximation

to a cylindrical coil. This was based on the assumption that the temperature dependence had no B-dependence.

In this chapter, however, the effect of heat transfer along the coil spirals will be

investigated. Once the model has been constructed, the first task will be to test whether the B-dependence of temperature is appreciable or not. A comparison will also be made between the time evolution of the hottest and coldest regions in the coil. This will be done under temperature conditions similar to the batch annealing process.

3.2 Finite Difference Model

This section deals with the discretization of both the cylindrical coil and the heat equation for theta-dependence. The grid structure of the coils is introduced prior to the discussion regarding the heat equation.

3.2.1 Numerical Grid

The cylindrical coil was discretized into a series of coil slices stacked in the z-direction. These coil slices were divided into an equal number of sectors (see figure 3.1).

z

Figure 3.1: The cylindrical coil is discretized into sectors and slices.

Each individual slice is described by the following equation in polar coordinates:

r=mB +c, (3.1)

where m and c are constants. The grid, in this case, was constructed using cartesian coordinates. This means that the x and y-values of each point on the coil grid are given by:

x = rcosB (3.2)

y= rsinB (3.3)

Note that the grid points of the coil are on the edges of each sector, such that if there are 16 sectors, then each coil spiral will have 16 elements. Thus, the entire grid can be viewed as a stack of matrices that specify each coil slice, as shown in figure 3.2.

z

X Points YPoints

Y

x

~

\COil3

J

Figure 3.2: The X and Y matrix stacks specify the stack of coils.

Note that the matrices are arranged in the manner shown in figure 3.3. Each row of theX

and Y matrices constitutes a spiral. In addition to that, each column of the X and Y matrices constitutes a sector.

Both X and Y matrices are structured in the following way: ~

-•

•

•

L_

•

_""

....

•

•

•

•

...

•

•

•

•

•

•

•

·

...

...1'

_l'

,\:'

Sector 1 _{Sector 2 • • • • Sector N}

Spirall Spiral2

Spiral N

Figure 3.3: The X and Y matrix structure.

Therefore each spiral, with all its sections, is packed into a separate row in the matrix structure and each matrix in a stack describes a coil.

3.2.2 Discretized Heat Equation

From (2.2), the heat equation in cylindrical coordinates is:

1

et

a2T 1 a2T a2T

et

Kr--+K"--2 +K_-2 --2 +K_--2 =-.

r ar ar _ r a () _ az at

Since the thermal conductivity is the same along the spirals as it is in the axial direction, it follows that:

(3.4)

_!_(

Til i+l,j,k - Til i-l,j,k

J

+ (Til i+l,j.k - 2TIl i,j,k +Til i-l,j,k

J

_1 (Til i.j+l,k - 2TIl i,j,k +Til i,j-l,k

J

Kr Kr 2

+

Ko 2 2

'i 2hr hr ri he

After some re-arranging, the result is:

TIl+I" - Til, , _!_[Tlli+l,j'k -Tlli_I,j,k

J

I,j,k - I,j,k

+

rKr

r.

2h,. (Tlli+l,j'k - 2T'';,j,k

+

Tlli_I,j,k

J

+rx , 2 " h; (3.5)

Since each row in the temperature matrix represents a spiral (figure 3.3), a central difference (appendix A) at a point with an index i =1 will use a point ati =N(or i=-1 is

equivalent to i

=

N), where N is the number of points in a row. In a similar manner, i

=

N+ 1 is equivalent toi

=

1. This is true for all values of indices j and k.

The stability conditions were obtained by substituting

Til _;:11 i(iIlJ,+ jllJn+kcv_)

i,j - '=' ev' (3.6)

-8A

(3.7)

T ~ (B2

+

16A2)' where A ~ ( - :'; - : " - : "

J

(3.8)

r () = and B=~.

(3.9)

hrrmin

3.3 Initial Conditions

Initially, the coil is at room temperature. Then, a heating cycle (figure 1.2) is applied to all the boundary surfaces of the coil. In terms of the grid arrangement, it means that the first and last spiral of each coil in the stack must also be heated directly (figure 3.4).

Top view:

first spiral

last spiral

Figure 3.4: Top view of a cylindrical coil, with the first and last spirals highlighted.

In addition to this, the first and last coil slices in the stack must also be heated directly (these coils form the top and bottom of the cylindrical coil). In terms of the matrix arrangement of the grid, one can use a stack of temperature matrices and correlate each grid point to its corresponding temperature (as shown in figure 3.5).

x

Sector Sector ••••

f

Sector Coil3

Figure 3.5: Structure of temperature matrix.

The initial conditions can now be implemented by allowing both the first and last rows of each temperature matrix to have the same value as the heating cycle at every instant. Once again, in addition to this, the first and last matrices in the stack must also have the same value as the heating cycle at every instant. The remaining matrix elements, corresponding to the interior of the cylindrical coil, are all at room temperature initially.

The heat equation was solved subject to the previously mentioned initial conditions. A computer program (called thermd.m) that solves equation (3.5) subject to the initial conditions (in the form of temperature matrices) was written in MATLAB (appendix C). This program also displays the calculated temperatures on the corresponding grid points. Note that the room temperature values for Kr, K= and p were used throughout the

calculation.

3.4 Numerical Results

The temperature distribution of the first sector of a coil with 16 sectors is shown in figure 3.6. In this particular case, the edges of the coil were heated at a constant temperature of

Time = 2 hours

spiral having a width of O.lm, which is 100 times the typical thickness (this was done for display purposes). The coil height is 1m.

lime =0 hours 700 600 600 ~ 400 ~ 500

i

400 200

i

~ 300 ~ 200₁₀₀ 1.5 0 1 0 ₀ (a) lime= 4hours 1.5 (c) 1.5 o 0 (b) (d)

Time=3 hours 700 540

I

~

li

~

_l\

~ _r-,

_V

~

v-~

-,

:~

.'\

r-~ ...

-/

680 680 640 ~ 620

i

600 ~ 560 560 520 500 o 0.5 1,5 (e)

Figure 3.6: (a) to (c) Temperature profiles for the first sector a cylindrical coil with dimensions as stated in the beginning of section 3.4. (d) A cross sectional profile of the first sector, in the z-direction. (e) A cross sectional profile of the first sector, in the

r-direction.

It can be seen, from figure 3.6, that the temperature profile follows a similar pattern as the solution of a solid cylinder with a concentric hole. There is symmetry in the

z-direction and anti-symmetry in the r-z-direction. The shape of the curve in the r-z-direction is not completely identical to the two-dimensional case, this difference can be attributed to heat flow along the spirals (remember that the spirals in this case are 100 times thicker than usual).

Another display of the coil during the heating process is shown in figure 3.7. A coil slice from the center is shown together with a side view of the first sector, at a time of2 hours.

lime= 2 hours -1 -0.5 0 0.5 1.5 (a) 100 600 Outside temperature= 670°C 1 500 0.81-+-+---+--+-+-+---4-+-+-4 400 0.4I-I-+-+-I-+-+-II-+-+---l 300 0.2 I-I-+-+--I-+-+-II-+-+---l 200 O~~~~~~~~~ 0.4 0.6 0.8 1.2 (b)

o

Figure 3.7: (a) A diagram of the temperature distribution in the center slice of a coil, with dimensions as stated at the beginning of section 3.4. (b) The temperature distribution of the first sector, in the z-direction. These distributions were taken after 2 hours.

Note that a cold region also develops in this case. From figure 3.6, one can see where it is situated. It appears to be a ring situated in the center of the cylindrical coil. When the coil is viewed from the top, the radius of the cold region is positioned at around half the distance between the inner and outer radius.

The time evolution of the coil edge and cold region, when the coil is subjected to a temperature ramp (see figure 1.2), is shown below (figure 3.8). These results were also generated from the program thermd.m in appendix C.

700 600 500 ~ 400

i

300 ~ 200

100 Hottest_Coldesl region in coil_{region in coil}

O~~--~--~~--~~~~--'_~--~

o 5 10 15 20 25 30 35 40 45 50

lime (hours)

Figure 3.8: The time evolution of the hottest region and coldest region during batch annealing.

3.5 Theta dependence

A possible way to examine whether the heat flow along the O-direction is appreciable or not is to compare the temperature distributions of the sectors. Note that the top and bottom slices have equal temperatures (all the sectors are at the same temperature) because of the effect of direct heating. The coil slices just below these two surfaces also experience this effect, but to a lesser degree. As a result, the middle coil slice should have the largest temperature differences between sectors, as a result of heat flow along the spirals. For this reason, the center coil slice will be used to test for O-dependence of temperature.

The extent of O-dependence can be measured by monitoring the maximum difference between sectors that are situated opposite to each other. For axial symmetry, the maximum difference would be zero.

The dimensions for all three cases are as follows:

amount of computation time necessary to complete the calculation. This problem was surmounted by performing a c9-dependence test for three simpler cases and using the resulting trend to make a projection for a realistic coil.

Coil height: 1m

Inner radius: 0.6m

Outer radius: l.4m

Number of sectors: 8

Number of slices: 10

The width wand number of spirals Nwere varied as follows:

Case 1: w

=

O.lm N=8 Case 2: w=0.05m N= 16 Case 3: w= 0.025m N=32

The number of spirals double as the spiral width is halved. The maximum temperature difference between sectors 1 and 4 are shown in tables 3.1, 3.2 and 3.3 for the three cases. This difference is shown with the corresponding number of iterations. The number of elapsed iterations are proportional to the elapsed time, with lIdt being the proportionality

Table 3.1: Case 1,dt =O.1993m, w =O.lm with 8 coils turns.

Number of Maximum Difference

Iterations

eq

1 0 2 6.917 3 12.97 4 18.27 5 22.91 10 33.03 20 29.37 30 20.82 40 13.99 50 9.233 60 6.051 70 3.949 80 2.570 90 1.669 100 l.081

Table 3.2: Case 2,dt = 0.1614m, w=O.05m with 16 coil turns.

Number of Iterations Maximum Difference

CC) 1 0 2 6.388 3 10.48 4 13.17 5 14.98 10 16.51 20 l3.06 30 10.53 40 7.910 50 5.776 60 4.167 70 2.985 80 2.147 90 1.546 100 1.110

Table 3.3: Case 3, dt

=

0.0917 hours, w

=

0.025m with 32 coil turns.

Number of Iterations Maximum Difference

(OC) 1 0 2 3.889 3 4.715 4 5.127 5 5.369 10 6.368 20 6.764 30 6.294 40 5.598 50 4.809 60 4.123 70 3.517 80 2.971 90 2.493 100 2.083

35 ~ 10

T

M 5 x _{x Case 1}

...

Case 2 x 0 Case 3 x x x

...

...

_x

..

...

...

...

x

...

•

0 0 0 x

?

0

...

0

..

~

_i

i

~

-30

o

o

10 20 30 40 50 60 70 80 90 100 Number of iterations

Figure 3.9: A graph of maximum change in temperature versus number of iterations, for cases 1, 2and 3.

According to figure 3.9, at the instant when the coil is heated initially, the maximum difference between sectors 1 and 4 is zero. This is because heat flow along the spirals has not occurred yet. But, as the number of iterations increase, the heat flow along the spirals increases the maximum difference between sectors 1 and 4 until the maximum difference curve reaches its highest value (at about 10 iterations). Beyond this point, the effect of heat transfer in the radial and z-direction increases, as the temperature moves closer to steady state. This causes the difference between sectors to decrease gradually until it eventually reaches zero, when the thermal gradients have leveled out.

Figure 3.9 also shows that the highest point in the maximum difference curve decreases as number of spirals increase and as the spiral width decreases. This observation can be used to determine a trend. Itis useful to determine how the highest point in the maximum difference curve changes as the division number increases. The division number is

division number is 1, for case 2, the division number is 2 and for case 3, the division number is 3.

The maximum values for the three cases were plotted against division number, and a curve was fitted through the relevant points, see figure 3.10.

35 P

I~

Fitted CUM

II

Numerical Experiment Data 30 ~--...:.,,--:::; ~ <= ;; 25

i

20 ~ Ijl I>

i

15 m ₁₀ u

i

M ₅ 0

\

0 10 20 30 40 50 60 70 80 90 100 DilAsion Number

Figure 3.10: A graph of maximum change in temperature versus division number.

The equation of the fitted curve is,

Y=Aexp(-Ax), (3.10)

where A

=

7S.4SoC and A

=

0.7929m-1. Itis expected that A and Awill depend on Kr, K=

maximal temperature difference at a division number of 100 is 2.77xlO-33 °C. This shows that for a realistic coil, the B-dependence is negligible and that the coil can be treated as a solid cylinder when modeling heat transfer.

",

#.~ -t, _" •

_-,·:·'.f<:~{>·~"

Chapter 4

ELASTICITY

THEORY

4.1 Introduction

The theory of elasticity basically deals with the study of stress, strain and deformation of materials in the elastic limit. In this field, the behaviour of the material as a whole 'is of importance (individual molecules are not considered). This is the continuum concept of matter. Stress is defined as the force experienced by a plane of a body per unit area

(stress=force/area). Strain, on the other hand, is the resulting distortion of the body divided by the original dimension tstrain=distortion/original dimension).

The main objective of the following sections is to derive some equations that will be instrumental in determining the thermal stress in a material that has a given temperature distribution. Pertinent concepts of elasticity will first be introduced in cartesian coordinates and thereafter in polar coordinates.

Brief discussions regarding stress, strain and the notion of a tensor have been included as precursors to the derivations. The equilibrium conditions for boundary traction, in

developed from section 4.6 to 4.7. These equations of strain are later used to arrive at the compatibility equations (section 4.8).

Hooke's law for isotropic media is introduced in section 4.9, and it provides the necessary background for the discussion on thermal stress that follows afterwards. The equilibrium conditions and the equations of strain are derived for polar coordinates in the subsequent sections.

4.2 Stress and Strain

Figure 4.1, which shows a bar of initial length La and cross-sectional area A, illustrates basic one-dimensional stress and strain. When an axial force F acts on area A, as shown in figure 4.1 b, the resulting stress and strain is given by the equations in figure 4.1b.

Area: A (a) F (force) F (force) Stress =FIA Strain

=

MILo (b)

Figure 4.1: (a) Initial bar. (b) Bar during the application of a force, with a force-induced lengthening of M.

displacement occurs in a direction that is perpendicular to the surface on which the stress is applied.

A shearing stress occurs when a force is applied in a parallel direction to a particular surface on a body. In a similar way to normal strain, shearing strain is also caused by a shearing stress. But in this case, the resulting displacement occurs in a direction that is parallel to the face upon which the stress has been applied (see figure 4.2).

(a)

F

(b) (c)

(d)

Figure 4.2: (a) Initial block element. (b) Deformation resulting from a horizontal normal stress. (c) Deformation resulting from a vertical normal stress. (c) Deformation resulting from a shearing stress.

It should be noted, however, that in realistic situations, materials are subjected to a mixture of normal and shearing stresses.

4.3 The Stress Tensor

The notion of a tensor will be used extensively throughout the treatment of elasticity in this text. A basic understanding of tensors can be gained by examining some of the stress tensor properties.

When dealing with a two or three-dimensional body, it is useful to express stresses and strains in tensor notation. The stress tensor 0i gives all the stresses that act on each point throughout the material. It is therefore a function of position. The stress tensor for a three-dimensional material, with axes pointing in theXI, X2 andX3 directions is,

0'13

1

0' 23 •

0'33

(4.1)

Figure 4.3 shows the stress tensor components on a material element. In figure 4.3, the xt face is the face that is perpendicular to the Xl axis. This convention also holds for the X2

andX3 faces.

The stress 0'11 represents the stress on thex, face in the xi direction. The stresses 0'22 and

0'33 can be explained in a similar manner. These are all normal stresses. Note, however,

that 0'12 is the shearing stress on the xt face in the X2 direction. Similarly, 0'23 and 0'13

represent the stress on theX2 face in the X3 direction and the stress on the Xl face in the X3

direction, respectively. Since each element in a continuum should be in equilibrium, it follows that:

and (4.2)

If (4.2) does not hold, the element will rotate.

4.4 Equilibrium Conditions for Boundary Traction

This section will seek to obtain some equilibrium conditions for a material element that experiences boundary traction [12]. The relationship between an external applied stress tn

and the internal stress tensor Oij is investigated here.

Consider a tetrahedron element of some continuum with one face on the surface, as shown in figure 4.4. Let ABC be the face on the surface of the continuum.

Area: n7dS

-

- -

--(b)

(c)

Figure 4.4: (a) Continuum with a tetrahedron element. (b) An enlargement of the element with axes drawn along the perpendicular edges of the tetrahedron. (c) An illustration of the external forces on the element.

Suppose that there are body forces acting on this element and b is the resultant body force per unit volume. Suppose also, that some external stress t, (also referred to as the traction)

is applied on the various faces of the element (tn on ABC, tI on BDC, t2 on ADC and t3 on

ADB), as shown in Fig 4.4c.

The traction on face BDC (which can also be considered as the xi face because it is perpendicular to thexi axis) of the element (tl) is given by:

tI = alli

+

anj

+

aJ3k, (4.3a)

or alternatively: tI is equal to the stress on face x, in the xi direction plus the stress on the

xi face in the X2 direction plus the stress on the x, face in the X3 direction. By using similar

reasoning, the following relations can be written:

and

As shown in figure 4.4b, dS is the area of face ABC. It thus follows that the area of face

xi is the projected area nldS. Similar results also hold for faces X2 and X3 which have areas

n2dS and n3dS respectively. In order for the element to be in equilibrium, the forces due to

tI, t2, t3, tn and b must add up to zero, that is:

(4.4)

In the limit as dS tends to 0, dV tends to

°

much faster, since dS is a square of a side while

dV is a cube of a side. Therefore, equation (4.4) can be approximated by

Note that the dS cancels out. If the values of tI, t2, t3 (equations 4.3a to 4.3c) are substituted into (4.5), the following expression is obtained:

=2:n, where 0'131

[n1-0'23 and n= n2 . 0'33 », (4.6) Therefore: tn

=L:n.

(4.7) This expression relates the external applied stress tn to the internal stress tensor

I

and it will be incorporated into the equilibrium conditions of the material as a whole in the following section.

4.5 Derivation of Equilibrium Conditions

In order for each element to remain in equilibrium, certain relationships between the various stress components and the body force must exist. These equations are known as the equilibrium conditions [12]. These conditions must be satisfied at all points throughout the interior of a body, regardless of the applied stress.

Consider an arbitrary volume Vof a continuum that is subjected to both surface and body forces (see figure 4.5).

n

Figure 4.5: A continuum of volume V that is subjected to surface and body forces.

In order for this body to be in equilibrium, both the surface and body forces must add up to zero at each point of the continuum. Itthus follows that,

Floral

=

ftn .

ndS

+

fpbdV

=

0,

S V

where n is a unit vector normal to the surface. According to Gauss's theorem:

fv

·ndS

=

fV. VdV,

s v

therefore

fv.tndV

+

fpbdV=O.

v v

Since this expression is valid for any shape, it follows that:

v .

(L)

+

pb

=

o.

In tensor notation (see appendix D), the above mentioned result can be expressed as:

(Jji,j

+

pbi = O. (4.9)

This is the set of equilibrium conditions for the interior points of a body. The expanded form is shown below:

4.6 Lagrangian and Eulerian Descriptions of Deformation

When a continuous medium is deformed, the particles that comprise the medium move along particular paths. The motion of these particles can be expressed by equations of the form

Xi=xi(X], X2, X3, t) =Xi (X,t) or x = x(X,t) (4.10)

where X]X2X3 and X]X2X3 are two sets of coordinate systems corresponding to the initial

u

X3

t=O

Figure 4.6: At t=O, the material is in an undeformed state. The material is deformed at t=t. There are two sets of coordinate systems corresponding to the initial and final states. Vector X points from the origin of the XlX2X3 axes to Po and vector x points from the

origin of theXlX2X3axes toP.

The particle that occupied the point (X], X2, X3) at t

=

0 now occupies (Xl, X2, X3) at t

=

t

(after deformation has occurred). So (4.10) is a mapping of the initial coordinates att

=

0 onto the coordinate system at t

=

t. A one-to-one mapping with continuous partial

derivatives is assumed. The deformation, as given by (4.10), is known as the Lagrangian description.

Ifthe deformation is given by equations of the form

X; =X;(Xl, X2, X3,t) or X

=

X(x,t), (4.11)

then it is known as the Eulerian description. Once again, a one-to-one mapping with continuous partial derivatives is assumed. As a result of this assumption, it follows that (4.10) and (4.11) are unique inverses of each other. This also implies that the Jacobian