CC BY
23
UDC 621.3: 537.8
doi: 10.20998/2074-272X.2016.4.03
M.G. Pantelyat, I. Dolezel
FINITE ELEMENT TECHNIQUE FOR SOLUTION OF THERMO-CONTACT PROBLEMS AND ITS APPLICATION IN NUMERICAL ANALYSIS OF DEVICES WORKING WITH INDUCTION HEATING
Purpose. To develop an effective approach for the numerical solution of transient thermo-contact problems and present a typical example of its utilization regarding devices working on the principle of thermoelasticity produced by induction heating and specific technological processes intended for assembly and disassembly of systems containing shrink fits. Methodology. A finite element technique for solution of 2D multiphysics (electromagnetic, thermal and structural) problems is developed, taking into account temperature dependences of material properties and continuous variations of the contact surfaces. Modeling of the contact interaction between two parts is based on the concept of a special contact finite element having no thickness. The functional for the temperature problem is supplemented with components corresponding to the thermal conductivity of this contact layer. The heat generated due to mutual sliding of both parts can also be taken into account, but the heat capacity (specific heat) of the contact layer is neglected. Using a special 1D 4-node finite elements a system of equations for the description of the thermo-contact problem is obtained. Originality. Relatively simple analytical formulae for calculation of the contact thermal resistances occurring in specific parts of electrical machines are known. The paper offers an alternative approach for the numerical solution of transient thermo-contact problems based on the concept of a special 1D contact finite element having no thickness. Results. The presented technique is applied for the computer simulation of assembly and disassembly of a shrink fit using induction heating. Conclusions regarding the choice of technological modes are made. Comparative computations for drills made from hard alloy and alloyed tool steel are carried out. References 8, figures 6.
Key words: induction heating, thermoelasticity, multiphysics problems, numerical analysis, finite element method.
Описана методика решения методом конечных элементов мультифизических (электромагнитных, тепловых и механических) задач с учетом зависимостей свойств материалов от температуры и изменения контактных поверхностей. Предложенный подход использован для численного анализа устройств, функционирующих на базе явления термоупругости в процессе индукционного нагрева и оригинальных технологических процессов, предназначенных для сборки и разборки конструкций с напряженными посадками. Применение разработанной методики проиллюстрировано на конкретном примере. Приведен анализ полученных результатов. Библ. 8, рис. 6.
Ключевые слова: индукционный нагрев, термоупругость, мультифизические задачи, численный анализ, метод конечных элементов.
Introduction. Investigation of behavior of numerous electrical devices in different operation regimes often requires considering relevant multiphysics phenomena of electromagnetic, thermal and structural origins. In many cases, multiphysics analysis of such devices must also include the influence of thin insulation layers and contact thermal resistances for obtaining more realistic results. Thin layers occur, for example, in the form of various bandages [1], contact resistances (of thermal origin) play a significant role in different shrink fits [2] and devices for a number of industrial purposes [3].
Relatively simple analytical formulae for calculation of the contact thermal resistances occurring in specific parts of electrical machines are presented in [1]. These can directly be used as material parameters during the solution of the temperature problem.
The paper offers an alternative approach for the numerical solution of transient thermo-contact problems. Currently the proposed technique is used for solving such problems in 2D Cartesian and axisymmetric systems. A similar technique for solving 3D problems is being developed nowadays. The aim of the paper is to describe this technique in detail and present a typical example of its utilization.
Formulation of technical problem. Many modern industrial technologies are based on the principle of induction heating. The paper, however, will focus on its application in production of shrink fits for specific purposes. Typical is, for example, setting the disks on shafts, fixing high-speed machine tools, or connecting pipes by
fixing sleeves, which always represent connecting of two metal parts with an interference whose value is decisive for the transferable mechanical force or torque.
The process of manufacturing shrink fits starts with induction heating of one of the parts, which leads to increase of its dimensions. Then it is connected with another part and the whole system is cooled. The shrink fit is obtained after cooling. A typical example of fixing a drilling tool in the chuck is depicted in Fig. 1.
I
shank
inductor
Fig. 1. Arrangement of drilling machine: I - heating of chuck and inserting of shank into it II - cooling of system and fixing of shank in chuck
© M.G. Pantelyat, I. Dolezel
From the physical viewpoint, the process represents a strongly nonlinear and nonstationary multiphysics problem characterized by mutual interaction of magnetic field, temperature field and field of thermoelastic displacements (and corresponding strains and stresses). Another aspect to be involved in the model is the contact problem that plays an important role when quantifying transfer of heat between both connected parts.
The numerical solution of the model should provide a sufficiently accurate idea about the whole process and mainly the steady-state parameters of the system under investigation.
Continuous mathematical model. The mathematical model of the process consists of three partial differential equations (PDEs) describing three involved physical fields and relations describing the contact problem.
Magnetic field. Distribution of magnetic field in the system may be described using several formulations [4]. When using, for example, the magnetic vector potential A, the above field obeys the equation
Curl(—Curl A) + = Jext' ß dt
(1)
wj =y
dA
dt
div(2 grad T ) = pCp—— wj,
dT dt
(4)
3 dT
-= a
dn
gen
(T - Text)
(5)
where n denotes the outward normal, agen is a coefficient generally quantifying both convection and radiation and Text stands for the temperature of sufficiently distant environment.
Thermoelastic problem. Mechanical status of structural parts of the considered devices is described by the following system of three-dimensional tensor equations [6]
Oji + fj =
1/
Sj =T(ui, j + uj,i )
2
(6)
E , ^ s? \ a, =-S +—-öje\
lj 1 + vlj
1 - 2v
where ^ is the magnetic permeability, y stands for the electric conductivity and Jext represents the current density applied to the inductor. The boundary condition along a sufficiently distant artificial boundary is of the Dirichlet type and reads A = 0.
The eddy currents produced by time-variable magnetic field in electrically conductive bodies (whose density is given by the second term on the left-hand side in (1)) give rise to the volumetric Joule losses wJ (the influence of magnetization losses being neglected)
2
(2)
where oj- is the tensor of mechanical stresses, sj is the tensor of mechanical strains, u is the vector of mechanical displacements, /is the vector of the external volumetric forces, E is the temperature-dependent modulus of elasticity of the material, v is the temperature-dependent Poisson ratio of the material, Sj is the Kronecker delta, and e = skk(k = i, j).
The first part of (6) is the system of equilibrium equations describing the correlation between the mechanical stress tensor
aij =
ai1CT21CT31 ^ a12a22a32
(7)
VCT13CT23CT33 y and given volumetric forces components f
The second part of (6) is the system of cinematic equations representing the correlation between the strain tensor
sj =
whose magnitude decreases roughly exponentially with the distance from the surface of the heated body.
In fact, the complete solution of (1) is practically unfeasible due to relatively long time of the heating process. That is why the model was simplified by considering the magnetic field harmonic. Now (1) can be rewritten in terms of the phasor A of the magnetic vector potential A in the form
curl (curl a) + j ■ coy/Li A = /J ext, (3)
where j = and c denotes the angular frequency.
The computations must be now, however, carried out iteratively, and at every step the permeability ^ in any element containing ferromagnetic material has to be adjusted to the real value of the local magnetic flux density.
Temperature field. The temperature field T is described by the equation [5]
S1S21S31 A
S12S22S32 VS13S23S33 y
(8)
and components ut of the mechanical displacement vector.
The third part of (6) is the constitutive equations presenting the correlation between the mechanical stress
tensor a, and strain tensor s,.
In the system of equations (6) the following tensor and vector operations are used
lj,i
O,
dxl
dul
du,
«u j
= , Ujii =, (9)
dx,
dxl
where i = 1, 2, 3 and j = 1, 2, 3 are the indices of the coordinate axes.
In the axisymmetric case, the mechanical deformed state of the device is described by the system of equations [6, 7]
d ar
d x
>ee
d z
d r
d x r z + d a z z + x r z + f = Q.
+fr = Q.
d r d z
d Ur
d uz
(1Q)
where 2 is the thermal conductivity, p denotes the mass density and cp stands for the specific heat at a constant pressure.
The boundary condition on the surface of system is given by the formula
rz
S = -. S = -
°r r - ' °z z -, d r d z
d ur d uz
I r z
See = -
' r z ~ ~ 1 _ ' °ee
d z d r r
where orr, ozz, oee, Trz are the radial, axial, circumferential and shear mechanical stresses, respectively, srr, szz, see, yrz
r r
+
r
r
u
r
are the radial, axial, azimuthal and tangential mechanical strains, respectively, ur and uz are the radial and axial mechanical displacements, and f, f are the given radial and axial forces, respectively.
In order to simulate the mechanical state of the structure under consideration, elastic, plastic and thermal deformations are generally considered. Therefore, the strain tensor e^ is presented as a sum of the elastic j plastic e/ and thermal sj components [7], thus
(11)
The elastic strains are described in the following way [6]
— e — P — T S ij — S ij + S ij + S ij ■
e
E [(1+у)ау -ôyvcrkk ]
where k = i, j.
The thermal strains are represented in the form
sij — aAT8ij,
(12)
(13)
where a is the temperature-dependent coefficient of the temperature expansion of the material, and AT is the difference of temperatures.
Irreversible and, therefore, unacceptable plastic deformations described by tensor ijp are not considered in this paper.
Final remark. Many material parameters are significantly dependent on temperature. Mentioned can be, for example, electric conductivity, magnetic permeability, thermal conductivity, specific heat, coefficient of linear expansion, etc. All these dependences (as far as they are known) are included in the computations.
Numerical solution. Magnetic field is solved in a classic way using second-order finite element method. The discretization mesh must cover not only the system, but also its relatively large neighborhood. It is calculated independently and takes into account only temperature variations of the electric conductivity y and magnetic permeability p. The mesh remains the same in the course of the calculations, the geometric changes due to thermoe-lastic dilatation are neglected.
Solution of the temperature field is carried out by an algorithm based on the generalized Crank-Nicolson method that respects the temperature variations of parameters 2 and pcp.
The solution of the thermoelastic problem by the finite element method uses at every time step a linearized Lagrange variational equation for increments [7] in the form
jjAaiJSAey + oijSArij-AFiSAui) r ■ dS -jAP18Aulr ■ dl + S0 l0
+jj (aijSAetj - FiSAu1) r ■ dS - J PiSAu1r ■ dl = 0 ;
S0 l0
jj (AolJSAelJ +oij8Arij -AFlSAul )rdS - jAPl8Au1rdl +
S0 L0
jj(olJSAelJ -FiSAui )rdS - jPlSAulrdL, (14)
S0 L0 where S0 and L0 are the surface and boundary of the meridian cross-section of the structure, d and Ad denote the components of the stress tensor and their increments,
Де- and Д^- stand for the increments of the linear and nonlinear parts of the strain tensor, Ди- are the increments of components of the displacement vector, f ' and Д/' are the components of the volume loads (for example gravitational) and their increments in one time step and, finally, P' and ДР' denote the components of the surface load and their increments.
Unfortunately, in the devices under consideration, from time to time it is necessary to take into account the exchange of heat inside the gap between individual structural parts of the device. This exchange is realized through the contact zone that varies in time and that must be determined in the course of solution of the thermoelas-tic problem. There are several possibilities to take into account the dependence of the contact thermal conductivity Kn on the contact pressure pc.
In many cases it is enough to use only two values of Kn. For pc > 0 (absence of the contact) we put
Kn (T, pc) = K+, for pc < 0 (presence of the contact)
Kn (T, pc) = k -.
Here, Kn+ is a relatively small value (or even zero) describing the heat conductivity of the contact layer
through the ambient environment, while КЩ~Л/h,
where 2 is the thermal conductivity of the layer, and h denotes its roughness.
For more accurate computations it is recommended to use more complicated empirical formulae such as that given in [8].
Even when the mechanical problem is considered to be linear (without plastic deformations), searching of the unknown contact domain (or domains) between both parts in each step is realized by means of an iterative process.
Contact problem. Modeling of the contact interaction between two parts is based on the concept of a special contact finite element having no thickness. The functional for the temperature problem is supplemented with components corresponding to the thermal conductivity of this contact layer. The heat generated due to mutual sliding of both parts can also be taken into account, but the heat capacity (specific heat) of the contact layer is neglected. Using a special 1D 4-node finite elements we obtain a system of equations for the description of the thermo-contact problem.
The functional of the temperature problem for the contact layer can be written in the form [7]
I = 2 ¡K{sM,Vn{ -Ti)2 + Q(Sp,VnЛ1 + eis^n,v)T2]rA-
(15)
- J)Tlr0dT- Jq2(S^,an)T2r0dr + Lqi Lqi
+ Ja^Sp^nIt -T»((,CTn)]rr<A +
La1
+ Ja2|t2- T»(S,CTn)r0dr,
where Kn(Sp,a„) is the contact thermal conductivity
(not the contact thermal resistance), which can be estimated as a ratio «thermal conductivity of material of the layer» / «thickness of the layer», Sp are the subdomains
(«macro-elements») made of different materials, a„ is the
normal contact pressure, T1,T2 are the surface temperatures of contacting parts (bodies), Q(SL,on ,vr) denotes
the internal heat source arising due to relative sliding of surfaces at the given friction, vr stands for the relative
velocity of sliding contacting surfaces, q1, q2 are the heat fluxes into the contacting bodies, respectively, Tx, a1, a2 denote the temperature of the medium and heat transfer coefficients for the first and second bodies, respectively and r0 is the distance from the contact point to the z-axis (for the axisymmetric problem) or the thickness of the subdomain (for the planar problem).
The internal heat sources due to mutual sliding of both parts are calculated as follows
Q = OnVrff, (16)
where ff is the coefficient of friction between both surfaces.
These sources can be distributed between both contacting parts (bodies) by using the corresponding velocities vr1 and vr 2 of sliding for both surfaces
Q1 = Onvr1 ff , Q2 = Onvr2 ff (17)
to take into account boundary conditions of the 2nd kind for both surfaces at the contact point.
Taking (17) into account, the variation of the functional (15) can be written in the form
S = j[(/On) -T1 )(-ST!)]-
- jq\(s/,on)r0dr+
Lq1
+ ja1 (s/,On
M,°n )]ST1r0 dT +
+ jOnv1ff ST1r0 dr- (18)
Lk
- jq2(SL,On^ST2r0dr +
Lq2
+ ja2 (s/,On)[T2
- T„(s
M,°n )]s2r0dT +
La2
+ jOnv2 ff ST2r0 dT = 0
The temperature distribution in the 1D 4-node finite element can be represented in the form
T (r) = Tp (r) + Tj^j (r) ;
T2T = T,2<pt (r) + Tj2Wj (r), (19)
where Ti1,Ti2,Tj1, Tj2 are the nodal temperatures on the
contact surfaces of both bodies, respectively, and < (r), Pj (r) are the coordinate functions
T] -T
9, (-t) =-; 9j t) =
Tj -T,
t-t,
Tj
(20)
Substituting (19) into (20) and collecting terms at the same variations, we can obtain a system of equations for the thermo-contact element as follows
A,, , T, + a, , ,+1T,+1 + A, jTj + A,, j+1Tj+1 = b,; Ai+1,jTi + Ai+1,i+1Ti+1 + Ai+1, jTj + Ai+1,j+1Tj+1 = bt+1; Aj,,T, + Aj,i+iTi+i + Aj,jTj + Aj,j+{Tj+1 = bj; (21)
Aj+1,iTi + Aj+1,i+1Ti+1 + Aj+1, jTj + Aj+1,j+1Tj+1 = bj+1; where
Ak,k = j Kn (S/,On )pl(T)ro(r)dT +
Lk
j a1(S/, °n )Pk(T)ro(T)dT;
La1
Ak,k+1 = Ak+1,k = -jKn (SL,On)pk2(т)r0(т)dт;
Lk
Ak+1,k+1 = jKn (SL,On')pl(T)ro(T)dr +
Lk
+ ja2(S/,°n )Pp2(T)ro(T)dr;
La2
A,, j = A j,, = -j Kn (S/On )p, (r)Pj (r)r0(r)dr +
Lk
+ ja1(S/,On)p, (r)pj (r)r0(r)dr;
La1
A,+1, j+1 = Aj+u+1 =-j Kn (S/,On )<, (r)pj (T)r()(r)dT +
Lk
+ ja2(S/,On)<,(t)pJ (T)r()(T)dr;
La2
Ai, j+1 = Aj+1,i = Ai+1, j = Aj,i+1 = = - j Kn (S/,On) p, (r)pj (r)r0 (r)dr;
Lk
bk = j a1(S/ ,On T (S/ ,On )Pk (r)r0 (r)dr +
+ j q1(S/,on )Pk (r)ro(r)dT + jon (r)v1ffPk(r) r0(r)dr;
Lq1 Lk
bk+1 = ja2(S/,On T (S/,On )Pk (т)r0(т)dт +
La2
j q2(S/,On )Pk (т)ro(т)dт+ (22)
Lq2
+ jon (r)v2 ff Pk (t) ro(r)dr;
Lk
where k = i, j;
After discretizing the model consisting of (15) - (22), it remains to cope with the displacements and mechanical strains and stresses along particular sliding surfaces (see Figs. 2 and 3).
At the points of the interpenetration of contact surfaces in the normal direction with condition
iiln - ul -Sn < 0; (23)
where u1n , un2 are the displacements of the contact surfaces and Sn is the width of the gap in the direction of the normal, we introduce sufficiently high contact stiffness Cn
K
in the normal direction in order to prevent the surfaces from penetrating each other. Should friction be taken into account, we further introduce an analogous tangential stiffness C .
r
T
z
Fig. 2. Contact length Lk and gap S„ between contacting bodies 1 and 2
n
T
'V., u 1.
"P * J
On the other hand, the process of disassembly is characterized by fast heating of the system that causes different displacements in the chuck (that are relatively high) and in the drill shank (much lower). In a short time the pressed joint is released and the shank can be drawn out of the hole.
chuck
inductor drill shank
10
0
50 /
m / V jf / / -1 ^
0
20
40
Fig. 3. Nodal displacements of the contact surfaces
The stresses in the contact layer are now given by the formulae
an = Cn (u\ - ul ~5n); °T= CT (u 1 - u2 ~8T). (24)
As far as the condition <?nff <|ctr | of slipping is
satisfied, the tangential stresses are expressed as
r = anff sign(u 1 - u2 -8T), (25)
where u1 , u2 r are the shifts of both surfaces and ST is their difference at the beginning of the contact. In the zone of slipping the tangential stiffness is equal to zero and the functional is supplemented with the work of the friction forces on the corresponding shifts.
Illustrative example. A typical example is solved concerning the problem of assembly and disassembly of high-speed machine tools. The basic axisymmetric arrangement of the system is depicted in Fig. 4, together with the principal dimensions. The interference between the two connected parts is 0.01 mm.
During the process of assembly the chuck spindle has to be heated by the inductor until the internal diameter of the bore exceeds the diameter of the shank. The drill shank is then put into the bore and the system is cooled until we obtain a shrink fit.
60 80 100 120 z AB = 40 mm
Fig. 4. Basic disposition of investigated system (dimensions in mm)
The chuck is made from steel. Comparative computations for drills made from hard alloy and alloyed tool steel are carried out, the inductor remaining the same.
The principal data of the problem follow (at the room temperature T = 20 °C):
• steel chuck: electric conductivity y20 = 3.2 • 106 S/m, thermal conductivity 120 = 47 W/mK, heat capacity pcp20 = 4-106 J/m3K, coefficient of linear thermal expansion cct20 = 2-10~5 / K, modulus of elasticity E20 = 2-1011 Pa, Poisson ratio v= 0.3.
• Hard alloy: electric conductivity y20 = 2.8-106 S/m, thermal conductivity 120 = 85 W/mK, heat capacity pcp20 = 2-106 J/m3K, coefficient of linear thermal expansion aT20 = 0.5-10~5 / K, modulus of elasticity E20 = 5.3-1011 Pa, Poisson ratio v= 0.25.
• Alloyed tool steel: electric conductivity y20 = 2.5-106 S/m, thermal conductivity = 37 W/mK, heat capacity pcp20 = 3.2-106 J/m3K, coefficient of linear thermal expansion aT20 = 1.1-10"5 / K, modulus of elasticity E20 = 1.91011 Pa, Poisson ratio v= 0.3.
• Field current density Jext = 28.275 A/m2, f = 1 kHz.
• Cooling by forced air: a = 200 W/m2K.
The contact mechanical problem is the main part of numerical analysis. Fig. 5 and 6 represent the calculated distributions of temperature and contact pressure at various time instants of the assembly and disassembly processes for drills made from hard alloy and alloyed tool steel. The maps were performed for the cross section z = 85 mm. Some preliminary results can be found in [2]. The most important conclusions follow:
• Regarding assembly for both of materials (hard alloy or alloyed tool steel) of the drill: the chuck has to be heated during 3.5 s. Maximal temperature in it is 228 °C, maximal difference through the thickness is 86 °C.
• Regarding disassembly:
- Drill made of hard alloy: complete disassembly of the joint takes place in about 4 s after switching on the inductor. Maximum temperature in the chuck is 231 °C (external edge of its face), in the drill is 85 °C.
- Drill made of alloyed tool steel: duration of treatment is about 4.25 s.
r
T
z
T (°C) 130
110 90 70 50 30 10
(
T (°C) 190
170
150
130
110
90
70
50
30
10
" 11 -----—___
32 to
- / ,3
■ ,
8 r (mm) 10
^y/
■ / 3.5 „ r 37
-£—" __ii
- 2.5_. — — 0
2 4 6 8 r (mm) 10
Fig. 5. Temperature distribution along joint radius (z = 85 mm) at various time levels (in s): up - assembly, bottom - disassembly. -drill made of hard alloy,
■ drill made of alloyed tool steel
p (MPa) 0
-50
-100
-150
-200
0
p (MPa)
0
10
20
30 z (mm)
40
-100
-150
-200
10
20
30 z (mm) 40
Fig. 6. Contact pressure distribution at various time levels (in s): up - assembly, bottom - disassembly
-drill made of hard alloy,
.....drill made of alloyed tool steel
Conclusions. A 2D finite element technique for solution of multiphysics (electromagnetic, thermal and structural) problems taking into account contact interaction between structural parts is proposed. The contact problem solution is based on a concept of a special 1D contact finite element (layer) having no thickness. The presented technique is applied for the computer simulation of assembly of a shrink fit.
REFERENCES
1. Driesen J., Belmans R.J.M., Hameyer K. Finite-element modeling of thermal contact resistances and insulation layers in electrical machines. IEEE Transactions on Industry Applications,, 2001, vol.37, no.1, pp. 15-20. doi: 10.1109/28.903121.
2. Shulzhenko N.G., Gontarowsky P.P., Matyukhin Yu.I., Pantelyat M.G., Dolezel I., Ulrych B., Benes K. Computer modeling of induction heating-based assembly and disassembly of shrink fits. Acta Technica CSAV, 2004, vol.49, no.2, pp. 169-183.
3. Dolezel I., Karban P., Ulrych B., Pantelyat M.G., Matyukhin Yu.I., Gontarowsky P.P., Shulzhenko N.G. Limit operation regimes of actuators working on principle of thermoelasticity. IEEE Transactions on Magnetics, 2008, vol.44, no.6, pp. 810813. doi: 10.1109/tmag.2007.916226.
4. Chari M.V.K., Salon S.J. Numerical Methods in Electro-magnetism. Academic Press, New York, 2000. ISBN 9780126157604.
5. Holman J.P. Heat Transfer, 10th edition. McGraw-Hill, New York, 2002. ISBN 978-0073529363.
6. Timoshenko S., Goodier J.N. Theory of Elasticity, 3rd Edition. McGraw Hill, New York, 1970. ISBN 978-0070701229.
7. Podgorny A.N., Gontarowsky P.P., Kirkatsch B.N., Matyukhin Yu.I., Khavin G.L. Zadachi kontaktnogo vzaimodeist-viia elementov konstruktsii [Tasks of contact interaction in construction elements]. Kiev, Naukova Dumka Publ., 1989. ISBN 5-12-000891-7. (Rus).
8. Schlykov Yu.P., Ganin E.A., Carevskij S.N. Kontaktnoe termicheskoe soprotivlenie [Contact thermal resistance]. Moscow, Energiya Publ., 1977. (Rus).
Received 15.03.2016
M.G. Pantelyat1, Candidate of Physics and Mathematics, Associate Professor, Ivo Dolezel2, Professor,
1 National Technical University «Kharkiv Polytechnic Institute», 21, Kyrpychova Str., Kharkiv, 61002, Ukraine,
e-mail: [email protected]
2 Czech Academy of Sciences, Institute of Thermomechanics, 5, Dolejskova, 182 00 Praha 8, Czech Republic,
e-mail: [email protected]
4
6
