- Email: [email protected]
Modified ballistic–diffusive equations for transient non-continuum heat conduction
International Journal of Heat and Mass Transfer 83 (2015) 51–63
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Modified ballistic–diffusive equations for transient non-continuum heat conduction Yujie Zhang, Wenjing Ye ⇑ Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology, Hong Kong
a r t i c l e
i n f o
Article history: Received 1 August 2014 Received in revised form 4 November 2014 Accepted 5 November 2014
Keywords: Non-continuum phonon transport Boltzmann transport equation Ballistic–diffusive approximation
a b s t r a c t For the modeling of non-continuum phonon transport, the phonon Boltzmann transport equation is regarded as an accurate model when the wave effect is negligible. However it is difficult and computational intensive to solve. The ballistic–diffusive approach (BDE) proposed by Chen (2001) greatly simplifies the solution procedure and reduces the computational cost by modeling the ballistic and diffusive parts separately and using the first-order spherical harmonic function to approximate the diffusive intensity. The accuracy of the BDE remains to be improved particularly near the boundary and at large time scales. In this paper, a deficiency in the BDE is identified, and a new physically sound boundary model is proposed, which leads to a new set of ballistic–diffusive equations. Several benchmark problems are employed to validate the performance of the modified ballistic–diffusive approach, and results indicate the new approach is much more accurate than the origin BDE and yet retains the same level of efficiency. Hence it is an effective method for the modeling of transient non-continuum heat conduction. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The advent of microfabrication technology has enabled the fast development of micro devices based on semiconductor materials. Many of these devices involve heat transfer [1,2] and the need for the accurate modeling of phonon transport, which is often the dominant transport mechanism, is of particularly importance to the development of micro systems. It is well known that both Fourier law and Cattaneo equation cannot satisfactorily predict heat transfer at small spatial and time scales due to the inherent defects such as the infinite propagating speed in Fourier law and the smaller speed predicted in Cattaneo equation [3]. The Boltzmann transport equation (BTE) [4], on the other hand, is regarded as an accurate description of phonon transport in the regime when the wave effect of phonon is negligible [3]. Unfortunately, the Boltzmann transport equation is very difficult to solve as it involves spatial, time and momentum variables. The complexity of the equation has greatly limited its applications in realistic devise/systems. By dividing the phonon distribution into ballistic and diffusive parts and formulating a transport equation for each part, Chen [5,6] derived the ballistic–diffuse equations (BDEs) based on the BTE under the relaxation time approximation. The ballistic transport equation is solved by the ray-tracing ⇑ Corresponding author. E-mail address: [email protected] (W. Ye). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.11.020 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.
method. For the diffusive component, spherical harmonic expansion is applied and the P1-approximation, first introduced by Cheng in the modeling of the radiation heat transferred [7,8], is employed, which leads to a flux equation with a form of Cattaneo equation, and a heat conduction equation that involves only the time and spatial variables. As such, BDEs are much simpler to solve compared to the BTE and can be easily incorporated into engineering software for the modeling and design of realistic systems. A few benchmark problems were used to demonstrate the accuracy of the BDEs. Compared to the Fourier law and Cattaneo equation, the BDEs predict much accurate results when benchmarked with the BTE solutions. The agreement between the results obtained from BDEs and the BTE is satisfactory inside the domain and within a relatively short time period. Near the boundary and at long periods of time, discrepancies in the results particularly the temperature profile can be observed in a few cases, which indicates some deficiency in the method. Indeed if one examines carefully the P1-approximation and the proposed boundary conditions, there seems to be an inconsistence between the two, which could be the reason that leads to inaccurate results near the boundary. Higher order sphere harmonic approximation (PN approximation) could be used to improve the results as was done for radiative heat transfer in [9]. However due to the lack of clear physical meaning for the higher moments of distribution function, it is generally difficult to define the relevant boundary conditions. In 2011, Mittal and Mazumder [10] proposed generalized ballistic–diffusive
52
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
equations (GBDEs). The basic equations for the ballistic and diffusive parts are the same as that of Chen’s BDEs. The main difference lies in the definition of temperature. In the BDEs, temperature is defined as a measure of the local internal energy, while in the GBDEs, the Bose–Einstein distribution is employed to link the internal energy with the temperature, hence local thermal equilibrium is assumed. Another difference is in the method used to solve the ballistic transport equation. In the GBDEs, to deal with the complex geometry, the ballistic transport equation was solved it by the discrete ordinate method (DOM) and the control angle discrete ordinate method (CADOM) instead of the ray-tracing method. The boundary conditions for the diffusive component in the GBDEs are the same as that in the BDEs. In this paper, a modified BDE (MBDE) method is presented in which the P1-approximation is applied piecewisely and a new boundary model for the diffusive part, which is consistent with the piecewise P1-approximation, is proposed. Based on this, a set of new ballistic–diffusion equations is derived. The accuracy of the MBDE method is tested on several benchmark problems including those used to validate previous BDE methods.
2.1. Boltzmann transport equation (BTE) Under the relaxation time approximation, the phonon Boltzmann transport equation [4] can be written in the phonon intensity form as [11]
@Ix I x I 0x þ v g~ ; s rI x ¼ @t sx
ð1Þ
where Ix ¼ 41p hxv g fDðxÞ, f is the distribution function for phonon, h is the Planck constant divided by 2p, x is the phonon frequency, D is the density of states, vg is the group velocity of phonon, sx is the phonon relaxation time, ~ s is phonon propagation direction and I0x is the local equilibrium phonon intensity. Macroscopic quantities can be obtained by proper integration of the phonon intensity, for example, temperature T and heat flux ~ q can be calculated as
T¼
Ix dXdx=C v g ;
@Ib Ib þ v g~ s rI b ¼ @t s
ð3Þ
@Im Im I 0 þ v g~ s rI m ¼ @t s
ð4Þ
where Ib is the ballistic part for phonon intensity and Im is the diffusive part. Eq. (3) can be solved analytically in 1-D problems or by ray-tracing method in 2-D and 3-D cases, similar to the method employed in Chen’s work [5,6]. Eq. (4) is solved by the numerical methods used to solve the original BTE. We denote this approach as the modified BTE approach in this paper. 2.3. ballistic–diffusive equations (BDE)
2. Methods
Z Z
the finite volume method for phonon transport. This idea has also been used to derive the modified discrete ordinate method and is proven to be very efficient in reducing false scattering and ray effect, especially when there is a wall temperature discontinuity [14–16]. The essential idea of the method is the same as that of the BDEs, that is, by dividing the phonon distribution function into a ballistic part and a diffusive part, and the transport equations for the two parts are formulated as [13]:
~ q¼
Z Z
Ix~ sdxdx
ð2Þ
where C is the specific heat capacity. The temperature in this work is defined as a measure of the local internal energy since the phonon distribution could deviate far from its local equilibrium distribution [4]. The focus of this work is on the development of the BDE method as an efficient alternative to the BTE. Hence for simplicity, a temperature-independent specific heat capacity value is used, and the gray model without sources is considered. In the subsequent sections, the subscript x is omitted and the relevant quantities refer to the average values over the frequency domain. 2.2. Modified Boltzmann transport equation (MBTE) Numerical simulation of the BTE faces two challenges. One is the ray effect caused by the discontinuity in wall temperature and resulted from the angular discretization. The other is the false scattering caused by the sharp gradients of the transient temperature profile, and resulted from space discretization. Much effort has been put forth to resolve these two challenges. In 1996, Ramankutty and Crosbie [12] developed a modified discrete ordinate solution of radiative transfer to resolve anomalies in the temperature and heat flux caused by the ray effects. In 2003, Murthy and Mathur [13] combined a ray-tracing technique with the finite volume method to improve the predictive accuracy of
The BDEs were proposed by Chen [5,6] as an efficient approximation of the phonon BTE under the relaxation approximation. The central idea of the ballistic–diffusive approximation is to divide the intensity into a ballistic term and a diffusive term as
I ¼ I b þ Im
ð5Þ
where Ib represents the intensity of phonons originated from the boundary and experiencing out-scattering only,
@Ib Ib þ v g~ s rI b ¼ @t s
ð6Þ
The remaining part, denoted as the diffusive term Im, represents the intensity of carries scattered or emitted from points inside the domain,
@Im Im I 0 þ v g~ s rI m ¼ @t s
ð7Þ
This part of the intensity is more isotropic due to multiple scattering and is then approximated by the first order spherical harmonic function (P1 approximation), that is,
Im ¼ g 0 þ ~ s ~ g1 where g0 is a scalar representing the average value of intensity over all propagation directions, and ~ g 1 is a vector that is related to the heat flux. This approximation together with the transport Eq. (7) leads to a constitutive relation for the diffusive component
s
@~ qm þ~ qm ¼ krT m @t
ð8Þ
qm and Tm are the heat flux and where k is the thermal conductivity, ~ temperature corresponding the diffusive part, defined as
~ qm ¼
Z
Im~ sdX;
Tm ¼
Z
Im dX=ðC v g Þ
ð9Þ
By utilizing the energy conservation equation, the governing equation for Tm is obtained as shown in Eq. (10)
@ 2 T m @T m C s þ @t @t2
! ¼ rðkrT m Þ r ~ qb
ð10Þ
The boundary conditions of the BDEs are prescribed as follows. At the boundary, the emitted phonon intensity contributes solely to the ballistic part. The diffusive heat flux at the boundary is,
53
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
~ qm ~ n¼
Z
s ~ ndX ¼ C v g T m =2 Im~
ð11Þ
~ s~ n<0
Table 1 Differences in the boundary conditions of the original BDE and the MBDE. Outgoing phonon intensity
Incoming phonon intensity
Original BDE
Ib = 0 Im – 0
Ib ¼
MBDE
Ib = 0 Im – 0
Ib ¼
Substituting Eq. (8) into the above equation yields
@T 2K s m þ T m ¼ rT m ~ n; @t 3
ð12Þ
where K is the mean free path for phonon and the kinetic relation k = CvgK/3 is used. The detailed procedure can be found in papers [5,6,17]. 2.4. Modified ballistic–diffusive equations (MBDE) In the ballistic–diffusive approach, all phonons emitted from the boundary and entering into the domain, that is, the incoming phonons as shown in Fig. 1, are ballistic. Hence the incoming diffusive component should be zero. However, according to the P1-approximation, which is applied over the entire solid angle in Chen’s approach, the diffusive phonon intensity of the outgoing phonons, that is, phonons leaving the domain, at the boundary should also be zero. This clearly is not consistent with the fact that generally the outgoing diffusive phonons are nonzero. Hence in the original ballistic–diffusive approach, the incoming diffusive phonon intensity is actually not zero, which conflicts with the fact that all incoming phonons are ballistic. In our new approach, that is, the modified ballistic–diffusive equations, the P1-approximation is not applied in the entire solid angle. Instead, the solid angle space is discretized and the P1-approximation is applied piecewisely within each section. Such a treatment enables a more natural boundary condition to be proposed, which therefore resolves the inconsistence between the P1-approximation and boundary conditions inherent in the original BDE. It can be seen later that for 1-D problems, only two sections of the solid angle are necessary, and for 2-D problems with rectangle domains, four sections are sufficient. Table 1 shows the differences in the boundary conditions of the original BDE and the MBDE. For a concise presentation, the 1-D modified ballistic– diffusive equations are presented first followed by the MBDEs for 2-D problems. 2.4.1. 1-D modified ballistic–diffusive equations (MBDE) In 1-D cases, the diffusive phonon intensity is further separated into two parts according to its propagation direction. Those propagated along the positive x-axis are grouped into Iþ m , while the rest is grouped into I m . The P1-approximation is then applied to each part as shown in Eqs. (13) and (14):
! þ g þ1 ^s
ð13Þ
! Im ¼ g 0 þ g 1 ^s
ð14Þ
Iþm
¼
g þ0
Consequently, the temperature and flux are also split and defined as follows,
R Cv g T þm ¼ ^s~i>0 Iþm dX; R Q þm ¼ ^s~i>0 Iþm^sdX;
R Cv g T m ¼ ^s~i<0 Im dX R Q m ¼ ^s~i<0 Im^sdX
ð15Þ
i is the unit vector along the positive x-axis direction. where ~
Fig. 1. Schematics of the incoming and outgoing phonons at the boundary.
Cv g T 4p Im Cv g T 4p Im
– 0 (incorrect) = 0 (correct)
The total temperature and heat flux can then be calculated by summing up all parts:
T ¼ T b þ T þm þ T m
ð16Þ
Q ¼ Q b þ Q þm þ Q m
ð17Þ
where Tb and Qb are the ballistic temperature and heat flux. The transport equation for Im is
@Im Im I 0 þ v g~ s rI m ¼ @t s
ð18Þ
Integrating the transport equation over the solid angle in each half space leads to
Z ^s~i>0
dX : Cv g
@ þ @ C v g T þm 2pI0 T þ T=2 T m þ v g Q þm ¼ ¼ C v g m @t @x s s ð19Þ
and
Z ^s~i<0
dX : Cv g
@ @ C v g T m 2pI0 T T=2 T m þ v g Q m ¼ ¼ C v g m @t @x s s ð20Þ
It is clear that the energy conversation is satisfied if one combines Eqs. (19) and (20). Multiplying the transport equation with ~ S and integrating the resulting equation over the solid angle within each half space leads to
Z
^sdX :
^s~i>0
@ þ @ 1 Q m þ vg Q þm C v g T þm @t @x 6
¼
Q þm pI0
s
¼
Q þm C v g T=4
s
ð21Þ
and
Z ^s~i<0
@ @ 1 Q m þ vg Q m C v g T m @t @x 6 Q m þ pI0 Q m þ C v g T=4 ¼ ¼
^sdX :
s
s
ð22Þ
Eqs. (19)–(22), together with the ballistic transport equation, form the complete set of the modified ballistic–diffusive equations. With the decoupling of Iþ m and Im , the boundary conditions for the MBDEs can be easily defined and are more natural. The emitted phonons from the boundary are ballistic and hence the incoming diffusive intensity should be zero at the boundary, which leads to the following conditions at two ends:
x ¼ 0 : T þm ¼ 0;
Q þm ¼ 0
ð23Þ
x ¼ L : T m ¼ 0;
Q m ¼ 0
ð24Þ
2.4.2. 2-D modified ballistic–diffusive equations (MBDE) For 2-D problems, the solid angle space can be divided into Nh Nu sections. Within each section, the P1 approximation is employed to approximate the diffusive intensity. For problems with rectangular domains, four sections, that is, X1 : 0 h
54
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
p2 ; 0 u p2 ; X2 : p2 h p; 0 u p2 ; X3 : p2 h p; p2 u p; X4 : 0 h p2 ; p2 u p , are usual sufficient. Note that the range of u is from 0 to p due to symmetry. Let 0 1 g0 i i i i Im ¼ g 0 þ g 1 cos h þ g 2 sin h cos u ¼ ð1; cos h; sin h cos uÞ @ g 1 A, g2 ð1; cos h; sin h cos uÞ Gi , i = 1, . . . , N = Nh Nu, and define
0
R
0
1i
Cv g T m X Im dX C B R i B Q C C ¼ Ai Gi : B Q ,@ mx A ¼ @ Xi Im coshdX A R Q my Xi I m sinhcosudX
Dt ð25Þ
Integrating the transport Eq. (7) over each solid angle section Xi yields
@ i i @ @ 1 I0 A G þ v g Bi Gi þ v g Ci Gi ¼ Ai Gi þ Fi ; @t @x @y s s
ð26Þ
where the detailed expressions for Ai, Bi, Ci, Fi are presented in the appendix. Substituting Eq. (25) into Eq. (26), one obtains the governing equations for the diffusive temperature and heat flux in the 2-D space.
þ ny sin
hi1 þhi 2
cos
hi1 þhi 2
< 0 and nx and ny are
components of the outward normal vector of the boundary. For cases with four sections, one has:
x ¼ 0 : Q 1 ¼ 0;
Q 4 ¼ 0;
x ¼ Lx : Q 2 ¼ 0;
Q3 ¼ 0
ð29Þ
y ¼ 0 : Q 1 ¼ 0;
Q 2 ¼ 0;
y ¼ Ly : Q 3 ¼ 0;
Q4 ¼ 0
ð30Þ
The total temperature and heat flux can be readily obtained after Qi is found.
Cv g T ¼ Cv g T b þ
N X
Q i ð1Þ
ð31Þ
i¼1
Q x ¼ Q xb þ
N X Q i ð2Þ
ð32Þ
i¼1
Q y ¼ Q yb þ
N X Q i ð3Þ
Dx
1
Inþ1 0i;j
s
s
m1
þhm
2
Q nþ1;m i;j
Dt
þ vg
e m Q nþ1;m B i;j
hm1 þhm 2
Q nþ1;m i;j
þ vg
ð33Þ
Dt
m1 þum
2
Q nþ1;m i1;j
s
s
s
2
2
s
þ vg
um1 þum 2
e m Q nþ1;m Q nþ1;m C i;j i;j1
Dy
<0
þ vg
Dx
s
>0
ð36Þ
e m Q nþ1;m Q nþ1;m B iþ1;j i;j
1 þ ¼ Q nþ1;m i;j
Dy
Fm
< 0 and cos
Inþ1 0i;j
e m Q nþ1;m Q nþ1;m C i;jþ1 i;j
ð35Þ um1 þum
þ vg
<0
þ vg
Dx Inþ1 0i;j
Dy
Fm
< 0 and cos
1
e m Q nþ1;m Q nþ1;m B iþ1;j i;j
hm1 þhm
Q nþ1;m i;j
Inþ1 0i;j
e m ðQ nþ1;m Q nþ1;m Þ C i;j i;j1
ð34Þ u
Dx
s
IF cos
þ vg
Fm
> 0 and cos
1 þ ¼ Q nþ1;m i;j
IF cos
ð28Þ 2
IF cos
i 1
The boundary conditions for a rectangular domain are given as follows: hi1 þhi
e m ðQ nþ1;m Q nþ1;m Þ B i;j i1;j
þ ¼ Q nþ1;m i;j
0 1 1 i@ A i i i i i e e where B ¼ B ðA Þ ; C ¼ C ðA Þ ; F ¼ A 0 : 0
if nx cos
þ vg
h
ð27Þ
Q i ¼ 0;
u hm1 þ hm m1 þ um > 0; and cos >0 2 2
þ ¼ Q nþ1;m i;j
Dt
@ i @ ei i @ ei i 1 I0 B Q þ vg G þ vg C Q ¼ Q i þ Fi @t @x @y s s i 1
If cos Q nþ1;m i;j
1
i
upstream scheme is adopted. Within each solid angle section Xm(hm 1 6 h < hm, um 1 6 u < um), the discretized equation is listed as follows
e m Q nþ1;m Q nþ1;m C i;jþ1 i;j
Fm
Dy ð37Þ
In Eqs. (34)–(37), the subscripts i, j are the nodal indexes in the x and y directions, n is the index for the time step and m is the angle index. As shown in Eqs. (34)–(37), time derivatives are discretized using an implicit time matching scheme, and hence these equations are solved iteratively. In our calculation, the stopping criterion is that the maximum relative error of Q between the two iteration steps is less than 1 106. The validation of the developed codes for the BTEs and the MBTEs is carried out by comparing the simulations results with the analytic 1-D Heaslet Warming solution for steady problems. For transient problems, our simulation results of both BTE and BDEs are compared with the results of a 1-D problem in reference [6] and a 2-D problem with the emitted temperature boundary condition in references [6,17]. Very good agreements have been obtained for all cases. 4. Results and discussions
i¼1
3. Numerical procedures In this work, the finite difference method and the upwind scheme are employed to solve the BTE and the Gauss–Legendre quadrature is used to calculate the equilibrium distribution function. To solve for the BDEs, we use the ray-tracing method to calculate the ballistic intensity, and the finite difference method with the central differential scheme to solve the diffusive equations. Detailed procedures can be found in reference [17]. Similar numerical schemes are used to solve the modified BTEs. For the MBDEs, the ballistic intensity is again calculated using the ray-tracing method. To solve Eq. (27), the finite difference
In this section, simulation results of a few test cases calculated from the modified ballistic–diffusive equations are presented and compared with results obtained from the original BDEs, the BTE and the MBTEs. These problems have been studied in the previous work [5,6,17] and are the model problems for a wide range of applications. Since the focus of this work is on the accuracy and efficiency of the proposed MBDEs, non-dimensional parameters defined in Eq. (38) are used in all calculations:
T T0 Tb ¼ ; T1 T0 t x t ¼ ; x ¼ ; L s
b ¼ Q
Q C v g ðT 1 T 0 Þ
y y ¼ ; L
Kn ¼
K vgs ¼ L L
ð38Þ
55
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63 1
0.3
+
BDE + + MBDE + + BTE + + + MBTE + + + + + 0.6 + + + + + + + + + + + + t*=100 + + + + + 0.4 + + + + + + + + + + t*=10 + + + + + + 0.2 * + t =1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ ++ ++ + ++ ++ ++ ++ ++ ++ 0 0.8 +
0
0.2
0.4
0.6
0.8
0.25
Non-dimensional Heat Flux
Non-dimensional Temperature
+ +
+
+ t*=1 + 0.1 +
++ * + t =10 + + + + 0.05 + + ++ + + + ++ 0
1
+
+
+
+
+
+
+
+
+
+
+
+
+
+
t =10
+
0.2
+
+
t*=1
+
+
+
+
+
+ t =0.1 *
+
+
+
+
+
0
+
0.2
+
+
+
+
0.4
+
+
+
+
+
+
0.6
+
+
+
+
+ + + + ++
0.8
0.2
0.4
0.6
0.8
1
Non-dimensional Coordinate
BDE MBDE BTE MBTE
+ +
+
++
+
+
+ +
+ +
+
+
+
+ +
+ + ++
+
t*=10
+ +
+
+
+
+
+
+
+
+
+
+
+
*
t =0.1
+
+ +
+
0
+
+ t =1 + *
0.1
0
1
+
0.2
+
+
+
+
0.4
+
+
0.6
+
+
+
+
0.8
+ + +
1
Non-dimensional Coordinate
0.3
+ + + +
+ +
+ +
+
*
+ +
+ +
+
t =1
+
+
+
+
+
+
+
+
+
+
+
+ t =0.01 + *
+
+
+
0
+
0.2
+
+
+
0.4
+
+
+
0.6
+
0.2
+
+
+
+
+
0.8
Non-dimensional Coordinate
+ + ++
++ + +
+ +
+
+
+
+
+
+
+
+
+
t =1
+ + t*=0.1
+
+ + + +
+
+ + +
+ +
*
t =0.01
+ +
+ 1
0
+
*
+
0.1
BDE MBDE BTE MBDE
+
+ +
+
0.05
+
+ +
+
t =0.1
+
+ ++ +
0.15
*
+ 0.2
0.25 + +
BDE MBDE BTE MBTE
0.8
Non-dimensional Heat Flux
Non-dimensional Temperature
+ + + + + + + + + + + + ++ ++ ++ ++ ++ ++
+
+
1
0
+
0.2
0.05
+
0.4
0
+
+
0.15
*
+
+
Non-dimensional Heat Flux
Non-dimensional Temperature
+
0.25
BDE MBDE BTE MBTE
+
0.6 +
0.6
+
0.3
+
0
+
*
t =100
+ +
Non-dimensional Coordinate
1
0.4
+
0.15 +
Non-dimensional Coordinate
0.8
BDE MBDE BTE MBTE
0.2 +
0
+
+
0.2
+
+
+
0.4
+
+
+
0.6
+
+
+
0.8
Non-dimensional Coordinate
+
+ + +
1
Fig. 2. Comparison of the temperature and the heat flux profiles obtained from the BDE, the MBDE, the BTE, and the MBTE at different time and phonon Knudsen number.
where L is the characteristic length, T0 and T1 are the temperatures at the boundary. 4.1. Case 1 Heat conduction inside a 1-D thin film of thickness L is simulated. The film is initially at a uniform temperature T0. At time t,
one of its surfaces is subjected to a phonon flux at a temperature of T1. This problem was studied using the original BDEs and results were presented in [5,6]. A fine spatial mesh with Nx = 800 is used for the discretization of the BTE to accurately capture the phenomenon and a relatively coarse mesh with Nx = 200 is used for the simulations with all the other methods in which the ballistic component can be calculated analytically. For the angular
56
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63 0 .2
0.0 6 *
Non-dimensional Temperature
K n= 1 , t = 1
0 .1 5
Non-dimensional Heat Flux
D iffusive P a rt *
K n= 1 , t = 1
BD E M BD E 0 .1
P ositive
0 .0 5
BD E M BD E
0.0 4 P ositive
0.0 2
0
D iffusive P a rt
-0.0 2
N e ga tive
N e ga tive
0
-0.0 4 0
0 .2
0.4
0 .6
0 .8
1
0
0.2
N on-dim e nsiona l C oordina te
0.4
0 .6
0 .8
1
N on-dim e nsiona l C oordina te
Fig. 3. Comparison of the diffusive temperature and heat flux obtained from the BDE and the MBDE; Positive indicates the temperature or the heat flux calculated from the phonon intensity along the positive x-axis, Iþ m ; Negative refers the temperature or the heat flux corresponding to Im .
0.52 0.51 0.5+
+ + 10
-2
10
-1
10
t
0.72 0.7
Temperature
0.68
1-D Kn=1
0.66
BDE MBDE BTE MBTE
0.64 0.62
+
0.6 0.58 0.56 0.54 0.52 0.5 -2 10
+
+ + + 10-1
+
t
+
0.24
0.235
0.23
0
10
+ + + + + + + + + + + + + + + +
0.74
1
10
+
0.7 0.65 0.6 0.55
+ + +
0.5 -1 10
0
t
+ + +
0.2
100
10
0.14 10-2
1
10
1
+
++
+ ++
+
10-1
+ ++ ++
0.2
0.15
+
+
10
-1
0.5+
10
10
1
-2
10
-1
10
t
0.516
+ + + + + + + + + + + + + + + + + + +
t
100
0
10
1
0.514
1-D Kn=1
0.512
BDE MBDE BTE MBTE
0.51
+
0.508 0.506 0.504 0.502
+
0.5
++
10-2
101
+ + + 10-1
+
t
+
+ +
+ + + + + + + + + + + + + +
++
100
101
0.515
+
+
+
+ + ++ + + + + + + + + + + + + + + + + + + + + + + + 1-D Kn=0.1 + + BDE +
+
0.05 2
0.505
0.518
+ + + + + + + 1-D Kn=0.1 + BDE + + MBDE + BTE + MBTE + + + + + +
0.1
10
0
0.51
++++
0.52
1-D Kn=1 BDE MBDE BTE MBTE +
0.18
10
t
0.16
+ + + 1-D Kn=0.1 + BDE + + MBDE + BTE + + MBTE + + + + + + + + + + + + + + 10
-1
0.22
Q
Temperature
0.75
+ + 1-D Kn=10 + + BDE + MBDE + BTE + MBTE + + + + + + + + + + + + + ++ +++++++ 10
0.24
0.25
0.8
-2
++
+ + +
0.9 0.85
+
0.245
++ + + + + + + 1-D Kn=10 + BDE + MBDE + + BTE + MBTE + + + + + + + + + + + ++ ++ + + ++
Q+T/2
0.53
+ + ++ +
Q+T/2
0.54
0.25+
++++
10
0
t
10
1
0.51
Q+T/2
Temperature
0.55
++
Q
0.56
+ + + + + + + + 1-D Kn=10 + + BDE + MBDE + BTE + MBTE + + + + + + + + + + ++ ++
Q
0.57
0.505
+ +
++
++
+ + +
10
+
+
+ +
+
MBDE BTE MBTE
0.5 2
10
-1
10
0
t
10
1
10
2
Fig. 4. Comparison of the surface temperature, heat flux and the sum of the heat flux and half of temperature at x = 0 as functions of time at different Knudsen numbers.
discretization, Nh = 160 is used for the BTE, the MBTE and the calculations of the ballistic component in the BDE and the MBDE. The time step is chosen to be 0.05, 0.005 and 0.0005 for Knudsen
number of 0.1, 1.0 and 10.0 respectively. Both the temperature and the heat flux profiles at different time and different Knudsen number are calculated using the modified BDEs and plotted in
57
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
Fig. 5. Schematic illustration of the problem domain, the boundary conditions and the local coordinate system of the 2-D problem.
0.6 +
0.25
+ + + + + + + + +
0.4
0.3
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+ 0.5 +
+
0.2
+ + +
0.1
0
+ +
0
+
+
0.2 +
+ 0.15
+
0.5
+ + + + + + + + + + + 1
1.5
+ + +
0.1
+ +
+ + +
0.05
+ + +
0
2
0
+ + + + + + + + + + +
0.5
1
1.5
2
Non-Dimentional Coordinate 0.2
0.8
+ + 0.6 + + + + + + 0.4 + + + + + + 0.2
0
+
+
+
+
++
BTE MBTE BDE MBDE-2*2 MBDE-4*4
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+ +
+ +
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
+ +
Non-Dimentional Coordinate
0
+
++ + + + + + + + + + ++ + + 2
4
Non-Dimentional Coordinate
+ + + 0.15 + + + + + + 0.1 + + + + + + 0.05 +
0
0
+
+
+
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+ + + + + + + + + + + + + + + 2
4
Non-Dimentional Coordinate
Fig. 6. Temperature and the y-direction heat flux along the centerline (x = Lx/2) at Kn = 1, t⁄ = 1, 10.
Fig. 2. The comparison with the reference solutions obtained from the BTE and the MBTEs indicates that the original BDEs predicts both the temperature and heat flux quite well when the Knudsen number is large and the time scale is small. However in both the intermediate and small Knudsen number regimes in which the diffusive component plays an important role, when the time scale is large, the temperature profile within the Knudsen layer deviates significantly from the reference solutions. Although the heat flux profile still matches well with the reference profiles, small discrepancies can be observed near the boundary. Overall, the prediction of the heat flux is more accurate than that of the temperature.
On the contrary, the results obtained from our modified BDEs agree with the reference solutions very well in all cases. The differences in the results from the BDEs and the MBDEs are largely caused by the boundary conditions employed in the two methods. As previously mentioned, the boundary conditions used in the BDE method are inconsistent with the P1-approximation for the diffusive component. The new boundary conditions employed in the MBDEs do not have such an inconsistence, which lead to a much improved accuracy particularly near the boundary of the problem domain. One can further understand this point by examining closely the differences in the two sets of results.
58
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.3
+
0.5 +
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
+ 0.4 +
+
+ + + +
+
0.3
+ +
+
+
0.2
+
+
+
t*=0.1
+
0.1
0
Non-Dimentional Heat Flux
Non-Dimentional Temperature
0.6
+
1
+
+ 2
+
3
+
+ 4
+ 5 /0
1
+
+ + +
+
+ +
+ +
+
+
0.1
+
+
+ t*=0.1
*
t =1.0
+
+
+ +
+
+
BTE BDE MBTE MBDE-2*2 MBDE-4*4
0.2 +
+
+
0
*
t =1.0
+
+
+ +
2
+
3
+
4
+
+ 0
5
0
+ 1
+
+
2
+ 3
+
+
4
5 /0
1
2
+
3
+ 4
+ 5
Non-Dimentional Coordinate
Non-Dimentional Coordinate
Fig. 7. Temperature and the y-direction heat flux along the centerline x = Lx/2 at Kn = 10.
0.3
0.6
t*=1
0.5
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
0.4
t*=1
+++ + + + + + + +
+++
Non-Dimentional Heat Flux
Non-Dimentional Temperature
0.7
+ + + + + +
0.3 0.2
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+ 0.2
+ + +++ + + +
++
+++
+ + +
0.1
0.1 0+ + + + + + + + + + 0 2 4
+ +
+
0+ + + + + + + + + +
+ + + + + + + + +
6
8
10
0
2
++ +
+
4
+ + + + + + + +
6
8
10
Non-Dimentional Coordinate
Non-Dimentional Coordinate 0.8 *
0.6
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
t*=10
+++ + + + + + + +
++ ++
Non-Dimentional Heat Flux
Non-Dimentional Temperature
t =10
+ + + + +
0.4
0.2
+ ++ + + + + + + + 0
1
2
3
4
+ +
+ + 5
0.2
BTE MBTE + + BDE + + + MBDE-2*2 MBDE-4*4
+
0+ + + + + + + + +
+ 6
+ + + + + + + + + 7
Non-Dimentional Coordinate
8
9
10
0
1
2
3
+ 4
+ +
+ + 5
+ + + + + + +
+ 6
++
+ + + + + + +
+ + + + + + + ++ 7
8
9
Non-Dimentional Coordinate
Fig. 8. Comparison of the temperature and the y-direction heat flux along the line y = 0 and at t⁄ = 1, 10. The Knudsen number is 1.0.
10
59
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.3
+
0.5
+
+
*
t =0.1 0.4
0.3
+
+
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
++
Non-Dimentional Heat Flux
Non-Dimentional Temperature
0.6
+
+
0.2
0.2
+ +
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
++++
+
0.1
0.1
0+ 0
+
+
+
2
+
+
+
+
4
6
+
+
+
8
+
+
0+ 0
+
10
+
+
+
+
+
2
+
+
+
+
4
+
6
+
+
8
+ 10
Non-Dimentional Coordinate
Non-Dimentional Coordinate 0.6
0.3
+
0.5
+
*
t =1 0.4
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
0.3
+
++
++++
+
+
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+
t*=0.1
+
0.2
+ t*=1 0.2
+
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
+
++++
+
0.1
0.1
0+ 0
+
+
+
2
+
+
+ 4
+
+
+
+
6
8
+
+
+
10
Non-Dimentional Coordinate
0+
+
+
0
+
+
2
+
+
4
+
+
+
+
+
6
+
8
+ 10
Non-Dimentional Coordinate
Fig. 9. Comparison of the temperature and the y-direction heat flux along the line y = 0 and at t⁄ = 0.11. The Knudsen number is 10.
Fig. 3 plots the temperature and heat flux calculated from the diffusive part of the intensity. Those predicted by the BDEs, particularly the results corresponding to Iþ m , exhibit weak but artificial wave fronts. As discussed in the reference [6], these wave fronts are non-physical. On the other hand, results calculated from the MBDEs do not have such non-physical wave fronts and are more accurate. Fig. 4 presents the boundary temperature and heat flux at x = 0 as functions of time. It is evident that the MBDEs predict both the transient temperature and heat flux at the boundary much more accurately than the BDEs, especially at small Knudsen numbers when the contribution from the diffusive part plays an important role. As shown in the third column of Fig. 4 where the quantity Q + T/2 is plotted, the MBDEs can capture the time dependent trend very well, while the results obtained from the BDEs remain constant. The constant value is due to the inherent assumption in the BDE approach that Q + T/2 is equal to half of the wall temperature. It is also noted that the results predicted by the BTE are not accurate at the beginning, and the reason is that the upstream scheme adopted in the BTE is unable to resolve the rapidly changed temperature at the beginning. By applying the MBTEs, such a shape change can be well captured by the ballistic part. The remaining part, that is the diffusive part, is smooth and can then be accurately simulated.
4.2. Case 2 A two-dimensional problem similar to the one in [6,17] is solved using the MBDEs. The problem domain, the boundary conditions and the local coordinate system are illustrated in Fig. 5, and the initial temperature everywhere inside the domain is T0. The dimensions of the domain are chosen to be:
Lh ¼ L;
Lx ¼ 10L;
Ly ¼ 5L
where L is the size of the heater region and served as the characteristic length based on which the Knudsen number is calculated. In this case, the same non-uniform grid system as described in [17] is used for the simulation of the BTE. Near the heated region, that is, 4.5 6 x 6 6.5, 0 6 y⁄ 6 1, a uniform mesh of Nh = 60 is used and outside the heated region, the grid is refined using a geometric series. A total mesh of Nx Ny = 228 140 is used. For the simulations with other methods, a uniform mesh is used in the entire simulation domain. For Kn = 1, Nx Ny = 280 200 is used for the MBTE and the MBDE, Nx Ny = 400 200 is used for the BDE and Nh Nu = 100 20 is used for all the other methods. For Kn = 10, Nx Ny = 200 140 is used in the simulations of the MBTE, the BDE and the MBDE; Nh Nu = 120 20 is used for the simulations of the BTE and the MBTE; Nh Nu = 120 40 is used for the calculations of the ballistic components of BDE and MBDE. The time
60
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.0008
0.0006 *
+
0.0006
t =0.05
MBTE BDE MBDE-2*2 MBDE-4*4
0.0004
++ + + + + + +
0.0002
+
+
+ +
+ + + +
+ 0+ + + + + + + + + + 0 1 2 3 4
Non-Dimentional Heat Flux
Non-Dimentional Temperature
*
t =0.05
5
+
MBTE BDE MBDE-2*2 MBDE-4*4
0.0004
++ + + + + + +
0.0002
+ + +
+ + + + + + + + +
6
7
8
9
8
9
10
0.2
+
BTE + + MBTE BDE + + MBDE-2*2 MBDE-4*4 +
+
0.15
+
+
+
+ + Non-Dimentional Heat Flux
Non-Dimentional Temperature
7
+ + + + +
t*=0.1
+
+
+
+ +
0.1
+ + + +
0.05
+ +
+
+
2
+
+
+
t*=0.1
0.15
+ BTE + MBTE BDE MBDE-2*2 + MBDE-4*4
+ 0.1
+
4
6
+
+
+
+
+ +
+ + +
0.05
+ +
+
+
+
+
+
+ +
+
+
Non-Dimentional Coordinate
8
+
+
0+ 0
10
+
+
+
2
+
+
+ +
4
+
6
8
+
+
Non-Dimentional Coordinate
+
10
0.3
0.3
+ + +
t*=1
BTE + MBTE BDE + MBDE-2*2 MBDE-4*4 +
+
+ ++ ++ + +
+ +
+
+
+
+
+
+
+ + 0.1
+ +
+ + +
2
BTE MBTE + BDE + MBDE-2*2 MBDE-4*4
+
0.2
+ + 0.1
+
+ + +
+
+
6
8
+ 10
Non-Dimentional Coordinate
0+ 0
2
+
+
+
+
+ + + +
+
+
+
+ + +
+
+ 4
+ ++ ++
+
+
+
+
+
+
+
+ + + + +
+ + +
t*=1
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+ + + + + + + + +
6
0.2
0.25
0+ 0
5
Non-Dimentional Coordinate
0.3
0.2
+ + + +
+
0+ + + + + + + + + + 0 1 2 3 4
10
Non-Dimentional Coordinate
0+ 0
+
+
+
+ 4
+ + +
6
8
+ 10
Non-Dimentional Coordinate
Fig. 10. Comparison of the temperature and the y-direction heat flux at y = 0.5, Kn = 10.
step is chosen to be 0.05 and 0.005 for Knudsen number of 1.0 and 10.0 respectively. Temperature and heat flux of the entire domain are calculated at various time instants and Knudsen numbers using the four
methods described in Section 2. In the modified ballistic–diffusive approach, the division of the solid angle space is done at two levels, 2 2 and 4 4 sections, to check the convergence. Fig. 6 shows the calculated five sets of the temperature and the y-direction heat flux
61
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
4.3. Case 3 The last example considered is the problem of radiative heat transfer in an isotropically scattering media which was studied previously [18]. The problem domain and the local coordinate system are the same as that of case 2. The media is initially at a uniform temperature T0, and then the temperature of the side y = 0 is raised to and maintained at T1. The dimensions and the Knudsen number are:
Lx ¼ L;
Ly ¼ 10; Lx
K Kn ¼ ¼ 10 L
0.6
0.5 X Crosbie and Schrenker C.Cai BTE MBDE-2*2 BDE MBTE
Temperature
X X
X
0.4
X 0.3
X X
0.2
X X X
0.1
0
X
0
0.5
1
y
1.5
X
2
Fig. 11. Steady-state temperature profile along the centerline x = Lx/2. Data indicated by dots and crosses are taken from the papers [18,19].
0.53 BTE MBDE-2*2 BDE MBTE
0.525
Temperature
along the centerline of the domain (x = Lx/2) at Kn = 1 and t⁄ = 1, 10. Similar to that of case 1, the results predicted by the BDEs are less accurate near the boundary, while the MBDEs can capture both the temperature and the heat flux quite well along the entire line. As the Knudsen number increases, the ballistic component becomes dominant and consequently the results predicted by the BDEs agree better with the reference solutions obtained from the BTE and MBTE as demonstrated in Fig. 7 where the temperature and the y-direction heat flux at Kn = 10 are plotted. The temperature and the heat flux profiles along the lines of y = 0 and y = 0.5 are also examined and plotted in Figs. 8–10. The Knudsen number in Fig. 8 is 1, while in Figs. 9 and 10 it is 10. Although the overall profiles obtained from the BDEs agree well with the reference solutions obtained from the MBTEs, discrepancies can be observed in the peaks of the profiles as shown in the insets of the plots particularly in Fig. 8 when Kn = 1, in this case the MBDE is much accurate than the original BDE. For the cases of Kn = 10, those results obtained from the MBDEs are generally more accurate than those from the original BDEs except the profiles along y = 0.5 and at t⁄ = 0.1 in which case the original BDEs perform better than the MBDEs. However, as shown in Fig. 10, at t⁄ = 0.05 and t⁄ = 0.1, the results from the BDEs are worse than that of the MBDEs, indicating that t⁄ = 0.1 is simply a special case for the BDEs. Hence even in high Knudsen number cases when the ballistic transport is dominant, the MBDEs perform better than the original BDEs. Also it should be pointed out that there is a small difference between the results obtained from solving the BTE and the MBTE in the region of 4.5 < x < 5.5. Such a difference is caused by the ray effect and false scattering in solving the BTE due to the discontinue wall temperatures [14] at x = 4.5 and x = 5.5. The MBTE approach eliminates these effects and produces more accurate results.
0.52
0.515
0.51
0.505
0.5
0
0.2
0.4
x
0.6
0.8
1
Fig. 12. Comparison of the steady state temperature profiles along the line y = 0 obtained from the BTE, the MBTE, the BDE and the MBDE at Kn = 10.
where L is chosen to be the characteristic length. Again, the steadystate temperature field is calculated using the four methods. In this case, Nx Ny = 40 200 and Nh Nu= 120 40 are used for all the methods. Results along the centerline (x = Lx/2) are presented in Fig. 11. Reference solutions extracted from reference [18] are also plotted for comparison. In this case, all four sets of temperature profiles agree well with the reference solutions from [18]. However, for the steady-state temperature profile along the bottom side, plotted in Fig. 12, the original BDE approach predicts a constant profile, which is significantly different from the curved profiles obtained from the BTE, the MBTEs and the MBDEs, which agree with each other quite well even quantitatively. This case further demonstrates the accuracy of the proposed MBDEs.
5. Efficiency To evaluate the efficiency of the MBDE approach, we compared the CPU time consumed by the four methods in solving case 2, which are listed in Table 2. A relatively coarse spatial mesh, Nx Ny = 100 100, is employed in the numerical simulation. To solve the BTEs, the MBTEs and the ballistic transport equation of the BDEs and the MBDEs, the solid angle is discretized with a mesh of Nh Nu = 80 20. To solve the diffusive equations of the MBDEs, the solid angle is divided into 2 2 and 4 4 sections. The time step in the simulation is 0.005 and the Knudsen number is 10. As shown in Table 2, both the original BDE and the MBDE are more efficient than the BTE and the MBTE. The CPU time required in the MBTE is only slightly larger than that of the BTE, indicating that the ballistic transport equation can be solved much more efficiently compared to the diffusive transport equation. Hence the major factor contributing to the efficiency of the ballistic–diffusive approaches lies in the diffusive transport equation. In the BTE and the MBTE, the phonon intensity is solved explicitly, while in the BDE and the MBDE, diffusive equations in terms of temperature and flux are solved and therefore are more efficient. It is also noticed that in the case when the solid angle is discretized into
Table 2 Comparison of computational time for a simulation duration of 200 time steps. Methods
BTE
MBTE
BDE
MBDE(2 2)
MBDE(4 4)
Time (s)
338.7
342.48
88.3
77.4
82.0
62
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.6 MBDE (2 X 2) : NθXNϕ=160 X 40
0.5
Temperature
and a modified ballistic–diffusive approach (MBDE) is proposed by introducing a new boundary model. A new set of ballistic–diffusive equations is derived. Several benchmark problems are employed to validate the accuracy of the new approach, and the results show that the new method is a much more accurate approximation of the BTE than the original BDE while retaining the same level of the efficiency of the original BDE. In particular, the MBDE eliminates the distortion of the temperature profiles predicted by the BDE at the boundary, and provides much more accurate steady-state solutions. In addition, similar to the MBTE, the MBDE is quite effective in eliminating the ray effects and false scattering due to the splitting of the ballistic and the diffusive intensity. Overall, the MBDE is an efficient method for the modeling of transient non-continuum heat conduction with a good accuracy and much less computational time and memory consumption. The extension of the MBDE for 3-D problems should be straightforward and is the subject of the future work.
MBDE-50 X 30 MBDE-100 X 60 MBDE-100 X 60 for Diffusive and 50 X 30 for Ballistic MBDE-200 X 120
0.4
0.3
0.2
0.1
0
0
1
2
y
3
4
5
Fig. 13. Temperature profiles along the centerline x = Lx/2 at Kn = 10.
Conflict of interest None declared.
Table 3 Mesh size and CPU time for the four meshes. Mesh
Parameters
CPU time (s)
Ballistic Part
1 2 3 4
Diffusive Part
Nx
Ny
Nx
Ny
50 100 50 200
30 60 30 120
50 100 100 200
30 60 60 120
Appendix A. Appendix The detailed expressions for Ai, Bi, Ci, Fi in Eq. (26) are listed as follows:
28.3 111.0 39.1 439.6
R
0 B Ai ¼ B @
R Xi
R Xi
four sections, that is, Nh Nu = 2 2, the time used by the MBDE is even less than that of the BDE. Such a reduction in the CPU time is due to the number of iterations of the MBDE, which is smaller than that of the BDE at each time step. In addition to the CPU time, the memory consumptions of the BDE and the MBDE are also much less than that of the BTE and the MBTE due to the fact that the detailed information of the solid angle, which is implicitly included in the macroscopic quantities, does not need to be stored in the ballistic–diffusive approaches. Hence they are more suitable for solving 3-D problems than the transport equation based approaches. It is realized that the computational time can be further reduced without losing accuracy by using a finer mesh for the diffusive part and a coarse mesh for the ballistic part. As the ballistic transport equation is solved in an analytical form, a coarse mesh is sufficient to obtain accurate results. And for the diffusive equations, a fine mesh is required in order to reduce the numerical errors caused by the finite difference scheme. Table 3 shows the mesh size and computational time for the four meshes employed in the simulation. The corresponding temperature profiles along the centerline x = Lx/2 are shown in Fig. 13. In all calculations, 160 40 sections of the solid angle are employed in solving the ballistic transport equation and 2 2 sections are used for solving the diffusive equations. As evident from Fig. 13 and Table 3, among the four sets of meshes, Mesh 3, which adopts a fine mesh for the diffusive part and a coarse mesh for the ballistic part, provides both efficient and accurate solution. 6. Summary In this paper, the inconsistence between the boundary conditions and the P1 approximation in the original BDE is identified
Xi
R
B Bi ¼ B @
R
R
coshdX
sinhcosudX
0
R
dX
Xi
R Xi
Xi Xi
coshdX cos2 hdX
sinhcoshcosudX R
coshdX
R
Xi
cos2 hdX
R R
Xi
sinhcosudX
1
C sinhcoshcosudX C A R 2 2 Xi sin hcos udX
Xi
R R
Xi
sinhcoshcosudX
1
C sinhcos2 hcosudX C A R R R 2 2 2 Xi sinhcoshcosudX Xi sinhcos hcosudX Xi sin hcoshcos udX cos2 hdX
Xi
Xi
cos3 hdX
Xi
0
1 R R R 2 2 sinhcosudX Xi sinhcoshcosudX Xi sin hcos udX B R Xi C R R B C 2 Ci ¼ B Xi sinhcoshcosudX Xi sinhcos2 hcosudX Xi sin hcoshcos2 udX C @ R A R R 2 2 3 2 2 3 Xi sin hcos udX Xi sin hcoshcos udX Xi sin hcos udX
R
0 B Fi ¼ B @
R Xi
R Xi
Xi
dX
coshdX
sinhcosudX
1 C C A
References [1] D.G. Cahill, W.K. Ford, K.E. Goodson, G.D. Mahan, A. Majumdar, H.J. Maris, R. Merlin, S.R. Phillpot, Nanoscale thermal transport, J. Appl. Phys. 93 (2) (2003) 793. [2] S.B. Riffat, X. Ma, Thermoelectrics: a review of present and potential applications, Appl. Therm. Eng. 23 (8) (2003) 913–935. [3] A.A. Joshi, A. Majumdar, Transient ballistic and diffusive phonon heat transport in thin films, J. Appl. Phys. 74 (1) (1993) 31. [4] J.M. Ziman, Electrons and Phonons: the Theory of Transport Phenomena in Solids, Oxford University Press, 2001. [5] G. Chen, Ballistic–diffusive heat-conduction equations, Phys. Rev. Lett. 86 (11) (2001) 2297–2300. [6] G. Chen, Ballistic–diffusive equations for transient heat conduction from nano to macroscales, J. Heat Transfer 124 (2) (2002) 320. [7] P. Cheng, Dynamics of a radiating gas with application to flow over a wavy wall, AIAA J. 4 (2) (1966) 238–245.
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63 [8] P. Cheng, Study of the flow of a radiating gas by a differential approximation, in: DTIC Document, 1965. [9] M.F. Modest, Radiative Heat Transfer, McGraw-Hill, New York, 1993. [10] A. Mittal, S. Mazumder, Generalized ballistic–diffusive formulation and hybrid SN-PN solution of the Boltzmann transport equation for phonons for nonequilibrium heat conduction, J. Heat Transfer 133 (9) (2011) 092402. [11] A. Majumdar, Microscale heat conduction in dielectric thin films, J. Heat Transfer 115 (1) (1993) 7–16. [12] M.A. Ramankutty, A.L. Crosbie, Modified discrete ordinates solution of radiative transfer in two-dimensional rectangular enclosures, J. Quant. Spectrosc. Radiat. Transfer 57 (1) (1997) 107–140. [13] J.Y. Murthy, S.R. Mathur, An improved computational procedure for submicron heat conduction, J. Heat Transfer 125 (5) (2003) 904. [14] P.J. Coelho, The role of ray effects and false scattering on the accuracy of the standard and modified discrete ordinates methods, J. Quant. Spectrosc. Radiat. Transfer 73 (2) (2002) 231–238.
63
[15] M. Sakami, A. El Kasmi, A. Charette, Analysis of radiative heat transfer in complex two-dimensional enclosures with obstacles using the modified discrete ordinates method, J. Heat Transfer 123 (5) (2001) 892–900. [16] M. Sakami, A. Charette, Application of a modified discrete ordinates method to two-dimensional enclosures of irregular geometry, J. Quant. Spectrosc. Radiat. Transfer 64 (3) (2000) 275–298. [17] R. Yang, G. Chen, M. Laroche, Y. Taur, Simulation of nanoscale multidimensional transient heat conduction problems using ballistic– diffusive equations and phonon Boltzmann equation, J. Heat Transfer 127 (3) (2005) 298. [18] J.C. Chai, H.S. Lee, S.V. Patankar, Finite volume method for radiation heat transfer, J. Thermophys. Heat Transfer 8 (3) (1994) 419–425. [19] A.L. Crosbie, R.G. Schrenker, Radiative transfer in a two-dimensional rectangular medium exposed to diffuse radiation, J. Quant. Spectrosc. Radiat. Transfer 31 (4) (1984) 339–372.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Modified ballistic–diffusive equations for transient non-continuum heat conduction Yujie Zhang, Wenjing Ye ⇑ Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology, Hong Kong
a r t i c l e
i n f o
Article history: Received 1 August 2014 Received in revised form 4 November 2014 Accepted 5 November 2014
Keywords: Non-continuum phonon transport Boltzmann transport equation Ballistic–diffusive approximation
a b s t r a c t For the modeling of non-continuum phonon transport, the phonon Boltzmann transport equation is regarded as an accurate model when the wave effect is negligible. However it is difficult and computational intensive to solve. The ballistic–diffusive approach (BDE) proposed by Chen (2001) greatly simplifies the solution procedure and reduces the computational cost by modeling the ballistic and diffusive parts separately and using the first-order spherical harmonic function to approximate the diffusive intensity. The accuracy of the BDE remains to be improved particularly near the boundary and at large time scales. In this paper, a deficiency in the BDE is identified, and a new physically sound boundary model is proposed, which leads to a new set of ballistic–diffusive equations. Several benchmark problems are employed to validate the performance of the modified ballistic–diffusive approach, and results indicate the new approach is much more accurate than the origin BDE and yet retains the same level of efficiency. Hence it is an effective method for the modeling of transient non-continuum heat conduction. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction The advent of microfabrication technology has enabled the fast development of micro devices based on semiconductor materials. Many of these devices involve heat transfer [1,2] and the need for the accurate modeling of phonon transport, which is often the dominant transport mechanism, is of particularly importance to the development of micro systems. It is well known that both Fourier law and Cattaneo equation cannot satisfactorily predict heat transfer at small spatial and time scales due to the inherent defects such as the infinite propagating speed in Fourier law and the smaller speed predicted in Cattaneo equation [3]. The Boltzmann transport equation (BTE) [4], on the other hand, is regarded as an accurate description of phonon transport in the regime when the wave effect of phonon is negligible [3]. Unfortunately, the Boltzmann transport equation is very difficult to solve as it involves spatial, time and momentum variables. The complexity of the equation has greatly limited its applications in realistic devise/systems. By dividing the phonon distribution into ballistic and diffusive parts and formulating a transport equation for each part, Chen [5,6] derived the ballistic–diffuse equations (BDEs) based on the BTE under the relaxation time approximation. The ballistic transport equation is solved by the ray-tracing ⇑ Corresponding author. E-mail address: [email protected] (W. Ye). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.11.020 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved.
method. For the diffusive component, spherical harmonic expansion is applied and the P1-approximation, first introduced by Cheng in the modeling of the radiation heat transferred [7,8], is employed, which leads to a flux equation with a form of Cattaneo equation, and a heat conduction equation that involves only the time and spatial variables. As such, BDEs are much simpler to solve compared to the BTE and can be easily incorporated into engineering software for the modeling and design of realistic systems. A few benchmark problems were used to demonstrate the accuracy of the BDEs. Compared to the Fourier law and Cattaneo equation, the BDEs predict much accurate results when benchmarked with the BTE solutions. The agreement between the results obtained from BDEs and the BTE is satisfactory inside the domain and within a relatively short time period. Near the boundary and at long periods of time, discrepancies in the results particularly the temperature profile can be observed in a few cases, which indicates some deficiency in the method. Indeed if one examines carefully the P1-approximation and the proposed boundary conditions, there seems to be an inconsistence between the two, which could be the reason that leads to inaccurate results near the boundary. Higher order sphere harmonic approximation (PN approximation) could be used to improve the results as was done for radiative heat transfer in [9]. However due to the lack of clear physical meaning for the higher moments of distribution function, it is generally difficult to define the relevant boundary conditions. In 2011, Mittal and Mazumder [10] proposed generalized ballistic–diffusive
52
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
equations (GBDEs). The basic equations for the ballistic and diffusive parts are the same as that of Chen’s BDEs. The main difference lies in the definition of temperature. In the BDEs, temperature is defined as a measure of the local internal energy, while in the GBDEs, the Bose–Einstein distribution is employed to link the internal energy with the temperature, hence local thermal equilibrium is assumed. Another difference is in the method used to solve the ballistic transport equation. In the GBDEs, to deal with the complex geometry, the ballistic transport equation was solved it by the discrete ordinate method (DOM) and the control angle discrete ordinate method (CADOM) instead of the ray-tracing method. The boundary conditions for the diffusive component in the GBDEs are the same as that in the BDEs. In this paper, a modified BDE (MBDE) method is presented in which the P1-approximation is applied piecewisely and a new boundary model for the diffusive part, which is consistent with the piecewise P1-approximation, is proposed. Based on this, a set of new ballistic–diffusion equations is derived. The accuracy of the MBDE method is tested on several benchmark problems including those used to validate previous BDE methods.
2.1. Boltzmann transport equation (BTE) Under the relaxation time approximation, the phonon Boltzmann transport equation [4] can be written in the phonon intensity form as [11]
@Ix I x I 0x þ v g~ ; s rI x ¼ @t sx
ð1Þ
where Ix ¼ 41p hxv g fDðxÞ, f is the distribution function for phonon, h is the Planck constant divided by 2p, x is the phonon frequency, D is the density of states, vg is the group velocity of phonon, sx is the phonon relaxation time, ~ s is phonon propagation direction and I0x is the local equilibrium phonon intensity. Macroscopic quantities can be obtained by proper integration of the phonon intensity, for example, temperature T and heat flux ~ q can be calculated as
T¼
Ix dXdx=C v g ;
@Ib Ib þ v g~ s rI b ¼ @t s
ð3Þ
@Im Im I 0 þ v g~ s rI m ¼ @t s
ð4Þ
where Ib is the ballistic part for phonon intensity and Im is the diffusive part. Eq. (3) can be solved analytically in 1-D problems or by ray-tracing method in 2-D and 3-D cases, similar to the method employed in Chen’s work [5,6]. Eq. (4) is solved by the numerical methods used to solve the original BTE. We denote this approach as the modified BTE approach in this paper. 2.3. ballistic–diffusive equations (BDE)
2. Methods
Z Z
the finite volume method for phonon transport. This idea has also been used to derive the modified discrete ordinate method and is proven to be very efficient in reducing false scattering and ray effect, especially when there is a wall temperature discontinuity [14–16]. The essential idea of the method is the same as that of the BDEs, that is, by dividing the phonon distribution function into a ballistic part and a diffusive part, and the transport equations for the two parts are formulated as [13]:
~ q¼
Z Z
Ix~ sdxdx
ð2Þ
where C is the specific heat capacity. The temperature in this work is defined as a measure of the local internal energy since the phonon distribution could deviate far from its local equilibrium distribution [4]. The focus of this work is on the development of the BDE method as an efficient alternative to the BTE. Hence for simplicity, a temperature-independent specific heat capacity value is used, and the gray model without sources is considered. In the subsequent sections, the subscript x is omitted and the relevant quantities refer to the average values over the frequency domain. 2.2. Modified Boltzmann transport equation (MBTE) Numerical simulation of the BTE faces two challenges. One is the ray effect caused by the discontinuity in wall temperature and resulted from the angular discretization. The other is the false scattering caused by the sharp gradients of the transient temperature profile, and resulted from space discretization. Much effort has been put forth to resolve these two challenges. In 1996, Ramankutty and Crosbie [12] developed a modified discrete ordinate solution of radiative transfer to resolve anomalies in the temperature and heat flux caused by the ray effects. In 2003, Murthy and Mathur [13] combined a ray-tracing technique with the finite volume method to improve the predictive accuracy of
The BDEs were proposed by Chen [5,6] as an efficient approximation of the phonon BTE under the relaxation approximation. The central idea of the ballistic–diffusive approximation is to divide the intensity into a ballistic term and a diffusive term as
I ¼ I b þ Im
ð5Þ
where Ib represents the intensity of phonons originated from the boundary and experiencing out-scattering only,
@Ib Ib þ v g~ s rI b ¼ @t s
ð6Þ
The remaining part, denoted as the diffusive term Im, represents the intensity of carries scattered or emitted from points inside the domain,
@Im Im I 0 þ v g~ s rI m ¼ @t s
ð7Þ
This part of the intensity is more isotropic due to multiple scattering and is then approximated by the first order spherical harmonic function (P1 approximation), that is,
Im ¼ g 0 þ ~ s ~ g1 where g0 is a scalar representing the average value of intensity over all propagation directions, and ~ g 1 is a vector that is related to the heat flux. This approximation together with the transport Eq. (7) leads to a constitutive relation for the diffusive component
s
@~ qm þ~ qm ¼ krT m @t
ð8Þ
qm and Tm are the heat flux and where k is the thermal conductivity, ~ temperature corresponding the diffusive part, defined as
~ qm ¼
Z
Im~ sdX;
Tm ¼
Z
Im dX=ðC v g Þ
ð9Þ
By utilizing the energy conservation equation, the governing equation for Tm is obtained as shown in Eq. (10)
@ 2 T m @T m C s þ @t @t2
! ¼ rðkrT m Þ r ~ qb
ð10Þ
The boundary conditions of the BDEs are prescribed as follows. At the boundary, the emitted phonon intensity contributes solely to the ballistic part. The diffusive heat flux at the boundary is,
53
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
~ qm ~ n¼
Z
s ~ ndX ¼ C v g T m =2 Im~
ð11Þ
~ s~ n<0
Table 1 Differences in the boundary conditions of the original BDE and the MBDE. Outgoing phonon intensity
Incoming phonon intensity
Original BDE
Ib = 0 Im – 0
Ib ¼
MBDE
Ib = 0 Im – 0
Ib ¼
Substituting Eq. (8) into the above equation yields
@T 2K s m þ T m ¼ rT m ~ n; @t 3
ð12Þ
where K is the mean free path for phonon and the kinetic relation k = CvgK/3 is used. The detailed procedure can be found in papers [5,6,17]. 2.4. Modified ballistic–diffusive equations (MBDE) In the ballistic–diffusive approach, all phonons emitted from the boundary and entering into the domain, that is, the incoming phonons as shown in Fig. 1, are ballistic. Hence the incoming diffusive component should be zero. However, according to the P1-approximation, which is applied over the entire solid angle in Chen’s approach, the diffusive phonon intensity of the outgoing phonons, that is, phonons leaving the domain, at the boundary should also be zero. This clearly is not consistent with the fact that generally the outgoing diffusive phonons are nonzero. Hence in the original ballistic–diffusive approach, the incoming diffusive phonon intensity is actually not zero, which conflicts with the fact that all incoming phonons are ballistic. In our new approach, that is, the modified ballistic–diffusive equations, the P1-approximation is not applied in the entire solid angle. Instead, the solid angle space is discretized and the P1-approximation is applied piecewisely within each section. Such a treatment enables a more natural boundary condition to be proposed, which therefore resolves the inconsistence between the P1-approximation and boundary conditions inherent in the original BDE. It can be seen later that for 1-D problems, only two sections of the solid angle are necessary, and for 2-D problems with rectangle domains, four sections are sufficient. Table 1 shows the differences in the boundary conditions of the original BDE and the MBDE. For a concise presentation, the 1-D modified ballistic– diffusive equations are presented first followed by the MBDEs for 2-D problems. 2.4.1. 1-D modified ballistic–diffusive equations (MBDE) In 1-D cases, the diffusive phonon intensity is further separated into two parts according to its propagation direction. Those propagated along the positive x-axis are grouped into Iþ m , while the rest is grouped into I m . The P1-approximation is then applied to each part as shown in Eqs. (13) and (14):
! þ g þ1 ^s
ð13Þ
! Im ¼ g 0 þ g 1 ^s
ð14Þ
Iþm
¼
g þ0
Consequently, the temperature and flux are also split and defined as follows,
R Cv g T þm ¼ ^s~i>0 Iþm dX; R Q þm ¼ ^s~i>0 Iþm^sdX;
R Cv g T m ¼ ^s~i<0 Im dX R Q m ¼ ^s~i<0 Im^sdX
ð15Þ
i is the unit vector along the positive x-axis direction. where ~
Fig. 1. Schematics of the incoming and outgoing phonons at the boundary.
Cv g T 4p Im Cv g T 4p Im
– 0 (incorrect) = 0 (correct)
The total temperature and heat flux can then be calculated by summing up all parts:
T ¼ T b þ T þm þ T m
ð16Þ
Q ¼ Q b þ Q þm þ Q m
ð17Þ
where Tb and Qb are the ballistic temperature and heat flux. The transport equation for Im is
@Im Im I 0 þ v g~ s rI m ¼ @t s
ð18Þ
Integrating the transport equation over the solid angle in each half space leads to
Z ^s~i>0
dX : Cv g
@ þ @ C v g T þm 2pI0 T þ T=2 T m þ v g Q þm ¼ ¼ C v g m @t @x s s ð19Þ
and
Z ^s~i<0
dX : Cv g
@ @ C v g T m 2pI0 T T=2 T m þ v g Q m ¼ ¼ C v g m @t @x s s ð20Þ
It is clear that the energy conversation is satisfied if one combines Eqs. (19) and (20). Multiplying the transport equation with ~ S and integrating the resulting equation over the solid angle within each half space leads to
Z
^sdX :
^s~i>0
@ þ @ 1 Q m þ vg Q þm C v g T þm @t @x 6
¼
Q þm pI0
s
¼
Q þm C v g T=4
s
ð21Þ
and
Z ^s~i<0
@ @ 1 Q m þ vg Q m C v g T m @t @x 6 Q m þ pI0 Q m þ C v g T=4 ¼ ¼
^sdX :
s
s
ð22Þ
Eqs. (19)–(22), together with the ballistic transport equation, form the complete set of the modified ballistic–diffusive equations. With the decoupling of Iþ m and Im , the boundary conditions for the MBDEs can be easily defined and are more natural. The emitted phonons from the boundary are ballistic and hence the incoming diffusive intensity should be zero at the boundary, which leads to the following conditions at two ends:
x ¼ 0 : T þm ¼ 0;
Q þm ¼ 0
ð23Þ
x ¼ L : T m ¼ 0;
Q m ¼ 0
ð24Þ
2.4.2. 2-D modified ballistic–diffusive equations (MBDE) For 2-D problems, the solid angle space can be divided into Nh Nu sections. Within each section, the P1 approximation is employed to approximate the diffusive intensity. For problems with rectangular domains, four sections, that is, X1 : 0 h
54
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
p2 ; 0 u p2 ; X2 : p2 h p; 0 u p2 ; X3 : p2 h p; p2 u p; X4 : 0 h p2 ; p2 u p , are usual sufficient. Note that the range of u is from 0 to p due to symmetry. Let 0 1 g0 i i i i Im ¼ g 0 þ g 1 cos h þ g 2 sin h cos u ¼ ð1; cos h; sin h cos uÞ @ g 1 A, g2 ð1; cos h; sin h cos uÞ Gi , i = 1, . . . , N = Nh Nu, and define
0
R
0
1i
Cv g T m X Im dX C B R i B Q C C ¼ Ai Gi : B Q ,@ mx A ¼ @ Xi Im coshdX A R Q my Xi I m sinhcosudX
Dt ð25Þ
Integrating the transport Eq. (7) over each solid angle section Xi yields
@ i i @ @ 1 I0 A G þ v g Bi Gi þ v g Ci Gi ¼ Ai Gi þ Fi ; @t @x @y s s
ð26Þ
where the detailed expressions for Ai, Bi, Ci, Fi are presented in the appendix. Substituting Eq. (25) into Eq. (26), one obtains the governing equations for the diffusive temperature and heat flux in the 2-D space.
þ ny sin
hi1 þhi 2
cos
hi1 þhi 2
< 0 and nx and ny are
components of the outward normal vector of the boundary. For cases with four sections, one has:
x ¼ 0 : Q 1 ¼ 0;
Q 4 ¼ 0;
x ¼ Lx : Q 2 ¼ 0;
Q3 ¼ 0
ð29Þ
y ¼ 0 : Q 1 ¼ 0;
Q 2 ¼ 0;
y ¼ Ly : Q 3 ¼ 0;
Q4 ¼ 0
ð30Þ
The total temperature and heat flux can be readily obtained after Qi is found.
Cv g T ¼ Cv g T b þ
N X
Q i ð1Þ
ð31Þ
i¼1
Q x ¼ Q xb þ
N X Q i ð2Þ
ð32Þ
i¼1
Q y ¼ Q yb þ
N X Q i ð3Þ
Dx
1
Inþ1 0i;j
s
s
m1
þhm
2
Q nþ1;m i;j
Dt
þ vg
e m Q nþ1;m B i;j
hm1 þhm 2
Q nþ1;m i;j
þ vg
ð33Þ
Dt
m1 þum
2
Q nþ1;m i1;j
s
s
s
2
2
s
þ vg
um1 þum 2
e m Q nþ1;m Q nþ1;m C i;j i;j1
Dy
<0
þ vg
Dx
s
>0
ð36Þ
e m Q nþ1;m Q nþ1;m B iþ1;j i;j
1 þ ¼ Q nþ1;m i;j
Dy
Fm
< 0 and cos
Inþ1 0i;j
e m Q nþ1;m Q nþ1;m C i;jþ1 i;j
ð35Þ um1 þum
þ vg
<0
þ vg
Dx Inþ1 0i;j
Dy
Fm
< 0 and cos
1
e m Q nþ1;m Q nþ1;m B iþ1;j i;j
hm1 þhm
Q nþ1;m i;j
Inþ1 0i;j
e m ðQ nþ1;m Q nþ1;m Þ C i;j i;j1
ð34Þ u
Dx
s
IF cos
þ vg
Fm
> 0 and cos
1 þ ¼ Q nþ1;m i;j
IF cos
ð28Þ 2
IF cos
i 1
The boundary conditions for a rectangular domain are given as follows: hi1 þhi
e m ðQ nþ1;m Q nþ1;m Þ B i;j i1;j
þ ¼ Q nþ1;m i;j
0 1 1 i@ A i i i i i e e where B ¼ B ðA Þ ; C ¼ C ðA Þ ; F ¼ A 0 : 0
if nx cos
þ vg
h
ð27Þ
Q i ¼ 0;
u hm1 þ hm m1 þ um > 0; and cos >0 2 2
þ ¼ Q nþ1;m i;j
Dt
@ i @ ei i @ ei i 1 I0 B Q þ vg G þ vg C Q ¼ Q i þ Fi @t @x @y s s i 1
If cos Q nþ1;m i;j
1
i
upstream scheme is adopted. Within each solid angle section Xm(hm 1 6 h < hm, um 1 6 u < um), the discretized equation is listed as follows
e m Q nþ1;m Q nþ1;m C i;jþ1 i;j
Fm
Dy ð37Þ
In Eqs. (34)–(37), the subscripts i, j are the nodal indexes in the x and y directions, n is the index for the time step and m is the angle index. As shown in Eqs. (34)–(37), time derivatives are discretized using an implicit time matching scheme, and hence these equations are solved iteratively. In our calculation, the stopping criterion is that the maximum relative error of Q between the two iteration steps is less than 1 106. The validation of the developed codes for the BTEs and the MBTEs is carried out by comparing the simulations results with the analytic 1-D Heaslet Warming solution for steady problems. For transient problems, our simulation results of both BTE and BDEs are compared with the results of a 1-D problem in reference [6] and a 2-D problem with the emitted temperature boundary condition in references [6,17]. Very good agreements have been obtained for all cases. 4. Results and discussions
i¼1
3. Numerical procedures In this work, the finite difference method and the upwind scheme are employed to solve the BTE and the Gauss–Legendre quadrature is used to calculate the equilibrium distribution function. To solve for the BDEs, we use the ray-tracing method to calculate the ballistic intensity, and the finite difference method with the central differential scheme to solve the diffusive equations. Detailed procedures can be found in reference [17]. Similar numerical schemes are used to solve the modified BTEs. For the MBDEs, the ballistic intensity is again calculated using the ray-tracing method. To solve Eq. (27), the finite difference
In this section, simulation results of a few test cases calculated from the modified ballistic–diffusive equations are presented and compared with results obtained from the original BDEs, the BTE and the MBTEs. These problems have been studied in the previous work [5,6,17] and are the model problems for a wide range of applications. Since the focus of this work is on the accuracy and efficiency of the proposed MBDEs, non-dimensional parameters defined in Eq. (38) are used in all calculations:
T T0 Tb ¼ ; T1 T0 t x t ¼ ; x ¼ ; L s
b ¼ Q
Q C v g ðT 1 T 0 Þ
y y ¼ ; L
Kn ¼
K vgs ¼ L L
ð38Þ
55
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63 1
0.3
+
BDE + + MBDE + + BTE + + + MBTE + + + + + 0.6 + + + + + + + + + + + + t*=100 + + + + + 0.4 + + + + + + + + + + t*=10 + + + + + + 0.2 * + t =1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++ ++ ++ + ++ ++ ++ ++ ++ ++ 0 0.8 +
0
0.2
0.4
0.6
0.8
0.25
Non-dimensional Heat Flux
Non-dimensional Temperature
+ +
+
+ t*=1 + 0.1 +
++ * + t =10 + + + + 0.05 + + ++ + + + ++ 0
1
+
+
+
+
+
+
+
+
+
+
+
+
+
+
t =10
+
0.2
+
+
t*=1
+
+
+
+
+
+ t =0.1 *
+
+
+
+
+
0
+
0.2
+
+
+
+
0.4
+
+
+
+
+
+
0.6
+
+
+
+
+ + + + ++
0.8
0.2
0.4
0.6
0.8
1
Non-dimensional Coordinate
BDE MBDE BTE MBTE
+ +
+
++
+
+
+ +
+ +
+
+
+
+ +
+ + ++
+
t*=10
+ +
+
+
+
+
+
+
+
+
+
+
+
*
t =0.1
+
+ +
+
0
+
+ t =1 + *
0.1
0
1
+
0.2
+
+
+
+
0.4
+
+
0.6
+
+
+
+
0.8
+ + +
1
Non-dimensional Coordinate
0.3
+ + + +
+ +
+ +
+
*
+ +
+ +
+
t =1
+
+
+
+
+
+
+
+
+
+
+
+ t =0.01 + *
+
+
+
0
+
0.2
+
+
+
0.4
+
+
+
0.6
+
0.2
+
+
+
+
+
0.8
Non-dimensional Coordinate
+ + ++
++ + +
+ +
+
+
+
+
+
+
+
+
+
t =1
+ + t*=0.1
+
+ + + +
+
+ + +
+ +
*
t =0.01
+ +
+ 1
0
+
*
+
0.1
BDE MBDE BTE MBDE
+
+ +
+
0.05
+
+ +
+
t =0.1
+
+ ++ +
0.15
*
+ 0.2
0.25 + +
BDE MBDE BTE MBTE
0.8
Non-dimensional Heat Flux
Non-dimensional Temperature
+ + + + + + + + + + + + ++ ++ ++ ++ ++ ++
+
+
1
0
+
0.2
0.05
+
0.4
0
+
+
0.15
*
+
+
Non-dimensional Heat Flux
Non-dimensional Temperature
+
0.25
BDE MBDE BTE MBTE
+
0.6 +
0.6
+
0.3
+
0
+
*
t =100
+ +
Non-dimensional Coordinate
1
0.4
+
0.15 +
Non-dimensional Coordinate
0.8
BDE MBDE BTE MBTE
0.2 +
0
+
+
0.2
+
+
+
0.4
+
+
+
0.6
+
+
+
0.8
Non-dimensional Coordinate
+
+ + +
1
Fig. 2. Comparison of the temperature and the heat flux profiles obtained from the BDE, the MBDE, the BTE, and the MBTE at different time and phonon Knudsen number.
where L is the characteristic length, T0 and T1 are the temperatures at the boundary. 4.1. Case 1 Heat conduction inside a 1-D thin film of thickness L is simulated. The film is initially at a uniform temperature T0. At time t,
one of its surfaces is subjected to a phonon flux at a temperature of T1. This problem was studied using the original BDEs and results were presented in [5,6]. A fine spatial mesh with Nx = 800 is used for the discretization of the BTE to accurately capture the phenomenon and a relatively coarse mesh with Nx = 200 is used for the simulations with all the other methods in which the ballistic component can be calculated analytically. For the angular
56
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63 0 .2
0.0 6 *
Non-dimensional Temperature
K n= 1 , t = 1
0 .1 5
Non-dimensional Heat Flux
D iffusive P a rt *
K n= 1 , t = 1
BD E M BD E 0 .1
P ositive
0 .0 5
BD E M BD E
0.0 4 P ositive
0.0 2
0
D iffusive P a rt
-0.0 2
N e ga tive
N e ga tive
0
-0.0 4 0
0 .2
0.4
0 .6
0 .8
1
0
0.2
N on-dim e nsiona l C oordina te
0.4
0 .6
0 .8
1
N on-dim e nsiona l C oordina te
Fig. 3. Comparison of the diffusive temperature and heat flux obtained from the BDE and the MBDE; Positive indicates the temperature or the heat flux calculated from the phonon intensity along the positive x-axis, Iþ m ; Negative refers the temperature or the heat flux corresponding to Im .
0.52 0.51 0.5+
+ + 10
-2
10
-1
10
t
0.72 0.7
Temperature
0.68
1-D Kn=1
0.66
BDE MBDE BTE MBTE
0.64 0.62
+
0.6 0.58 0.56 0.54 0.52 0.5 -2 10
+
+ + + 10-1
+
t
+
0.24
0.235
0.23
0
10
+ + + + + + + + + + + + + + + +
0.74
1
10
+
0.7 0.65 0.6 0.55
+ + +
0.5 -1 10
0
t
+ + +
0.2
100
10
0.14 10-2
1
10
1
+
++
+ ++
+
10-1
+ ++ ++
0.2
0.15
+
+
10
-1
0.5+
10
10
1
-2
10
-1
10
t
0.516
+ + + + + + + + + + + + + + + + + + +
t
100
0
10
1
0.514
1-D Kn=1
0.512
BDE MBDE BTE MBTE
0.51
+
0.508 0.506 0.504 0.502
+
0.5
++
10-2
101
+ + + 10-1
+
t
+
+ +
+ + + + + + + + + + + + + +
++
100
101
0.515
+
+
+
+ + ++ + + + + + + + + + + + + + + + + + + + + + + + 1-D Kn=0.1 + + BDE +
+
0.05 2
0.505
0.518
+ + + + + + + 1-D Kn=0.1 + BDE + + MBDE + BTE + MBTE + + + + + +
0.1
10
0
0.51
++++
0.52
1-D Kn=1 BDE MBDE BTE MBTE +
0.18
10
t
0.16
+ + + 1-D Kn=0.1 + BDE + + MBDE + BTE + + MBTE + + + + + + + + + + + + + + 10
-1
0.22
Q
Temperature
0.75
+ + 1-D Kn=10 + + BDE + MBDE + BTE + MBTE + + + + + + + + + + + + + ++ +++++++ 10
0.24
0.25
0.8
-2
++
+ + +
0.9 0.85
+
0.245
++ + + + + + + 1-D Kn=10 + BDE + MBDE + + BTE + MBTE + + + + + + + + + + + ++ ++ + + ++
Q+T/2
0.53
+ + ++ +
Q+T/2
0.54
0.25+
++++
10
0
t
10
1
0.51
Q+T/2
Temperature
0.55
++
Q
0.56
+ + + + + + + + 1-D Kn=10 + + BDE + MBDE + BTE + MBTE + + + + + + + + + + ++ ++
Q
0.57
0.505
+ +
++
++
+ + +
10
+
+
+ +
+
MBDE BTE MBTE
0.5 2
10
-1
10
0
t
10
1
10
2
Fig. 4. Comparison of the surface temperature, heat flux and the sum of the heat flux and half of temperature at x = 0 as functions of time at different Knudsen numbers.
discretization, Nh = 160 is used for the BTE, the MBTE and the calculations of the ballistic component in the BDE and the MBDE. The time step is chosen to be 0.05, 0.005 and 0.0005 for Knudsen
number of 0.1, 1.0 and 10.0 respectively. Both the temperature and the heat flux profiles at different time and different Knudsen number are calculated using the modified BDEs and plotted in
57
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
Fig. 5. Schematic illustration of the problem domain, the boundary conditions and the local coordinate system of the 2-D problem.
0.6 +
0.25
+ + + + + + + + +
0.4
0.3
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+ 0.5 +
+
0.2
+ + +
0.1
0
+ +
0
+
+
0.2 +
+ 0.15
+
0.5
+ + + + + + + + + + + 1
1.5
+ + +
0.1
+ +
+ + +
0.05
+ + +
0
2
0
+ + + + + + + + + + +
0.5
1
1.5
2
Non-Dimentional Coordinate 0.2
0.8
+ + 0.6 + + + + + + 0.4 + + + + + + 0.2
0
+
+
+
+
++
BTE MBTE BDE MBDE-2*2 MBDE-4*4
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+ +
+ +
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
+ +
Non-Dimentional Coordinate
0
+
++ + + + + + + + + + ++ + + 2
4
Non-Dimentional Coordinate
+ + + 0.15 + + + + + + 0.1 + + + + + + 0.05 +
0
0
+
+
+
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+ + + + + + + + + + + + + + + 2
4
Non-Dimentional Coordinate
Fig. 6. Temperature and the y-direction heat flux along the centerline (x = Lx/2) at Kn = 1, t⁄ = 1, 10.
Fig. 2. The comparison with the reference solutions obtained from the BTE and the MBTEs indicates that the original BDEs predicts both the temperature and heat flux quite well when the Knudsen number is large and the time scale is small. However in both the intermediate and small Knudsen number regimes in which the diffusive component plays an important role, when the time scale is large, the temperature profile within the Knudsen layer deviates significantly from the reference solutions. Although the heat flux profile still matches well with the reference profiles, small discrepancies can be observed near the boundary. Overall, the prediction of the heat flux is more accurate than that of the temperature.
On the contrary, the results obtained from our modified BDEs agree with the reference solutions very well in all cases. The differences in the results from the BDEs and the MBDEs are largely caused by the boundary conditions employed in the two methods. As previously mentioned, the boundary conditions used in the BDE method are inconsistent with the P1-approximation for the diffusive component. The new boundary conditions employed in the MBDEs do not have such an inconsistence, which lead to a much improved accuracy particularly near the boundary of the problem domain. One can further understand this point by examining closely the differences in the two sets of results.
58
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.3
+
0.5 +
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
+ 0.4 +
+
+ + + +
+
0.3
+ +
+
+
0.2
+
+
+
t*=0.1
+
0.1
0
Non-Dimentional Heat Flux
Non-Dimentional Temperature
0.6
+
1
+
+ 2
+
3
+
+ 4
+ 5 /0
1
+
+ + +
+
+ +
+ +
+
+
0.1
+
+
+ t*=0.1
*
t =1.0
+
+
+ +
+
+
BTE BDE MBTE MBDE-2*2 MBDE-4*4
0.2 +
+
+
0
*
t =1.0
+
+
+ +
2
+
3
+
4
+
+ 0
5
0
+ 1
+
+
2
+ 3
+
+
4
5 /0
1
2
+
3
+ 4
+ 5
Non-Dimentional Coordinate
Non-Dimentional Coordinate
Fig. 7. Temperature and the y-direction heat flux along the centerline x = Lx/2 at Kn = 10.
0.3
0.6
t*=1
0.5
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
0.4
t*=1
+++ + + + + + + +
+++
Non-Dimentional Heat Flux
Non-Dimentional Temperature
0.7
+ + + + + +
0.3 0.2
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+ 0.2
+ + +++ + + +
++
+++
+ + +
0.1
0.1 0+ + + + + + + + + + 0 2 4
+ +
+
0+ + + + + + + + + +
+ + + + + + + + +
6
8
10
0
2
++ +
+
4
+ + + + + + + +
6
8
10
Non-Dimentional Coordinate
Non-Dimentional Coordinate 0.8 *
0.6
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
t*=10
+++ + + + + + + +
++ ++
Non-Dimentional Heat Flux
Non-Dimentional Temperature
t =10
+ + + + +
0.4
0.2
+ ++ + + + + + + + 0
1
2
3
4
+ +
+ + 5
0.2
BTE MBTE + + BDE + + + MBDE-2*2 MBDE-4*4
+
0+ + + + + + + + +
+ 6
+ + + + + + + + + 7
Non-Dimentional Coordinate
8
9
10
0
1
2
3
+ 4
+ +
+ + 5
+ + + + + + +
+ 6
++
+ + + + + + +
+ + + + + + + ++ 7
8
9
Non-Dimentional Coordinate
Fig. 8. Comparison of the temperature and the y-direction heat flux along the line y = 0 and at t⁄ = 1, 10. The Knudsen number is 1.0.
10
59
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.3
+
0.5
+
+
*
t =0.1 0.4
0.3
+
+
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
++
Non-Dimentional Heat Flux
Non-Dimentional Temperature
0.6
+
+
0.2
0.2
+ +
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
++++
+
0.1
0.1
0+ 0
+
+
+
2
+
+
+
+
4
6
+
+
+
8
+
+
0+ 0
+
10
+
+
+
+
+
2
+
+
+
+
4
+
6
+
+
8
+ 10
Non-Dimentional Coordinate
Non-Dimentional Coordinate 0.6
0.3
+
0.5
+
*
t =1 0.4
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
0.3
+
++
++++
+
+
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+
t*=0.1
+
0.2
+ t*=1 0.2
+
+
BTE MBTE BDE MBDE-2*2 MBDE-4*4
+
+
++++
+
0.1
0.1
0+ 0
+
+
+
2
+
+
+ 4
+
+
+
+
6
8
+
+
+
10
Non-Dimentional Coordinate
0+
+
+
0
+
+
2
+
+
4
+
+
+
+
+
6
+
8
+ 10
Non-Dimentional Coordinate
Fig. 9. Comparison of the temperature and the y-direction heat flux along the line y = 0 and at t⁄ = 0.11. The Knudsen number is 10.
Fig. 3 plots the temperature and heat flux calculated from the diffusive part of the intensity. Those predicted by the BDEs, particularly the results corresponding to Iþ m , exhibit weak but artificial wave fronts. As discussed in the reference [6], these wave fronts are non-physical. On the other hand, results calculated from the MBDEs do not have such non-physical wave fronts and are more accurate. Fig. 4 presents the boundary temperature and heat flux at x = 0 as functions of time. It is evident that the MBDEs predict both the transient temperature and heat flux at the boundary much more accurately than the BDEs, especially at small Knudsen numbers when the contribution from the diffusive part plays an important role. As shown in the third column of Fig. 4 where the quantity Q + T/2 is plotted, the MBDEs can capture the time dependent trend very well, while the results obtained from the BDEs remain constant. The constant value is due to the inherent assumption in the BDE approach that Q + T/2 is equal to half of the wall temperature. It is also noted that the results predicted by the BTE are not accurate at the beginning, and the reason is that the upstream scheme adopted in the BTE is unable to resolve the rapidly changed temperature at the beginning. By applying the MBTEs, such a shape change can be well captured by the ballistic part. The remaining part, that is the diffusive part, is smooth and can then be accurately simulated.
4.2. Case 2 A two-dimensional problem similar to the one in [6,17] is solved using the MBDEs. The problem domain, the boundary conditions and the local coordinate system are illustrated in Fig. 5, and the initial temperature everywhere inside the domain is T0. The dimensions of the domain are chosen to be:
Lh ¼ L;
Lx ¼ 10L;
Ly ¼ 5L
where L is the size of the heater region and served as the characteristic length based on which the Knudsen number is calculated. In this case, the same non-uniform grid system as described in [17] is used for the simulation of the BTE. Near the heated region, that is, 4.5 6 x 6 6.5, 0 6 y⁄ 6 1, a uniform mesh of Nh = 60 is used and outside the heated region, the grid is refined using a geometric series. A total mesh of Nx Ny = 228 140 is used. For the simulations with other methods, a uniform mesh is used in the entire simulation domain. For Kn = 1, Nx Ny = 280 200 is used for the MBTE and the MBDE, Nx Ny = 400 200 is used for the BDE and Nh Nu = 100 20 is used for all the other methods. For Kn = 10, Nx Ny = 200 140 is used in the simulations of the MBTE, the BDE and the MBDE; Nh Nu = 120 20 is used for the simulations of the BTE and the MBTE; Nh Nu = 120 40 is used for the calculations of the ballistic components of BDE and MBDE. The time
60
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.0008
0.0006 *
+
0.0006
t =0.05
MBTE BDE MBDE-2*2 MBDE-4*4
0.0004
++ + + + + + +
0.0002
+
+
+ +
+ + + +
+ 0+ + + + + + + + + + 0 1 2 3 4
Non-Dimentional Heat Flux
Non-Dimentional Temperature
*
t =0.05
5
+
MBTE BDE MBDE-2*2 MBDE-4*4
0.0004
++ + + + + + +
0.0002
+ + +
+ + + + + + + + +
6
7
8
9
8
9
10
0.2
+
BTE + + MBTE BDE + + MBDE-2*2 MBDE-4*4 +
+
0.15
+
+
+
+ + Non-Dimentional Heat Flux
Non-Dimentional Temperature
7
+ + + + +
t*=0.1
+
+
+
+ +
0.1
+ + + +
0.05
+ +
+
+
2
+
+
+
t*=0.1
0.15
+ BTE + MBTE BDE MBDE-2*2 + MBDE-4*4
+ 0.1
+
4
6
+
+
+
+
+ +
+ + +
0.05
+ +
+
+
+
+
+
+ +
+
+
Non-Dimentional Coordinate
8
+
+
0+ 0
10
+
+
+
2
+
+
+ +
4
+
6
8
+
+
Non-Dimentional Coordinate
+
10
0.3
0.3
+ + +
t*=1
BTE + MBTE BDE + MBDE-2*2 MBDE-4*4 +
+
+ ++ ++ + +
+ +
+
+
+
+
+
+
+ + 0.1
+ +
+ + +
2
BTE MBTE + BDE + MBDE-2*2 MBDE-4*4
+
0.2
+ + 0.1
+
+ + +
+
+
6
8
+ 10
Non-Dimentional Coordinate
0+ 0
2
+
+
+
+
+ + + +
+
+
+
+ + +
+
+ 4
+ ++ ++
+
+
+
+
+
+
+
+ + + + +
+ + +
t*=1
Non-Dimentional Heat Flux
Non-Dimentional Temperature
+ + + + + + + + +
6
0.2
0.25
0+ 0
5
Non-Dimentional Coordinate
0.3
0.2
+ + + +
+
0+ + + + + + + + + + 0 1 2 3 4
10
Non-Dimentional Coordinate
0+ 0
+
+
+
+ 4
+ + +
6
8
+ 10
Non-Dimentional Coordinate
Fig. 10. Comparison of the temperature and the y-direction heat flux at y = 0.5, Kn = 10.
step is chosen to be 0.05 and 0.005 for Knudsen number of 1.0 and 10.0 respectively. Temperature and heat flux of the entire domain are calculated at various time instants and Knudsen numbers using the four
methods described in Section 2. In the modified ballistic–diffusive approach, the division of the solid angle space is done at two levels, 2 2 and 4 4 sections, to check the convergence. Fig. 6 shows the calculated five sets of the temperature and the y-direction heat flux
61
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
4.3. Case 3 The last example considered is the problem of radiative heat transfer in an isotropically scattering media which was studied previously [18]. The problem domain and the local coordinate system are the same as that of case 2. The media is initially at a uniform temperature T0, and then the temperature of the side y = 0 is raised to and maintained at T1. The dimensions and the Knudsen number are:
Lx ¼ L;
Ly ¼ 10; Lx
K Kn ¼ ¼ 10 L
0.6
0.5 X Crosbie and Schrenker C.Cai BTE MBDE-2*2 BDE MBTE
Temperature
X X
X
0.4
X 0.3
X X
0.2
X X X
0.1
0
X
0
0.5
1
y
1.5
X
2
Fig. 11. Steady-state temperature profile along the centerline x = Lx/2. Data indicated by dots and crosses are taken from the papers [18,19].
0.53 BTE MBDE-2*2 BDE MBTE
0.525
Temperature
along the centerline of the domain (x = Lx/2) at Kn = 1 and t⁄ = 1, 10. Similar to that of case 1, the results predicted by the BDEs are less accurate near the boundary, while the MBDEs can capture both the temperature and the heat flux quite well along the entire line. As the Knudsen number increases, the ballistic component becomes dominant and consequently the results predicted by the BDEs agree better with the reference solutions obtained from the BTE and MBTE as demonstrated in Fig. 7 where the temperature and the y-direction heat flux at Kn = 10 are plotted. The temperature and the heat flux profiles along the lines of y = 0 and y = 0.5 are also examined and plotted in Figs. 8–10. The Knudsen number in Fig. 8 is 1, while in Figs. 9 and 10 it is 10. Although the overall profiles obtained from the BDEs agree well with the reference solutions obtained from the MBTEs, discrepancies can be observed in the peaks of the profiles as shown in the insets of the plots particularly in Fig. 8 when Kn = 1, in this case the MBDE is much accurate than the original BDE. For the cases of Kn = 10, those results obtained from the MBDEs are generally more accurate than those from the original BDEs except the profiles along y = 0.5 and at t⁄ = 0.1 in which case the original BDEs perform better than the MBDEs. However, as shown in Fig. 10, at t⁄ = 0.05 and t⁄ = 0.1, the results from the BDEs are worse than that of the MBDEs, indicating that t⁄ = 0.1 is simply a special case for the BDEs. Hence even in high Knudsen number cases when the ballistic transport is dominant, the MBDEs perform better than the original BDEs. Also it should be pointed out that there is a small difference between the results obtained from solving the BTE and the MBTE in the region of 4.5 < x < 5.5. Such a difference is caused by the ray effect and false scattering in solving the BTE due to the discontinue wall temperatures [14] at x = 4.5 and x = 5.5. The MBTE approach eliminates these effects and produces more accurate results.
0.52
0.515
0.51
0.505
0.5
0
0.2
0.4
x
0.6
0.8
1
Fig. 12. Comparison of the steady state temperature profiles along the line y = 0 obtained from the BTE, the MBTE, the BDE and the MBDE at Kn = 10.
where L is chosen to be the characteristic length. Again, the steadystate temperature field is calculated using the four methods. In this case, Nx Ny = 40 200 and Nh Nu= 120 40 are used for all the methods. Results along the centerline (x = Lx/2) are presented in Fig. 11. Reference solutions extracted from reference [18] are also plotted for comparison. In this case, all four sets of temperature profiles agree well with the reference solutions from [18]. However, for the steady-state temperature profile along the bottom side, plotted in Fig. 12, the original BDE approach predicts a constant profile, which is significantly different from the curved profiles obtained from the BTE, the MBTEs and the MBDEs, which agree with each other quite well even quantitatively. This case further demonstrates the accuracy of the proposed MBDEs.
5. Efficiency To evaluate the efficiency of the MBDE approach, we compared the CPU time consumed by the four methods in solving case 2, which are listed in Table 2. A relatively coarse spatial mesh, Nx Ny = 100 100, is employed in the numerical simulation. To solve the BTEs, the MBTEs and the ballistic transport equation of the BDEs and the MBDEs, the solid angle is discretized with a mesh of Nh Nu = 80 20. To solve the diffusive equations of the MBDEs, the solid angle is divided into 2 2 and 4 4 sections. The time step in the simulation is 0.005 and the Knudsen number is 10. As shown in Table 2, both the original BDE and the MBDE are more efficient than the BTE and the MBTE. The CPU time required in the MBTE is only slightly larger than that of the BTE, indicating that the ballistic transport equation can be solved much more efficiently compared to the diffusive transport equation. Hence the major factor contributing to the efficiency of the ballistic–diffusive approaches lies in the diffusive transport equation. In the BTE and the MBTE, the phonon intensity is solved explicitly, while in the BDE and the MBDE, diffusive equations in terms of temperature and flux are solved and therefore are more efficient. It is also noticed that in the case when the solid angle is discretized into
Table 2 Comparison of computational time for a simulation duration of 200 time steps. Methods
BTE
MBTE
BDE
MBDE(2 2)
MBDE(4 4)
Time (s)
338.7
342.48
88.3
77.4
82.0
62
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63
0.6 MBDE (2 X 2) : NθXNϕ=160 X 40
0.5
Temperature
and a modified ballistic–diffusive approach (MBDE) is proposed by introducing a new boundary model. A new set of ballistic–diffusive equations is derived. Several benchmark problems are employed to validate the accuracy of the new approach, and the results show that the new method is a much more accurate approximation of the BTE than the original BDE while retaining the same level of the efficiency of the original BDE. In particular, the MBDE eliminates the distortion of the temperature profiles predicted by the BDE at the boundary, and provides much more accurate steady-state solutions. In addition, similar to the MBTE, the MBDE is quite effective in eliminating the ray effects and false scattering due to the splitting of the ballistic and the diffusive intensity. Overall, the MBDE is an efficient method for the modeling of transient non-continuum heat conduction with a good accuracy and much less computational time and memory consumption. The extension of the MBDE for 3-D problems should be straightforward and is the subject of the future work.
MBDE-50 X 30 MBDE-100 X 60 MBDE-100 X 60 for Diffusive and 50 X 30 for Ballistic MBDE-200 X 120
0.4
0.3
0.2
0.1
0
0
1
2
y
3
4
5
Fig. 13. Temperature profiles along the centerline x = Lx/2 at Kn = 10.
Conflict of interest None declared.
Table 3 Mesh size and CPU time for the four meshes. Mesh
Parameters
CPU time (s)
Ballistic Part
1 2 3 4
Diffusive Part
Nx
Ny
Nx
Ny
50 100 50 200
30 60 30 120
50 100 100 200
30 60 60 120
Appendix A. Appendix The detailed expressions for Ai, Bi, Ci, Fi in Eq. (26) are listed as follows:
28.3 111.0 39.1 439.6
R
0 B Ai ¼ B @
R Xi
R Xi
four sections, that is, Nh Nu = 2 2, the time used by the MBDE is even less than that of the BDE. Such a reduction in the CPU time is due to the number of iterations of the MBDE, which is smaller than that of the BDE at each time step. In addition to the CPU time, the memory consumptions of the BDE and the MBDE are also much less than that of the BTE and the MBTE due to the fact that the detailed information of the solid angle, which is implicitly included in the macroscopic quantities, does not need to be stored in the ballistic–diffusive approaches. Hence they are more suitable for solving 3-D problems than the transport equation based approaches. It is realized that the computational time can be further reduced without losing accuracy by using a finer mesh for the diffusive part and a coarse mesh for the ballistic part. As the ballistic transport equation is solved in an analytical form, a coarse mesh is sufficient to obtain accurate results. And for the diffusive equations, a fine mesh is required in order to reduce the numerical errors caused by the finite difference scheme. Table 3 shows the mesh size and computational time for the four meshes employed in the simulation. The corresponding temperature profiles along the centerline x = Lx/2 are shown in Fig. 13. In all calculations, 160 40 sections of the solid angle are employed in solving the ballistic transport equation and 2 2 sections are used for solving the diffusive equations. As evident from Fig. 13 and Table 3, among the four sets of meshes, Mesh 3, which adopts a fine mesh for the diffusive part and a coarse mesh for the ballistic part, provides both efficient and accurate solution. 6. Summary In this paper, the inconsistence between the boundary conditions and the P1 approximation in the original BDE is identified
Xi
R
B Bi ¼ B @
R
R
coshdX
sinhcosudX
0
R
dX
Xi
R Xi
Xi Xi
coshdX cos2 hdX
sinhcoshcosudX R
coshdX
R
Xi
cos2 hdX
R R
Xi
sinhcosudX
1
C sinhcoshcosudX C A R 2 2 Xi sin hcos udX
Xi
R R
Xi
sinhcoshcosudX
1
C sinhcos2 hcosudX C A R R R 2 2 2 Xi sinhcoshcosudX Xi sinhcos hcosudX Xi sin hcoshcos udX cos2 hdX
Xi
Xi
cos3 hdX
Xi
0
1 R R R 2 2 sinhcosudX Xi sinhcoshcosudX Xi sin hcos udX B R Xi C R R B C 2 Ci ¼ B Xi sinhcoshcosudX Xi sinhcos2 hcosudX Xi sin hcoshcos2 udX C @ R A R R 2 2 3 2 2 3 Xi sin hcos udX Xi sin hcoshcos udX Xi sin hcos udX
R
0 B Fi ¼ B @
R Xi
R Xi
Xi
dX
coshdX
sinhcosudX
1 C C A
References [1] D.G. Cahill, W.K. Ford, K.E. Goodson, G.D. Mahan, A. Majumdar, H.J. Maris, R. Merlin, S.R. Phillpot, Nanoscale thermal transport, J. Appl. Phys. 93 (2) (2003) 793. [2] S.B. Riffat, X. Ma, Thermoelectrics: a review of present and potential applications, Appl. Therm. Eng. 23 (8) (2003) 913–935. [3] A.A. Joshi, A. Majumdar, Transient ballistic and diffusive phonon heat transport in thin films, J. Appl. Phys. 74 (1) (1993) 31. [4] J.M. Ziman, Electrons and Phonons: the Theory of Transport Phenomena in Solids, Oxford University Press, 2001. [5] G. Chen, Ballistic–diffusive heat-conduction equations, Phys. Rev. Lett. 86 (11) (2001) 2297–2300. [6] G. Chen, Ballistic–diffusive equations for transient heat conduction from nano to macroscales, J. Heat Transfer 124 (2) (2002) 320. [7] P. Cheng, Dynamics of a radiating gas with application to flow over a wavy wall, AIAA J. 4 (2) (1966) 238–245.
Y. Zhang, W. Ye / International Journal of Heat and Mass Transfer 83 (2015) 51–63 [8] P. Cheng, Study of the flow of a radiating gas by a differential approximation, in: DTIC Document, 1965. [9] M.F. Modest, Radiative Heat Transfer, McGraw-Hill, New York, 1993. [10] A. Mittal, S. Mazumder, Generalized ballistic–diffusive formulation and hybrid SN-PN solution of the Boltzmann transport equation for phonons for nonequilibrium heat conduction, J. Heat Transfer 133 (9) (2011) 092402. [11] A. Majumdar, Microscale heat conduction in dielectric thin films, J. Heat Transfer 115 (1) (1993) 7–16. [12] M.A. Ramankutty, A.L. Crosbie, Modified discrete ordinates solution of radiative transfer in two-dimensional rectangular enclosures, J. Quant. Spectrosc. Radiat. Transfer 57 (1) (1997) 107–140. [13] J.Y. Murthy, S.R. Mathur, An improved computational procedure for submicron heat conduction, J. Heat Transfer 125 (5) (2003) 904. [14] P.J. Coelho, The role of ray effects and false scattering on the accuracy of the standard and modified discrete ordinates methods, J. Quant. Spectrosc. Radiat. Transfer 73 (2) (2002) 231–238.
63
[15] M. Sakami, A. El Kasmi, A. Charette, Analysis of radiative heat transfer in complex two-dimensional enclosures with obstacles using the modified discrete ordinates method, J. Heat Transfer 123 (5) (2001) 892–900. [16] M. Sakami, A. Charette, Application of a modified discrete ordinates method to two-dimensional enclosures of irregular geometry, J. Quant. Spectrosc. Radiat. Transfer 64 (3) (2000) 275–298. [17] R. Yang, G. Chen, M. Laroche, Y. Taur, Simulation of nanoscale multidimensional transient heat conduction problems using ballistic– diffusive equations and phonon Boltzmann equation, J. Heat Transfer 127 (3) (2005) 298. [18] J.C. Chai, H.S. Lee, S.V. Patankar, Finite volume method for radiation heat transfer, J. Thermophys. Heat Transfer 8 (3) (1994) 419–425. [19] A.L. Crosbie, R.G. Schrenker, Radiative transfer in a two-dimensional rectangular medium exposed to diffuse radiation, J. Quant. Spectrosc. Radiat. Transfer 31 (4) (1984) 339–372.