- Email: [email protected]
finite element method for transient thermoelastic problem of composite hollow cylinder
Compurers & Strucrures Vol. 36, No. 5. pp. .353-860, 1990 Printed in Great Britain.
0045-7949/90 13.00 + o&l Q 1990 Pqamon Press plc
HYBRID LAPLACE TRANSFORM/FINITE ELEMENT METHOD FOR TRANSIENT THERMOELASTIC PROBLEM OF COMPOSITE HOLLOW CYLINDER L. S. CHEN and H. S. CHU Department
of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan 30049, Republic of China (Received 24 April 1989)
Abstract-The transient thermal stresses and temperature distribution of a composite hollow cylinder are investigated theoretically. The cylinder is composed of two different materials and heated by a moving line heat source on its inner boundary and cooled convectively on the exterior surface. The heat source is assumed to be axisymmetric, moving along the axis of the hollow cylinder with constant velocity, and decaying exponentially with time. The hybrid Laplace transform/finite element method is adopted for the solution of this problem. Finally, a numerical scheme, Fourier series technique, is utilized to calculate the inverse transform. The results show that, with the exception of the tangential thermal stress, the positions of the peak values of the temperature and thermal stresses are at exactly the same position as that at which the line heat source is applied. The effects of a number of parameters are also investigated.
%.Uz
NOTATION
PI, V%,l gradient matrices for temperature and stress
KI, K2
field, respectively specific heat c/J d exponential decaying coefficient of line heat source D dimensionless exponential decaying coefficient of line heat source material property matrices for temperature and PI, RI stress field, respectively element number e E total number of elements IFI, {fm} force vectors for temperature and stress field, respectively h heat transfer coefficient thermal conductivities of inner and outer k,,k, cylinder, respectively K thermal conductivity ratio of thP compo?te hollow cylinder stiffness matrices for temperature Jnd strc .s WI, Kl field, respectively total number of nodes interpolation function for both temperature and IMNI stress field number of nodes for each element P heat flux of the line moving heat source e cylindrical coordinates dimensionless inner radius of composite hello :z cylinder dimensionless interface radius of composite R? hollow cylinder dimensionless outer radius of composite hollow R, cylinder dimensionless cylindrical coordinates R, Z s the variable used in the Laplace transform t time T temperature ambient temperature T, 1’ velocity of heat source V dimensionless velocity of heat source
K e ? 111.& p “1, “2
thermal expansion coefficients of inner and outer cylinder, respectively thermal diffusivities of inner and outer cylinder, respectively thermal diffusivity ratio of composite hollow cylinder dimensionless temperature dimensionless time shear moduli of inner and outer cylinder, respectively shear modulus ratio of composite hollow cylinder Poisson’s ratios of inner and outer cylinder, respectively stress field mass density
INTRODUCIION
The problem of temperature and thermal stress distribution in a cylinder due to a time dependent heat source has many engineering applications, for example solid propellants, rocket motors, gun barrels, nuclear reactor elements, etc. and has thus been the subject of numerous investigations [l]. The unsteady state thermoelastic problems of a cylinder have also been studied intensively, for example the one dimensional quasi-static coupled thermoelastic ,>roblem of an axisymmetric infinitely long cylinder [2,3]. The transient uncoupled thermal stress distriTVrtion has also been investigated in the axisymmetric caoling case [4], in the cooling problem during tran.ient convection heating [S] and in a cylinder subJetted to the conductive heat transfer with thermal boundary conductance [6]. The transient thermoelastic contact problem of an infinitely long solid cylinder
853
L. S. C&EN and H. S. CHU
854
perfectly bonded by a heated rigid circular ring has been analyzed with homogeneous materials [7,8]. The above studies were all based on the analytical models. The unsteady-state problems have also been solved by finite element methods; for instance, the temperature distribution due to a moving heat source has been solved by the moving mesh finite element method [9], and the thermal stress of sintering has been studied by the moving grid method [lo]. Finally, since the use of the finite element method is time consuming for a large problem, the hybrid Laplace transform/finite element method has been used [ll]. However, most of the previous studies deal only with cylinders with either infinite length or homogenous material, while few of them deal with cylinders with finite length and composite materials. In order to study the behavior of transient thermal stresses and the temperature distribution of a composite hollow cylinder with finite length, the hybrid Laplace transform/finite element method is adopted. The cylinder is composed of two different materials and is heated by a moving line heat source on its inner boundary and cooled convectively on the exterior surface. The heat source is assumed to be axisymmetric, moving along the axis of the hollow cylinder with a constant velocity and decaying exponentially with time. In this study, the coupling of temperature and strain rate as well as the inertia effect in the stress field are ignored. The temperature and thermal stress fields are obtained by applying Laplace transform with respect to the time variable and then the finite element method is used to solve the problem in the transformed region. Then, the inverse Laplace transform is accomplished by a numerical method, the Fourier series technique [12, 131. Numerical results for both temperature and thermal stresses are presented graphically with respect to a number of parameters.
r
FLFig. I. Physical model and coordinate system of composite hollow cylinder. axisymmetrical
a2Ti
(
a2Ti
1 ar,
F+;yjy’g
)
1 aTi =K_dtv f
rl
o
(I)
where i = 1,2, which refer to the regions I and 2, respectively. The initial and boundary conditions for the temperature field are Ti=T,, aT. -=O,
z=O
az
t=O
(2) z=L
and
-k,$!=Q=q,e-d’h(z-vt),
-k2$
=h(T,-
(3)
r=r,
T,),
r =r3,
(4)
(5)
where 6 is the Kronecker delta function. The matching conditions at the interface are T,=T,,
MATHEMATICAL
field in cylindrical coor-
temperature
dinates is
r=r2
(6)
ANALYSIS
kaT,=kaT, I ar 2 ar
Basic equations A finite-length composite hollow cylinder, which is composed of two homogeneous and isotropic materials with different properties, is considered. The physical model and coordinate system are shown in Fig. 1. The inner part of the cylinder is denoted as region 1 and the outer part as region 2. The material properties are assumed to be constant and perfect contact is maintained at the interface. Initially, the cylinder is assumed to be in a stress-free state and to have a uniform temperature T,. At time t = 0, an axisymmetrical line heat source, Q, moving with constant velocity v in the axial direction and decaying exponentially with time, acts at the inner surface of the cylinder. The heat is convectively transferred to the surroundings from the outer surface. Both ends of the cylinder are thermally insulated and fixed at the axial direction. Under these assumptions, the heat conduction equation for the
r=r2,
The generalized field equations of thermal stress distribution can be written in terms of displacement as 1 rz rr
v2uri--n.+--_ p2u,i+--
1 hi 1 - 2vi &
2(1 + vi)qdT, =o 1 - 2vi ar
I ae, 2(1 + ~,)a~ aT, = o ---I - 2vi az I - 2vi az ’
where v2=
a2 I a dr2+;ar+aZ’
e=!2+;+~
a2
(8)
(9)
Hybrid Laplace transform/finite
element method
1 -v, k, g z---s,] ” 1 + v, a,q,L
and i = 1,2. The initial, boundary and matching conditions of the thermal stress field are *wi
=
t=0
a,:, = ObU = 0,
urri =
855
_ 1 -v, ;G, =afl 1 +v, 2p,a,q0L
(10)
=---Q.,
z=O
u;, = u,* = 0,
and
z=L
‘I
(11) Laplace transform /finite element formulation
O,,l = a,, = 0, ~w2 =
%I =
g,,, =
42,
~,,z,
0,
t7,:2 =
%I =
~,I
=
r = r,
(12)
r =
(13)
42,
~,ZZ?
r,
r = r2
(14)
r =r,.
(1%
In order to solve this initial and boundary value problem, the method of Laplace transform is utilized. Writing the governing equations and conditions in terms of dimensionless quantities and taking the Laplace transform with respect to time, eqns (I)--(7) become
In order to simplify the problem, the following dimensionless variables and parameters are defined: R,
O
(16)
R,
O
(17)
Z=l)
(18)
R, = r,/L R,=r,/L R,=r,/L R = r/L
ag. -z&=0,
OnS,(Z=O
and
Z=z/L Bi = hL /k,
on S,(R = R,)
v = VL/K,
atT2 -dR = -B/B*, on S,(R = R,)
K=k,/k, ‘5=
8~8 I
K, t/L2
%=K!% on
2'aR
dR'
S,(R = R2).
(19)
(20)
(21)
D = dL=/K, The temperature distribution &(R, Z, s) can be determined inside the composite hollow cylinder which minimizes the functional
K=K,/K2
a = a, Ia2
7’cc)k,lq~L
ei=(Ti-
iz = Q/e P=
&dS+$
@dS,
(22)
P,lPz 1 -v,
BI=-I +v, 1 -v, “‘(1
+
fv,)a
“R,,LLLLLU
1 + v, a, q0L2 ‘.’
where w = 1 for region 1 and o = K for region 2. The above functional will satisfy the governing equations and boundary and matching conditions of eqns (16x21). The functional I can be written as the sum of the individual functionals defined for all elements of the assemblage, i.e.
(23)
L. S. CHENand H. S. CHU
856
where the superscript T represents the transpose of a matrix. Substituting eqn (27) into eqns (26) yields
[KlI& = IFI, +
@+‘dS s S2
+ 2
&)* dS.
(24)
(30)
where
s%
Suppose that the domain V is divided into E finite elements of p nodes each and that the distribution of 6 in element e is assumed to be &‘(R, 2, S) = f: N:” By’ = [N”‘]{&‘}, j-l
(25)
where 8y’ is the temperature of node j and NY’ the interpolation function corresponding to node j of element e. For the minimization of the functional Z, the necessary conditions given in eqn (26) should be satisfied. !&i,&
j=1,2...M,
(26)
where M is the total number of nodal temperature unknowns. Substituting eqn (25) into eqn (24), and then differentiating the functional Z with respect to 4 of each element, the following equations expressed by matrices can be obtained:
!!g=[K”‘]{&‘}
+ {F”‘},
The matrix [K] is an (M x M) band matrix with complex number, {&} is an (M x 1) vector representing the unknown temperatures and (F} is an (M x I) vector representing the unknown forcing terms. For calculating the thermal stresses distribution, the strain induced by the thermal expansion of the temperature variation can be regarded as the initial force to result in the stress distribution. The general finite element formulation for solving the thermal stress distribution due to temperature variations can be found from many references [14,15], and thus the functional which will satisfy the governing equation and boundary and match conditions of the thermal stress field for each element can be written as
(27)
.,=A2
J
where
{b”‘}~D$]({@}
- {@)})dV,
(31)
” where
[B(“]~D”‘][B”‘] dV
[K”‘] = bw
wRs [N'e)]Z[N(e'] dV
+ ml
+
,
s[N'p)]~N'p'] dS
(28)
ok’ = [N"']{ it$‘}’
(2%
1 - vj vi vi 1 - vi
and dNb” aR ahry,
az
[N”‘] = [NI”, N’,“.
. N”‘] P
’
(1 -2Vj) [@I = ~ ri
1
r,= 1, ’ {p,
1 1 fori=l fori=
vi vi
0 0 0
vi
v,
l-vi
0
0
0
1
1 - 2vi J
Hybrid Laplace transform/finite element method By using the same procedure as that for the temperature field, the formula of the finite element method for the thermal stress field can be summarized as (32)
857
where the chosen positive arbitrary constant, y, should be greater than the real parts of all sinaularities of J(s~. The numerical inversion of the Laplace transform can be written in the development of a Fourier transform by letting s = y + io; then eqn (34) becomes
where cb(r)=;
; rj(r + iw)e’““do. s m
(35)
The above integral can be approximated by taking a general one-dimensional Riemann sum and letting Ao = z/p. We then have f#J(t)1;
{fml=j,Jj
[B$J’&‘]{S6”}
Considering
dI’
Y(r)
k_g
cc
f#& + ikn/p)e’~(X”)‘.
(36)
the relation
f$(r + ikz/p)eik(n’p)r+ $(y - ikn/p)e-‘k(n’P)’
and
= 2Re[Q(y + ika/p)e’k(r’~“l (37)
[B;‘] =
84
dR’ N,
-1
R 0,
w
0,
Z’
0
N2
3
..
*
0,
-, R
0,
...
-5 aN,
0,
-g’ aN2
..*
aN,
aN,
alv,
alv,
XT’
dR’
dZ’
iii?’
..*
alv, 1
aR
N
2, R
0,
and truncating the infinite series by N terms, the numerical inversion formula can be expressed as
(3 1 0 aN
T$
’
aNpaNp az’ aR
where [K,,,] is a (2M x 2M) band matrix, {a,,,} is a (2M x 1) vector representing the unknown displacements and {fm} is a (2M x 1) vector representing the unknown forcing terms. Once the displacements {&,I are obtained, the stresses can be calculated using the following formula for every position within each element: {a(e)} = [D$]({P}
- {$)}),
(33)
From this equation, the numerical value of 4(t) at any time t in the interval 0 & t < 2p can be calculated. It is obvious from eqn (38) that the round-off and truncation errors can be made arbitrarily small if the product it is sufficiently large. However, if the value of yt is quite large, N should be increased accordingly in order to obtain a good numerical result. Therefore, an optimal value of the constant y should be chosen when t is within the interval 0-2~. At f =p, eqn (38) becomes
4@)
z:
_
~~(y)+Re~~,g(~+ikrlp)(--l)”
1,
(39)
where
from which the numerical result at the point p can be calculated. The present method can efficiently compute the temperatures and thermal stresses of the total nodes at a specific time by directly using the optimal value of yp.
{P’} = [B$‘]{&‘}. METHOD OF INVERSION
RESULTS AND DISCUSSION
After a Laplace transform function, d(r), has been obtained, the formula for inversion of the Laplace transform is defined as
4(t) =
s
&, y,‘% &s)es’ ds, i’ ICC
(34)
The finite element technique with a four-node rectangular element has been used to discretize the domain and the isoparametric interpolation function [NJ is used for both the temperature and thermal stress fields. The three-point Gauss quadrature integration formula is utilized to calculate the volume
858
L. S.
CHEN
and H. S. CHU
Table 1. Values for different element meshes in the Z-direction
10 20 30 50 70
1.171418 I .256332 1.344753 1.389504
1.388238
-0.108645 -0.130223 -0.143938 -0.166078 -0.164262
- 0.209348 -0.280972 -0.421900 -0.411316 -0.401433
-0.483090 -0.458018 -0.496581 -0.457327 -0.450191
-0.013552 -0.015889 -0.030492 -0.035005 -0.035388
237.304 468.916 700.340 1162.261 1625.148
The material properties are selected as K = V = Bi = 5.0 p = 2.0 at r = 0.1. integral in forming the stiffness matrices [K] and [K,,,] and the surface integral in forming the forcing matrices {F} and {f,). In this study, the dimensionless radii of the composite hollow cylinder, R, , R2 and R,, are selected as 0.04, 0.06 and 0.1, respectively. The thermal diffusivity ratio is chosen as K = 5.0, the thermal expansion coefficient ratio a is 2.0 and the Poisson ratios v, and v2 are selected as 0.15 and 0.25, respectively. The dimensionless decaying constant, D, is 0.01. The numerical results of the dimensionless temperature and thermal stresses for a number of different values of the dimensionless parameters are presented. In this study, the finite element meshes of the domain in the radical direction are chosen as six elements. In order to eveluate the accuracy and the convergence of this method, five different values of element meshes in the axial direction are used to carry out the calculation. The results of dimensionless temperature and thermal stresses at the mid-point of the interface are shown in Table 1. From this table, it can be seen that the computer time increases rapidly as the number of element meshes increase. Since both the values of temperature and thermal stresses converged completely at 50 element meshes in an axial direction, therefore only 50 element meshes in the axial direction are chosen in the following calculation. Figures 2 and 3 display the distribution of dimensionless temperature and hoop stress at the interface and outer surface for different values of thermal conductivity ratio. From these figures, it can be seen
that the highest values of the temperature and hoop stress are at exactly the same position as that where the line heat source is applied. It can also be seen that the temperature and hoop stress are increased when the thermal conductivity ratio, K, increases. This phenomenon is due to the fact that a large value of K means high values of the internal resistance of conduction heat flow in region 2, and thus the heat conduction from region 1 to the surface is impeded. The above situations cause large quantities of heat to be stored in the cylinder and consequently the cylinder temperature and hoop stress will be higher. Another interesting phenomenon is that the compressive hoop stress occurs at the interface and the tensile hoop stress occurs at the outer surface. The reason for this situation is that the thermal expansion of the inner part of the cylinder is constrained by the outer part and both ends of the cylinder, but that of the outer surface is constrained only by both ends of the cylinder. When the heat is applied at the inner surface, the hoop stress is compressive at the inner part of the cylinder but tensile at the outer surface. Figure 4 shows the radial stress distribution on the interface of the cylinder for Biot numbers, Bi equal to 0.1, 1.O, 5.0, and 50.0. The radial stress is compressive in the composite cylinder but on the inner and outer surface is 0, because both the surfaces are in a traction-free state. The effect of the Biot number on the axial stress is also shown in Fig. 5. The axial stress is compressive at the interface and tensile on the outer surface. This phenomenon is similar to the hoop stress distribution shown in Fig. 3. Examining both
e
Fig. 2. Temperature distributions for various values of thermal conductivity ratio, K.
Fig. 3. Hoop stress distribution for various values of thermal conductivity ratio, K.
Hybrid Laplace transform/finite element method
859
-0.08
aim -0.111
-0.15
-0.20
L
Fig. 4. Radial stress distribution for various values of Biot number, Bi.
Fig. 6. Axial stress distribution for various values of shear modulus ratio, p.
Figs 4 and 5, it can be Seen that the radial and axial stresses are increased as the Biot number, Bi, is
This is due to the insulation of the ends. The same situation for the hoop stress is shown in Fig. 9. The hoop stress distribution is also compressive at the interface and tensile on the outer surface.
decreased. Figures 6 and 7 depict variations of axial and shear stresses for different shear modulus ratios, )I. Both the axial and shear stresses will be higher when the shear modulus ratio is decreased. This means that the larger shear modulus at the inner region will prohibit the stress propagation and, therefore, the stresses at the interface and the outer surface will become smaller. The fact that the axial stress is compressive at the interface and is tensile on the outer surface can also be Seen from Fig. 6. Figure 7 shows that the position of the peak value of the shear stress is not the same as the position of the applied line heat source; instead, it occurs at both sides of the applied heat source. Figures 8 and 9 display the dimensionless temperature and hoop stress as a function of dimensionless axial position for several values of the dimensionless velocity of the moving heat source. The phenomenon that can be observed from Fig. 8 is that the peak value of temperature distribution is decreased as the velocity is increased, but at both ends of the hollow cylinder the higher peak value of temperature occurs.
Fig. 5. Axial stress distribution for various values of Biot number, Bi.
CONCLUSIONS
In this paper, a general Laplace transform/finite element technique is utilized to analyze the transient thermoelastic problem of a composite hollow cylinder. To illustrate the numerical accuracy, five different element meshes in the Z-direction have been analyzed. Based on the results, the conclusions drawn from this study can be summarized as follows. The positions of the peak values of temperature and thermal stresses (with the exception of shear stress) are at the same position as the line heat source. The axial and hoop thermal stresses are compressive in the inner part of the composite cylinder, but become tensile at the outer surface. The increase of thermal conductivity ratio, K,will cause the temperature and thermal stresses to increase.
Fig. 7. Shear stress distribution for various values of shear modulus ratio, 1.
860
L. S. CHENand H. S. CHU REFERENCES
2.6 2.0 6 l.6
0.6 0.0 0.0
0.2
0.4
0.6 a-#
0.6
LO
Fig. 8. Temperature distribution for various values of heat source velocity, V.
0.4
0.0 s -0.4
-0.6
1. M. P. Singh, R. A. Heller and S. Thangjitham, Thermal stresses in concentric cylinders due to asymmetric and time dependent temperature inputs. J. fherm. Stresses 7, 183-195 (1984). 2. Y. C. Yang, T. S. Wang and C. K. Chen, Thennoelastic transient response of an infinitely long annular cylinder. J. therm. Stresses 9, 19-30 (1986). 3. A. H. Ghosen and M. Sabbaghia, Quasi-static coupled problems of thermoelasticity for cylindrical regions. J. therm. Stresses 5, 299-311 (1982). 4. E. Pardo, G. S. Sarmiento, P. A. A. Laura and R. H. Gutierrez, Analytical solution for unsteady thermal stresses in an infinite cylinder composed of two materials. J. therm. Stresses 10, 29-43 (1987). 5. J. R. Thomas Jr, J. P. Singh, H. Tawil, L. Powers and D. P. H. Hasselman, Thermal stresses in a long circular cylinder subjected to sudden cooling during transient convection heating. J. therm. Stresses 8,249-260 (1985). 6. M. Gundappa, J. R. Thomas Jr and D. P. H. Hasselman, Analysis of thermal stresses in a solid circular cylinder subjected to conductive heat transfer with thermal boundary conductance. J. therm. Stresses 8, 265-276 (1985). 7. N. Noda, Transient thermoelastic contact problem in a long, circular cylinder. J. therm. Stresses 7, 135-147 (1984). 8. N. Noda, Transient thermoelastic constant stress in a short length circular cylinder. J. therm. Stresses 8, 413424 (1985). 9. Z.-B. Kuang and S. N. Atluri, Temperature field due to a moving heat source: a moving mesh finite element analysis. J. appl. Mech. 52, 274-280 (1985). 10. B. Dorri, V. Kadambi and F. W. Staub, Thermal stress analysis of sintering using a moving grid. Int. J. Numer. Meth. Engng 24, 47-57 (1987).
-1.2 0.0
0.2
0.4
0.6
0.e
1.0
z
Fig. 9. Hoop stress distribution for various values of heat
source velocity, V. 4. The temperature
and thermal stresses distribution is increased as the Biot number, Bi, is decreased. 5. The thermal stress distribution increases when the shear modulus ratio, g, decreases.
Il. T. C. Chen and C. I. Weng, Generalized coupled transient thermoelastic plane problems by Laplace transform/finite element method. J. appl. Mech. 55, 377-382 (1988). 12. H. S. Chu, C. K. Chen and C. I. Weng, Applications of Fourier series technique to transient heat transfer problem. Chem. Engng Commun. 16, 215-227 (1982). 13. H. S. Chu, C. I. Weng and C. K. Chen, Transient response of a composite straight fin. Trans. ASME J. Heat Transf. 105, 307-311 (1983). 14. S. S. Rao, The Finite Element Method in Engineering. Pergamon Press, New York (1982). 15. J. N. Reddy, An Inlroducrion IO the Finite Element Method. McGraw-Hill, New York (1984).
0045-7949/90 13.00 + o&l Q 1990 Pqamon Press plc
HYBRID LAPLACE TRANSFORM/FINITE ELEMENT METHOD FOR TRANSIENT THERMOELASTIC PROBLEM OF COMPOSITE HOLLOW CYLINDER L. S. CHEN and H. S. CHU Department
of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan 30049, Republic of China (Received 24 April 1989)
Abstract-The transient thermal stresses and temperature distribution of a composite hollow cylinder are investigated theoretically. The cylinder is composed of two different materials and heated by a moving line heat source on its inner boundary and cooled convectively on the exterior surface. The heat source is assumed to be axisymmetric, moving along the axis of the hollow cylinder with constant velocity, and decaying exponentially with time. The hybrid Laplace transform/finite element method is adopted for the solution of this problem. Finally, a numerical scheme, Fourier series technique, is utilized to calculate the inverse transform. The results show that, with the exception of the tangential thermal stress, the positions of the peak values of the temperature and thermal stresses are at exactly the same position as that at which the line heat source is applied. The effects of a number of parameters are also investigated.
%.Uz
NOTATION
PI, V%,l gradient matrices for temperature and stress
KI, K2
field, respectively specific heat c/J d exponential decaying coefficient of line heat source D dimensionless exponential decaying coefficient of line heat source material property matrices for temperature and PI, RI stress field, respectively element number e E total number of elements IFI, {fm} force vectors for temperature and stress field, respectively h heat transfer coefficient thermal conductivities of inner and outer k,,k, cylinder, respectively K thermal conductivity ratio of thP compo?te hollow cylinder stiffness matrices for temperature Jnd strc .s WI, Kl field, respectively total number of nodes interpolation function for both temperature and IMNI stress field number of nodes for each element P heat flux of the line moving heat source e cylindrical coordinates dimensionless inner radius of composite hello :z cylinder dimensionless interface radius of composite R? hollow cylinder dimensionless outer radius of composite hollow R, cylinder dimensionless cylindrical coordinates R, Z s the variable used in the Laplace transform t time T temperature ambient temperature T, 1’ velocity of heat source V dimensionless velocity of heat source
K e ? 111.& p “1, “2
thermal expansion coefficients of inner and outer cylinder, respectively thermal diffusivities of inner and outer cylinder, respectively thermal diffusivity ratio of composite hollow cylinder dimensionless temperature dimensionless time shear moduli of inner and outer cylinder, respectively shear modulus ratio of composite hollow cylinder Poisson’s ratios of inner and outer cylinder, respectively stress field mass density
INTRODUCIION
The problem of temperature and thermal stress distribution in a cylinder due to a time dependent heat source has many engineering applications, for example solid propellants, rocket motors, gun barrels, nuclear reactor elements, etc. and has thus been the subject of numerous investigations [l]. The unsteady state thermoelastic problems of a cylinder have also been studied intensively, for example the one dimensional quasi-static coupled thermoelastic ,>roblem of an axisymmetric infinitely long cylinder [2,3]. The transient uncoupled thermal stress distriTVrtion has also been investigated in the axisymmetric caoling case [4], in the cooling problem during tran.ient convection heating [S] and in a cylinder subJetted to the conductive heat transfer with thermal boundary conductance [6]. The transient thermoelastic contact problem of an infinitely long solid cylinder
853
L. S. C&EN and H. S. CHU
854
perfectly bonded by a heated rigid circular ring has been analyzed with homogeneous materials [7,8]. The above studies were all based on the analytical models. The unsteady-state problems have also been solved by finite element methods; for instance, the temperature distribution due to a moving heat source has been solved by the moving mesh finite element method [9], and the thermal stress of sintering has been studied by the moving grid method [lo]. Finally, since the use of the finite element method is time consuming for a large problem, the hybrid Laplace transform/finite element method has been used [ll]. However, most of the previous studies deal only with cylinders with either infinite length or homogenous material, while few of them deal with cylinders with finite length and composite materials. In order to study the behavior of transient thermal stresses and the temperature distribution of a composite hollow cylinder with finite length, the hybrid Laplace transform/finite element method is adopted. The cylinder is composed of two different materials and is heated by a moving line heat source on its inner boundary and cooled convectively on the exterior surface. The heat source is assumed to be axisymmetric, moving along the axis of the hollow cylinder with a constant velocity and decaying exponentially with time. In this study, the coupling of temperature and strain rate as well as the inertia effect in the stress field are ignored. The temperature and thermal stress fields are obtained by applying Laplace transform with respect to the time variable and then the finite element method is used to solve the problem in the transformed region. Then, the inverse Laplace transform is accomplished by a numerical method, the Fourier series technique [12, 131. Numerical results for both temperature and thermal stresses are presented graphically with respect to a number of parameters.
r
FLFig. I. Physical model and coordinate system of composite hollow cylinder. axisymmetrical
a2Ti
(
a2Ti
1 ar,
F+;yjy’g
)
1 aTi =K_dtv f
rl
o
(I)
where i = 1,2, which refer to the regions I and 2, respectively. The initial and boundary conditions for the temperature field are Ti=T,, aT. -=O,
z=O
az
t=O
(2) z=L
and
-k,$!=Q=q,e-d’h(z-vt),
-k2$
=h(T,-
(3)
r=r,
T,),
r =r3,
(4)
(5)
where 6 is the Kronecker delta function. The matching conditions at the interface are T,=T,,
MATHEMATICAL
field in cylindrical coor-
temperature
dinates is
r=r2
(6)
ANALYSIS
kaT,=kaT, I ar 2 ar
Basic equations A finite-length composite hollow cylinder, which is composed of two homogeneous and isotropic materials with different properties, is considered. The physical model and coordinate system are shown in Fig. 1. The inner part of the cylinder is denoted as region 1 and the outer part as region 2. The material properties are assumed to be constant and perfect contact is maintained at the interface. Initially, the cylinder is assumed to be in a stress-free state and to have a uniform temperature T,. At time t = 0, an axisymmetrical line heat source, Q, moving with constant velocity v in the axial direction and decaying exponentially with time, acts at the inner surface of the cylinder. The heat is convectively transferred to the surroundings from the outer surface. Both ends of the cylinder are thermally insulated and fixed at the axial direction. Under these assumptions, the heat conduction equation for the
r=r2,
The generalized field equations of thermal stress distribution can be written in terms of displacement as 1 rz rr
v2uri--n.+--_ p2u,i+--
1 hi 1 - 2vi &
2(1 + vi)qdT, =o 1 - 2vi ar
I ae, 2(1 + ~,)a~ aT, = o ---I - 2vi az I - 2vi az ’
where v2=
a2 I a dr2+;ar+aZ’
e=!2+;+~
a2
(8)
(9)
Hybrid Laplace transform/finite
element method
1 -v, k, g z---s,] ” 1 + v, a,q,L
and i = 1,2. The initial, boundary and matching conditions of the thermal stress field are *wi
=
t=0
a,:, = ObU = 0,
urri =
855
_ 1 -v, ;G, =afl 1 +v, 2p,a,q0L
(10)
=---Q.,
z=O
u;, = u,* = 0,
and
z=L
‘I
(11) Laplace transform /finite element formulation
O,,l = a,, = 0, ~w2 =
%I =
g,,, =
42,
~,,z,
0,
t7,:2 =
%I =
~,I
=
r = r,
(12)
r =
(13)
42,
~,ZZ?
r,
r = r2
(14)
r =r,.
(1%
In order to solve this initial and boundary value problem, the method of Laplace transform is utilized. Writing the governing equations and conditions in terms of dimensionless quantities and taking the Laplace transform with respect to time, eqns (I)--(7) become
In order to simplify the problem, the following dimensionless variables and parameters are defined: R,
O
(16)
R,
O
(17)
Z=l)
(18)
R, = r,/L R,=r,/L R,=r,/L R = r/L
ag. -z&=0,
OnS,(Z=O
and
Z=z/L Bi = hL /k,
on S,(R = R,)
v = VL/K,
atT2 -dR = -B/B*, on S,(R = R,)
K=k,/k, ‘5=
8~8 I
K, t/L2
%=K!% on
2'aR
dR'
S,(R = R2).
(19)
(20)
(21)
D = dL=/K, The temperature distribution &(R, Z, s) can be determined inside the composite hollow cylinder which minimizes the functional
K=K,/K2
a = a, Ia2
7’cc)k,lq~L
ei=(Ti-
iz = Q/e P=
&dS+$
@dS,
(22)
P,lPz 1 -v,
BI=-I +v, 1 -v, “‘(1
+
fv,)a
“R,,LLLLLU
1 + v, a, q0L2 ‘.’
where w = 1 for region 1 and o = K for region 2. The above functional will satisfy the governing equations and boundary and matching conditions of eqns (16x21). The functional I can be written as the sum of the individual functionals defined for all elements of the assemblage, i.e.
(23)
L. S. CHENand H. S. CHU
856
where the superscript T represents the transpose of a matrix. Substituting eqn (27) into eqns (26) yields
[KlI& = IFI, +
@+‘dS s S2
+ 2
&)* dS.
(24)
(30)
where
s%
Suppose that the domain V is divided into E finite elements of p nodes each and that the distribution of 6 in element e is assumed to be &‘(R, 2, S) = f: N:” By’ = [N”‘]{&‘}, j-l
(25)
where 8y’ is the temperature of node j and NY’ the interpolation function corresponding to node j of element e. For the minimization of the functional Z, the necessary conditions given in eqn (26) should be satisfied. !&i,&
j=1,2...M,
(26)
where M is the total number of nodal temperature unknowns. Substituting eqn (25) into eqn (24), and then differentiating the functional Z with respect to 4 of each element, the following equations expressed by matrices can be obtained:
!!g=[K”‘]{&‘}
+ {F”‘},
The matrix [K] is an (M x M) band matrix with complex number, {&} is an (M x 1) vector representing the unknown temperatures and (F} is an (M x I) vector representing the unknown forcing terms. For calculating the thermal stresses distribution, the strain induced by the thermal expansion of the temperature variation can be regarded as the initial force to result in the stress distribution. The general finite element formulation for solving the thermal stress distribution due to temperature variations can be found from many references [14,15], and thus the functional which will satisfy the governing equation and boundary and match conditions of the thermal stress field for each element can be written as
(27)
.,=A2
J
where
{b”‘}~D$]({@}
- {@)})dV,
(31)
” where
[B(“]~D”‘][B”‘] dV
[K”‘] = bw
wRs [N'e)]Z[N(e'] dV
+ ml
+
,
s[N'p)]~N'p'] dS
(28)
ok’ = [N"']{ it$‘}’
(2%
1 - vj vi vi 1 - vi
and dNb” aR ahry,
az
[N”‘] = [NI”, N’,“.
. N”‘] P
’
(1 -2Vj) [@I = ~ ri
1
r,= 1, ’ {p,
1 1 fori=l fori=
vi vi
0 0 0
vi
v,
l-vi
0
0
0
1
1 - 2vi J
Hybrid Laplace transform/finite element method By using the same procedure as that for the temperature field, the formula of the finite element method for the thermal stress field can be summarized as (32)
857
where the chosen positive arbitrary constant, y, should be greater than the real parts of all sinaularities of J(s~. The numerical inversion of the Laplace transform can be written in the development of a Fourier transform by letting s = y + io; then eqn (34) becomes
where cb(r)=;
; rj(r + iw)e’““do. s m
(35)
The above integral can be approximated by taking a general one-dimensional Riemann sum and letting Ao = z/p. We then have f#J(t)1;
{fml=j,Jj
[B$J’&‘]{S6”}
Considering
dI’
Y(r)
k_g
cc
f#& + ikn/p)e’~(X”)‘.
(36)
the relation
f$(r + ikz/p)eik(n’p)r+ $(y - ikn/p)e-‘k(n’P)’
and
= 2Re[Q(y + ika/p)e’k(r’~“l (37)
[B;‘] =
84
dR’ N,
-1
R 0,
w
0,
Z’
0
N2
3
..
*
0,
-, R
0,
...
-5 aN,
0,
-g’ aN2
..*
aN,
aN,
alv,
alv,
XT’
dR’
dZ’
iii?’
..*
alv, 1
aR
N
2, R
0,
and truncating the infinite series by N terms, the numerical inversion formula can be expressed as
(3 1 0 aN
T$
’
aNpaNp az’ aR
where [K,,,] is a (2M x 2M) band matrix, {a,,,} is a (2M x 1) vector representing the unknown displacements and {fm} is a (2M x 1) vector representing the unknown forcing terms. Once the displacements {&,I are obtained, the stresses can be calculated using the following formula for every position within each element: {a(e)} = [D$]({P}
- {$)}),
(33)
From this equation, the numerical value of 4(t) at any time t in the interval 0 & t < 2p can be calculated. It is obvious from eqn (38) that the round-off and truncation errors can be made arbitrarily small if the product it is sufficiently large. However, if the value of yt is quite large, N should be increased accordingly in order to obtain a good numerical result. Therefore, an optimal value of the constant y should be chosen when t is within the interval 0-2~. At f =p, eqn (38) becomes
4@)
z:
_
~~(y)+Re~~,g(~+ikrlp)(--l)”
1,
(39)
where
from which the numerical result at the point p can be calculated. The present method can efficiently compute the temperatures and thermal stresses of the total nodes at a specific time by directly using the optimal value of yp.
{P’} = [B$‘]{&‘}. METHOD OF INVERSION
RESULTS AND DISCUSSION
After a Laplace transform function, d(r), has been obtained, the formula for inversion of the Laplace transform is defined as
4(t) =
s
&, y,‘% &s)es’ ds, i’ ICC
(34)
The finite element technique with a four-node rectangular element has been used to discretize the domain and the isoparametric interpolation function [NJ is used for both the temperature and thermal stress fields. The three-point Gauss quadrature integration formula is utilized to calculate the volume
858
L. S.
CHEN
and H. S. CHU
Table 1. Values for different element meshes in the Z-direction
10 20 30 50 70
1.171418 I .256332 1.344753 1.389504
1.388238
-0.108645 -0.130223 -0.143938 -0.166078 -0.164262
- 0.209348 -0.280972 -0.421900 -0.411316 -0.401433
-0.483090 -0.458018 -0.496581 -0.457327 -0.450191
-0.013552 -0.015889 -0.030492 -0.035005 -0.035388
237.304 468.916 700.340 1162.261 1625.148
The material properties are selected as K = V = Bi = 5.0 p = 2.0 at r = 0.1. integral in forming the stiffness matrices [K] and [K,,,] and the surface integral in forming the forcing matrices {F} and {f,). In this study, the dimensionless radii of the composite hollow cylinder, R, , R2 and R,, are selected as 0.04, 0.06 and 0.1, respectively. The thermal diffusivity ratio is chosen as K = 5.0, the thermal expansion coefficient ratio a is 2.0 and the Poisson ratios v, and v2 are selected as 0.15 and 0.25, respectively. The dimensionless decaying constant, D, is 0.01. The numerical results of the dimensionless temperature and thermal stresses for a number of different values of the dimensionless parameters are presented. In this study, the finite element meshes of the domain in the radical direction are chosen as six elements. In order to eveluate the accuracy and the convergence of this method, five different values of element meshes in the axial direction are used to carry out the calculation. The results of dimensionless temperature and thermal stresses at the mid-point of the interface are shown in Table 1. From this table, it can be seen that the computer time increases rapidly as the number of element meshes increase. Since both the values of temperature and thermal stresses converged completely at 50 element meshes in an axial direction, therefore only 50 element meshes in the axial direction are chosen in the following calculation. Figures 2 and 3 display the distribution of dimensionless temperature and hoop stress at the interface and outer surface for different values of thermal conductivity ratio. From these figures, it can be seen
that the highest values of the temperature and hoop stress are at exactly the same position as that where the line heat source is applied. It can also be seen that the temperature and hoop stress are increased when the thermal conductivity ratio, K, increases. This phenomenon is due to the fact that a large value of K means high values of the internal resistance of conduction heat flow in region 2, and thus the heat conduction from region 1 to the surface is impeded. The above situations cause large quantities of heat to be stored in the cylinder and consequently the cylinder temperature and hoop stress will be higher. Another interesting phenomenon is that the compressive hoop stress occurs at the interface and the tensile hoop stress occurs at the outer surface. The reason for this situation is that the thermal expansion of the inner part of the cylinder is constrained by the outer part and both ends of the cylinder, but that of the outer surface is constrained only by both ends of the cylinder. When the heat is applied at the inner surface, the hoop stress is compressive at the inner part of the cylinder but tensile at the outer surface. Figure 4 shows the radial stress distribution on the interface of the cylinder for Biot numbers, Bi equal to 0.1, 1.O, 5.0, and 50.0. The radial stress is compressive in the composite cylinder but on the inner and outer surface is 0, because both the surfaces are in a traction-free state. The effect of the Biot number on the axial stress is also shown in Fig. 5. The axial stress is compressive at the interface and tensile on the outer surface. This phenomenon is similar to the hoop stress distribution shown in Fig. 3. Examining both
e
Fig. 2. Temperature distributions for various values of thermal conductivity ratio, K.
Fig. 3. Hoop stress distribution for various values of thermal conductivity ratio, K.
Hybrid Laplace transform/finite element method
859
-0.08
aim -0.111
-0.15
-0.20
L
Fig. 4. Radial stress distribution for various values of Biot number, Bi.
Fig. 6. Axial stress distribution for various values of shear modulus ratio, p.
Figs 4 and 5, it can be Seen that the radial and axial stresses are increased as the Biot number, Bi, is
This is due to the insulation of the ends. The same situation for the hoop stress is shown in Fig. 9. The hoop stress distribution is also compressive at the interface and tensile on the outer surface.
decreased. Figures 6 and 7 depict variations of axial and shear stresses for different shear modulus ratios, )I. Both the axial and shear stresses will be higher when the shear modulus ratio is decreased. This means that the larger shear modulus at the inner region will prohibit the stress propagation and, therefore, the stresses at the interface and the outer surface will become smaller. The fact that the axial stress is compressive at the interface and is tensile on the outer surface can also be Seen from Fig. 6. Figure 7 shows that the position of the peak value of the shear stress is not the same as the position of the applied line heat source; instead, it occurs at both sides of the applied heat source. Figures 8 and 9 display the dimensionless temperature and hoop stress as a function of dimensionless axial position for several values of the dimensionless velocity of the moving heat source. The phenomenon that can be observed from Fig. 8 is that the peak value of temperature distribution is decreased as the velocity is increased, but at both ends of the hollow cylinder the higher peak value of temperature occurs.
Fig. 5. Axial stress distribution for various values of Biot number, Bi.
CONCLUSIONS
In this paper, a general Laplace transform/finite element technique is utilized to analyze the transient thermoelastic problem of a composite hollow cylinder. To illustrate the numerical accuracy, five different element meshes in the Z-direction have been analyzed. Based on the results, the conclusions drawn from this study can be summarized as follows. The positions of the peak values of temperature and thermal stresses (with the exception of shear stress) are at the same position as the line heat source. The axial and hoop thermal stresses are compressive in the inner part of the composite cylinder, but become tensile at the outer surface. The increase of thermal conductivity ratio, K,will cause the temperature and thermal stresses to increase.
Fig. 7. Shear stress distribution for various values of shear modulus ratio, 1.
860
L. S. CHENand H. S. CHU REFERENCES
2.6 2.0 6 l.6
0.6 0.0 0.0
0.2
0.4
0.6 a-#
0.6
LO
Fig. 8. Temperature distribution for various values of heat source velocity, V.
0.4
0.0 s -0.4
-0.6
1. M. P. Singh, R. A. Heller and S. Thangjitham, Thermal stresses in concentric cylinders due to asymmetric and time dependent temperature inputs. J. fherm. Stresses 7, 183-195 (1984). 2. Y. C. Yang, T. S. Wang and C. K. Chen, Thennoelastic transient response of an infinitely long annular cylinder. J. therm. Stresses 9, 19-30 (1986). 3. A. H. Ghosen and M. Sabbaghia, Quasi-static coupled problems of thermoelasticity for cylindrical regions. J. therm. Stresses 5, 299-311 (1982). 4. E. Pardo, G. S. Sarmiento, P. A. A. Laura and R. H. Gutierrez, Analytical solution for unsteady thermal stresses in an infinite cylinder composed of two materials. J. therm. Stresses 10, 29-43 (1987). 5. J. R. Thomas Jr, J. P. Singh, H. Tawil, L. Powers and D. P. H. Hasselman, Thermal stresses in a long circular cylinder subjected to sudden cooling during transient convection heating. J. therm. Stresses 8,249-260 (1985). 6. M. Gundappa, J. R. Thomas Jr and D. P. H. Hasselman, Analysis of thermal stresses in a solid circular cylinder subjected to conductive heat transfer with thermal boundary conductance. J. therm. Stresses 8, 265-276 (1985). 7. N. Noda, Transient thermoelastic contact problem in a long, circular cylinder. J. therm. Stresses 7, 135-147 (1984). 8. N. Noda, Transient thermoelastic constant stress in a short length circular cylinder. J. therm. Stresses 8, 413424 (1985). 9. Z.-B. Kuang and S. N. Atluri, Temperature field due to a moving heat source: a moving mesh finite element analysis. J. appl. Mech. 52, 274-280 (1985). 10. B. Dorri, V. Kadambi and F. W. Staub, Thermal stress analysis of sintering using a moving grid. Int. J. Numer. Meth. Engng 24, 47-57 (1987).
-1.2 0.0
0.2
0.4
0.6
0.e
1.0
z
Fig. 9. Hoop stress distribution for various values of heat
source velocity, V. 4. The temperature
and thermal stresses distribution is increased as the Biot number, Bi, is decreased. 5. The thermal stress distribution increases when the shear modulus ratio, g, decreases.
Il. T. C. Chen and C. I. Weng, Generalized coupled transient thermoelastic plane problems by Laplace transform/finite element method. J. appl. Mech. 55, 377-382 (1988). 12. H. S. Chu, C. K. Chen and C. I. Weng, Applications of Fourier series technique to transient heat transfer problem. Chem. Engng Commun. 16, 215-227 (1982). 13. H. S. Chu, C. I. Weng and C. K. Chen, Transient response of a composite straight fin. Trans. ASME J. Heat Transf. 105, 307-311 (1983). 14. S. S. Rao, The Finite Element Method in Engineering. Pergamon Press, New York (1982). 15. J. N. Reddy, An Inlroducrion IO the Finite Element Method. McGraw-Hill, New York (1984).