- Email: [email protected]
Numerical solution of the heat conduction equation with the electro-thermal analogy and the code PSPICE
Applied Mathematics and Computation 162 (2005) 103–113 www.elsevier.com/locate/amc
Numerical solution of the heat conduction equation with the electro-thermal analogy and the code PSPICE Francisco Alhama a, Antonio Campo a
b,*
, Joaquın Zueco
c
Departmento de Fısica Aplicada, E.T.S. Ingenieros Industriales, Universidad Politecnica de Cartagena, 30203 Cartagena, Murcia, Spain b Department of Mechanical Engineering, The University of Vermont, Burlington, VT 05405, USA c Departmento de Ingenierıa, Termica y de Fluidos, E.T.S Ingenieros Industriales, Universidad Politecnica de Cartagena, 30203 Cartagena, Murcia, Spain
Abstract The central objective of this paper is to provide numerical analysts with a new procedure named the network simulation method (NSM) for solving the heat conduction equation in bodies of regular shape. In principle, NSM rests on the electro-thermal analogy (loosely called the resistance–capacitance analogy or the RC analogy) that exists between the unsteady, unidirectional conduction of heat and the unsteady flow of electric current. Once the electric network model has been set up for the heat conduction equation, the numerical treatment of the analog electric circuit equation can be easily done with the computer code PSPICE. As a conceptual example, the allied numerical solution of the heat conduction equation for a large plate with symmetric surface temperatures has been carried out, demonstrating that the temperature–time histories and the heat flux–time histories can be obtained simultaneously, quickly, and accurately for the entire time domain. Ó 2004 Elsevier Inc. All rights reserved.
*
Corresponding author. E-mail address: [email protected] (A. Campo).
0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.123
104
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
Nomenclature A CV Ci j k L N q RiD t T Ti Tw x X
surface area heat capacity at constant volume capacitance of cell i heat flux density i electric current thermal conductivity semi-thickness number of cells in the electric network heat flux resistance of cell i time temperature or voltage initial temperature surface temperature space variable dimensionless x, Eq. (11)
Greek letters a thermal diffusivity, k=CV h dimensionless T , Eq. (11) s dimensionless t or Fourier number, Eq. (11)
1. Introduction The finite-difference method (FDM) is undoubtedly the most popular numerical technique for solving partial differential equations, such as the heat conduction equation [1]. In general terms, FDM is classified as explicit or implicit, depending on the kind of finite-difference analog that is employed for the transformation of the time derivative in the heat conduction equation. The underlying theory and a vast number of applications of the FDM to the heat conduction equation may be found in [1]. Alternatively, in this paper a new procedure for the numerical solution of the heat conduction equation termed the network simulation method (NSM) is developed, tested and compared utilizing the computer code PSPICE. Historically, the genesis of the NSM lies in the electro-thermal analogy that exists between unsteady, unidirectional conduction of heat and unsteady flow of electric current. In the platform of the electro-thermal analogy, there are two characteristic variables in the heat conduction equation: (1) temperature, T , (analogous to the voltage, E) and (2) the heat flux density, j, (analogous to the electric current, i).
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
105
From the standpoint of numerical analysis, NSM falls under the category of a hybrid technique possessing differential components of analytic form as well as finite-difference components of approximate form. The SPICE family of codes belongs to the electrical engineering arena and was constructed by electrical engineers in order to speed up the computer-aided design of electric circuits. The numerical solution of the heat conduction equation in a large plate with symmetrical prescribed surface temperature has been presented as a conceptual example. Extensions to other more involved boundary conditions is a straight forward procedure.
2. Mathematical model A physical system consists of unsteady and unidirectional heat conduction in a large plate of thickness 2L maintained at a uniform temperature Ti . For t > 0, a prescribed temperature Tw is applied at the two exposed surfaces of the plate. For a plate made from a solid material with constant heat capacity at constant volume CV and constant thermal conductivity k, the mathematical model is expressed by the heat conduction equation CV
oT oj ¼ ot ox
ð1Þ
subject to the initial condition T ¼ Ti ;
t¼0
ð2Þ
and the boundary conditions oT ¼ 0; x ¼ 0; ox T ¼ Tw ; x ¼ L:
j¼k
ð3aÞ ð3bÞ
The primary dependent variables in the preceding model are the temperature T and the heat flux density j. The rest of the symbols is defined in Nomenclature.
3. Brief outline of the finite-difference method (FDM) and the method of lines (MOL) Without any doubt, the most popular computational procedure for the numerical solution of the heat conduction equation (1) is the FDM [1]. FDM is classified as explicit or implicit, depending on how the finite-difference analog is employed for the treatment of the time derivative in Eq. (1). The solution of the explicit scheme amounts to the evaluation of the dependent variable T at each location and at each time interval in terms of the value of T at the same
106
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
location and the neighboring values of T taken from the preceding time interval. Starting with the initial condition as given in Eq. (2), the numerical procedure marches from time equal to zero for a pre-determined sequence of time intervals. Although the explicit scheme is simple, its major drawback is linked to the size of the permissible time step which is severely restricted by stability considerations [1]. Conversely, the numerical solution of Eq. (1) via the implicit scheme reduces to an algebraic equation at each location and each time interval involving the dependent variable T at the same location and the neighboring T expressed in terms of the following time interval. Thus, this scheme gives rise to a set of algebraic equations to be solved numerically at each future time step. Logically, a solution with the implicit scheme is more involved computationally than that of the explicit scheme. Nevertheless, the implicit scheme is advantageous in that there is no restrictions on the size of the time step which is connected to stability limitations [1]. In addition to the FDM, another semi-numerical method exists for solving Eq. (1) which is less popular than the FDM. The method is called the method of lines (MOL) and possesses hybrid characteristics of differential and finitedifference nature [2]. MOL is essentially a transformation technique that contains two viable options: (1) the longitudinal method of lines (LMOL) and (2) the transversal method of lines (TMOL). The implementation of LMOL in Eq. (1) begins with the replacement of the derivative with respect to the space variable with a finite-difference analog while the first-order derivative with respect to time remains in continuous form. Consequently, this elementary operation generates a system of ordinary differential equations of first-order, i.e., an initial value problem (IVP) which depending of its structure may be easily solved analytically or numerically. If the central finite-difference analog is employed in Eq. (1), the LMOL-based solutions are second-order accurate. The numerical calculations via LMOL may be envisioned as long time asymptotic solutions. The other viable option of MOL is the TMOL. When TMOL is applied to Eq. (1), the first-order derivative with respect to time is discretized while leaving the second-order derivative with respect to the space variable in continuous form. Correspondingly, Eq. (1) is replaced by a sequence of nonhomogeneous ordinary differential equations of second-order (two-point value problem (BVP)) that are created at successive preselected time levels. As expected, depending on the form of the ode’s and the number of time levels involved, analytical or numerical solutions may be obtained methodically at each time level. As far as accuracy is concerned, the primary limitation inherent to the TMOL when applied to Eq. (1) is the intrinsic first-order error which is attached to the backward finite-difference analog used for the discretization of the time derivative. The one-time step TMOL solution (say 1-TMOL for short) delivers short-time asymptotic solutions.
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
107
4. Network simulation method (NSM) The network simulation method (also known by the acronym NSM) is essentially a numerical method, but distinct to the finite-difference methods described in the preceding section. NSM has its origins in the theory of electric circuits, specifically in the electro-thermal analogy (the resistance–capacitance analogy or RC analogy) that prevails between the time-dependent flow of heat by conduction in one dimension and the time-dependent flow of electric current. 1 Electrical and thermal systems are said to be analogous if they are formulated in the same domain with similar equations and identical initial and boundary conditions. Accordingly, the equations describing the behavior of thermal systems can be transformed into the equations that represent electrical systems by simply changing the variables. The basic equations that sustain the electro-thermal analogy for applications of unsteady heat conduction have been unraveled in Jakob [3]. A list of thermal and electric variables, their units and the proper one-to-one electro-thermal analogies are given in Table 1. Historically speaking, the electro-thermal analogy was first attempted by Paschkis and Baker [4] for solving the unsteady unidirectional heat conduction equation in a large plate with an analog computer which was called by the authors the Heat and Mass Flow Analyzer. As witnessed in Fig. 1, the analyzer was a massive apparatus composed by a large array of resistors and capacitors for the molding of the electric circuit networks that mimic the heat conduction processes. As a consequence of the ensuing one-to-one correspondence between heat flow and electric current (both time-dependent), the instantaneous voltage differences (indicative of temperature differences) and the instantaneous electric currents (indicative of heat flux densities) were measured simultaneously in [4]. The implementation of the NSM commences with the partition of the large plate comprised by 0 < x < L into a finite number of cells i, i ¼ 1; 2 . . . ; N of equal thickness Dx forming the computational domain. An energy balance in a typical cell i yields the generic relation ðCV Þi
dTi Dj jiD jiþD ¼ ¼ ; Dx dt Dx
ð4Þ
where jiD , and jiþD , signify the respective heat flux densities entering and leaving the cell i as depicted in Fig. 2a. In turn, D and jiD may be defined by Fourier’s law in the following manner jiD ¼ k
1
Ti TiD : Dx=2
ð5Þ
A simplified version of the electro-thermal analogy deals with the unidirectional conduction of heat and the steady flow of electric current and is called Ohm’s law [3].
108
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
Table 1 Analogy between the thermal and electrical quantities Thermal quantity Parameter T j Rth CV sth
Temperature Heat flux density Thermal resistance Heat capacity Time constant
Electrical quantity Unit
Parameter
C W/m2 C m2 /W J/Kg C s
E i R C s
Unit Voltage Current Resistance Capacitance Time constant
V A X F s
Fig. 1. A picture of the Heat and Mass Flow Thermal Analyzer taken from Ref. [4].
In the trio of temperatures TiD , Ti and TiþD participating in the preceding equation, TiD denotes the temperatures at the left extreme, Ti represents the temperature at the center and TiþD , signifies the temperature at the right extreme of the cell i, respectively. Hence, for a typical cell i, i ¼ 1; 2 . . . ; N , the combination of Eqs. (4) and (5) yields DxðCV Þi
dTi TiD Ti Ti TiþD ¼k k : dt Dx=2 Dx=2
ð6Þ
At this point, the reader may recognize the differential-difference structure of Eq. (6) which is emblematic of the electro-thermal analogy within the NSM platform. A distinguishable feature of the NSM approach stems from the fact that the analogy demands discrete intervals of the space variable x and real, continuous time for the independent variable t. Thereby, it is easy to perceive a
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
109
Fig. 2. Electronic network model for a large plate: (a) a typical interior cell i and (b) a cell i ¼ 0 at the midplane of symmetry ðX ¼ 0Þ.
parallelism between the physically-intensive NSM and the mathematicallyoriented LMOL [2]. Next, each term in Eq. (6) is appropriately redefined as an electric current, resulting in the following relation jiD jiþD ji;c ¼ 0;
ð7Þ
where dTi : ð8Þ dt In the context of the prevalent equivalence, Eq. (7) or Eq. (8) may be interpreted as the Kirchhoff’s current law (implying energy conservation) where the temperature T is a continuous single-valued variable that satisfies Kirchhoff’s ji;c ¼ DxðCV Þi
110
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
voltage law while time is the continuous independent variable. In synthesis, the numerical simulation in the sense of an electric network relies on two characteristic variables usable in a typical cell i: one is temperature, Ti , (analogous to the voltage Ei ) and the other is the heat flux density ji (analogous to the electric current i). The independent variable is the time, t. Eq. (6) envisions the presence of two equal resistors with resistances Dx Dx ; RiþD ¼ ð9Þ RiþD ¼ 2k 2k and one capacitor whose capacitance is Ci ¼ DxðCV Þi : ð10Þ In view of this, in a typical cell i, the interconnection between the resistor RiD at the neighbor cell i 1, the capacitor Ci at the cell i and the resistor RiþD at the neighbor cell i þ 1 is illustrated graphically in the portion of the electric network in Fig. 2a. The three vital elements RiD , Ci and RiþD affixed to a typical cell i, may be easily generalized to the total number of N cells making up the complete electric network model of the heat conduction equation in the plate. For the completion of the electric network model, the final requirement is pertinent to incorporation of the initial and boundary conditions. The initial condition is easily adjusted by charging the capacitors Ci of the electric network model to the initial temperature specified in Eq. (2). The prescribed temperature boundary condition, Eq. (8b) requires the application of a voltage source to the ends of cell N as indicated in Fig. 2. 5. Numerical computations Once the electric network model of the heat conduction equation (1) in the large plate has been set up, the temperature field T ðx; tÞ and the heat flux density field jðx; tÞ of the plate may be simulated rapidly and efficiently with a computer code named Simulation Program with Integrated Circuit Emphasis (SPICE). This code was assembled by a team of electrical engineers headed by Nagel [5] with affiliation to the Electronic Research Laboratory of the University of California at Berkeley. Owing to its simplicity, functionability and popularity, the SPICE code has become the staple for the computer-aided design of electric circuits in industry worldwide. A version of SPICE specialized to personal computers is PSPICE [6]. 6. Discussion of results For convenience, let us adopt the customary set of dimensionless variables h¼
T Tw ; Ti Tw
s¼
t ; L2 =a
x X ¼ ; L
ð11Þ
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
111
Fig. 3. The dimensionless temperature history in a large plate parameterized by the dimensionless space variable.
where the scales for temperature, time and space variable are: ðTi T1 Þ, L, and L2 =a 2, respectively. The temperature field T ðx; tÞ and the heat flux density field jðx; tÞ in a large plate modeled by Eqs. (1)–(3) have been accurately determined by the code PSPICE. Certainly, a sensitivity analysis of the number of cells forming the computational domain is a mandatory step in numerical calculus. The numerical experiments revealed that to achieve satisfactory results the optimal number of cells must be N ¼ 100. This choice insures that the intrinsic errors are nearly zero in all cells occupying the computational domain. Fig. 3 shows the trend exhibited by the dimensionless temperature ðT Tw Þ= ðTi Tw Þ as a function of the dimensionless time T ¼ t=ðL2 =aÞ at different locations inside the large plate. The family of curves consists of ten curves, the uppermost curve represents the temperature at the centerplane of the plate ðX ¼ 0Þ and the lowermost curve corresponds to the temperature at a location near the surface of the plate, X ¼ 0:9. Fig. 4 is the companion figure to Fig. 3 which displays the response of the dimensionless instantaneous heat flow rate qL=kAðTw Ti Þ to continuous changes in the dimensionless time h ¼ t=ðL2 =aÞ. In the family of 11 curves, the uppermost curve identifies the instantaneous heat flow at the surface of the plate ðX ¼ 1Þ, whereas the lowermost curve symbolizes the instantaneous heat flow at the midplane of the plate ðX ¼ 0Þ.
2
R2 =as , has been identified as the heat diffusion time in [7].
112
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
Fig. 4. The dimensionless instantaneous heat flux history as a in a large plate parameterized by the dimensionless space variable.
Note that the union of NSM with the code PSPICE may be extended to unsteady heat conduction in long cylinders and spheres requires the specification of two unequal resistors in the group of cells. This modification is basically due to the differences in surface areas of each half of the cells. Likewise, the numerical solution of the heat conduction equation in multidimensional regular bodies do not pose any obstacle to the duo NSM and the code PSPICE either. Additional information about the use of the electro-thermal analogy, the NSM and PSPICE for solving the heat conduction equation may be found in Alhama [8]. 7. Conclusions The coupling of the electro-thermal analogy and the NSM has proved to be a versatile numerical procedure for solving the heat conduction equation with the help of the computer code PSPICE. The spatio-temporal temperatures and the spatio-temporal heat fluxes in large plates may be predicted with good accuracy simultaneously in the entire time domain ð0 < s < 1Þ. References [1] W.F. Ames, Numerical Methods for Partial Differential Equations, third ed., Academic Press, San Diego, CA, 1992.
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
113
[2] O.A. Liskovets, The method of lines (review), Different. Equat. 1 (1965) 1308–1323. [3] M. Jakob, in: Heat Transfer, vol. 1, McGraw-Hill, New York, 1949. [4] V. Paschkis, H.D. Baker, A method for determining unsteady-state heat transfer by means of an electrical analogy, Trans. ASME 64 (1942) 105–112. [5] L.W. Nagel, SPICE: a computer program to simulate semiconductor circuits Memo No. UCB/ ERL M520, Electronic Research Laboratory, University of California, Berkeley, CA, 1975. [6] PSPICE, version 6.0, 1994. Distributor: Microsim Corp., Irvine, CA. [7] P.S. Lykoudis, Non-dimensional numbers as ratios of characteristic times, Int. J. Heat Mass Transfer 33 (1990) 1560–1570. [8] F. Alhama, Estudio de las Respuestas Termicas Transitorias en Procesos no Lineales de Conducci on de Calor Mediante el Metodo de Simulaci on por Redes, Doctoral Thesis (in Spanish), Universidad de Murcia, Cartagena, Murcia, Spain, 1998.
Numerical solution of the heat conduction equation with the electro-thermal analogy and the code PSPICE Francisco Alhama a, Antonio Campo a
b,*
, Joaquın Zueco
c
Departmento de Fısica Aplicada, E.T.S. Ingenieros Industriales, Universidad Politecnica de Cartagena, 30203 Cartagena, Murcia, Spain b Department of Mechanical Engineering, The University of Vermont, Burlington, VT 05405, USA c Departmento de Ingenierıa, Termica y de Fluidos, E.T.S Ingenieros Industriales, Universidad Politecnica de Cartagena, 30203 Cartagena, Murcia, Spain
Abstract The central objective of this paper is to provide numerical analysts with a new procedure named the network simulation method (NSM) for solving the heat conduction equation in bodies of regular shape. In principle, NSM rests on the electro-thermal analogy (loosely called the resistance–capacitance analogy or the RC analogy) that exists between the unsteady, unidirectional conduction of heat and the unsteady flow of electric current. Once the electric network model has been set up for the heat conduction equation, the numerical treatment of the analog electric circuit equation can be easily done with the computer code PSPICE. As a conceptual example, the allied numerical solution of the heat conduction equation for a large plate with symmetric surface temperatures has been carried out, demonstrating that the temperature–time histories and the heat flux–time histories can be obtained simultaneously, quickly, and accurately for the entire time domain. Ó 2004 Elsevier Inc. All rights reserved.
*
Corresponding author. E-mail address: [email protected] (A. Campo).
0096-3003/$ - see front matter Ó 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.123
104
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
Nomenclature A CV Ci j k L N q RiD t T Ti Tw x X
surface area heat capacity at constant volume capacitance of cell i heat flux density i electric current thermal conductivity semi-thickness number of cells in the electric network heat flux resistance of cell i time temperature or voltage initial temperature surface temperature space variable dimensionless x, Eq. (11)
Greek letters a thermal diffusivity, k=CV h dimensionless T , Eq. (11) s dimensionless t or Fourier number, Eq. (11)
1. Introduction The finite-difference method (FDM) is undoubtedly the most popular numerical technique for solving partial differential equations, such as the heat conduction equation [1]. In general terms, FDM is classified as explicit or implicit, depending on the kind of finite-difference analog that is employed for the transformation of the time derivative in the heat conduction equation. The underlying theory and a vast number of applications of the FDM to the heat conduction equation may be found in [1]. Alternatively, in this paper a new procedure for the numerical solution of the heat conduction equation termed the network simulation method (NSM) is developed, tested and compared utilizing the computer code PSPICE. Historically, the genesis of the NSM lies in the electro-thermal analogy that exists between unsteady, unidirectional conduction of heat and unsteady flow of electric current. In the platform of the electro-thermal analogy, there are two characteristic variables in the heat conduction equation: (1) temperature, T , (analogous to the voltage, E) and (2) the heat flux density, j, (analogous to the electric current, i).
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
105
From the standpoint of numerical analysis, NSM falls under the category of a hybrid technique possessing differential components of analytic form as well as finite-difference components of approximate form. The SPICE family of codes belongs to the electrical engineering arena and was constructed by electrical engineers in order to speed up the computer-aided design of electric circuits. The numerical solution of the heat conduction equation in a large plate with symmetrical prescribed surface temperature has been presented as a conceptual example. Extensions to other more involved boundary conditions is a straight forward procedure.
2. Mathematical model A physical system consists of unsteady and unidirectional heat conduction in a large plate of thickness 2L maintained at a uniform temperature Ti . For t > 0, a prescribed temperature Tw is applied at the two exposed surfaces of the plate. For a plate made from a solid material with constant heat capacity at constant volume CV and constant thermal conductivity k, the mathematical model is expressed by the heat conduction equation CV
oT oj ¼ ot ox
ð1Þ
subject to the initial condition T ¼ Ti ;
t¼0
ð2Þ
and the boundary conditions oT ¼ 0; x ¼ 0; ox T ¼ Tw ; x ¼ L:
j¼k
ð3aÞ ð3bÞ
The primary dependent variables in the preceding model are the temperature T and the heat flux density j. The rest of the symbols is defined in Nomenclature.
3. Brief outline of the finite-difference method (FDM) and the method of lines (MOL) Without any doubt, the most popular computational procedure for the numerical solution of the heat conduction equation (1) is the FDM [1]. FDM is classified as explicit or implicit, depending on how the finite-difference analog is employed for the treatment of the time derivative in Eq. (1). The solution of the explicit scheme amounts to the evaluation of the dependent variable T at each location and at each time interval in terms of the value of T at the same
106
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
location and the neighboring values of T taken from the preceding time interval. Starting with the initial condition as given in Eq. (2), the numerical procedure marches from time equal to zero for a pre-determined sequence of time intervals. Although the explicit scheme is simple, its major drawback is linked to the size of the permissible time step which is severely restricted by stability considerations [1]. Conversely, the numerical solution of Eq. (1) via the implicit scheme reduces to an algebraic equation at each location and each time interval involving the dependent variable T at the same location and the neighboring T expressed in terms of the following time interval. Thus, this scheme gives rise to a set of algebraic equations to be solved numerically at each future time step. Logically, a solution with the implicit scheme is more involved computationally than that of the explicit scheme. Nevertheless, the implicit scheme is advantageous in that there is no restrictions on the size of the time step which is connected to stability limitations [1]. In addition to the FDM, another semi-numerical method exists for solving Eq. (1) which is less popular than the FDM. The method is called the method of lines (MOL) and possesses hybrid characteristics of differential and finitedifference nature [2]. MOL is essentially a transformation technique that contains two viable options: (1) the longitudinal method of lines (LMOL) and (2) the transversal method of lines (TMOL). The implementation of LMOL in Eq. (1) begins with the replacement of the derivative with respect to the space variable with a finite-difference analog while the first-order derivative with respect to time remains in continuous form. Consequently, this elementary operation generates a system of ordinary differential equations of first-order, i.e., an initial value problem (IVP) which depending of its structure may be easily solved analytically or numerically. If the central finite-difference analog is employed in Eq. (1), the LMOL-based solutions are second-order accurate. The numerical calculations via LMOL may be envisioned as long time asymptotic solutions. The other viable option of MOL is the TMOL. When TMOL is applied to Eq. (1), the first-order derivative with respect to time is discretized while leaving the second-order derivative with respect to the space variable in continuous form. Correspondingly, Eq. (1) is replaced by a sequence of nonhomogeneous ordinary differential equations of second-order (two-point value problem (BVP)) that are created at successive preselected time levels. As expected, depending on the form of the ode’s and the number of time levels involved, analytical or numerical solutions may be obtained methodically at each time level. As far as accuracy is concerned, the primary limitation inherent to the TMOL when applied to Eq. (1) is the intrinsic first-order error which is attached to the backward finite-difference analog used for the discretization of the time derivative. The one-time step TMOL solution (say 1-TMOL for short) delivers short-time asymptotic solutions.
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
107
4. Network simulation method (NSM) The network simulation method (also known by the acronym NSM) is essentially a numerical method, but distinct to the finite-difference methods described in the preceding section. NSM has its origins in the theory of electric circuits, specifically in the electro-thermal analogy (the resistance–capacitance analogy or RC analogy) that prevails between the time-dependent flow of heat by conduction in one dimension and the time-dependent flow of electric current. 1 Electrical and thermal systems are said to be analogous if they are formulated in the same domain with similar equations and identical initial and boundary conditions. Accordingly, the equations describing the behavior of thermal systems can be transformed into the equations that represent electrical systems by simply changing the variables. The basic equations that sustain the electro-thermal analogy for applications of unsteady heat conduction have been unraveled in Jakob [3]. A list of thermal and electric variables, their units and the proper one-to-one electro-thermal analogies are given in Table 1. Historically speaking, the electro-thermal analogy was first attempted by Paschkis and Baker [4] for solving the unsteady unidirectional heat conduction equation in a large plate with an analog computer which was called by the authors the Heat and Mass Flow Analyzer. As witnessed in Fig. 1, the analyzer was a massive apparatus composed by a large array of resistors and capacitors for the molding of the electric circuit networks that mimic the heat conduction processes. As a consequence of the ensuing one-to-one correspondence between heat flow and electric current (both time-dependent), the instantaneous voltage differences (indicative of temperature differences) and the instantaneous electric currents (indicative of heat flux densities) were measured simultaneously in [4]. The implementation of the NSM commences with the partition of the large plate comprised by 0 < x < L into a finite number of cells i, i ¼ 1; 2 . . . ; N of equal thickness Dx forming the computational domain. An energy balance in a typical cell i yields the generic relation ðCV Þi
dTi Dj jiD jiþD ¼ ¼ ; Dx dt Dx
ð4Þ
where jiD , and jiþD , signify the respective heat flux densities entering and leaving the cell i as depicted in Fig. 2a. In turn, D and jiD may be defined by Fourier’s law in the following manner jiD ¼ k
1
Ti TiD : Dx=2
ð5Þ
A simplified version of the electro-thermal analogy deals with the unidirectional conduction of heat and the steady flow of electric current and is called Ohm’s law [3].
108
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
Table 1 Analogy between the thermal and electrical quantities Thermal quantity Parameter T j Rth CV sth
Temperature Heat flux density Thermal resistance Heat capacity Time constant
Electrical quantity Unit
Parameter
C W/m2 C m2 /W J/Kg C s
E i R C s
Unit Voltage Current Resistance Capacitance Time constant
V A X F s
Fig. 1. A picture of the Heat and Mass Flow Thermal Analyzer taken from Ref. [4].
In the trio of temperatures TiD , Ti and TiþD participating in the preceding equation, TiD denotes the temperatures at the left extreme, Ti represents the temperature at the center and TiþD , signifies the temperature at the right extreme of the cell i, respectively. Hence, for a typical cell i, i ¼ 1; 2 . . . ; N , the combination of Eqs. (4) and (5) yields DxðCV Þi
dTi TiD Ti Ti TiþD ¼k k : dt Dx=2 Dx=2
ð6Þ
At this point, the reader may recognize the differential-difference structure of Eq. (6) which is emblematic of the electro-thermal analogy within the NSM platform. A distinguishable feature of the NSM approach stems from the fact that the analogy demands discrete intervals of the space variable x and real, continuous time for the independent variable t. Thereby, it is easy to perceive a
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
109
Fig. 2. Electronic network model for a large plate: (a) a typical interior cell i and (b) a cell i ¼ 0 at the midplane of symmetry ðX ¼ 0Þ.
parallelism between the physically-intensive NSM and the mathematicallyoriented LMOL [2]. Next, each term in Eq. (6) is appropriately redefined as an electric current, resulting in the following relation jiD jiþD ji;c ¼ 0;
ð7Þ
where dTi : ð8Þ dt In the context of the prevalent equivalence, Eq. (7) or Eq. (8) may be interpreted as the Kirchhoff’s current law (implying energy conservation) where the temperature T is a continuous single-valued variable that satisfies Kirchhoff’s ji;c ¼ DxðCV Þi
110
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
voltage law while time is the continuous independent variable. In synthesis, the numerical simulation in the sense of an electric network relies on two characteristic variables usable in a typical cell i: one is temperature, Ti , (analogous to the voltage Ei ) and the other is the heat flux density ji (analogous to the electric current i). The independent variable is the time, t. Eq. (6) envisions the presence of two equal resistors with resistances Dx Dx ; RiþD ¼ ð9Þ RiþD ¼ 2k 2k and one capacitor whose capacitance is Ci ¼ DxðCV Þi : ð10Þ In view of this, in a typical cell i, the interconnection between the resistor RiD at the neighbor cell i 1, the capacitor Ci at the cell i and the resistor RiþD at the neighbor cell i þ 1 is illustrated graphically in the portion of the electric network in Fig. 2a. The three vital elements RiD , Ci and RiþD affixed to a typical cell i, may be easily generalized to the total number of N cells making up the complete electric network model of the heat conduction equation in the plate. For the completion of the electric network model, the final requirement is pertinent to incorporation of the initial and boundary conditions. The initial condition is easily adjusted by charging the capacitors Ci of the electric network model to the initial temperature specified in Eq. (2). The prescribed temperature boundary condition, Eq. (8b) requires the application of a voltage source to the ends of cell N as indicated in Fig. 2. 5. Numerical computations Once the electric network model of the heat conduction equation (1) in the large plate has been set up, the temperature field T ðx; tÞ and the heat flux density field jðx; tÞ of the plate may be simulated rapidly and efficiently with a computer code named Simulation Program with Integrated Circuit Emphasis (SPICE). This code was assembled by a team of electrical engineers headed by Nagel [5] with affiliation to the Electronic Research Laboratory of the University of California at Berkeley. Owing to its simplicity, functionability and popularity, the SPICE code has become the staple for the computer-aided design of electric circuits in industry worldwide. A version of SPICE specialized to personal computers is PSPICE [6]. 6. Discussion of results For convenience, let us adopt the customary set of dimensionless variables h¼
T Tw ; Ti Tw
s¼
t ; L2 =a
x X ¼ ; L
ð11Þ
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
111
Fig. 3. The dimensionless temperature history in a large plate parameterized by the dimensionless space variable.
where the scales for temperature, time and space variable are: ðTi T1 Þ, L, and L2 =a 2, respectively. The temperature field T ðx; tÞ and the heat flux density field jðx; tÞ in a large plate modeled by Eqs. (1)–(3) have been accurately determined by the code PSPICE. Certainly, a sensitivity analysis of the number of cells forming the computational domain is a mandatory step in numerical calculus. The numerical experiments revealed that to achieve satisfactory results the optimal number of cells must be N ¼ 100. This choice insures that the intrinsic errors are nearly zero in all cells occupying the computational domain. Fig. 3 shows the trend exhibited by the dimensionless temperature ðT Tw Þ= ðTi Tw Þ as a function of the dimensionless time T ¼ t=ðL2 =aÞ at different locations inside the large plate. The family of curves consists of ten curves, the uppermost curve represents the temperature at the centerplane of the plate ðX ¼ 0Þ and the lowermost curve corresponds to the temperature at a location near the surface of the plate, X ¼ 0:9. Fig. 4 is the companion figure to Fig. 3 which displays the response of the dimensionless instantaneous heat flow rate qL=kAðTw Ti Þ to continuous changes in the dimensionless time h ¼ t=ðL2 =aÞ. In the family of 11 curves, the uppermost curve identifies the instantaneous heat flow at the surface of the plate ðX ¼ 1Þ, whereas the lowermost curve symbolizes the instantaneous heat flow at the midplane of the plate ðX ¼ 0Þ.
2
R2 =as , has been identified as the heat diffusion time in [7].
112
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
Fig. 4. The dimensionless instantaneous heat flux history as a in a large plate parameterized by the dimensionless space variable.
Note that the union of NSM with the code PSPICE may be extended to unsteady heat conduction in long cylinders and spheres requires the specification of two unequal resistors in the group of cells. This modification is basically due to the differences in surface areas of each half of the cells. Likewise, the numerical solution of the heat conduction equation in multidimensional regular bodies do not pose any obstacle to the duo NSM and the code PSPICE either. Additional information about the use of the electro-thermal analogy, the NSM and PSPICE for solving the heat conduction equation may be found in Alhama [8]. 7. Conclusions The coupling of the electro-thermal analogy and the NSM has proved to be a versatile numerical procedure for solving the heat conduction equation with the help of the computer code PSPICE. The spatio-temporal temperatures and the spatio-temporal heat fluxes in large plates may be predicted with good accuracy simultaneously in the entire time domain ð0 < s < 1Þ. References [1] W.F. Ames, Numerical Methods for Partial Differential Equations, third ed., Academic Press, San Diego, CA, 1992.
F. Alhama et al. / Appl. Math. Comput. 162 (2005) 103–113
113
[2] O.A. Liskovets, The method of lines (review), Different. Equat. 1 (1965) 1308–1323. [3] M. Jakob, in: Heat Transfer, vol. 1, McGraw-Hill, New York, 1949. [4] V. Paschkis, H.D. Baker, A method for determining unsteady-state heat transfer by means of an electrical analogy, Trans. ASME 64 (1942) 105–112. [5] L.W. Nagel, SPICE: a computer program to simulate semiconductor circuits Memo No. UCB/ ERL M520, Electronic Research Laboratory, University of California, Berkeley, CA, 1975. [6] PSPICE, version 6.0, 1994. Distributor: Microsim Corp., Irvine, CA. [7] P.S. Lykoudis, Non-dimensional numbers as ratios of characteristic times, Int. J. Heat Mass Transfer 33 (1990) 1560–1570. [8] F. Alhama, Estudio de las Respuestas Termicas Transitorias en Procesos no Lineales de Conducci on de Calor Mediante el Metodo de Simulaci on por Redes, Doctoral Thesis (in Spanish), Universidad de Murcia, Cartagena, Murcia, Spain, 1998.