- Email: [email protected]
Transient heat flow analysis using the fundamental collocation method
Applied Thermal Engineering 23 (2003) 893–904 www.elsevier.com/locate/apthermeng
Transient heat flow analysis using the fundamental collocation method Gary Burgess b
a,*
, Enayat Mahajerin
b
a School of Packaging, Michigan State University, East Lansing, Michigan, MI 48824, USA Department of Mechanical Engineering, Saginaw Valley State University, University Center, Michigan, MI 48710, USA
Received 2 March 2001; accepted 13 January 2003
Abstract The fundamental collocation method for general two dimensional transient heat flow problems in solids is presented. The numerical procedure is applied directly to heat conduction in the time domain and is illustrated using several example problems. The method is capable of handling problems of arbitrary shapes subjected to arbitrary initial conditions and mixed time-dependent boundary conditions. Its simplicity and minimal computational effort combined with accuracy and stability of the results make the method an excellent alternate approach to domain methods such as the finite difference and finite element methods. Ó 2003 Elsevier Science Ltd. All rights reserved. Keywords: Fundamental collocation method; Finite difference; Finite element; Transient heat flow
1. Introduction The general problem of transient heat flow in two dimensions involves conduction heat transfer in regions of arbitrary shape where the temperature is a function of both time and position. Such problems generally occur before steady state operating conditions are reached. The time duration of the transient part can be of importance in the design of machine parts such as those in turbines and auto engines. Analytical solutions may be found for simple geometries and initial/boundary conditions. For the general case, a numerical solution is preferred. Domain type numerical methods such as finite difference and finite element methods are not necessarily the best approaches due to setup difficulties for arbitrary regions, accuracy, and the fact that the solution is
*
Corresponding author. Tel.: +1-517-355-9504; fax: +1-517-535-8999. E-mail address: [email protected] (G. Burgess).
1359-4311/03/$ - see front matter Ó 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S1359-4311(03)00026-7
894
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
available only at the grid or the element points. In this paper, a much simpler approach based on the fundamental collocation method (FCM) [1] is presented. It can handle regions of arbitrary shapes subjected to general initial and boundary conditions and gives the solution anywhere inside the region and on the boundary. The FCM algorithm presented here is applied directly in the time domain where a GreenÕs function solution to the governing partial differential equation is readily available.
2. Problem statement Consider a general two-dimensional region R with its boundary oR (Fig. 1). We are interested in the time-dependent temperature field, T ðx; y; tÞ. The governing equation for unsteady Fourier heat conduction in a homogeneous isotropic material with no internal heat generation is well-known [2]: oT ¼ ar2 T ðx; y; tÞ ot
ð1Þ
for x; y in R
The initial temperature distribution and the boundary conditions are taken to be T ðx; y; 0Þ ¼ gðx; yÞ aT ðxB ; yB ; tÞ þ b
ð2Þ
in R
oT ðxB ; yB ; tÞ ¼c on
on oR
ð3Þ
where a, b, c, specified functions of position and time; a, thermal diffusivity of the material, assumed to be constant; r2 , o2 =ox2 þ o2 =oy 2 is the Laplace operator; gðx; yÞ, specified initial temperature distribution; xB , yB , boundary point; n, outward normal to the boundary; oT =on ¼ ðoT =oxÞ cos h þ ðoT =oyÞ sin h; h, the angle that outward normal to the boundary makes with the +x axis.
Fig. 1. The problem of interest.
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
895
3. Numerical analysis Steady state heat conduction problems have been handled successfully using the FCM with excellent results [3–5]. The addition of a new variable, time, does not add any real complexity to the FCM algorithm as it does with most other numerical methods. Analytical methods such as separation of variables tend to be so complicated that closed form (including series) solutions are nearly impossible, especially when the region under consideration is arbitrary and the boundary conditions are nonhomogeneous and possibly time-dependent. Lumped capacity problems are about the only problems that can be handled analytically for general regions. Graphical methods like the Schmidt Plots and Heisler Charts [6] are based on the finite difference approach and are only rough approximations. Numerical methods for the unsteady case usually produce sequential solutions over time in which the next generation result depends on the entire history of previous generation results. Numerical stability of the solution is often questionable due to the accumulation of round off error. Domain approaches such as the finite difference and finite element methods have been successfully applied to unsteady cases [7,8]. However, these methods tend to be complicated and the accuracy of the results for regions of arbitrary shapes is questionable. Boundary approaches such as the boundary integral equation method and the boundary element method can significantly reduce the time and effort spent on the computation because they reduce the dimensionality of the problem by one [9]. However, these methods require special treatment when the solution point is near the boundary of the region [10,11]. A much simpler boundary-oriented method called the FCM [1,3,4] significantly reduces computational effort and usually eliminates difficulties with the solution point being close to the boundary. The method consists of embedding the region of interest inside an infinite plane of the same material. A finite number of ‘‘heat sources’’ are placed anywhere outside the boundary of the region, oR, as shown in Fig. 2. The simplest way to do this is to locate the source points on a congruent boundary, oB, at a distance S away from the actual boundary. In this case, NB boundary points are selected for the satisfaction of prescribed boundary conditions by collocation and NB sources are placed along the normals to these points a distance S away. At regular time intervals Dt, each source j, j ¼ 1, NB, deposits a momentary blast of heat over an arbitrary but small patch of area DAj centered on the source point. This blast of heat is sufficient to raise the temperature of the patch from zero to Tj , where Tj is the jth ‘‘source temperature’’ to be determined by the FCM. If this is allowed to continue for ‘‘M’’ time steps, then the temperature due to all outside sources at location ðx; yÞ inside a region initially at zero temperature at time t ¼ M Dt is 3 2 ðxxj Þ2 þðyyj Þ2 NB X M1 exp 4aðt‘DtÞ X 5DAj Tj ð‘DtÞ4 ð4Þ T ðx; y; tÞ ¼ 4paðt ‘DtÞ j¼1 ‘¼0 The quantity in square brackets in Eq. (4) is a modified version of the GreenÕs function for a temperature spike at a point in an infinite medium [2]. This quantity is dimensionless and can be thought of as the fraction of the source temperature Tj imposed at time. ‘ Dt at source location ðxj ; yj Þ which is produced at the field location ðx; yÞ at time t ¼ M Dt, M > ‘. This fraction
896
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
Fig. 2. The numerical setup for FCM.
decreases with distance from the source and with the lag time ðM ‘Þ Dt. Eq. (4) satisfies Eq. (1) exactly for all ðx; y; tÞ. For simplicity, each of the patch areas DAj can be taken as unity. It can be shown [2] that an initial temperature distribution is equivalent to an infinite number of infinitesimal internal heat sources which release momentary blasts of heat at time t ¼ 0 only, and in amounts sufficient to raise their temperatures to the specified initial conditions. Subdividing the region into internal patches of area DAk centered on field points ðxk ; yk Þ, k ¼ 1, NF, and applying the same reasoning as in Eq. (4), the FCM solution is the sum of contributions from NF known internal sources and NB unknown external ones: 3 2 þðyyk Þ2 exp ðxxk Þ 4at 5DAk T ðx; y; tÞ ¼ gðxk ; yk Þ4 4pat k¼1 3 2 ðxxj Þ2 þðyyj Þ2 NB X M1 exp 4aðt‘DtÞ X 5 þ Tj ð‘DtÞ4 4paðt ‘DtÞ j¼1 ‘¼0 NF X
2
ð5Þ
where t ¼ M Dt, NB is the number of heat sources outside the boundary and NF is the number of heat sources inside the region representing the initial temperature distribution. Eq. (5) satisfies Eq. (1) everywhere inside the region, and Eq. (2) at the ðxk ; yk Þ. All that remains is to satisfy Eq. (3) at the NB collocation points on the boundary. Substituting Eq. (5) into Eq. (3) and splitting off the last term in the sum over time requires that
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904 NB X
897
Dij Tj ðt DtÞ
j¼1
2 3 exp ðxi xk Þ2 þðyi yk Þ2 4at ðxi xk Þ cos hi þ ðyi yk Þ sin hi 4 ¼ ci DAk 5 gðxk ; yk Þ ai bi 2at 4pat k¼1 3 2 exp ðxi xj Þ2 þðyi yj Þ2 NB X M 2 X 4aðt‘DtÞ ðxi xj Þ cos hi þ ðyi yj Þ sin hi 4 5 ð6Þ Tj ð‘DtÞ ai bi 2aðt ‘DtÞ 4paðt ‘DtÞ j¼1 ‘¼0 NF X
for each i ¼ 1 to NB at each time t ¼ M Dt, M ¼ 1; 2; 3; . . ., where D is the NBNB ‘‘influence matrix’’ with 3 2 exp ðxi xj Þ2 þðyi yj Þ2 4aDt ðxi xj Þ cos hi þ ðyi yj Þ sin hi 4 5 Dij ¼ ai bi ð7Þ 2aDt 4paDt If the current time at which boundary conditions are to be satisfied is t ¼ M Dt, then Eq. (6) gives an NBNB system of linear equations in the NB outside source temperatures Tj at the previous time step t Dt in terms of the time steps. ‘ Dt, ‘ ¼ 0 to M 2 leading up to this. Beginning with M ¼ 1, the source temperatures Tj at time ðM 1Þ Dt are found from Eq. (6) and stored in the Mth column of a solution matrix which contains the whole history of the source temperatures up to the present time M Dt. Since an NBNB system of equations must be solved at each time step, and since the influence matrix is independent of time, this process is best handled by inverting D. Once the source temperatures are known at a particular time, Eq. (5) and its derivatives are used to find the temperatures and fluxes at any point on the boundary or inside the region.
4. Stability and steady state For the case where temperature is specified everywhere on the boundary, ai ¼ 1, bi ¼ 0, and ci ¼ specified temperature. The influence matrix in Eq. (7) consists of terms of the form expðk=DtÞ=Dt, which can be shown to reach a maximum value when Dt ¼ k. In order for the numerical inversion of the influence matrix D to be reasonably accurate, the diagonal terms of D should be strong. Because the diagonal terms represent the influence of an outside source on its closest boundary point a distance S away, the diagonal elements can be made to be larger than any other element in the same row by choosing S2 ð8Þ 4a Eq. (8) gives relationship between the chosen time step and the distance S between the source and the region boundaries which should generate a stable solution. In practice, fairly large deviations from this can be tolerated as long as Dt is larger than the value in Eq. (8). Smaller DtÕs usually result in inaccurate transient responses. Physically, Eq. (8) says that the time step must be sufficiently long to allow a blast of heat at distance S away to reach the boundary. For general Dt ¼
898
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
boundary conditions ai 6¼ 0, bi 6¼ 0, it is not clear how to choose the time step in Eq. (7) to make the diagonal of D strong. In this case, it is best to choose Dt as in Eq. (8), since the physical interpretation still applies to general boundary conditions. The choice of the distance S is somewhat arbitrary. In practice, S is usually chosen to be somewhere between the smallest and the largest dimensions of the region being investigated. The largest dimension is defined to be the diameter of the smallest circle into which the region can be embedded. As S increases, the determinant of the influence matrix rapidly decreases, so there is a practical limit to the size of S. This limit is largely dependent on machine precision. Experience shows that the best way to choose S is to start with a small value and gradually increase it until the sum of the squares of the errors (SSE) between the prescribed boundary conditions and the delivered boundary conditions is a minimum, 2 NB
X oT cj aj Tj bj ¼ minimum ð9Þ SSE ¼ onj j¼1 where aj , bj and cj are the specified parameters in Eq. (3) and Tj and oT =onj are the values delivered by the FCM solution (Eq. (5) and its normal derivative). In the SSE, it is a good idea to include boundary points in-between the ones chosen for the collocation process. The resulting SSE will then include both roundoff errors from the collocation process and oscillation errors, if any, due to the method itself. As time gets large, the effect of the initial temperature distribution in Eq. (5) dies off. If the boundary conditions do not vary over time, a steady state solution exists and the source temperatures Tj ð‘DtÞ are expected to approach constant values Tj . When this happens, Eq. (5) takes the form 3 2 ðxxj Þ2 þðyyj Þ2 t=Dt exp NB 4aðt‘DtÞ X X 5 as t ! 1 Tj 4 T ðx; y; tÞ ¼ ð10Þ 4paðt ‘DtÞ j¼1 ‘¼0 The second sum in Eq. (10) can be evaluated as an integral having the form Z 1 expðk=tÞ dt IðkÞ ¼ t 0 To perform the integration, first differentiate with respect to k, Z 1 Z dI expðk=tÞ 1 1 d 1 ¼ ðexpðk=tÞÞ dt ¼ dt ¼ 2 dk t k 0 dt k 0
ð11Þ
ð12Þ
and integrate to get I back, I ¼ ‘nk þ C
ð13Þ
When applied to Eq. (10), the steady state solution has the following FCM representation: T ðx; yÞ ¼
NB X
Tj ‘n½ðx xj Þ2 þ ðy yj Þ2
j¼1
which is identical to the FCM solution of the Laplace equation used in [1].
ð14Þ
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
899
5. Examples The method was applied to several examples. Examples were chosen to illustrate performance for different geometries, initial conditions, and boundary conditions. In order to demonstrate the accuracy of the method, a comparison to the exact solution was made when available. Example 1. Consider a rectangular region AB made of bronze (a ¼ 14:1 mm2 /s) subjected to the initial and boundary conditions shown in Fig. 3. The exact solution can be obtained by separation of variables,
T ðx; y; tÞ ¼
2 2
1 X 1 X 200 ipx jpy ip j2 p2 þ 2 t Cij sin xy þ sin exp a AB A B A2 B i¼1 j¼1
ð15Þ
where Cij ¼
4 ½100 ð1Þiþj 200 2 ijp
ð16Þ
The FCM was applied to this problem using A ¼ 200 mm, B ¼ 100 mm. The boundary was represented by 40 equally spaced collocation points, 10 per side, and the external heat sources used to satisfy boundary conditions were placed on a congruent boundary located at a distance S ¼ 10 mm away. The interior of the region was divided into 100 rectangular patches and the temperature distribution was represented by internal heat sources located at the centers of these patches. A time increment of 10 s was used to compute the temperatures at five equally spaced points along a diagonal including the corners. The minimum recommended time step from Eq. (8) is Dt ¼ 1:8 s. The results for the temperature are shown in Table 1. Changing the distance S between the source boundary and the boundary of the region affects these results very little. The practice of changing S is recommended as a check on stability; the results should be the same over a wide range of values. In rare cases, very large or very small S will affect the solution at certain points, usually on the boundary. In these cases, the problem can usually be corrected by clustering collocation points at these problem points.
Fig. 3. Example 1.
900
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
Table 1 Results for Example 1 Time, s 10 20 30 40 50 100
Diagonal points ðx; yÞ (0,0)
(50,25)
(100,50)
(150,75)
(200,100)
0.25/0 )0.20/0 )0.11/0 )0.08/0 )0.05/0 )0.01/0
46.42/45.83 32.70/32.78 24.01/24.40 18.80/19.28 15.67/16.10 11.44/11.63
25.08/25.14 26.79/26.76 29.70/29.30 32.68/31.93 35.42/34.39 44.51/43.03
26.81/27.09 56.59/53.31 75.51/70.67 87.88/82.13 95.34/98.91 110.67/106.08
195.47/200 196.10/200 196.48/200 196.71/200 196.88/200 197.28/200
Table entries are FCM/exact temperatures.
Example 2. Consider a quarter circular plate of radius R ¼ 50 mm made of a material with a ¼ 10 mm2 /s, subjected to the initial and the mixed boundary conditions shown in Fig. 4. Since an exact solution for the transient case is not available, only the steady state results of the method will be compared to the following exact steady state solution (obtained using the method of Separation of Variables [12]):
1 200 X 1 r 2k kp sin cosð2khÞ ð17Þ T ðr; hÞ ¼ 50 þ p k¼1 k R 2 where ðr; hÞ are polar coordinates. Table 2 shows the results for T at selected points. The angles listed are in degrees and the column corresponding to t ¼ 1000 s is taken to be the steady state numerical solution. The FCM was applied using 16 collocation points, five on each straight segment and six on the arc, with 100 internal heat sources at the centers of a polar grid with Dr ¼ 1 and Dh ¼ 9°. The external sources were located S ¼ 20 mm away and the time step was 5 s.
Fig. 4. Example 2.
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
901
Table 2 Results for Example 2 r=R; h
t¼5s
t ¼ 10 s
t ¼ 20 s
t ¼ 100 s
t ¼ 1000 s
T Eq. (17)
0.1, 0 0.5, 0 0.9, 0
10.69 13.84 39.83
13.80 19.36 18.65
17.21 34.80 87.2
44.50 61.59 93.0
50.41 65.59 93.88
50.64 65.60 93.34
0.1, 22.5 0.5, 22.5 0.9, 22.5
10.76 16.76 44.11
10.80 19.40 86.88
17.13 32.31 87.28
44.31 57.50 90.64
50.23 61.49 91.43
50.45 61.48 90.72
0.1, 45 0.5, 45 0.9, 45
10.85 18.41 27.16
13.80 17.92 45.89
16.93 24.53 44.28
43.86 45.90 49.20
49.78 49.89 49.98
50 50 50
0.1, 67.5 0.5, 67.5 0.9, 67.5
10.83 16.68 9.19
13.81 15.54 3.99
16.73 16.68 4.65
43.42 34.29 7.24
49.33 38.27 8.52
49.55 38.52 9.28
0.1, 90 0.5, 90 0.9, 90
10.80 13.80 6.45
13.81 14.85 3.29
16.65 14.07 2.32
43.24 30.19 5.12
49.15 34.17 5.90
49.36 34.40 6.65
NB ¼ 16, NF ¼ 30, Dt ¼ 5 s, DS ¼ 20.
Example 3. To check the performance of FCM for irregular regions, consider the L-shaped region in Fig. 5. The source boundary was located at a distance S ¼ 1 m away from the region boundary and a total of 40 evenly spaced boundary were used for collocation: 10 on each long side and five on each of the remaining four short sides. The thermal diffusivity was a ¼ 0:01 m2 /s and the time step was 25 s in agreement with Eq. (8). The FCM results (Table 3) are compared with the finite
Fig. 5. Example 3.
902
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
Table 3 Results for Example 3 t ¼ 250 s
x, y 0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.2, 0.2, 0.3,
0.6 0.6 0.6 0.7 0.7 0.7 0.8 0.8 0.8 0.2 0.9 0.9
t ¼ 2500 s
TFCM
TFDM
TFCM
TFDM
651.93 357.54 150.51 634.93 336.66 131.62 586.42 293.08 112.80 440.50 189.98 66.66
656 370 162 639 353 150 590 307 126 444 202 78
747.99 504.50 274.01 717.50 454.22 214.25 644.18 375.27 170.05 470.72 233.30 95.81
738 492 259 706 448 220 637 374 174 468 236 103
difference method (FDM) reported in [12]. The FDM results are based on an explicit FDM version of Eq. (1) in which a forward difference in time, and a centered difference in both x and y are used. Using Dx ¼ Dy ¼ D, the largest stable time step is chosen such that Dt2 6 D2 =4a. For Dx ¼ Dy ¼ 0:1 m and a ¼ 0:01 m2 /s, we obtain Dt ¼ 0:25 s. The temperature solution after 100 and 1000 time steps are shown in Table 3. Example 4. To check the performance of FCM for time-dependent boundary conditions, we consider a quarter ellipse made of a material with a ¼ 10 mm2 /s as shown in Fig. 6. The initial temperature distribution is nonuniform and the boundary conditions consist of insulated straight edges and a time-dependent boundary condition on the curved part of the boundary. The FCM was applied to this problem using 30 collocation points (eight points on each side and 14 points on
Fig. 6. Example 4.
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
903
Table 4 Results for Example 4 Time, s 10 20 50 100 200
Radial points ðx; yÞ (0,0)
(10,10)
(20,20)
(30,30)
(40,40)
(44.72,44.72)
2.97/5 7.31/10 21.85/25 47.72/50 98.64/100
5.54/7.5 9.38/12.5 24.58/27.5 50.45/52.5 101.35/102
12.55/15 16.20/20 32.6/35 58.36/60 109.14/110
21.02/27.5 31.06/32.5 44.82/42.5 71.75/72.5 121.97/122
35.10/45 54.24/50 61.32/65 90.62/90 139.87/140
55.08/55 59.90/60 75.05/75 99.99/100 145/145
Table entries are FCM/exact.
the arc) and 100 internal field points centered on a grids system generated by 10 angular and 10 radial sectors. Using S ¼ 20 mm and Dt ¼ 10 s, temperatures were computed at six different field points located on a 45° radial line. The results (shown in Table 4) agree with the exact solution: T ðx; y; tÞ ¼
x2 y2 þ þ 0:5t 200 50
ð18Þ
6. Conclusions This paper presents an efficient application of the FCM to transient heat conduction problems in two-dimensional regions. The method utilizes GreenÕs function for the pde and can be used reliably in a wide variety of heat transfer problems. The results in the four example problems demonstrate that FCM: 1. is applicable to regions of arbitrary plan form; 2. can handle time-dependent boundary conditions; 3. utilizes a small matrix size which is independent of the number of solution points and needs to be inverted only once; 4. does not have problems with boundary and corner points; 5. can be used to find the solution anywhere in the region or on the boundary; 6. is accurate and compact. In most cases we found that the solution was stable for a wide range of the chosen distance S. As a result, S may be regarded as an optimization parameter. Experiments with the FCM algorithm on many problems show that this value should be taken to be between the smallest and largest diameter that characterizes the region. This gives a minimum value for the sum of the squares of the collocation errors in Eq. (9).
References [1] G. Burgess, E. Mahajerin, On the numerical solution of LaplaceÕs equation using personal computers, International Journal of Mechanical Engineering Education 13 (1) (1984) 45–54.
904
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
[2] M.N. Ozisik, Heat Conduction, Wiley, NY, 1980. [3] M.A. Golberg, Method of fundamental solution for PoisonÕs equation, Engineering Analysis with Boundary Elements 16 (3) (1995) 205–213. [4] R. Mathon, R.L. Johnson, The approximate solution of elliptic boundary value problems by fundamental solutions, SIAM Journal of Numerical Analysis 14 (1977) 638–650. [5] R.T. Fenner, Source field superposition analysis of two-dimensional potential problems, International Journal for Numerical Methods in Engineering 32 (1991) 1079–1091. [6] F.P. Incropera, D.P. Dewitt, Fundamentals of Heat Transfer, John Wiley, 1996. [7] G.D. Smith, Numerical Solution of Partial differential Equations, Oxford University Press, 1965. [8] G. Comini, S. Del Giudice, C. Nonino, Finite Element Analysis in Heat Transfer, Taylor & Francis, 1994. [9] M.A. Golberg, C.S. Chen, The method of fundamental solution for potential, Helmholtz, and diffusion problems, in: M.A. Golberg (Ed.), Boundary Integral Methods: Numerical and Mathematical Aspects, Chapter 4, WIT Press & Computational Mechanics Publications, Boston, Southampton, 1999, pp. 105–176. [10] S.P. Walker, R.T. Fenner, Treatment of corners in BIE: analysis of potential problems, International Journal for Numerical Methods in Engineering 28 (1989) 2569–2581. [11] A.J. Kassab, R.S. Nordlund, Addressing the corner problem in the BEM solution of heat conduction problems, Communications in Numerical Methods in Engineering 10 (1994) 385–392. [12] R. Haberman, Elementary Applied Partial Differential Equations, Prentice Hall, 1983.
Transient heat flow analysis using the fundamental collocation method Gary Burgess b
a,*
, Enayat Mahajerin
b
a School of Packaging, Michigan State University, East Lansing, Michigan, MI 48824, USA Department of Mechanical Engineering, Saginaw Valley State University, University Center, Michigan, MI 48710, USA
Received 2 March 2001; accepted 13 January 2003
Abstract The fundamental collocation method for general two dimensional transient heat flow problems in solids is presented. The numerical procedure is applied directly to heat conduction in the time domain and is illustrated using several example problems. The method is capable of handling problems of arbitrary shapes subjected to arbitrary initial conditions and mixed time-dependent boundary conditions. Its simplicity and minimal computational effort combined with accuracy and stability of the results make the method an excellent alternate approach to domain methods such as the finite difference and finite element methods. Ó 2003 Elsevier Science Ltd. All rights reserved. Keywords: Fundamental collocation method; Finite difference; Finite element; Transient heat flow
1. Introduction The general problem of transient heat flow in two dimensions involves conduction heat transfer in regions of arbitrary shape where the temperature is a function of both time and position. Such problems generally occur before steady state operating conditions are reached. The time duration of the transient part can be of importance in the design of machine parts such as those in turbines and auto engines. Analytical solutions may be found for simple geometries and initial/boundary conditions. For the general case, a numerical solution is preferred. Domain type numerical methods such as finite difference and finite element methods are not necessarily the best approaches due to setup difficulties for arbitrary regions, accuracy, and the fact that the solution is
*
Corresponding author. Tel.: +1-517-355-9504; fax: +1-517-535-8999. E-mail address: [email protected] (G. Burgess).
1359-4311/03/$ - see front matter Ó 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S1359-4311(03)00026-7
894
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
available only at the grid or the element points. In this paper, a much simpler approach based on the fundamental collocation method (FCM) [1] is presented. It can handle regions of arbitrary shapes subjected to general initial and boundary conditions and gives the solution anywhere inside the region and on the boundary. The FCM algorithm presented here is applied directly in the time domain where a GreenÕs function solution to the governing partial differential equation is readily available.
2. Problem statement Consider a general two-dimensional region R with its boundary oR (Fig. 1). We are interested in the time-dependent temperature field, T ðx; y; tÞ. The governing equation for unsteady Fourier heat conduction in a homogeneous isotropic material with no internal heat generation is well-known [2]: oT ¼ ar2 T ðx; y; tÞ ot
ð1Þ
for x; y in R
The initial temperature distribution and the boundary conditions are taken to be T ðx; y; 0Þ ¼ gðx; yÞ aT ðxB ; yB ; tÞ þ b
ð2Þ
in R
oT ðxB ; yB ; tÞ ¼c on
on oR
ð3Þ
where a, b, c, specified functions of position and time; a, thermal diffusivity of the material, assumed to be constant; r2 , o2 =ox2 þ o2 =oy 2 is the Laplace operator; gðx; yÞ, specified initial temperature distribution; xB , yB , boundary point; n, outward normal to the boundary; oT =on ¼ ðoT =oxÞ cos h þ ðoT =oyÞ sin h; h, the angle that outward normal to the boundary makes with the +x axis.
Fig. 1. The problem of interest.
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
895
3. Numerical analysis Steady state heat conduction problems have been handled successfully using the FCM with excellent results [3–5]. The addition of a new variable, time, does not add any real complexity to the FCM algorithm as it does with most other numerical methods. Analytical methods such as separation of variables tend to be so complicated that closed form (including series) solutions are nearly impossible, especially when the region under consideration is arbitrary and the boundary conditions are nonhomogeneous and possibly time-dependent. Lumped capacity problems are about the only problems that can be handled analytically for general regions. Graphical methods like the Schmidt Plots and Heisler Charts [6] are based on the finite difference approach and are only rough approximations. Numerical methods for the unsteady case usually produce sequential solutions over time in which the next generation result depends on the entire history of previous generation results. Numerical stability of the solution is often questionable due to the accumulation of round off error. Domain approaches such as the finite difference and finite element methods have been successfully applied to unsteady cases [7,8]. However, these methods tend to be complicated and the accuracy of the results for regions of arbitrary shapes is questionable. Boundary approaches such as the boundary integral equation method and the boundary element method can significantly reduce the time and effort spent on the computation because they reduce the dimensionality of the problem by one [9]. However, these methods require special treatment when the solution point is near the boundary of the region [10,11]. A much simpler boundary-oriented method called the FCM [1,3,4] significantly reduces computational effort and usually eliminates difficulties with the solution point being close to the boundary. The method consists of embedding the region of interest inside an infinite plane of the same material. A finite number of ‘‘heat sources’’ are placed anywhere outside the boundary of the region, oR, as shown in Fig. 2. The simplest way to do this is to locate the source points on a congruent boundary, oB, at a distance S away from the actual boundary. In this case, NB boundary points are selected for the satisfaction of prescribed boundary conditions by collocation and NB sources are placed along the normals to these points a distance S away. At regular time intervals Dt, each source j, j ¼ 1, NB, deposits a momentary blast of heat over an arbitrary but small patch of area DAj centered on the source point. This blast of heat is sufficient to raise the temperature of the patch from zero to Tj , where Tj is the jth ‘‘source temperature’’ to be determined by the FCM. If this is allowed to continue for ‘‘M’’ time steps, then the temperature due to all outside sources at location ðx; yÞ inside a region initially at zero temperature at time t ¼ M Dt is 3 2 ðxxj Þ2 þðyyj Þ2 NB X M1 exp 4aðt‘DtÞ X 5DAj Tj ð‘DtÞ4 ð4Þ T ðx; y; tÞ ¼ 4paðt ‘DtÞ j¼1 ‘¼0 The quantity in square brackets in Eq. (4) is a modified version of the GreenÕs function for a temperature spike at a point in an infinite medium [2]. This quantity is dimensionless and can be thought of as the fraction of the source temperature Tj imposed at time. ‘ Dt at source location ðxj ; yj Þ which is produced at the field location ðx; yÞ at time t ¼ M Dt, M > ‘. This fraction
896
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
Fig. 2. The numerical setup for FCM.
decreases with distance from the source and with the lag time ðM ‘Þ Dt. Eq. (4) satisfies Eq. (1) exactly for all ðx; y; tÞ. For simplicity, each of the patch areas DAj can be taken as unity. It can be shown [2] that an initial temperature distribution is equivalent to an infinite number of infinitesimal internal heat sources which release momentary blasts of heat at time t ¼ 0 only, and in amounts sufficient to raise their temperatures to the specified initial conditions. Subdividing the region into internal patches of area DAk centered on field points ðxk ; yk Þ, k ¼ 1, NF, and applying the same reasoning as in Eq. (4), the FCM solution is the sum of contributions from NF known internal sources and NB unknown external ones: 3 2 þðyyk Þ2 exp ðxxk Þ 4at 5DAk T ðx; y; tÞ ¼ gðxk ; yk Þ4 4pat k¼1 3 2 ðxxj Þ2 þðyyj Þ2 NB X M1 exp 4aðt‘DtÞ X 5 þ Tj ð‘DtÞ4 4paðt ‘DtÞ j¼1 ‘¼0 NF X
2
ð5Þ
where t ¼ M Dt, NB is the number of heat sources outside the boundary and NF is the number of heat sources inside the region representing the initial temperature distribution. Eq. (5) satisfies Eq. (1) everywhere inside the region, and Eq. (2) at the ðxk ; yk Þ. All that remains is to satisfy Eq. (3) at the NB collocation points on the boundary. Substituting Eq. (5) into Eq. (3) and splitting off the last term in the sum over time requires that
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904 NB X
897
Dij Tj ðt DtÞ
j¼1
2 3 exp ðxi xk Þ2 þðyi yk Þ2 4at ðxi xk Þ cos hi þ ðyi yk Þ sin hi 4 ¼ ci DAk 5 gðxk ; yk Þ ai bi 2at 4pat k¼1 3 2 exp ðxi xj Þ2 þðyi yj Þ2 NB X M 2 X 4aðt‘DtÞ ðxi xj Þ cos hi þ ðyi yj Þ sin hi 4 5 ð6Þ Tj ð‘DtÞ ai bi 2aðt ‘DtÞ 4paðt ‘DtÞ j¼1 ‘¼0 NF X
for each i ¼ 1 to NB at each time t ¼ M Dt, M ¼ 1; 2; 3; . . ., where D is the NBNB ‘‘influence matrix’’ with 3 2 exp ðxi xj Þ2 þðyi yj Þ2 4aDt ðxi xj Þ cos hi þ ðyi yj Þ sin hi 4 5 Dij ¼ ai bi ð7Þ 2aDt 4paDt If the current time at which boundary conditions are to be satisfied is t ¼ M Dt, then Eq. (6) gives an NBNB system of linear equations in the NB outside source temperatures Tj at the previous time step t Dt in terms of the time steps. ‘ Dt, ‘ ¼ 0 to M 2 leading up to this. Beginning with M ¼ 1, the source temperatures Tj at time ðM 1Þ Dt are found from Eq. (6) and stored in the Mth column of a solution matrix which contains the whole history of the source temperatures up to the present time M Dt. Since an NBNB system of equations must be solved at each time step, and since the influence matrix is independent of time, this process is best handled by inverting D. Once the source temperatures are known at a particular time, Eq. (5) and its derivatives are used to find the temperatures and fluxes at any point on the boundary or inside the region.
4. Stability and steady state For the case where temperature is specified everywhere on the boundary, ai ¼ 1, bi ¼ 0, and ci ¼ specified temperature. The influence matrix in Eq. (7) consists of terms of the form expðk=DtÞ=Dt, which can be shown to reach a maximum value when Dt ¼ k. In order for the numerical inversion of the influence matrix D to be reasonably accurate, the diagonal terms of D should be strong. Because the diagonal terms represent the influence of an outside source on its closest boundary point a distance S away, the diagonal elements can be made to be larger than any other element in the same row by choosing S2 ð8Þ 4a Eq. (8) gives relationship between the chosen time step and the distance S between the source and the region boundaries which should generate a stable solution. In practice, fairly large deviations from this can be tolerated as long as Dt is larger than the value in Eq. (8). Smaller DtÕs usually result in inaccurate transient responses. Physically, Eq. (8) says that the time step must be sufficiently long to allow a blast of heat at distance S away to reach the boundary. For general Dt ¼
898
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
boundary conditions ai 6¼ 0, bi 6¼ 0, it is not clear how to choose the time step in Eq. (7) to make the diagonal of D strong. In this case, it is best to choose Dt as in Eq. (8), since the physical interpretation still applies to general boundary conditions. The choice of the distance S is somewhat arbitrary. In practice, S is usually chosen to be somewhere between the smallest and the largest dimensions of the region being investigated. The largest dimension is defined to be the diameter of the smallest circle into which the region can be embedded. As S increases, the determinant of the influence matrix rapidly decreases, so there is a practical limit to the size of S. This limit is largely dependent on machine precision. Experience shows that the best way to choose S is to start with a small value and gradually increase it until the sum of the squares of the errors (SSE) between the prescribed boundary conditions and the delivered boundary conditions is a minimum, 2 NB
X oT cj aj Tj bj ¼ minimum ð9Þ SSE ¼ onj j¼1 where aj , bj and cj are the specified parameters in Eq. (3) and Tj and oT =onj are the values delivered by the FCM solution (Eq. (5) and its normal derivative). In the SSE, it is a good idea to include boundary points in-between the ones chosen for the collocation process. The resulting SSE will then include both roundoff errors from the collocation process and oscillation errors, if any, due to the method itself. As time gets large, the effect of the initial temperature distribution in Eq. (5) dies off. If the boundary conditions do not vary over time, a steady state solution exists and the source temperatures Tj ð‘DtÞ are expected to approach constant values Tj . When this happens, Eq. (5) takes the form 3 2 ðxxj Þ2 þðyyj Þ2 t=Dt exp NB 4aðt‘DtÞ X X 5 as t ! 1 Tj 4 T ðx; y; tÞ ¼ ð10Þ 4paðt ‘DtÞ j¼1 ‘¼0 The second sum in Eq. (10) can be evaluated as an integral having the form Z 1 expðk=tÞ dt IðkÞ ¼ t 0 To perform the integration, first differentiate with respect to k, Z 1 Z dI expðk=tÞ 1 1 d 1 ¼ ðexpðk=tÞÞ dt ¼ dt ¼ 2 dk t k 0 dt k 0
ð11Þ
ð12Þ
and integrate to get I back, I ¼ ‘nk þ C
ð13Þ
When applied to Eq. (10), the steady state solution has the following FCM representation: T ðx; yÞ ¼
NB X
Tj ‘n½ðx xj Þ2 þ ðy yj Þ2
j¼1
which is identical to the FCM solution of the Laplace equation used in [1].
ð14Þ
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
899
5. Examples The method was applied to several examples. Examples were chosen to illustrate performance for different geometries, initial conditions, and boundary conditions. In order to demonstrate the accuracy of the method, a comparison to the exact solution was made when available. Example 1. Consider a rectangular region AB made of bronze (a ¼ 14:1 mm2 /s) subjected to the initial and boundary conditions shown in Fig. 3. The exact solution can be obtained by separation of variables,
T ðx; y; tÞ ¼
2 2
1 X 1 X 200 ipx jpy ip j2 p2 þ 2 t Cij sin xy þ sin exp a AB A B A2 B i¼1 j¼1
ð15Þ
where Cij ¼
4 ½100 ð1Þiþj 200 2 ijp
ð16Þ
The FCM was applied to this problem using A ¼ 200 mm, B ¼ 100 mm. The boundary was represented by 40 equally spaced collocation points, 10 per side, and the external heat sources used to satisfy boundary conditions were placed on a congruent boundary located at a distance S ¼ 10 mm away. The interior of the region was divided into 100 rectangular patches and the temperature distribution was represented by internal heat sources located at the centers of these patches. A time increment of 10 s was used to compute the temperatures at five equally spaced points along a diagonal including the corners. The minimum recommended time step from Eq. (8) is Dt ¼ 1:8 s. The results for the temperature are shown in Table 1. Changing the distance S between the source boundary and the boundary of the region affects these results very little. The practice of changing S is recommended as a check on stability; the results should be the same over a wide range of values. In rare cases, very large or very small S will affect the solution at certain points, usually on the boundary. In these cases, the problem can usually be corrected by clustering collocation points at these problem points.
Fig. 3. Example 1.
900
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
Table 1 Results for Example 1 Time, s 10 20 30 40 50 100
Diagonal points ðx; yÞ (0,0)
(50,25)
(100,50)
(150,75)
(200,100)
0.25/0 )0.20/0 )0.11/0 )0.08/0 )0.05/0 )0.01/0
46.42/45.83 32.70/32.78 24.01/24.40 18.80/19.28 15.67/16.10 11.44/11.63
25.08/25.14 26.79/26.76 29.70/29.30 32.68/31.93 35.42/34.39 44.51/43.03
26.81/27.09 56.59/53.31 75.51/70.67 87.88/82.13 95.34/98.91 110.67/106.08
195.47/200 196.10/200 196.48/200 196.71/200 196.88/200 197.28/200
Table entries are FCM/exact temperatures.
Example 2. Consider a quarter circular plate of radius R ¼ 50 mm made of a material with a ¼ 10 mm2 /s, subjected to the initial and the mixed boundary conditions shown in Fig. 4. Since an exact solution for the transient case is not available, only the steady state results of the method will be compared to the following exact steady state solution (obtained using the method of Separation of Variables [12]):
1 200 X 1 r 2k kp sin cosð2khÞ ð17Þ T ðr; hÞ ¼ 50 þ p k¼1 k R 2 where ðr; hÞ are polar coordinates. Table 2 shows the results for T at selected points. The angles listed are in degrees and the column corresponding to t ¼ 1000 s is taken to be the steady state numerical solution. The FCM was applied using 16 collocation points, five on each straight segment and six on the arc, with 100 internal heat sources at the centers of a polar grid with Dr ¼ 1 and Dh ¼ 9°. The external sources were located S ¼ 20 mm away and the time step was 5 s.
Fig. 4. Example 2.
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
901
Table 2 Results for Example 2 r=R; h
t¼5s
t ¼ 10 s
t ¼ 20 s
t ¼ 100 s
t ¼ 1000 s
T Eq. (17)
0.1, 0 0.5, 0 0.9, 0
10.69 13.84 39.83
13.80 19.36 18.65
17.21 34.80 87.2
44.50 61.59 93.0
50.41 65.59 93.88
50.64 65.60 93.34
0.1, 22.5 0.5, 22.5 0.9, 22.5
10.76 16.76 44.11
10.80 19.40 86.88
17.13 32.31 87.28
44.31 57.50 90.64
50.23 61.49 91.43
50.45 61.48 90.72
0.1, 45 0.5, 45 0.9, 45
10.85 18.41 27.16
13.80 17.92 45.89
16.93 24.53 44.28
43.86 45.90 49.20
49.78 49.89 49.98
50 50 50
0.1, 67.5 0.5, 67.5 0.9, 67.5
10.83 16.68 9.19
13.81 15.54 3.99
16.73 16.68 4.65
43.42 34.29 7.24
49.33 38.27 8.52
49.55 38.52 9.28
0.1, 90 0.5, 90 0.9, 90
10.80 13.80 6.45
13.81 14.85 3.29
16.65 14.07 2.32
43.24 30.19 5.12
49.15 34.17 5.90
49.36 34.40 6.65
NB ¼ 16, NF ¼ 30, Dt ¼ 5 s, DS ¼ 20.
Example 3. To check the performance of FCM for irregular regions, consider the L-shaped region in Fig. 5. The source boundary was located at a distance S ¼ 1 m away from the region boundary and a total of 40 evenly spaced boundary were used for collocation: 10 on each long side and five on each of the remaining four short sides. The thermal diffusivity was a ¼ 0:01 m2 /s and the time step was 25 s in agreement with Eq. (8). The FCM results (Table 3) are compared with the finite
Fig. 5. Example 3.
902
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
Table 3 Results for Example 3 t ¼ 250 s
x, y 0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3, 0.2, 0.2, 0.3,
0.6 0.6 0.6 0.7 0.7 0.7 0.8 0.8 0.8 0.2 0.9 0.9
t ¼ 2500 s
TFCM
TFDM
TFCM
TFDM
651.93 357.54 150.51 634.93 336.66 131.62 586.42 293.08 112.80 440.50 189.98 66.66
656 370 162 639 353 150 590 307 126 444 202 78
747.99 504.50 274.01 717.50 454.22 214.25 644.18 375.27 170.05 470.72 233.30 95.81
738 492 259 706 448 220 637 374 174 468 236 103
difference method (FDM) reported in [12]. The FDM results are based on an explicit FDM version of Eq. (1) in which a forward difference in time, and a centered difference in both x and y are used. Using Dx ¼ Dy ¼ D, the largest stable time step is chosen such that Dt2 6 D2 =4a. For Dx ¼ Dy ¼ 0:1 m and a ¼ 0:01 m2 /s, we obtain Dt ¼ 0:25 s. The temperature solution after 100 and 1000 time steps are shown in Table 3. Example 4. To check the performance of FCM for time-dependent boundary conditions, we consider a quarter ellipse made of a material with a ¼ 10 mm2 /s as shown in Fig. 6. The initial temperature distribution is nonuniform and the boundary conditions consist of insulated straight edges and a time-dependent boundary condition on the curved part of the boundary. The FCM was applied to this problem using 30 collocation points (eight points on each side and 14 points on
Fig. 6. Example 4.
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
903
Table 4 Results for Example 4 Time, s 10 20 50 100 200
Radial points ðx; yÞ (0,0)
(10,10)
(20,20)
(30,30)
(40,40)
(44.72,44.72)
2.97/5 7.31/10 21.85/25 47.72/50 98.64/100
5.54/7.5 9.38/12.5 24.58/27.5 50.45/52.5 101.35/102
12.55/15 16.20/20 32.6/35 58.36/60 109.14/110
21.02/27.5 31.06/32.5 44.82/42.5 71.75/72.5 121.97/122
35.10/45 54.24/50 61.32/65 90.62/90 139.87/140
55.08/55 59.90/60 75.05/75 99.99/100 145/145
Table entries are FCM/exact.
the arc) and 100 internal field points centered on a grids system generated by 10 angular and 10 radial sectors. Using S ¼ 20 mm and Dt ¼ 10 s, temperatures were computed at six different field points located on a 45° radial line. The results (shown in Table 4) agree with the exact solution: T ðx; y; tÞ ¼
x2 y2 þ þ 0:5t 200 50
ð18Þ
6. Conclusions This paper presents an efficient application of the FCM to transient heat conduction problems in two-dimensional regions. The method utilizes GreenÕs function for the pde and can be used reliably in a wide variety of heat transfer problems. The results in the four example problems demonstrate that FCM: 1. is applicable to regions of arbitrary plan form; 2. can handle time-dependent boundary conditions; 3. utilizes a small matrix size which is independent of the number of solution points and needs to be inverted only once; 4. does not have problems with boundary and corner points; 5. can be used to find the solution anywhere in the region or on the boundary; 6. is accurate and compact. In most cases we found that the solution was stable for a wide range of the chosen distance S. As a result, S may be regarded as an optimization parameter. Experiments with the FCM algorithm on many problems show that this value should be taken to be between the smallest and largest diameter that characterizes the region. This gives a minimum value for the sum of the squares of the collocation errors in Eq. (9).
References [1] G. Burgess, E. Mahajerin, On the numerical solution of LaplaceÕs equation using personal computers, International Journal of Mechanical Engineering Education 13 (1) (1984) 45–54.
904
G. Burgess, E. Mahajerin / Applied Thermal Engineering 23 (2003) 893–904
[2] M.N. Ozisik, Heat Conduction, Wiley, NY, 1980. [3] M.A. Golberg, Method of fundamental solution for PoisonÕs equation, Engineering Analysis with Boundary Elements 16 (3) (1995) 205–213. [4] R. Mathon, R.L. Johnson, The approximate solution of elliptic boundary value problems by fundamental solutions, SIAM Journal of Numerical Analysis 14 (1977) 638–650. [5] R.T. Fenner, Source field superposition analysis of two-dimensional potential problems, International Journal for Numerical Methods in Engineering 32 (1991) 1079–1091. [6] F.P. Incropera, D.P. Dewitt, Fundamentals of Heat Transfer, John Wiley, 1996. [7] G.D. Smith, Numerical Solution of Partial differential Equations, Oxford University Press, 1965. [8] G. Comini, S. Del Giudice, C. Nonino, Finite Element Analysis in Heat Transfer, Taylor & Francis, 1994. [9] M.A. Golberg, C.S. Chen, The method of fundamental solution for potential, Helmholtz, and diffusion problems, in: M.A. Golberg (Ed.), Boundary Integral Methods: Numerical and Mathematical Aspects, Chapter 4, WIT Press & Computational Mechanics Publications, Boston, Southampton, 1999, pp. 105–176. [10] S.P. Walker, R.T. Fenner, Treatment of corners in BIE: analysis of potential problems, International Journal for Numerical Methods in Engineering 28 (1989) 2569–2581. [11] A.J. Kassab, R.S. Nordlund, Addressing the corner problem in the BEM solution of heat conduction problems, Communications in Numerical Methods in Engineering 10 (1994) 385–392. [12] R. Haberman, Elementary Applied Partial Differential Equations, Prentice Hall, 1983.