- Email: [email protected]
Application of bem to solve two-dimensional thermoelastic contact problems with convection and radiation conditions
Compufers 0
Pergamon
PII: s0045-7949(97)00052-7
& S~rucrures Vol. 66. No. I, pp. 115-125, 1998 1997 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0045-7949/98 $19.00 + 0.00
APPLXATION OF BEM TO SOLVE TWO-DIMENSIONAL THERMOELASTIC CONTACT PROBLEMS WITH CONVECTION AND RADIATION CONDITIONS M. P. Alonso, J. A. Garrido and A. Faces Department of Strength of Materials and Structures, University of Valladolid, Valladolid, Spain (Received 3 December 1996; accepted 27 March 1997)
Abstract-In this paper a technique, based on the boundary element method (BEM), for the treatment of two-dimen.sional thermoelastic contact problems without friction is developed. Thermal and elastic deformations are considered simultaneously. Additionally, the effects of the thermal resistance to the heat flow between both solids are included in the formulation. This magnitude depends on the local pressure at the contact zone and on the gap at the separated ones (interstitial regions). On the other hand, simultaneous convection and radiation conditions at the free boundary of the solids and at interstitial regions are admitted in the proposed technique of analysis. 0 1997 Elsevier Science Ltd Key words-Boundary
elements, Conduction, Contact resistance, Convection, Radiation, Thermoelastic
INTRODUCTION
ary. Since the nonlinearities and the imposed boundary conditions are located at the boundaries, the boundary element method (BEM) appears to be a suitable technique for the analysis of these kind of problems. In this paper an iterative schema of resolution is developed. In it the boundary integral equations corresponding to the linear elastoestatic problem and the potential one are operated simultaneously.
As is well known, the analysis of thermoelastic contact problems is a very important item for mechanical and civil engineering. This is because in any physical system, actions are transmitted throughout contacting solids, which have different temperatures in most practical cases. Unfortunately, previous study results are quite complex due to the existence of a great number of factors which have a decisive influence over the elastic and thermal fields. Thus, apart from the macroscopic characteristics of the solids, it is necessary to take into account the microscopic ones. These have a special interest at the contact zone. For that reason, the microscopic analysis of the interface constrictions has been the aim of numerous investigations [ 1,2]. On the other hand, as these systems are not isolated in practice, it is not possible to suppose adiabatic conditions. Hence it is imperative to evaluate the environmental influence which imposes a variation on the thermal and elastic conditions of the contact structure. In the same manner, it is necessary to consider the possible radiating heat interchange between the contacting solids in the regions surrounding the interface (interstitial zones). Furthermore, in these regions a heat transference by means of conduction and convection processes takes place between the contact structure and the interstitial fluids. All these factors, together with the nonlinearities inherent in thermoelastic contact problems, make it impossible to attain an analytical solution for the greater part of the practical cases. This is why a numerical treatment of the problem becomes necess-
THERMOELASTIC CONTACTPROBLEM Problem definition
Situations of plane contact without friction between elastic solids are considered, The solids are subject to small strains and small deformations due to the simultaneous action of any system of static loads applied at the boundary and any system of stationary thermal loads transmitted by conduction. The problem is depicted in Fig. 1. At any point y belonging to the boundary of each body k(k = A,B) a stress vector is associated, t$) (j = 1,2), as well as a displacement vector, u:(y), a temperature, t?b), and the derivative of temperature, qk(y), along the inner normal to the surface. The boundary of each solid can be divided in different zones in terms of the kind of conditions imposed on it (see Garrido et al. [3]). Thus, a first differentiation can be made between free regions, in which contact is not possible, and a potential contact zone,r$)A,B, in which contact can occur. In this last region it can be distinguished between the natural contact zone, I$@A*B (Ffo)AsB C TLp)AyB),corresponding to the interface of the structure with II5
M. P. Alonso et al.
116
r-tB$/L/iy_!p-d l&$c problem. Undeformed configuration.
riJ
E&k problem. Deformed configuration.
lhr’
Thermal problem. Undeformed configuration.
Thermal problem. Deformed configuration.
Fig. 1. Regionalization of the thermoelastic problem.
neither thermal nor static loads, and the final contact zone, I$B (I$” c Trp)*qB), obtained in the deformed configuration. The partitions of the free boundaries of both solids for the elastic and thermal problems are different. The regionalization corresponding to the elastic problem is depicted in Fig. l(a). Here, If: represents the free zone of the boundary of body k in which the displacement vector (utj) is determined and the stress vector is unknown (cj). On the contrary, in If the stress is prescribed in both directions (t$) and the displacement vector (I$) is unknown. In rfU there are mixed conditions: either the first component of the displacement vector and the second one of the stress vector (I$,,,, c,,,) are prescribed or, on the contrary, the second component of the displacement vector and the first one of the stress vector (&, &) are determined while the other two are unknown. The
deformed configuration of the system is shown in Fig. l(b), even though the partition is referred to the initial shape (this is depicted with a dashed line in the figure). The undeformed configuration and the final one for the thermal problem are depicted in Fig. l(c) and (d), respectively. The boundary conditions imposed at the free zones of both solids are a known temperature (@b) in I’: and a prescribed derivative (g(y) in I’:. It,, are the regions belonging to the potential contact zone that remain separated in the final configuration (I$ u r$ = rip’“). In these zones a heat interchange can take place by means of conduction (rEA*B) or forced convection (I$, jASB)through an interstitial fluid. Additionally, an energy transference by means of radiation can occur in these zones. Finally, in I;, there is a heat interchange by convection processes, natural or forced, between the surfaces of
.4pplication of BEM to solve two-dimensional thermoelastic contact problems
117
ratio of solid k, respectively). For more details about the dependence between thermal resistance and contact pressure see Refs [ 1,2]. Then, the thermal conditions at the contact zone are: (5) O,B(crB) = @(aA) + R[&aA)].
AA. d(~r*).
(6)
Fig. 2. Modified local system. the solids and the surrounding fluid and a radiant energy transference between the solids and an isothermal environment. So, the heat flow in each point y of these regions can be obtained as follows:
lk . AAY)= hk(!4. p:c $4-
e;,c y,]
+a.BYY) [ [h4]” - [gr(y)]4),(1) where L!$,is the ab’solute temperature of body k at point y, I!$(y) and hkh) represent the absolute temperature and the convective coefficient of the surrounded fluid, ik is the thermal conductivity, o is the Stefan-Boltzma:nn constant and Bkb) represents the factor of heat interchange by means of radiation between solid k anld its environment. Finally, e(y) represents the absolute temperature of the environment of body k. Elastic and thermal conditions in the potential contact zone
Each pair of contacting nodes c( (formed by aA and cr’) of rzP)A,B which will belong to I-tsB in the final configuration has to satisfy a set of elastic equations (the contact direction is assumed to be the straight line between the nodes, see Fig. 2):
4 (aA)-- u,“I(aB)= g@)(aA,aB), 6 (XA)-
t;1(aB)= 0,
Gc~~A, = t,B,(aB)=
0,
(2)
uz,@“) + ~&,(a~) = g@(aA,aB) - g(aA,aB),
(4)
(7)
&,(z”) = &(aB) = 0,
(8)
&(a”)
(9)
= t,Bh2(aB) = 0,
where g(aA,aB) is the final separation between aA and tlB. With regard to the thermal conditions in these zones, it will be considered that an energy transference can take place through an interstitial fluid. This heat interchange may be due to either conduction (rgAsB ) or forced convection (r$ jAsB).In this last case, the temperature of the fluid, or, is known. Consequently, it must be verified:
(3)
where g(“)(tlA,cxB)l:he initial gap between CI* and xB. It is possible to eliminate some variables related to body B through the use of the previous equations, and then only four elastic unknowns will remain each pair CI of contacting points: &(a”), &(a”), ~,“,(a~) and e,@“). In relation with the thermal conditions in the contact zone, it will be considered, in accordance with the analyses deduced by Barber et al. [4,5] and Comninou et al. [61, a value for the thermal resistance R[tt,(aA)] that varies inversely with the contact pressure (imperfect contact) when heat flows towards the solid with the lesser distortivity, 6, [sk = tik/lk(l + L”)], being gk and vk the coefficient of linear thermal expansion and Poisson’s CAS 6611-E
Again, it is possible to eliminate the variables of solid B. Consequently, only two thermal unknowns remain in each pair of contacting points: @‘(a”) and d(~r*). It is important to note that the direct imposition of Equation (6) is only possible in those cases in which thermal resistance is not considered. If the thermal resistance is a known function of the pressure, it is necessary to determine the local compression [t$l(aA)] for establishing Equation (6) between the temperatures of both solids. For that reason, a coupling between thermal and elastic problems takes place. On the other hand, the pairs of points belonging to ry in the final configuration must satisfy the following:
q$ jk(ak) + $.
[CJ$ jk(ak) - e,]= 0,
(10)
where I!?$,jk(ak) re p resents the temperature of point LX~ belonging to r$, )k, qi; jk(ak) is the derivative of the temperature at this point and h represents the convective coefficient of the interstitial fluid. If the heat interchange in r,AdBis due to a conduction process, the thermal jump in these zones depends on the separation between the surfaces of both solids: &)a(aa) _ &)*(aA) = g(aA,aB) ch
ch
If
. LA
. q”)A(aA) ch
’
(11)
where 1’ is the fluid conductivity. As g(aA,aB) is very small, in the previous formulation one-dimensional heat conduction has been assumed. Furthermore, the continuity of the flow in the normal direction to the surfaces must be verified:
M. P. Alonso et al.
118
It is advisable to emphasize that three equations
(12) are available in each point of the noncontacting It is important to notice that a stagnant fluid in F!$ produces the same effect as a thermal resistance of value:
g(aA,aB)
R[daA,aB)] = If
For convenience in the exposition of the formulation, it will be assumed that the size of the contact zone is known and also the values of R[$(cz*)] and R[g(aA,aB)]. Therefore, it is possible to disconnect the elastic problem from the thermal one. Later in this paper, the technique for coupling by means of an iterative process until adjusting the resistance law to the contact pressure in F,“,” and to the gap in rzA*B will be exposed.
NUMERICAL
FORMULATION
As is well known, in the absence of body forces, the application of Somigliana’s identity to an arbitrary boundary point of an elastic solid under small deformation conditions due to static loads and stationary thermal conduction ones implies the following boundary integral equation (see, for example, Brebbia et al. [7]):
sUi(53Y). s(k,Y) .
tj(Y) ’
-
dr( Y)
I-
=
O(Y) dU Y) A%’ + Afad x4’ = bO.
I-
J
- rQ;K~)
.q( y) . dU Y)>
(14)
where cd<) is the generally called free-term of the integral equation, whose value depends on the sh*arpness of *the boundary at point 5. Tensors T, (&Y) and Uii (
=
zones: Equation (15) and the two derived from Equation (14) (ij = 1,2). This number agrees with the total of unknowns associated to 5. In the same way, six unknowns for each pair of contacting points (u~,u,“,,u,“,,t,4,6l~,q~) remain in F$” after the application of the contact conditions and six integral equations are available, three for each contacting point. We will now refer to the numerical treatment of the problem. The boundary of each body is discretized in boundary elements. The discretization of the potential contact zone is accomplished so that the possible contacting elements are of the same size. Although several types of boundary elements can be employed in order to approximate the variables of the problem (straight or curved, constant, continuous or discontinuous ones, linear or parabolic, etc.), in this work discontinuous straight elements with linear approximation of temperatures, thermal gradients, displacements and stresses have been used. The nodes are placed inside the elements at a distance from the corner equal to a quarter of its length (for more details on this type of elements see Paris and Garrido [S]). Supposing that the size of r$B, the gap between both solids and the contact pressures are known, applying Equation (15) to each of the 2NA+ 2NB nodes of the discretization (Nk represents the number of boundary elements employed to discretize body k), after the integration over the elements and the imposition of the boundary thermal conditions, the following system of 2NA+2NB nonlinear equations is obtained:
q*b4 sI-
sre*(T,Y)
. Q(Y) . dU Y)
.4(Y) . dU Y)
(15)
Here, c(t) represents the free-term of the integral equation, 0*(&y) and q*(t,y) constitute the fundamental solution of the potential problem.
(16)
A superscript 0 denotes that the thermal problem is being referred to. A’ is a square matrix of 2NA+ 2NB rows containing the coefficients that result from the integrations indicated in Equation (15) except for the terms corresponding to the electromagnetic radiation. These coefficients have been organized taking the type of boundary conditions as a basis. Vector x0 includes the thermal unknowns of the problem. Afad is also a square matrix of 2NA + 2NB rows, and it contains the coefficients that result from integrating the terms corresponding to the electromagnetic radiation and zero in the rest of its elements. x4’ is a vector with the nonlinear unknowns of the problem (the fourth powers of the absolute temperatures of the nodes of rk). Finally, be is the vector containing the freeterms of the system, and it results from a proper operation of the integration coefficients and the prescribed boundary conditions at the noncontacting zones. The value of the free-term implied in Equation (15) is 0.5, as a discontinuous element has
Application of BEM to solve two-dimensional thermoelastic contact problems
a smooth boundary at every node (see, for example, Brebbia and Dominguez [9]). An iterative resolution of the previous system does not guarantee the convergence. This fact has been confirmed by Bialecki and Nowak [lo]. This problem has been solved by means of the application of a Newto:n-Rhapson method of adaptive step to solve Equa.tion (16). Using this technique convergence has been achieved in the whole of the analyzed cases. In the same way, with Equation (14) [if the size of IAsB is known and after solving Equation (16)], a set Gf 4NA + 4NB linear equations is obtained once the prescribed elastic boundary conditions have been applied: A“x” = b” + f’.
(17)
Superscript u refers to the contact problem with deformations due to static and thermal loads. The matrix A” contains the 4NA + 4NB integration coefficients that affect the unknown displacements and stresses at the nodes of the boundary elements. These displacements and stresses are included in the vector x’. b” contends the known values of the integral equation that results from operating the integration coefficients with the prescribed boundary conditions at the noncontacting zones. Finally, I@is a known vector derived from the integrations that appear in the second member of Equation (14). It is important to note that, the same as in the potential problem, the free-term of the elastic equation is c~=O.SS, (hii= 1 if i = j or a,=0 in another case), because of the continuity of the normal at the nodes of these elements. Solving Equation (17) the local components of the displacement and stress vectors that are unknown at the nodes of the boundary elements are obtained. ITER.QTIVE
PROCEDURE
In the previous formulation it was supposed that the length of TC A3Biand the value of the thermal resistance at each pair of nodes belonging to Irp)A,B were known. For that reason, it was possible to uncouple the thermal and elastic equations and to avoid the nonlinearities associated to the unknown of the size of the contact zone. Next, the strategy for determining these magnitudes is exposed. The solution of a generic problem must verify the following conditiotrs: (a) Normal stresses must be compressions over the whole of the contact zone, i.e. &@A)IO, vu.* E r$.
(18)
(b) No overlapping between the pairs of nodes of the potential contact zone that finally will not belong to T$B:
119
g(“)(aA,uB) - [l& (UA) + $1(cc”,] = g(aA,aB)20, VU*,U~ E rth.
(1%
(c) In each pair of contacting points, the resistance to the heat flow that was supposed in the resolution of the thermal problem (~*[tfq(c~*)]) has to agree with the value obtained by replacing the local compression determined after solving Equation (17) in the analytical expression of the resistance: R[$(c~*)] = ~*[&(a*)],
v~r E r;.
(20)
(d) If conduction through an interstitial fluid takes place, in each pair of nodes belonging to r(c)A3B,the resistance assumed in the resolution oFhthe thermal problem (R*[g(aA,aB)]) has to agree with the separation between both solids obtained after solving Equation (17): R[g(ctA,aB)] = R*[g(aA,aB)].
(21)
To make the convergence of the method possible a tolerance must be permitted in Equations (20) and (21). The designed process is always initiated supposing that the contact zone is f$J’)*,’ and considering perfect contact (without thermal resistance) in all pairs. Then, the system of thermal Equation (16) is solved and the thermal loads vector determined (8). Next, the system of elastic Equation (17) is solved and we proceed to check the compression in I$” [Equation (18)]. If this condition is not satisfied, the size of the contact zone is redefined and the pairs of nodes with tensile traction are excluded from this zone. The process begins again by solving the thermal equations under perfect contact assumption and considering a thermal resistance in the interstitial zones. This resistance is based on the gap obtained after solving the elastic problem in the previous iteration. On the contrary, if Equation (18) is fulfilled, either a check on overlaps [Equation (19)] or a convergence test [Equations (20) and (21)] follows. If Equation (19) is not satisfied it is necessary to include the overlapping nodes in lYt*B and to redefine it. The whole process begins again. In any case, the thermal resistance tests of Equations (20) and (21) are always performed by a contact zone with no tensile tractions and without overlapping at its vicinity. A simple convergence criterion for the thermal resistance in rLJ’)A,Bhas been implemented. This criterion lies in taking the value corresponding to the contact pressure distribution that was obtained in the previous iteration for determining the contact resistance. In the same manner, in the interstitial regions, the value of the gap obtained after solving the elastic problem is
M. P. Alonso et al.
120
emitted radiation in the subelement node a: can be obtained as follows: &(a:) =E (a:).
associated to
us [&(a:)]“,
(23)
where E (a:) and &(a:) represent the emissivity and the absolute temperature, respectively, at node a: of solid k. In the same manner, the quantity of reflected energy in each subelement is given by: &(a;)
where I
t
I
[p(af)
V B
9 Fig. 3. Geometric relations between the elements of r$B.
for calculating the resistance in the following step. The solution is reached when Equations (20) and (21) are verified at every pair of nodes belonging to rLp)A,B. employed
(24)
reflectivity
coefficient
E (a:) because the considered materials
are gray and opaque] and G(af) the quantity of energy that arrives to the subelement associated to a:. The sum of the quantities given by Equations (23) and (24) represents the loss of radiant energy in the considered node. This magnitude is called radiosity [J”(6)]. It is not determined until calculating G(a:) which depends fundamentally on the radiosity of the rest of the nodes and on the geometric view factor between the considered node and the other ones:
G(a”) =L$jL(ap). L($) j=l ENERGY INTERCHANGE BY MEANS OF RADIATION THE INTERSTITIAL REGIONS
p(af ) ’ G(aF),
the
P(& = l-
=
F(aT,af)
’
.J(ay),
(25)
IN
Next, an iterative technique for determining the interchange of radiant energy between the interstitial surfaces of both solids is exposed. The evaluation of these processes is complex because not only is it necessary to know the geometric view factors for each pair of elements of this zone, but also the quantity of heat that arrives to each of them proceeding from the rest is required. This heat quantity depends on the element temperatures, and so it will be referred to the nodes by employing nodal magnitudes in the analysis. Because of that, each element of riB has been divided in two subelements of equal size corresponding to one node each. Thus, the geometric view factor associated to the nodes c$ and a: is given by [l 1,121:
G(aB) =i~,(a~).F(~,aB).~(a~), L(ay) j=, J
(26)
where Equation (25) determines the quantity of radiant energy per unit of time that arrives to each node of l-2. Equation (26) is analogous for body B. In the previous analysis it has been assumed that the regions r$” adopt a convex shape, that is, a subelement of these zones is inside the shadow zone of any other subelement belonging to the same solid. Considering Equations (25) and (26), the radiosity at every node of the interstitial regions can be obtained as follows: J”(af)=[l-E(af)].&j
I F(aT,af)
=
cosO(aA 1 ’a!). I
cosO(aB a+)
2. r(af,a,p)
”
’
L(af).
(22) .$jL(a,B).
The different parameters that appear in the previous expression are represented in Fig. 3. L(af ) is the length of the subelement associated to node a?. The quantity of radiant energy falling upon each element of l-2” depends on the rate of energy that is absorbed and reflected by the rest of the elements. For determining this energetic interchange it has been considered that the materials are gray and opaque. For that reason, it is only necessary to calculate the quantity of heat that is emitted or reflected for each element of these regions. The
F(aT,$) . J’(g(jB)
j=l
+
E
(a:)
. u. [&
JB(aje) = [l-
.~.~L’“4)-F(o*,as)-J*(a*) 3 r-l
E (a;)]
)I’,
(27)
Application of BEM to solve two-dimensional thermoelastic contact problems
121
EXAMPLE
+ E
b,? ). Q. [e&,B)]4*
(28)
Equations
(27) and (28) constitute a system of equatisons (since N 2, the number of N a + N 5, nodes belonging to the interstitial region of solid k). For solving this system it is vital to know beforehand the temperatures of the N th + N :h nodes of IcAB. So, it is necessary to design an iterative procedure, beginning with the assumption that no radiant energy interchange exits between the nodes of the potential conta.ct zone that remain separated in the final configuration [lkqf(c$) = 0 1. Then, a first approximation for the value of the temperature in these nodes is obtained. This permits the calculation of each J’(c$) by solving Equations (27) and (28). The quantity of radiant energy that arrives to a node of I$, per unit of time and surface due to the heat interchange with every node of F,“h is given by:
The contacting solids depicted in Fig. 4 are considered. Solid A is made of aluminum, and solid B of steel. Thermal and elastic characteristics of both materials are indicated in the figure. Only half of the contact structure has been represented because it presents symmetry with regard to the vertical axis. As heat flows towards the solid with the lesser distortivity coefficient (sA= 1.49.10-‘m W-‘, dB= 3.12. lo-’ m W-l) thermal contact resistance is considered. The model proposed by Cooper et al. [2] has been employed. This model determines the coupling between resistance and pressure through a potential relation: $, @A)
R[&(aA)] =
___r 1 . Itan@ . 1.45 . [ H
1
-0.985
’
(32)
in which
and 111 -=2 2
.f($, A similar equation node of FF,: lB. q;(a;)
- JA(cq). can be obtained
(29) for every
= ‘GB(a;) - .?(ai”)
(30)
These values of the thermal flux must be considered in the integ,ral equation of the potential problem [Equation (1$1. Then, the heat flux at any point LX!of Fit, is given by Ik . qf(c$) + Ik. Q(af), since J ,,,($) is the heat flux due to either conduction or free convection at these regions. After this, a set of nonlinear equations analogous to Equation (16) is obtained. By solving this system a new approximation for the temperature at each node of the boundary is reached. The convergence is achieved when at each node of the boundary of the solids is verified: lek( j@ - @( Jo)“-“l
(31)
where ek(v)“’ represents the temperature of node y obtained in the n iteration and TOL is a prescribed tolerance.
(
>
’
where Hk represents the microhardness of body k, Q is the standard deviation of profile heights and ltanel represents the mean of absolute slope of a profile. The hardness and the parameters of the roughness distribution of each solid are: HA=930.8 MPa, HB=2413 MPa, Q = 9.510” m and [tan61 =0.08. For that reason, thermal resistance at the interface is given by the following expression: R[ti(aA)]
.JA(cY~)-JB(oL;).
1 FfF
= 1000. [ffj(~~~)]-~‘~*~m2 “C W-l.
(33)
First, the influence of this resistance over the different variables (temperature, heat flow and contact pressure) will be evaluated. The interstitial regions are assumed to be isolated. The temperature distributions at I’LpjAand F$‘)” obtained when perfect contact is assumed (R = 0), together with the ones reached when thermal resistance (R) is taken into account, are represented in Fig. 5. The thermal jump at the contact zone can clearly be appreciated when imperfect contact is considered. Moreover, a softer temperature transition from the interface to the separation region takes place in this case. Besides, when contact is perfect the temperatures of both solids are more similar at the interstitial regions (logically, the solids are at the same temperature at the contact zone). In both situations (perfect and imperfect contact) a strong increment of the heat flow at the right extreme of the interface appears (see Fig. 6). This is because, as the interstitial regions are isolated, a part of the flow is deviated to that zone. Nevertheless, this maximum decreases when thermal
122
M. P. Alonso
et at.
?A_?04 -& A.*=50 -& v*=O.3 v*=O.3 E*43
GPa
Et210 GPa aA2.34 .l$oc
B I
&B4*2
I p90>
w
)
0’0’0
JpJ)
($473.15
‘Y
“K
’
- p*=O.9 Fig. 4. Problem
__
0
definition.
--U--R=O, B
1
2
3 x ,W
Fig. 5. Temperature
distribution
at r$)A.B.
4
5
.1()-SC-1
-l
Application of BEM to solve two-dimensional thermoelastic contact problems
? 8; :: ::
--+--R=O,A
--t--RA --+--R=O,A
l-
o----0
123
1
2
3
4
0
5
i 2
1
:
‘\
3
d 4
5
x ,@w
x ,(cm)
Fig. 6. Heat flow distribution at rip)*.
Fig. 7. Pressure distribution at r:B.
resistance is considlered. This is related to the decrement of pressure that takes place at the last contacting nodes (see Fig. 7) which implies a rise in the resistance to the heat flow. The mentioned drop in flux is compensated by the increase of the energetic interchange at the rest of the contacting nodes. It can be appreciated, as it was previously indicated, that the heat flow is incoming in solid A (positive value of flow). The contact pressure distribution is represented in Fig. 7. In this case, the variations of contact pressure for considlering perfect or imperfect contact are not appreciable. Next, the influence of the different heat transference processes at the interstitial regions over the variables in the study is evaluated. For that reason, a heat interchange by means of conduction through an interstitial oil (If = 0.1 W m-’ ‘C-l) is considered. The results obtained for this case are compared with the ones reached when the interstitial regions are supposed to be isolated and with the ones obtained when a radiating energy transference takes place in these zones. The emissivities of both solids are: E* =0.2 and eB =0.4. Imperfect contact has been considered in any case.
The temperature distributions obtained at FfP)*,’ for the three cases are depicted in Fig. 8. As expected, a greater proximity between the temperatures of both solids exits when a heat interchange at the interstitial regions is possible. In this example, the influence of conduction through an interstitial oil (curves indicated with “cond”) over temperatures has a greater importance than radiation in these regions has (curves corresponding to “rad”). Logically, if another fluid with lesser conductivity (for example, air) were stagnant in the interstitial regions the heat transference would have lesser importance and the temperatures of both solids would be more approximate to the ones obtained when the separation region is isolated (curves indicated by “isol”). The heat flow at FLr)ASBis represented in Fig. 9. Practically, there is no difference between the fluxes at the contact zone in the cases of isolation and interstitial radiation. Nevertheless, when interstitial conduction is considered, the energetic interchange at the interface drops while an important increment of this takes place at the separation region. Here, the heat flow decreases slowly from a significant value at the vicinity of the interface to nearly zero
41.8 41.6 41,4 41.2 41 4’3,s
cl 0
1
2
3
4
x ,Mo
Fig. 8. Temperature distribution at I$r)As.
5
0
1
2
3
4
x ,(cno
Fig. 9. Heat flow distribution at I$@*.
5
124
M. P. Alonso et al.
35 ~_ --O-- conv,B ___&
&4
30 --
conv,A ---i-- RJ3 --+--
25 -20 --
Fig. 10. Temperature distribution at r$P)A.B.
Fig. 11. Heat flow distribution at l$‘jA.
at the last node of the potential contact zone. The radiating energy flow at the interstitial region is not very important and it has an approximately uniform value. The contact pressure distributions for the three cases are practically the same as the ones represented in Fig. 7. There is no important variation in compressions in the cases analyzed. Finally, the effects of the forced convection of oil (h = 1000 W m-* ‘C-‘, Br= 15X) in the interstitial regions has been analyzed. The obtained curves together with those ones corresponding to the problem with isolated interstitial zones are represented in Figs 10, 11 and 12. A strong cooling of both solids is appreciated since, as it can be seen in Fig. 11 (only solid A is represented), the bodies give up thermal energy to the fluid. Nevertheless, the flow is incoming through the contact zone of A, and this is reduced when oil is present because the totality of the flow proceeding from B does not arrive to the upper body. Important changes are also appreciated in the pressure distribution (see Fig. 12). It must be pointed out that a reduction of the length of r$” and a decrement of compression at the contact zone occur due to the greater contractions of both solids since they are at a lesser temperature when oil is present.
‘,
‘,
0 0
1
2
3
4
x ,@I
Fig. 12. Contact pressure distribution at l-2.
5
CONCLUSIONS
An integral formulation, based on the boundary element technique, for solving two-dimensional thermoelastic contact problems has been developed. This method is particularly appropriate for the treatment of this type of problem since the nonlinearities are associated at the boundary of the solids. The effects of thermal resistance at the interface, which varies inversely with contact pressure, have been included in the formulation. Then, thermal and elastic equations are coupled throughout the dependence resistance-pressure. In the developed technique both systems of equations are solved independently and the value of thermal resistance is obtained via an iterative process. Moreover, the possible heat interchange by means of conduction through a fluid contained at the interstitial regions has also been taken into account. This is equivalent to considering a thermal resistance at the interstitial regions which depends on the gap between the solids. Again, an iterative process has been employed for determining the heat flow rate in these zones. Besides, convection and radiation processes have been considered at the noncontatting regions of both solids. The radiation conditions introduce another nonlinearity in the problem. The Newton-Rhapson technique appears to be a suitable method for solving the thermal system of equations. Finally, an iterative process has been necessary for determining the radiating energy interchange at the interstitial regions. For future works, the inclusion of a suitable friction model and the extension of the method to the treatment of three-dimensional problems will be considered. Additionally, the consideration of the dependence between thermal and elastic characteristics of materials and temperature could be another important factor to develop.
.4pplication of BEM to solve two-dimensional thermoelastic contact problems REFERENCES
1. Madhusudana, C.V. and Fletcher, L.S. Contact heat transfer-the lasl decade. AZAA Journal, 1986, 24, 510-523. 2. Cooper,
M.G., Mikic, B.B. and Thermal contact conductance.
Yovanovich,
ZnrernarionuZ Journal of Heat and Mass Transfer, 1969, 12, 279-300.
M.M.
3. Garrido, J.A., Alonso, P., de1 Cairo, J.C. and Paris,
F. A boundary integral formulation for frictionless contact problems involving steady heat conduction. European Journal of Mechanics and Solids, 1995, 14, l-20. 4. Barber,
J.R. Contact problems involving a cooled punch. Journal of Elasticity, 1978, 8, 409-423. 5. Barber, J.R., Dundurs, J. and Comninou, M. Stability considerations in thermoelastic contact. Journal of Applied Mechunks,
1980, 47, 871-874.
6. Comninou,
M., Barber, J.R. and Dundurs, D. Heat conduction through a flat punch. Journal of Applied Mechanics, 1981, Qs, 871-875. 7. Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C., Boundary Element Methods: Techniques. Theory and Applicafions in Engineering. Springer, Berlin, 1984.
125
8. Paris, F. and Garrido, J. A., On the use of discontinuous elements in two-dimensional contact problemsIn Boundary Elements VZZ,Vol. 2. Springer, Berlin, 1985, pp. 13-27. 9. Brebbia, C. A. and Dominguez, J., Boundary Elements-An Introductory Course. Computational Mechanics Publications, McGraw-Hill, New York, 1988. 10. Bialecki, R. and Nowak, A. Boundary value problems in heat conduction with nonlinear material and nonlinear boundary conditions. Applied Mathematical Modelling, 1981, 5, 417-421.
11. Nowak, A. J., Solving coupled problems involving conduction, convection and thermal radiation. In Boundary
Elemeni
Method
in
Heat
Transfer.
Computational Mechanics Publications, Elsevier Applied Science, New York, 1992, pp. 145-173. 12. Bialecki, R., Solving Heat Radiation Problems using the Boundary Element Method. Computational Mechanics Publications, Elsevier Applied Science, New York, 1993.
Pergamon
PII: s0045-7949(97)00052-7
& S~rucrures Vol. 66. No. I, pp. 115-125, 1998 1997 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0045-7949/98 $19.00 + 0.00
APPLXATION OF BEM TO SOLVE TWO-DIMENSIONAL THERMOELASTIC CONTACT PROBLEMS WITH CONVECTION AND RADIATION CONDITIONS M. P. Alonso, J. A. Garrido and A. Faces Department of Strength of Materials and Structures, University of Valladolid, Valladolid, Spain (Received 3 December 1996; accepted 27 March 1997)
Abstract-In this paper a technique, based on the boundary element method (BEM), for the treatment of two-dimen.sional thermoelastic contact problems without friction is developed. Thermal and elastic deformations are considered simultaneously. Additionally, the effects of the thermal resistance to the heat flow between both solids are included in the formulation. This magnitude depends on the local pressure at the contact zone and on the gap at the separated ones (interstitial regions). On the other hand, simultaneous convection and radiation conditions at the free boundary of the solids and at interstitial regions are admitted in the proposed technique of analysis. 0 1997 Elsevier Science Ltd Key words-Boundary
elements, Conduction, Contact resistance, Convection, Radiation, Thermoelastic
INTRODUCTION
ary. Since the nonlinearities and the imposed boundary conditions are located at the boundaries, the boundary element method (BEM) appears to be a suitable technique for the analysis of these kind of problems. In this paper an iterative schema of resolution is developed. In it the boundary integral equations corresponding to the linear elastoestatic problem and the potential one are operated simultaneously.
As is well known, the analysis of thermoelastic contact problems is a very important item for mechanical and civil engineering. This is because in any physical system, actions are transmitted throughout contacting solids, which have different temperatures in most practical cases. Unfortunately, previous study results are quite complex due to the existence of a great number of factors which have a decisive influence over the elastic and thermal fields. Thus, apart from the macroscopic characteristics of the solids, it is necessary to take into account the microscopic ones. These have a special interest at the contact zone. For that reason, the microscopic analysis of the interface constrictions has been the aim of numerous investigations [ 1,2]. On the other hand, as these systems are not isolated in practice, it is not possible to suppose adiabatic conditions. Hence it is imperative to evaluate the environmental influence which imposes a variation on the thermal and elastic conditions of the contact structure. In the same manner, it is necessary to consider the possible radiating heat interchange between the contacting solids in the regions surrounding the interface (interstitial zones). Furthermore, in these regions a heat transference by means of conduction and convection processes takes place between the contact structure and the interstitial fluids. All these factors, together with the nonlinearities inherent in thermoelastic contact problems, make it impossible to attain an analytical solution for the greater part of the practical cases. This is why a numerical treatment of the problem becomes necess-
THERMOELASTIC CONTACTPROBLEM Problem definition
Situations of plane contact without friction between elastic solids are considered, The solids are subject to small strains and small deformations due to the simultaneous action of any system of static loads applied at the boundary and any system of stationary thermal loads transmitted by conduction. The problem is depicted in Fig. 1. At any point y belonging to the boundary of each body k(k = A,B) a stress vector is associated, t$) (j = 1,2), as well as a displacement vector, u:(y), a temperature, t?b), and the derivative of temperature, qk(y), along the inner normal to the surface. The boundary of each solid can be divided in different zones in terms of the kind of conditions imposed on it (see Garrido et al. [3]). Thus, a first differentiation can be made between free regions, in which contact is not possible, and a potential contact zone,r$)A,B, in which contact can occur. In this last region it can be distinguished between the natural contact zone, I$@A*B (Ffo)AsB C TLp)AyB),corresponding to the interface of the structure with II5
M. P. Alonso et al.
116
r-tB$/L/iy_!p-d l&$c problem. Undeformed configuration.
riJ
E&k problem. Deformed configuration.
lhr’
Thermal problem. Undeformed configuration.
Thermal problem. Deformed configuration.
Fig. 1. Regionalization of the thermoelastic problem.
neither thermal nor static loads, and the final contact zone, I$B (I$” c Trp)*qB), obtained in the deformed configuration. The partitions of the free boundaries of both solids for the elastic and thermal problems are different. The regionalization corresponding to the elastic problem is depicted in Fig. l(a). Here, If: represents the free zone of the boundary of body k in which the displacement vector (utj) is determined and the stress vector is unknown (cj). On the contrary, in If the stress is prescribed in both directions (t$) and the displacement vector (I$) is unknown. In rfU there are mixed conditions: either the first component of the displacement vector and the second one of the stress vector (I$,,,, c,,,) are prescribed or, on the contrary, the second component of the displacement vector and the first one of the stress vector (&, &) are determined while the other two are unknown. The
deformed configuration of the system is shown in Fig. l(b), even though the partition is referred to the initial shape (this is depicted with a dashed line in the figure). The undeformed configuration and the final one for the thermal problem are depicted in Fig. l(c) and (d), respectively. The boundary conditions imposed at the free zones of both solids are a known temperature (@b) in I’: and a prescribed derivative (g(y) in I’:. It,, are the regions belonging to the potential contact zone that remain separated in the final configuration (I$ u r$ = rip’“). In these zones a heat interchange can take place by means of conduction (rEA*B) or forced convection (I$, jASB)through an interstitial fluid. Additionally, an energy transference by means of radiation can occur in these zones. Finally, in I;, there is a heat interchange by convection processes, natural or forced, between the surfaces of
.4pplication of BEM to solve two-dimensional thermoelastic contact problems
117
ratio of solid k, respectively). For more details about the dependence between thermal resistance and contact pressure see Refs [ 1,2]. Then, the thermal conditions at the contact zone are: (5) O,B(crB) = @(aA) + R[&aA)].
AA. d(~r*).
(6)
Fig. 2. Modified local system. the solids and the surrounding fluid and a radiant energy transference between the solids and an isothermal environment. So, the heat flow in each point y of these regions can be obtained as follows:
lk . AAY)= hk(!4. p:c $4-
e;,c y,]
+a.BYY) [ [h4]” - [gr(y)]4),(1) where L!$,is the ab’solute temperature of body k at point y, I!$(y) and hkh) represent the absolute temperature and the convective coefficient of the surrounded fluid, ik is the thermal conductivity, o is the Stefan-Boltzma:nn constant and Bkb) represents the factor of heat interchange by means of radiation between solid k anld its environment. Finally, e(y) represents the absolute temperature of the environment of body k. Elastic and thermal conditions in the potential contact zone
Each pair of contacting nodes c( (formed by aA and cr’) of rzP)A,B which will belong to I-tsB in the final configuration has to satisfy a set of elastic equations (the contact direction is assumed to be the straight line between the nodes, see Fig. 2):
4 (aA)-- u,“I(aB)= g@)(aA,aB), 6 (XA)-
t;1(aB)= 0,
Gc~~A, = t,B,(aB)=
0,
(2)
uz,@“) + ~&,(a~) = g@(aA,aB) - g(aA,aB),
(4)
(7)
&,(z”) = &(aB) = 0,
(8)
&(a”)
(9)
= t,Bh2(aB) = 0,
where g(aA,aB) is the final separation between aA and tlB. With regard to the thermal conditions in these zones, it will be considered that an energy transference can take place through an interstitial fluid. This heat interchange may be due to either conduction (rgAsB ) or forced convection (r$ jAsB).In this last case, the temperature of the fluid, or, is known. Consequently, it must be verified:
(3)
where g(“)(tlA,cxB)l:he initial gap between CI* and xB. It is possible to eliminate some variables related to body B through the use of the previous equations, and then only four elastic unknowns will remain each pair CI of contacting points: &(a”), &(a”), ~,“,(a~) and e,@“). In relation with the thermal conditions in the contact zone, it will be considered, in accordance with the analyses deduced by Barber et al. [4,5] and Comninou et al. [61, a value for the thermal resistance R[tt,(aA)] that varies inversely with the contact pressure (imperfect contact) when heat flows towards the solid with the lesser distortivity, 6, [sk = tik/lk(l + L”)], being gk and vk the coefficient of linear thermal expansion and Poisson’s CAS 6611-E
Again, it is possible to eliminate the variables of solid B. Consequently, only two thermal unknowns remain in each pair of contacting points: @‘(a”) and d(~r*). It is important to note that the direct imposition of Equation (6) is only possible in those cases in which thermal resistance is not considered. If the thermal resistance is a known function of the pressure, it is necessary to determine the local compression [t$l(aA)] for establishing Equation (6) between the temperatures of both solids. For that reason, a coupling between thermal and elastic problems takes place. On the other hand, the pairs of points belonging to ry in the final configuration must satisfy the following:
q$ jk(ak) + $.
[CJ$ jk(ak) - e,]= 0,
(10)
where I!?$,jk(ak) re p resents the temperature of point LX~ belonging to r$, )k, qi; jk(ak) is the derivative of the temperature at this point and h represents the convective coefficient of the interstitial fluid. If the heat interchange in r,AdBis due to a conduction process, the thermal jump in these zones depends on the separation between the surfaces of both solids: &)a(aa) _ &)*(aA) = g(aA,aB) ch
ch
If
. LA
. q”)A(aA) ch
’
(11)
where 1’ is the fluid conductivity. As g(aA,aB) is very small, in the previous formulation one-dimensional heat conduction has been assumed. Furthermore, the continuity of the flow in the normal direction to the surfaces must be verified:
M. P. Alonso et al.
118
It is advisable to emphasize that three equations
(12) are available in each point of the noncontacting It is important to notice that a stagnant fluid in F!$ produces the same effect as a thermal resistance of value:
g(aA,aB)
R[daA,aB)] = If
For convenience in the exposition of the formulation, it will be assumed that the size of the contact zone is known and also the values of R[$(cz*)] and R[g(aA,aB)]. Therefore, it is possible to disconnect the elastic problem from the thermal one. Later in this paper, the technique for coupling by means of an iterative process until adjusting the resistance law to the contact pressure in F,“,” and to the gap in rzA*B will be exposed.
NUMERICAL
FORMULATION
As is well known, in the absence of body forces, the application of Somigliana’s identity to an arbitrary boundary point of an elastic solid under small deformation conditions due to static loads and stationary thermal conduction ones implies the following boundary integral equation (see, for example, Brebbia et al. [7]):
sUi(53Y). s(k,Y) .
tj(Y) ’
-
dr( Y)
I-
=
O(Y) dU Y) A%’ + Afad x4’ = bO.
I-
J
- rQ;K~)
.q( y) . dU Y)>
(14)
where cd<) is the generally called free-term of the integral equation, whose value depends on the sh*arpness of *the boundary at point 5. Tensors T, (&Y) and Uii (
=
zones: Equation (15) and the two derived from Equation (14) (ij = 1,2). This number agrees with the total of unknowns associated to 5. In the same way, six unknowns for each pair of contacting points (u~,u,“,,u,“,,t,4,6l~,q~) remain in F$” after the application of the contact conditions and six integral equations are available, three for each contacting point. We will now refer to the numerical treatment of the problem. The boundary of each body is discretized in boundary elements. The discretization of the potential contact zone is accomplished so that the possible contacting elements are of the same size. Although several types of boundary elements can be employed in order to approximate the variables of the problem (straight or curved, constant, continuous or discontinuous ones, linear or parabolic, etc.), in this work discontinuous straight elements with linear approximation of temperatures, thermal gradients, displacements and stresses have been used. The nodes are placed inside the elements at a distance from the corner equal to a quarter of its length (for more details on this type of elements see Paris and Garrido [S]). Supposing that the size of r$B, the gap between both solids and the contact pressures are known, applying Equation (15) to each of the 2NA+ 2NB nodes of the discretization (Nk represents the number of boundary elements employed to discretize body k), after the integration over the elements and the imposition of the boundary thermal conditions, the following system of 2NA+2NB nonlinear equations is obtained:
q*b4 sI-
sre*(T,Y)
. Q(Y) . dU Y)
.4(Y) . dU Y)
(15)
Here, c(t) represents the free-term of the integral equation, 0*(&y) and q*(t,y) constitute the fundamental solution of the potential problem.
(16)
A superscript 0 denotes that the thermal problem is being referred to. A’ is a square matrix of 2NA+ 2NB rows containing the coefficients that result from the integrations indicated in Equation (15) except for the terms corresponding to the electromagnetic radiation. These coefficients have been organized taking the type of boundary conditions as a basis. Vector x0 includes the thermal unknowns of the problem. Afad is also a square matrix of 2NA + 2NB rows, and it contains the coefficients that result from integrating the terms corresponding to the electromagnetic radiation and zero in the rest of its elements. x4’ is a vector with the nonlinear unknowns of the problem (the fourth powers of the absolute temperatures of the nodes of rk). Finally, be is the vector containing the freeterms of the system, and it results from a proper operation of the integration coefficients and the prescribed boundary conditions at the noncontacting zones. The value of the free-term implied in Equation (15) is 0.5, as a discontinuous element has
Application of BEM to solve two-dimensional thermoelastic contact problems
a smooth boundary at every node (see, for example, Brebbia and Dominguez [9]). An iterative resolution of the previous system does not guarantee the convergence. This fact has been confirmed by Bialecki and Nowak [lo]. This problem has been solved by means of the application of a Newto:n-Rhapson method of adaptive step to solve Equa.tion (16). Using this technique convergence has been achieved in the whole of the analyzed cases. In the same way, with Equation (14) [if the size of IAsB is known and after solving Equation (16)], a set Gf 4NA + 4NB linear equations is obtained once the prescribed elastic boundary conditions have been applied: A“x” = b” + f’.
(17)
Superscript u refers to the contact problem with deformations due to static and thermal loads. The matrix A” contains the 4NA + 4NB integration coefficients that affect the unknown displacements and stresses at the nodes of the boundary elements. These displacements and stresses are included in the vector x’. b” contends the known values of the integral equation that results from operating the integration coefficients with the prescribed boundary conditions at the noncontacting zones. Finally, I@is a known vector derived from the integrations that appear in the second member of Equation (14). It is important to note that, the same as in the potential problem, the free-term of the elastic equation is c~=O.SS, (hii= 1 if i = j or a,=0 in another case), because of the continuity of the normal at the nodes of these elements. Solving Equation (17) the local components of the displacement and stress vectors that are unknown at the nodes of the boundary elements are obtained. ITER.QTIVE
PROCEDURE
In the previous formulation it was supposed that the length of TC A3Biand the value of the thermal resistance at each pair of nodes belonging to Irp)A,B were known. For that reason, it was possible to uncouple the thermal and elastic equations and to avoid the nonlinearities associated to the unknown of the size of the contact zone. Next, the strategy for determining these magnitudes is exposed. The solution of a generic problem must verify the following conditiotrs: (a) Normal stresses must be compressions over the whole of the contact zone, i.e. &@A)IO, vu.* E r$.
(18)
(b) No overlapping between the pairs of nodes of the potential contact zone that finally will not belong to T$B:
119
g(“)(aA,uB) - [l& (UA) + $1(cc”,] = g(aA,aB)20, VU*,U~ E rth.
(1%
(c) In each pair of contacting points, the resistance to the heat flow that was supposed in the resolution of the thermal problem (~*[tfq(c~*)]) has to agree with the value obtained by replacing the local compression determined after solving Equation (17) in the analytical expression of the resistance: R[$(c~*)] = ~*[&(a*)],
v~r E r;.
(20)
(d) If conduction through an interstitial fluid takes place, in each pair of nodes belonging to r(c)A3B,the resistance assumed in the resolution oFhthe thermal problem (R*[g(aA,aB)]) has to agree with the separation between both solids obtained after solving Equation (17): R[g(ctA,aB)] = R*[g(aA,aB)].
(21)
To make the convergence of the method possible a tolerance must be permitted in Equations (20) and (21). The designed process is always initiated supposing that the contact zone is f$J’)*,’ and considering perfect contact (without thermal resistance) in all pairs. Then, the system of thermal Equation (16) is solved and the thermal loads vector determined (8). Next, the system of elastic Equation (17) is solved and we proceed to check the compression in I$” [Equation (18)]. If this condition is not satisfied, the size of the contact zone is redefined and the pairs of nodes with tensile traction are excluded from this zone. The process begins again by solving the thermal equations under perfect contact assumption and considering a thermal resistance in the interstitial zones. This resistance is based on the gap obtained after solving the elastic problem in the previous iteration. On the contrary, if Equation (18) is fulfilled, either a check on overlaps [Equation (19)] or a convergence test [Equations (20) and (21)] follows. If Equation (19) is not satisfied it is necessary to include the overlapping nodes in lYt*B and to redefine it. The whole process begins again. In any case, the thermal resistance tests of Equations (20) and (21) are always performed by a contact zone with no tensile tractions and without overlapping at its vicinity. A simple convergence criterion for the thermal resistance in rLJ’)A,Bhas been implemented. This criterion lies in taking the value corresponding to the contact pressure distribution that was obtained in the previous iteration for determining the contact resistance. In the same manner, in the interstitial regions, the value of the gap obtained after solving the elastic problem is
M. P. Alonso et al.
120
emitted radiation in the subelement node a: can be obtained as follows: &(a:) =E (a:).
associated to
us [&(a:)]“,
(23)
where E (a:) and &(a:) represent the emissivity and the absolute temperature, respectively, at node a: of solid k. In the same manner, the quantity of reflected energy in each subelement is given by: &(a;)
where I
t
I
[p(af)
V B
9 Fig. 3. Geometric relations between the elements of r$B.
for calculating the resistance in the following step. The solution is reached when Equations (20) and (21) are verified at every pair of nodes belonging to rLp)A,B. employed
(24)
reflectivity
coefficient
E (a:) because the considered materials
are gray and opaque] and G(af) the quantity of energy that arrives to the subelement associated to a:. The sum of the quantities given by Equations (23) and (24) represents the loss of radiant energy in the considered node. This magnitude is called radiosity [J”(6)]. It is not determined until calculating G(a:) which depends fundamentally on the radiosity of the rest of the nodes and on the geometric view factor between the considered node and the other ones:
G(a”) =L$jL(ap). L($) j=l ENERGY INTERCHANGE BY MEANS OF RADIATION THE INTERSTITIAL REGIONS
p(af ) ’ G(aF),
the
P(& = l-
=
F(aT,af)
’
.J(ay),
(25)
IN
Next, an iterative technique for determining the interchange of radiant energy between the interstitial surfaces of both solids is exposed. The evaluation of these processes is complex because not only is it necessary to know the geometric view factors for each pair of elements of this zone, but also the quantity of heat that arrives to each of them proceeding from the rest is required. This heat quantity depends on the element temperatures, and so it will be referred to the nodes by employing nodal magnitudes in the analysis. Because of that, each element of riB has been divided in two subelements of equal size corresponding to one node each. Thus, the geometric view factor associated to the nodes c$ and a: is given by [l 1,121:
G(aB) =i~,(a~).F(~,aB).~(a~), L(ay) j=, J
(26)
where Equation (25) determines the quantity of radiant energy per unit of time that arrives to each node of l-2. Equation (26) is analogous for body B. In the previous analysis it has been assumed that the regions r$” adopt a convex shape, that is, a subelement of these zones is inside the shadow zone of any other subelement belonging to the same solid. Considering Equations (25) and (26), the radiosity at every node of the interstitial regions can be obtained as follows: J”(af)=[l-E(af)].&j
I F(aT,af)
=
cosO(aA 1 ’a!). I
cosO(aB a+)
2. r(af,a,p)
”
’
L(af).
(22) .$jL(a,B).
The different parameters that appear in the previous expression are represented in Fig. 3. L(af ) is the length of the subelement associated to node a?. The quantity of radiant energy falling upon each element of l-2” depends on the rate of energy that is absorbed and reflected by the rest of the elements. For determining this energetic interchange it has been considered that the materials are gray and opaque. For that reason, it is only necessary to calculate the quantity of heat that is emitted or reflected for each element of these regions. The
F(aT,$) . J’(g(jB)
j=l
+
E
(a:)
. u. [&
JB(aje) = [l-
.~.~L’“4)-F(o*,as)-J*(a*) 3 r-l
E (a;)]
)I’,
(27)
Application of BEM to solve two-dimensional thermoelastic contact problems
121
EXAMPLE
+ E
b,? ). Q. [e&,B)]4*
(28)
Equations
(27) and (28) constitute a system of equatisons (since N 2, the number of N a + N 5, nodes belonging to the interstitial region of solid k). For solving this system it is vital to know beforehand the temperatures of the N th + N :h nodes of IcAB. So, it is necessary to design an iterative procedure, beginning with the assumption that no radiant energy interchange exits between the nodes of the potential conta.ct zone that remain separated in the final configuration [lkqf(c$) = 0 1. Then, a first approximation for the value of the temperature in these nodes is obtained. This permits the calculation of each J’(c$) by solving Equations (27) and (28). The quantity of radiant energy that arrives to a node of I$, per unit of time and surface due to the heat interchange with every node of F,“h is given by:
The contacting solids depicted in Fig. 4 are considered. Solid A is made of aluminum, and solid B of steel. Thermal and elastic characteristics of both materials are indicated in the figure. Only half of the contact structure has been represented because it presents symmetry with regard to the vertical axis. As heat flows towards the solid with the lesser distortivity coefficient (sA= 1.49.10-‘m W-‘, dB= 3.12. lo-’ m W-l) thermal contact resistance is considered. The model proposed by Cooper et al. [2] has been employed. This model determines the coupling between resistance and pressure through a potential relation: $, @A)
R[&(aA)] =
___r 1 . Itan@ . 1.45 . [ H
1
-0.985
’
(32)
in which
and 111 -=2 2
.f($, A similar equation node of FF,: lB. q;(a;)
- JA(cq). can be obtained
(29) for every
= ‘GB(a;) - .?(ai”)
(30)
These values of the thermal flux must be considered in the integ,ral equation of the potential problem [Equation (1$1. Then, the heat flux at any point LX!of Fit, is given by Ik . qf(c$) + Ik. Q(af), since J ,,,($) is the heat flux due to either conduction or free convection at these regions. After this, a set of nonlinear equations analogous to Equation (16) is obtained. By solving this system a new approximation for the temperature at each node of the boundary is reached. The convergence is achieved when at each node of the boundary of the solids is verified: lek( j@ - @( Jo)“-“l
(31)
where ek(v)“’ represents the temperature of node y obtained in the n iteration and TOL is a prescribed tolerance.
(
>
’
where Hk represents the microhardness of body k, Q is the standard deviation of profile heights and ltanel represents the mean of absolute slope of a profile. The hardness and the parameters of the roughness distribution of each solid are: HA=930.8 MPa, HB=2413 MPa, Q = 9.510” m and [tan61 =0.08. For that reason, thermal resistance at the interface is given by the following expression: R[ti(aA)]
.JA(cY~)-JB(oL;).
1 FfF
= 1000. [ffj(~~~)]-~‘~*~m2 “C W-l.
(33)
First, the influence of this resistance over the different variables (temperature, heat flow and contact pressure) will be evaluated. The interstitial regions are assumed to be isolated. The temperature distributions at I’LpjAand F$‘)” obtained when perfect contact is assumed (R = 0), together with the ones reached when thermal resistance (R) is taken into account, are represented in Fig. 5. The thermal jump at the contact zone can clearly be appreciated when imperfect contact is considered. Moreover, a softer temperature transition from the interface to the separation region takes place in this case. Besides, when contact is perfect the temperatures of both solids are more similar at the interstitial regions (logically, the solids are at the same temperature at the contact zone). In both situations (perfect and imperfect contact) a strong increment of the heat flow at the right extreme of the interface appears (see Fig. 6). This is because, as the interstitial regions are isolated, a part of the flow is deviated to that zone. Nevertheless, this maximum decreases when thermal
122
M. P. Alonso
et at.
?A_?04 -& A.*=50 -& v*=O.3 v*=O.3 E*43
GPa
Et210 GPa aA2.34 .l$oc
B I
&B4*2
I p90>
w
)
0’0’0
JpJ)
($473.15
‘Y
“K
’
- p*=O.9 Fig. 4. Problem
__
0
definition.
--U--R=O, B
1
2
3 x ,W
Fig. 5. Temperature
distribution
at r$)A.B.
4
5
.1()-SC-1
-l
Application of BEM to solve two-dimensional thermoelastic contact problems
? 8; :: ::
--+--R=O,A
--t--RA --+--R=O,A
l-
o----0
123
1
2
3
4
0
5
i 2
1
:
‘\
3
d 4
5
x ,@w
x ,(cm)
Fig. 6. Heat flow distribution at rip)*.
Fig. 7. Pressure distribution at r:B.
resistance is considlered. This is related to the decrement of pressure that takes place at the last contacting nodes (see Fig. 7) which implies a rise in the resistance to the heat flow. The mentioned drop in flux is compensated by the increase of the energetic interchange at the rest of the contacting nodes. It can be appreciated, as it was previously indicated, that the heat flow is incoming in solid A (positive value of flow). The contact pressure distribution is represented in Fig. 7. In this case, the variations of contact pressure for considlering perfect or imperfect contact are not appreciable. Next, the influence of the different heat transference processes at the interstitial regions over the variables in the study is evaluated. For that reason, a heat interchange by means of conduction through an interstitial oil (If = 0.1 W m-’ ‘C-l) is considered. The results obtained for this case are compared with the ones reached when the interstitial regions are supposed to be isolated and with the ones obtained when a radiating energy transference takes place in these zones. The emissivities of both solids are: E* =0.2 and eB =0.4. Imperfect contact has been considered in any case.
The temperature distributions obtained at FfP)*,’ for the three cases are depicted in Fig. 8. As expected, a greater proximity between the temperatures of both solids exits when a heat interchange at the interstitial regions is possible. In this example, the influence of conduction through an interstitial oil (curves indicated with “cond”) over temperatures has a greater importance than radiation in these regions has (curves corresponding to “rad”). Logically, if another fluid with lesser conductivity (for example, air) were stagnant in the interstitial regions the heat transference would have lesser importance and the temperatures of both solids would be more approximate to the ones obtained when the separation region is isolated (curves indicated by “isol”). The heat flow at FLr)ASBis represented in Fig. 9. Practically, there is no difference between the fluxes at the contact zone in the cases of isolation and interstitial radiation. Nevertheless, when interstitial conduction is considered, the energetic interchange at the interface drops while an important increment of this takes place at the separation region. Here, the heat flow decreases slowly from a significant value at the vicinity of the interface to nearly zero
41.8 41.6 41,4 41.2 41 4’3,s
cl 0
1
2
3
4
x ,Mo
Fig. 8. Temperature distribution at I$r)As.
5
0
1
2
3
4
x ,(cno
Fig. 9. Heat flow distribution at I$@*.
5
124
M. P. Alonso et al.
35 ~_ --O-- conv,B ___&
&4
30 --
conv,A ---i-- RJ3 --+--
25 -20 --
Fig. 10. Temperature distribution at r$P)A.B.
Fig. 11. Heat flow distribution at l$‘jA.
at the last node of the potential contact zone. The radiating energy flow at the interstitial region is not very important and it has an approximately uniform value. The contact pressure distributions for the three cases are practically the same as the ones represented in Fig. 7. There is no important variation in compressions in the cases analyzed. Finally, the effects of the forced convection of oil (h = 1000 W m-* ‘C-‘, Br= 15X) in the interstitial regions has been analyzed. The obtained curves together with those ones corresponding to the problem with isolated interstitial zones are represented in Figs 10, 11 and 12. A strong cooling of both solids is appreciated since, as it can be seen in Fig. 11 (only solid A is represented), the bodies give up thermal energy to the fluid. Nevertheless, the flow is incoming through the contact zone of A, and this is reduced when oil is present because the totality of the flow proceeding from B does not arrive to the upper body. Important changes are also appreciated in the pressure distribution (see Fig. 12). It must be pointed out that a reduction of the length of r$” and a decrement of compression at the contact zone occur due to the greater contractions of both solids since they are at a lesser temperature when oil is present.
‘,
‘,
0 0
1
2
3
4
x ,@I
Fig. 12. Contact pressure distribution at l-2.
5
CONCLUSIONS
An integral formulation, based on the boundary element technique, for solving two-dimensional thermoelastic contact problems has been developed. This method is particularly appropriate for the treatment of this type of problem since the nonlinearities are associated at the boundary of the solids. The effects of thermal resistance at the interface, which varies inversely with contact pressure, have been included in the formulation. Then, thermal and elastic equations are coupled throughout the dependence resistance-pressure. In the developed technique both systems of equations are solved independently and the value of thermal resistance is obtained via an iterative process. Moreover, the possible heat interchange by means of conduction through a fluid contained at the interstitial regions has also been taken into account. This is equivalent to considering a thermal resistance at the interstitial regions which depends on the gap between the solids. Again, an iterative process has been employed for determining the heat flow rate in these zones. Besides, convection and radiation processes have been considered at the noncontatting regions of both solids. The radiation conditions introduce another nonlinearity in the problem. The Newton-Rhapson technique appears to be a suitable method for solving the thermal system of equations. Finally, an iterative process has been necessary for determining the radiating energy interchange at the interstitial regions. For future works, the inclusion of a suitable friction model and the extension of the method to the treatment of three-dimensional problems will be considered. Additionally, the consideration of the dependence between thermal and elastic characteristics of materials and temperature could be another important factor to develop.
.4pplication of BEM to solve two-dimensional thermoelastic contact problems REFERENCES
1. Madhusudana, C.V. and Fletcher, L.S. Contact heat transfer-the lasl decade. AZAA Journal, 1986, 24, 510-523. 2. Cooper,
M.G., Mikic, B.B. and Thermal contact conductance.
Yovanovich,
ZnrernarionuZ Journal of Heat and Mass Transfer, 1969, 12, 279-300.
M.M.
3. Garrido, J.A., Alonso, P., de1 Cairo, J.C. and Paris,
F. A boundary integral formulation for frictionless contact problems involving steady heat conduction. European Journal of Mechanics and Solids, 1995, 14, l-20. 4. Barber,
J.R. Contact problems involving a cooled punch. Journal of Elasticity, 1978, 8, 409-423. 5. Barber, J.R., Dundurs, J. and Comninou, M. Stability considerations in thermoelastic contact. Journal of Applied Mechunks,
1980, 47, 871-874.
6. Comninou,
M., Barber, J.R. and Dundurs, D. Heat conduction through a flat punch. Journal of Applied Mechanics, 1981, Qs, 871-875. 7. Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C., Boundary Element Methods: Techniques. Theory and Applicafions in Engineering. Springer, Berlin, 1984.
125
8. Paris, F. and Garrido, J. A., On the use of discontinuous elements in two-dimensional contact problemsIn Boundary Elements VZZ,Vol. 2. Springer, Berlin, 1985, pp. 13-27. 9. Brebbia, C. A. and Dominguez, J., Boundary Elements-An Introductory Course. Computational Mechanics Publications, McGraw-Hill, New York, 1988. 10. Bialecki, R. and Nowak, A. Boundary value problems in heat conduction with nonlinear material and nonlinear boundary conditions. Applied Mathematical Modelling, 1981, 5, 417-421.
11. Nowak, A. J., Solving coupled problems involving conduction, convection and thermal radiation. In Boundary
Elemeni
Method
in
Heat
Transfer.
Computational Mechanics Publications, Elsevier Applied Science, New York, 1992, pp. 145-173. 12. Bialecki, R., Solving Heat Radiation Problems using the Boundary Element Method. Computational Mechanics Publications, Elsevier Applied Science, New York, 1993.