- Email: [email protected]
A boundary-domain integral equation method for solving convective heat transfer problems
International Journal of Heat and Mass Transfer 63 (2013) 183–190
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A boundary-domain integral equation method for solving convective heat transfer problems Xiao-Wei Gao ⇑, Hai-Feng Peng, Jian Liu State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China
a r t i c l e
i n f o
Article history: Received 12 December 2012 Received in revised form 20 March 2013 Accepted 29 March 2013 Available online 27 April 2013 Keywords: Convective heat transfer Navier–Stokes equations Energy equation Boundary element method Steady incompressible viscous flow
a b s t r a c t In this paper, a new boundary-domain integral equation for solving energy equation in convective heat transfer problems is derived. The derived integral equation accounts for fluid velocity which is computed using the velocity boundary-domain integral equations presented in [2,3] (X.W. Gao, 2004, 2005). Numerical implementation of these integral equations is discussed for steady incompressible fluid flows, in which the velocity and pressure integral equations are uncoupled from the temperature equation, and thus can be solved separately. By substituting the velocity, velocity gradient, and pressure integral equations into the energy integral equation, temperatures can be computed without iterative processes. Two numerical examples for 2D problems are given to validate the proposed method. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The governing equations for convective heat transfer problems are the continuity, momentum, and energy equations. In general, they are coupled nonlinear partial differential equations which are not amenable to analytical solutions except for a few simple problems [1]. Therefore, in the past half century, many works have been done to develop efficient numerical techniques to solve these equations. Currently, the most popular methods for solving these equations are the finite difference method, finite element method, and finite volume method. Their common feature is that they are based on domain variable representations and local interpolation schemes, resulting in systems of equations that are highly sparse matrices. Recently, the boundary element method (BEM) has been used to analyze convective heat transfer problems. Compared to other numerical methods, the distinct advantages of BEM are that the velocity and temperature gradient formulations can be explicitly derived from the velocity and temperature integral equations so that the pressure and heat flux have the same computational accuracy as velocity and temperature themselves (e.g. [2–4]), and the boundary conditions at infinity are automatically satisfied in BEM formulations. During the past three decades, BEM formulations for convective heat transfer problems have been extensively studied [5–15]. Part of the works focus on the governing convective diffusion equations [5–9]. The first integral equation analysis for steady convective ⇑ Corresponding author. Tel.: +86 411 84706332. E-mail address: [email protected] (X.-W. Gao). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.03.071
diffusion equations was carried out by Ikeuchi and Onishi [5] using the fundamental solutions of the diffusion-convection and Laplacian operator to develop direct and iterative boundary element methods respectively. Shi and Banerjee [6] developed a boundary element formulation for linearized convective heat transfer problems based on unsteady fundamental solutions. In a different way, free-space time-dependent convective diffusion fundamental solutions were used by Grigoriev and Dargush [7] to obtain higherorder BEM integral equations. Then the equations are extended to solve two-dimensional problems [8], from which numerical examples are given to investigate the performance of the BEM formulations [9]. In these mentioned methods, the viscous dissipation is not taken into consideration in the convective heat transfer problems since Navier–Stokes equations are not employed in the fluid flow computation. A BEM formulation was directly derived from the governing equations of two-dimensional incompressible viscous fluids by Dargush and Banerjee [10] for solving convective heat transfer problems. The used energy equation includes the terms of viscous dissipations as a portion of pseudo sources. The resulting integral equations are written exclusively in terms of velocities and temperatures, and the computation of velocity and temperature gradients is removed through integration by parts. The same technique was also employed by Banerjee and Honkala [11] to obtain a set of integral equations incorporating the buoyancy effect for the analysis of the natural convection heat transfer. Different BEM formulations were also established to solve natural convection problems by other authors, e.g. Onishi et al. [12], and Tosaka and Fukushima [13]. In addition, the dual reciprocity boundary element method
184
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
was used by Rahaim and Kassab [14] to solve incompressible laminar viscous fluid flows and heat transfer, in which all nonlinear terms were lumped into the forcing term. As to the transient convective heat transfer, Vajravelu et al. [15] developed a one-dimensional boundary integral equation to analyze the free heat transfer in a viscous, electrically conducting and heat-generating fluid past a vertical porous plate in the presence of free stream oscillations with the use of time-dependent fundamental solutions. Recently, a set of general boundary-domain integral equations for solving full Navier–Stokes equations using primitive variables was proposed by Gao [2–4]. The BEM formulations are expressed in terms of velocity, traction, and pressure, and are valid for steady, unsteady, compressible, and incompressible flows. Although these integral equations are derived under isothermal conditions, they can also be applied to analyze non-isothermal viscous fluids. However, as done in this paper, an integral equation from the energy equation is required to close the integral equation system. This paper presents a new boundary-domain integral equation for computing temperatures, which is derived from the full energy equation without heat source rates. In this equation, velocities, velocity gradients and pressures are included, which can be computed using the BEM formulations presented in Refs. [2–4]. Together with these integral equations, a complete set of general boundary-domain integral equations for solving convective heat transfer problems is established. These equations are valid for steady, unsteady, compressible, and incompressible viscous fluid flows. For compressible fluids, the number of unknown variables is one more than those provided by the integral equations. Therefore, an equation of state is required to close the equation system. However, for incompressible viscous fluid flows, the boundary-domain integral equation set given in this paper is closed, and the temperature integral equation is uncoupled from the velocity and pressure integral equations. After velocities, velocity gradients, and pressures are computed through evaluating the correspondent integral equations, the temperatures can be solved by using the derived temperature integral equation. Details of the numerical implementation involved are provided for steady incompressible viscous flows. Finally two numerical examples are given to validate the formulations derived in this paper.
wi ¼ rij uj qEui
ð4bÞ
where e is the internal energy per unit mass, and for the calorically perfect gas (constant specific heats)
e ¼ cv T
ð5Þ
where cv is the specific heat at constant volume. This equation is sometimes labeled the caloric equation of state. For incompressible viscous flows, cv = cp (here cp is the specific heat at constant pressure). For Newtonian fluids, the constitutive relationship between the stresses and velocities based on Stokes’ hypothesis can be expressed as:
rij ¼ pdij þ l
@ui @uj 2 @uk l þ dij 3 @xk @xj @xi
ð6Þ
in which, p is the pressure, l the dynamics coefficient of viscosity (constant here), and dij the Kronecker delta symbol. On the fluid surface with outward normal nj, the relationship between the stress and the traction ti (force per unit area) can be written as:
ti ¼ rij nj
ð7Þ
Based on these equations, a complete system of boundary integral equations for viscous fluid mechanics can be well established.
3. Review of velocity boundary-domain integral equations for viscous flows The velocity boundary-domain integrals for 2D and 3D fluid flows formulated with respect to primary variables (velocities, tractions and pressure) can be written as [2,3]:
ui ðxÞ ¼
2. Governing differential equations in viscous fluid mechanics
Z
Z uij ðx; yÞt j ðyÞdCðyÞ t ij ðx; yÞuj ðyÞdCðyÞ C C Z uij ðx; yÞnk ðyÞqðyÞuj ðyÞuk ðyÞdCðyÞ ZC þ uij;k ðx; yÞqðyÞuj ðyÞuk ðyÞdXðyÞ ZX Z þ uij;j ðx; yÞpðyÞdXðyÞ þ uij ðx; yÞqðyÞbj ðyÞdXðyÞ X ZX @ quj uij ðx; yÞ dXðyÞ @t X
ð8Þ
In order to analyze the fluid flow and heat transfer in fluid mechanics, one needs to consider the conservation laws of the mass, momentum and energy. After some operations on these conservation laws, the governing differential equations in viscous fluid mechanics can be expressed as:
where X denotes the domain of the problem and C the boundary of the domain X. x is referred to as the source point and y as the field point, and ð Þ;k ¼ @ð Þ=yk . The fundamental solutions appearing in integral Eq. (8) can be expressed as:
@ q @ qui þ ¼0 @t @xi
ð1Þ
uij ðx; yÞ
@ qui @ qui uj @ rij þ ¼ þ qbi @t @xj @xj
ð2Þ
tij ðx; yÞ ¼
@ qE @ @T @wi þ ¼ k þ qbi ui @t @xi @xi @xi
ð3Þ
In the above equations, q is the fluid density, t the time, ui the ith velocity component, xi the ith component of spatial coordinates at point x, bi the body force per unit mass (e.g. the gravity force), rij the stress tensor, E the total energy per unit mass, T the temperature, k the thermal conductivity, and the repeated subscripts stand for summation. Other quantities included in Eq. (3) can be determined using
E ¼ e þ ui ui =2
ð4aÞ
( ¼
1 f7dij 16apl
lnð1r Þ þ r;i r;j g for 2D
1 f7dij 16aplr
þ r;i r;j g
for 3D
ð9Þ
1 3ðni r;j nj r;i Þ þ ðbr;i r;j þ 3dij Þnk r;k 8apr a
ð10Þ
uij;k ðx; yÞ ¼
1 7dij r ;k dik r ;j djk r;i þ br;i r ;j r ;k 16aplr a
ð11Þ
uij;j ðx; yÞ ¼
3r ;i 8aplr a
ð12Þ
where a = b 1 with b = 2 for 2D and b = 3 for 3D problems, r is the distance from the source point x to the field point y, i.e. r = ky xk, and
r;i ¼
@r y xi ¼ i @yi r
ð13Þ
185
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Differentiating Eq. (8) with respect to the source points x and using the singular integral separation technique, the velocity gradient boundary-domain integral equations can be derived as follows [4]: Z Z @ui ðxÞ ¼ uij;l ðx; yÞtj ðyÞdCðyÞ þ t ij;l ðx; yÞuj ðyÞdCðyÞ @xl C C Z þ uij;l ðx; yÞnk ðyÞqðyÞuj ðyÞuk ðyÞdCðyÞ ZC Z @ quj uij;l ðx; yÞqðyÞbj ðyÞdXðyÞ þ uij;l ðx; yÞ dXðyÞ @t X X Z Z uij;kl ðx; yÞqðyÞuj ðyÞuk ðyÞdXðyÞ uij;jl ðx; yÞpðyÞdXðyÞ X
þ
X
1 3 fðb þ 7cÞdij dkl 2ðdik djl þ djk dil ÞgqðxÞuj ðxÞuk ðxÞ þ dil pðxÞ 8bcl 4bl ð14Þ
in which c = b + 2 and the used fundamental solutions are
1 3ðni djl þ nl dij nj dil Þ þ b½3ðnj r;i ni r;j Þr;l 8par b þnl r;i r;j þ br;m nm ½dil r;j þ djl r;i 3dij r;l cr;i r;j r;l
t ij;l ðx; yÞ ¼
uij;kl ðx;yÞ ¼
ð15Þ
1 dil djk þdik djl 7dij dkl 16aplrb bðdik r;j r;l þdjk r;i r;l 7dij r;k r;l þdjl r;i r;k þdkl r;i r;j þdil r;j r;k Þ þbcr;i r;j r;k r;l ð16Þ 3 fdil br;i r;l g 8aplr b
uij;jl ðx; yÞ ¼
ð17Þ
Incorporating the indices i and l in Eq. (14) and substituting the resulting velocity divergence boundary-domain integral equation into continuity Eq. (1), the integral equation for evaluating the pressure p can written as [4]:
3 pðxÞ ¼ 4l
Z
Z uij;i ðx; yÞtj ðyÞdCðyÞ t ij;i ðx; yÞuj ðyÞdCðyÞ C C Z uij;i ðx; yÞnk ðyÞqðyÞuj ðyÞuk ðyÞdCðyÞ ZC þ uij;ki ðx; yÞqðyÞuj ðyÞuk ðyÞdXðyÞ ZX Z @ quj þ uij;i ðx; yÞqðyÞbj ðyÞdXðyÞ uij;i ðx; yÞ dXðyÞ @t X X 3 1 @q @q qðxÞui ðxÞui ðxÞ u ðxÞ þ 4bl qðxÞ i @xi @t
other integrals included in Eqs.(8), (14), and (18) are strongly or weakly singular and are interpreted in the Cauchy principal value sense. The singularity separation technique [16] is applied to evaluate the domain integrals appearing in Eq. (14) for kernels uij;kl and uij;jl , and in Eq. (18) for kernel uij;ki , since they are strongly singular even for internal collocation points. For the weakly singular integrals, the element/cell sub-division techniques [16] are used to remove the integral singularities. Eqs. (8), (14), and (18) are general boundary-domain integral equations for the analysis of viscous fluid flows. For incompressible flows, the pressures appearing in the discretized algebraic equations from Eq. (8) can be eliminated by employing Eq. (18), such that the final system of equations only contains the velocities/tractions as unknowns. Using the results obtained by Eqs. (8) and (18), the velocity gradients can be solved through the evaluation of integral Eq. (14). They are only applied to analyze isothermal viscous flows. In order to solve convective heat transfer problems, the energy equation needs to be considered. For this purpose, a new boundary-domain integral equation for evaluating the temperature will be derived in the next section. 4. Temperature boundary-domain integral equation In order to make the symbols appearing in this section consistent with those used in Eqs. (8), (14), and (18), let yi represent xi in Eq. (3), and then the energy equation can be expressed as follows:
@ qE @ @T @wi þ ¼ k þ qbi ui @t @yi @yi @yi
Multiplying both sides of the above equation by a weight function T⁄ and integrating over the domain O bounded with the boundary C, it is obtained that
Z
T
X
uij;i ðx; yÞ ¼
ð19Þ
3r ;j 8aplra
ð20Þ
uij;ki ðx; yÞ ¼
3 8aplr
djk br;j r;k b
Z
T
X
Z X
ð21Þ
When the field points y approaches the source points x, the fundamental solutions in Eqs. (9)–(12), (15)–(17), and (19)–(21) are singular. From the kernel functions t ij;l included in Eq. (15) and tij;i in Eq. (19), it can be seen that the second boundary integrals in Eqs. (14) and (18) are hyper-singular and only exist in the sense of the Hadamard finite-part. Therefore, integral Eqs (14) and (18) are valid only for internal collocation points. For boundary points, the traction-recovery technique as done in literatures [3,16] is used to evaluate the velocity gradient and pressure. Except for the above hyper-singular integrals, all the
Z @ @T @wi dX þ T T k dX @y @y @yi X X i i Z þ T qbi ui dX Z
ð23Þ
Integrating the former two domain integrals on the right-hand side of Eq. (23) by parts and using Gauss’s divergence theorem, it follows that
where
3 fdjl br;j r;l gnl 4apr b
@ qE dX ¼ @t
X
ð18Þ
t ij;i ðx; yÞ ¼
ð22Þ
T
Z Z Z @ @T @T @T dX ¼ k T k dC Tk dC þ Tk @yi @yi @n @n C C X @ @T dX @yi @yi @wi dX ¼ @yi
Z
T wi ni dC
Z
C
X
wi
@T dX @yi
ð24aÞ
ð24bÞ
If the weight function T⁄ is the Green’s function T⁄(x, y) which satisfies the following equation:
@ @T ðx; yÞ þ dðy xÞ ¼ 0 @yi @yi
ð25Þ
where d(y x) is the Dirac delta function of y at the point x, then in terms of the integration property of the Dirac delta function [17], the last domain integral of Eq. (24a) can be written as:
Z X
TðyÞk
Z @ @T dXðyÞ ¼ TðyÞkðyÞdðy xÞdXðyÞ ¼ kTðxÞ @yi @yi X ð26Þ
Substituting this equation into Eq. (24a) and the result together with Eq. (24b) into Eq. (23), it follows that
186
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Z
T ;n ðx; yÞkTðyÞdCðyÞ Z Z þ T ðx; yÞni ðyÞwi ðyÞdCðyÞ T ;i ðx; yÞwi ðyÞdXðyÞ X ZC Z @ qE þ T ðx; yÞqðyÞbi ðyÞui ðyÞdXðyÞ T ðx; yÞ dXðyÞ @t X X ð27Þ
kTðxÞ ¼
T ðx; yÞqðyÞdCðyÞ
Z
C
C
where T ;n ðx; yÞ ¼ @T ðx; yÞ=@n, T ;i ðx; yÞ ¼ @T ðx; yÞ=@yi and q(y) is the heat flux, i.e.,
qðyÞ ¼ k
@TðyÞ @n
ð28Þ
The fundamental solution T⁄ and its derivatives in Eq. (25) can be written as:
(
T ðx; yÞ ¼
1 2p
1 ln rðx;yÞ for 2D
1 4prðx;yÞ
for 3D
1 T ;i ðx; yÞ ¼ r;i 2par a T ;n ðx; yÞ
1 ¼ r ;i ni 2par a
ð29Þ
@ui @uj 2 @uk uj l þ ui ðp þ qui ui =2Þui þ qcv ui T 3 xk @xj @xi
5. Numerical implementation for steady incompressible flows Eqs. (8), (14), (18), and (27) are general integral equations for solving convective heat transfer problems. To be considered in this section is the numerical implementation of these equations in steady incompressible flows. In this case, the density q is constant and the time-related domain integrals in Eqs. (8), (14), (18), and (27)) disappear. The details of the numerical implementation of integral Eqs. (8), (14), and (18) can be found in Refs. [2,3]. To numerically implement Eq. (27), the boundary C and the domain X of the problem should be discretized into boundary elements and internal cells, respectively. For each boundary element or internal cell, the temperature, heat flux and wi can be expressed in terms of their nodal values through shape functions as follows:
T¼
X g g N T
ð33aÞ
g
ð30Þ q¼
X g N qg
ð33bÞ
g
ð31Þ wi ¼
The thermal conductivity k is assumed to be constant in the derivation process of integral Eq. (27). If the thermal conductivity to be considered is the spatial function of coordinates y, a domain integral will appear in the resulting integral equation [18]. As a result, the complexity of the solution of convective heat transfer problems will be raised. The varying thermal conductivity problems will be treated with in the future research. Eq. (27) is a general energy boundary-domain integral equation valid for steady, unsteady, compressible, and incompressible viscous fluids. It is bounded only for internal collocation points, since the second integral associated with the kernel function T ;n ðx; yÞ in Eq. (27) is strongly singular when the source point x approaches the field point y. For boundary collocation points, the rigid body motion strategy can be used to treat the strong singularity as in conventional BEM [16]. Except for the second integral, all other integrals in Eq. (27) are weakly singular and therefore can be evaluated accurately by using the element/cell sub-division technique [16]. Eq. (27) together with Eqs. (8) and (18) described in the above section constitute a boundary-domain integral equation system for analyzing convective heat transfer problems using the primitive variables. For compressible viscous flows, since the density is not constant, the number of unknowns being velocity/traction, pressure, temperature, and density is one more than the number of integral Eqs. (8), (18), and (27). Therefore, an equation of state, such as p = qRT for perfect gases (here R is the specific gas constant), is required to close the final system of equations. For incompressible viscous flows, the system of Eqs. (18), (18), and (27) is closed. Since l is assumed constant as done in Section 2, integral Eqs. (8) and (18) are uncoupled from the temperature integral Eq. (27) and thus the velocity and pressure can be solved independently on the temperature, after which the temperature can be solved from the integral Eq. (27). It is noted that the temperature integral equation itself is not independent of the velocity and pressure, embodied in Eq. (27) through the term wi. Substituting Eq. (5) into Eq. (4a), and the result together with Eq. (6) into Eq. (4b) yields:
wi ¼ l
From Eq. (32), it can be seen that the quantity wi appearing in integral Eq. (27) includes the velocity, velocity gradient, pressure, and temperature.
ð32Þ
X g g N wi
ð33cÞ
g
where Ng are the shape functions for boundary elements and internal cells [16], and Tg represents the value of temperature T at node g. Substituting Eqs. (33) into the discretized form of Eq. (27) and collocating x for all boundary nodes (the rigid body motion is used for determination of diagonal terms) yield the following algebraic matrix equation:
½HfTg ¼ ½Gfqg þ fbg þ ½Eb fwg
ð34Þ
where {T} and {q} are vectors consisting of temperatures and heat fluxes at all boundary nodes, respectively, and {b} is a constant vector from body forces. The vector {w} is given through the evaluation of Eq. (32) or Eqs. (4) and (5) at all boundary and internal nodes. [Eb] is the coefficient matrix determined by both the third boundary integral and the fourth domain integral on the right-hand side of Eq. (27). At each boundary node, either temperature or heat flux is specified as a boundary condition. So after applying boundary conditions to Eq. (34) and rearranging the equation, it follows that
½Ab fXg ¼ fY b g þ ½Eb fwg
ð35Þ
where {X} is a vector consisting of unknown temperatures and unknown heat fluxes, and {Yb} is a known vector. The square matrix [Ab] is invertible and thereby
fXg ¼ fY b g þ ½Eb fwg
ð36Þ
where
fY b g ¼ ½Ab 1 fY b g ½Eb ¼ ½Ab 1 ½Eb
ð37Þ
Similarly, for internal nodes, Eq. (27) gives
fT I g ¼ ½AI fXg þ fY I g þ ½EI fwg
ð38Þ
where {TI} is the vector consisting of all temperatures at internal nodes. Substituting Eq. (36) into Eq. (38) yields:
fT I g ¼ fY I g þ ½EI fwg where
ð39Þ
187
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Table 1 The number of boundary elements, internal cells, boundary nodes, internal nodes and total nodes in four BEM models.
Cell240 Cell320 Cell480 Cell640
Boundary element
Internal cell
Boundary node
Internal node
Total node
76 96 92 112
240 320 480 640
76 96 92 112
203 273 435 585
279 369 527 697
Fig. 1. Geometry and boundary conditions for Poiseuille flow.
" I
b
fY g ¼ fY I g þ ½AI fY g
ð40Þ
½EI ¼ ½EI þ ½AI ½Eb
Upon close inspection of Eq. (32), it can be found that the vector {w} appearing in Eqs. (36) and (39) consists of the constant values and unknown temperatures over all nodes after the velocities, velocity gradients, and pressures are computed using integral equations (8), (14), and (18). Therefore, substituting the computational results into Eq. (32), the vector fwg can be written as:
fT b g fwg ¼ fcg þ sB; ½Dt fT I g
ð41Þ
where {c} is a constant vector consisting of both the term qcvuiT with known boundary temperatures and all other terms in Eq. (32), and {Tb} is the vector consisting of unknown temperatures at boundary nodes. The coefficient matrices [B] and [D] corresponding to the vectors {Tb} and {TI}, respectively, are determined by the term – qcvui. It is noted that after using boundary conditions, the columns of matrix [B] corresponding to the known boundary temperature nodes should be zero. Substituting Eq. (41) into Eqs. (36) and (39), it follows that
e bg ½Ab fXg þ ½ e E b fT I g ¼ f Y
ð42Þ
e Ig E I fT I g ¼ f Y ½AI fXg þ ½ e
ð43Þ
where
8 b b > > < ½A ¼ ½I ½E ½B ½e E b ¼ ½Eb ½D > > : eb f Y g ¼ fY b g þ ½Eb fcg
ð44Þ
8 I I > < ½A ¼ ½E ½B ½e E I ¼ ½I ½EI ½D > : I e g ¼ fY I g þ ½EI fcg fY
ð45Þ
Eb ½Ab ½ e I ½A ½ e EI
fXg fT I g
(
¼
e bg fY e Ig fY
) ð46Þ
The above equation is the final system of equations with boundary unknowns {X} and internal temperatures {TI} as unknowns. Since Eq. (46) is a linear equation about {X} and {TI}, it can be solved without iterative computation. 6. Numerical examples To verify the correctness of the proposed method in this paper, two numerical examples for two-dimensional steady incompressible flows are presented in this section. The first one is the wellknown Poiseuille flow which has analytical solutions to verify, and the second one is a 2D pipe flow with curvature and the comparison is made to FLUENT CFD results. 6.1. Poiseuille flow It is assumed that the fluid flows between two parallel plates without body forces (Fig. 1). The ‘no-slip’ condition is applied to the upper and lower sides. The left and right ends are insulated, while the upper and lower sides are imposed with the temperatures of 40 and 10, respectively. Other boundary conditions applied are t x ¼ p ¼ 50; uy ¼ 0 on left end and tx ¼ p ¼ 0; uy ¼ 0 on right end. Therefore the viscous fluid moves with a uniform x-component of velocity ux, and the velocity and temperature vary only in the y-direction. Referring to literature [19], the analytic solutions for the horizontal velocity and temperature can be expressed as:
ux ¼
in above equations, [I] is the identity matrix. Combining Eqs. (42) and (43) yields
#
T¼
p0 yðH yÞ 2l
p’2 ð2y HÞ4 15 p’2 H4 þ ð2y HÞ þ 25 þ H 192lk 192lk
ð47aÞ
ð47bÞ
where H = 1 is the distance between the upper and lower sides, p0 = dp/dx = 10 the gradient of pressure, k = 1 is the thermal conductivity and l = 1.
Fig. 2. BEM models with the different number of internal cells: (a) 240; (b) 320; (c) 480; (d) 640.
188
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Table 2 Computed temperatures along y-direction at x = 2.5. y
0.125
0.25
0.375
0.5
0.625
0.75
0.875
Cell240 Cell320 Cell480 Cell640 Analytical
14.118778 14.119476 14.108114 14.108603 14.106038
18.010649 18.011593 17.992470 17.993120 17.988281
21.797131 21.798141 21.774462 21.775145 21.768799
25.551184 25.552202 25.527060 25.527737 25.520833
29.297131 29.298143 29.274578 29.275251 29.268799
33.010651 33.011597 32.992710 32.993340 32.988281
36.618771 36.619469 36.608498 36.608940 36.606038
Table 3 Relative errors (%) between the numerical results and analytical solutions along ydirection at x = 2.5. y
0.125
0.25
0.375
0.5
0.625
0.75
0.875
Cell240 Cell320 Cell480 Cell640
0.0903 0.0953 0.0147 0.0182
0.1243 0.1296 0.0233 0.0269
0.1301 0.1348 0.0260 0.0292
0.1189 0.1229 0.0244 0.0271
0.0968 0.1003 0.0197 0.0220
0.0678 0.0707 0.0134 0.0153
0.0347 0.0367 0.0067 0.0079
Fig. 5. BEM model for 2D pipe with curved part.
35
Temperature
Fig. 3. Computed temperatures along y-direction at x = 2.5.
Current
30
Fluent
25
20
15
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
Length of centreline Fig. 6. Computational temperature along the center line of the pipe from the inlet to the outlet.
Fig. 4. Geometry and boundary conditions for 2D pipe flow.
To use the formulations presented in the paper, the computational domain and its boundary are required to be discretized into the internal cells and boundary elements, respectively. As shown in Fig. 2, four BEM models with the different number of internal cells (240, 320, 480 and 640) are considered to examine the influences of the number of internal cells on the computational results. The number of boundary elements, internal cells, boundary nodes, internal nodes, and total nodes for each BEM model is given in Table 1. Table 2 lists the computed temperatures and analytical solutions along y-direction of x = 2.5, the relative errors of which are shown in Table 3. Fig. 3 is a plot of the temperature profile. The
189
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
0.3
0.4
0.5
Y-coordinate
2426 22
Fig. 7. Computational temperature on the outlet of 2D pipe.
30
20
Temperature
Current 25
12
Fluent
34 32 30
14
28 26 24 22 16
20
18
16
14 12
0.2
20
0.1
18
10 0.0
2 8 30 3 2
15
12
28 30
20
2426 22 18
34
25
36
Temperature
Fluent
20
32
Current 30
14 1 6
35
16 14
12
20
Fig. 9. Contour plot of temperature computed by the current method.
15
The exact solutions for the heat flux q on the upper and lower sides of the fluid are 25.8333 and 34.1667, respectively, and the computed results using Cell640 BEM model are 25.8657 and 34.1359. The relative errors are 0.13% and 0.09%, respectively.
10 0.0
0.1
0.2
0.3
0.4
0.5
Length of line L Fig. 8. Computational temperature along the line L.
6.2. 2D angled pipe flow
computational results for velocity are not given in this paper as detailed description can be found in Ref. [2]. It is noted that the problem is independent of the density q and the specific heat at constant volume cv (see the analytical solution (47)). From Tables 2 and 3 and Fig. 3, it can be seen that the computational results with the use of 240, 320, 480, and 640 internal cells are in good agreement with the analytical solutions. The comparison of the results denoted by Cell240 and Cell480 to the analytical solutions shows that the results of 480 internal cells is closer to the analytical solutions than that of 240 internal cells. However, from the results denoted by Cell240 and Cell320, we can find that the relative errors between the results of 240 cells and analytical solutions are smaller than that of Cell320 results. These mean that the fining of the internal cells along y-direction (Fig. 2(a) and (c)) can make the numerical results more accurate, while the mesh refinement along x-direction improves computational accuracy very little (see Fig. 2(a) and (b)). The phenomenon is mainly due to the fluid velocity and temperature varying along y-direction, not along x-direction.
The second example deals with a convective heat transfer over 2D pipe with a curved part (Fig. 4). The curvature of the curved part is determined by a radius of R = 0.5. The fluid with the property of l = 1, q = 1 and cv = 0.4930966 is subjected to a pressure of p = 50 on the inlet of the pipe and is pressure free on the outlet. The pipe is approximated using 140 linear boundary elements with 140 boundary nodes and 600 internal cells with 531 internal nodes (Fig. 5). The boundary conditions applied are tx = 0, ty = 50, T = 30 on the inlet of the pipe, tx = ty = 0, q = 0 on the outlet, and ux = uy = 0, q = 50 on the upper wall and ux = uy = 0, T = 10 on the lower wall. Fig. 6 shows the computed temperatures along the center line of the pipe from the inlet to the outlet. Figs. 7 and 8 are the distributions of temperature along the y-axis on the outlet of the pipe and line L (see Fig. 4), respectively, and specific results are listed in Table 4, along with the relative errors. For the purpose of verification, this problem is also computed using the commercial software FLUENT, where a fine mesh with 2000 quadrilateral
Table 4 Temperatures on the outlet of pipe and the line L (see Fig. 4). Distance
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Outlet
Line L
Current
FLUENT
Error (%)
Current
FLUENT
Error (%)
12.4586 14.9283 17.4042 19.8831 22.3658 24.8533 27.3456 29.8417 32.3413 34.7998
12.5502 15.0785 1735955 20.1103 22.6288 25.1536 27.6841 30.1266 32.7439 34.9444
0.73 0.80 1.09 1.13 1.16 1.19 1.22 1.24 1.20 0.41
11.5245 13.1173 14.7776 16.5162 18.3448 20.2757 22.3215 24.4935 26.7990 29.2373
11.5275 13.1123 14.7654 16.4977 18.3206 20.2463 22.2871 24.4546 26.7568 29.1906
0.026 0.038 0.083 0.112 0.132 0.145 0.154 0.159 0.158 0.160
30 26 24 22 18 32
30 28
14
16
26 24 22 20
12
20
28
2426 22
18
14 16
20
28
34
12
14 16 2 8
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
32
190
18
14
16
12
Fig. 10. Contour plot of temperature computed by FLUENT.
elements and 2121 nodes is used and pressure boundary conditions are applied on the inlet and outlet. The computational results using FLUENT are also shown in Figs. 6–8 and Table 4. Figs. 9 and 10 are the contour plots of the temperature obtained by BEM and FLUENT, respectively. From Figs. 6–8 and Table 4, it can be seen that the BEM results are close to the results by FLUENT with a good accuracy. However, there is a small discrepancy between these two results near the outlet of the pipe. Table 4 shows the maximum error on the outlet of 2D pipe being 1.24%. 7. Conclusions A novel boundary-domain integral equation for computing the temperature in viscous fluids is presented based on the conservation form of the energy equation. The derived formulation is general and applicable to steady, unsteady, compressible, and incompressible flows. Together with the formulations for the velocity, pressure, and velocity gradient proposed in Refs. [2–4], a boundary-domain integral equation system for solving convective heat transfer problems is established. For incompressible flows, the system of equations is closed and the velocity and pressure integral equations are uncoupled from the temperature integral equation. However, the temperature integral equation involves the velocity, pressure, and velocity gradient, which need
to be computed in advance. Two numerical examples for convective heat transfer in steady incompressible fluids have been provided to demonstrate the correctness of the derived formulations. Acknowledgements The authors gratefully acknowledge the National Natural Science Foundation of China under Grant NSFC Nos. 10872050 and 11172055 for financial supports to this work.
References [1] A. Bourchtein, Exact solution of the generalized Navier–Stokes equations for benchmarking, Int. J. Numer. Methods Fluids 39 (2002) 1053–1071. [2] X.W. Gao, A boundary-domain integral equation method in viscous fluid flow, Int. J. Numer. Methods Fluids 45 (2004) 463–484. [3] X.W. Gao, A promising boundary element formulation for three-dimensional viscous flow, Int. J. Numer. Methods Fluids 47 (2005) 19–43. [4] X.W. Gao, Explicit formulations for evaluation of velocity gradients using boundary-domain integral equations in 2D and 3D viscous flows, Int. J. Numer. Methods Fluids 54 (2007) 1351–1368. [5] M. Ikeuchi, K. Onishi, Boundary element solutions to steady convective diffusion equations, Appl. Math. Model. 7 (1983) 115–118. [6] Y. Shi, P.K. Banerjee, Boundary element methods for convective heat transfer, Comput. Methods Appl. Mech. Eng. 105 (1993) 261–284. [7] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part I: General formulations and 1D implementation, Comput. Methods Appl. Mech. Eng. 192 (2003) 4281–4298. [8] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part II. 2D implementation, Comput. Methods Appl. Mech. Eng. 192 (2003) 4299–4312. [9] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part III: Numerical examples, Comput. Methods Appl. Mech. Eng. 192 (2003) 4313–4335. [10] G.F. Dargush, P.K. Banerjee, A boundary element method for steady incompressible thermoviscous flow, Int. J. Numer. Methods Eng. 31 (1991) 1605–1626. [11] P.K. Banerjee, K.A. Honkala, Development of BEM for thermoviscous flow, Comput. Methods Appl. Mech. Eng. 151 (1998) 43–62. [12] K. Onishi, T. Kuroki, M. Tanaka, An application of a boundary element method to natural convection, Appl. Math. Model. 8 (1984) 383–390. [13] N. Tosaka, N. Fukushima, Integral Equation Analysis of Laminar Natural Convection Problems, BEM VIII, Computation Mechanics Publication, SpringerVerlag, Southampton, Berlin, 1986. [14] C.P. Rahaim, A.J. Kassab, Pressure correction DRBEM solution for heat transfer and fluid flow in incompressible viscous fluids, Eng. Anal. Bound. Elem. 18 (1996) 265–272. [15] K. Vajravelu, A. Kassab, A. Hadjinicolaou, BEM solution to transient free convection heat transfer in a viscous, electrically conducting, and heat generating fluid, Numer. Heat Transfer Part A: Appl. 30 (1996) 555–574. [16] X.W. Gao, T.G. Davies, Boundary Element Programming in Mechanics, Cambridge University Press, Cambridge, 2002. [17] M.D. Greenberg, Application of Green’s Functions in Science and Engineering, Prentice-Hall, Englewood Cliffs, New Jersey, 1971. [18] X.W. Gao, A meshless BEM for isotropic heat conduction problems with heat generation and spatially varying conductivity, Int. J. Numer. Methods Eng. 66 (2006) 1411–1431. [19] R.A. Granger, Fluid Mechanics, Dover, New York, 1995.
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A boundary-domain integral equation method for solving convective heat transfer problems Xiao-Wei Gao ⇑, Hai-Feng Peng, Jian Liu State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China
a r t i c l e
i n f o
Article history: Received 12 December 2012 Received in revised form 20 March 2013 Accepted 29 March 2013 Available online 27 April 2013 Keywords: Convective heat transfer Navier–Stokes equations Energy equation Boundary element method Steady incompressible viscous flow
a b s t r a c t In this paper, a new boundary-domain integral equation for solving energy equation in convective heat transfer problems is derived. The derived integral equation accounts for fluid velocity which is computed using the velocity boundary-domain integral equations presented in [2,3] (X.W. Gao, 2004, 2005). Numerical implementation of these integral equations is discussed for steady incompressible fluid flows, in which the velocity and pressure integral equations are uncoupled from the temperature equation, and thus can be solved separately. By substituting the velocity, velocity gradient, and pressure integral equations into the energy integral equation, temperatures can be computed without iterative processes. Two numerical examples for 2D problems are given to validate the proposed method. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The governing equations for convective heat transfer problems are the continuity, momentum, and energy equations. In general, they are coupled nonlinear partial differential equations which are not amenable to analytical solutions except for a few simple problems [1]. Therefore, in the past half century, many works have been done to develop efficient numerical techniques to solve these equations. Currently, the most popular methods for solving these equations are the finite difference method, finite element method, and finite volume method. Their common feature is that they are based on domain variable representations and local interpolation schemes, resulting in systems of equations that are highly sparse matrices. Recently, the boundary element method (BEM) has been used to analyze convective heat transfer problems. Compared to other numerical methods, the distinct advantages of BEM are that the velocity and temperature gradient formulations can be explicitly derived from the velocity and temperature integral equations so that the pressure and heat flux have the same computational accuracy as velocity and temperature themselves (e.g. [2–4]), and the boundary conditions at infinity are automatically satisfied in BEM formulations. During the past three decades, BEM formulations for convective heat transfer problems have been extensively studied [5–15]. Part of the works focus on the governing convective diffusion equations [5–9]. The first integral equation analysis for steady convective ⇑ Corresponding author. Tel.: +86 411 84706332. E-mail address: [email protected] (X.-W. Gao). 0017-9310/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.03.071
diffusion equations was carried out by Ikeuchi and Onishi [5] using the fundamental solutions of the diffusion-convection and Laplacian operator to develop direct and iterative boundary element methods respectively. Shi and Banerjee [6] developed a boundary element formulation for linearized convective heat transfer problems based on unsteady fundamental solutions. In a different way, free-space time-dependent convective diffusion fundamental solutions were used by Grigoriev and Dargush [7] to obtain higherorder BEM integral equations. Then the equations are extended to solve two-dimensional problems [8], from which numerical examples are given to investigate the performance of the BEM formulations [9]. In these mentioned methods, the viscous dissipation is not taken into consideration in the convective heat transfer problems since Navier–Stokes equations are not employed in the fluid flow computation. A BEM formulation was directly derived from the governing equations of two-dimensional incompressible viscous fluids by Dargush and Banerjee [10] for solving convective heat transfer problems. The used energy equation includes the terms of viscous dissipations as a portion of pseudo sources. The resulting integral equations are written exclusively in terms of velocities and temperatures, and the computation of velocity and temperature gradients is removed through integration by parts. The same technique was also employed by Banerjee and Honkala [11] to obtain a set of integral equations incorporating the buoyancy effect for the analysis of the natural convection heat transfer. Different BEM formulations were also established to solve natural convection problems by other authors, e.g. Onishi et al. [12], and Tosaka and Fukushima [13]. In addition, the dual reciprocity boundary element method
184
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
was used by Rahaim and Kassab [14] to solve incompressible laminar viscous fluid flows and heat transfer, in which all nonlinear terms were lumped into the forcing term. As to the transient convective heat transfer, Vajravelu et al. [15] developed a one-dimensional boundary integral equation to analyze the free heat transfer in a viscous, electrically conducting and heat-generating fluid past a vertical porous plate in the presence of free stream oscillations with the use of time-dependent fundamental solutions. Recently, a set of general boundary-domain integral equations for solving full Navier–Stokes equations using primitive variables was proposed by Gao [2–4]. The BEM formulations are expressed in terms of velocity, traction, and pressure, and are valid for steady, unsteady, compressible, and incompressible flows. Although these integral equations are derived under isothermal conditions, they can also be applied to analyze non-isothermal viscous fluids. However, as done in this paper, an integral equation from the energy equation is required to close the integral equation system. This paper presents a new boundary-domain integral equation for computing temperatures, which is derived from the full energy equation without heat source rates. In this equation, velocities, velocity gradients and pressures are included, which can be computed using the BEM formulations presented in Refs. [2–4]. Together with these integral equations, a complete set of general boundary-domain integral equations for solving convective heat transfer problems is established. These equations are valid for steady, unsteady, compressible, and incompressible viscous fluid flows. For compressible fluids, the number of unknown variables is one more than those provided by the integral equations. Therefore, an equation of state is required to close the equation system. However, for incompressible viscous fluid flows, the boundary-domain integral equation set given in this paper is closed, and the temperature integral equation is uncoupled from the velocity and pressure integral equations. After velocities, velocity gradients, and pressures are computed through evaluating the correspondent integral equations, the temperatures can be solved by using the derived temperature integral equation. Details of the numerical implementation involved are provided for steady incompressible viscous flows. Finally two numerical examples are given to validate the formulations derived in this paper.
wi ¼ rij uj qEui
ð4bÞ
where e is the internal energy per unit mass, and for the calorically perfect gas (constant specific heats)
e ¼ cv T
ð5Þ
where cv is the specific heat at constant volume. This equation is sometimes labeled the caloric equation of state. For incompressible viscous flows, cv = cp (here cp is the specific heat at constant pressure). For Newtonian fluids, the constitutive relationship between the stresses and velocities based on Stokes’ hypothesis can be expressed as:
rij ¼ pdij þ l
@ui @uj 2 @uk l þ dij 3 @xk @xj @xi
ð6Þ
in which, p is the pressure, l the dynamics coefficient of viscosity (constant here), and dij the Kronecker delta symbol. On the fluid surface with outward normal nj, the relationship between the stress and the traction ti (force per unit area) can be written as:
ti ¼ rij nj
ð7Þ
Based on these equations, a complete system of boundary integral equations for viscous fluid mechanics can be well established.
3. Review of velocity boundary-domain integral equations for viscous flows The velocity boundary-domain integrals for 2D and 3D fluid flows formulated with respect to primary variables (velocities, tractions and pressure) can be written as [2,3]:
ui ðxÞ ¼
2. Governing differential equations in viscous fluid mechanics
Z
Z uij ðx; yÞt j ðyÞdCðyÞ t ij ðx; yÞuj ðyÞdCðyÞ C C Z uij ðx; yÞnk ðyÞqðyÞuj ðyÞuk ðyÞdCðyÞ ZC þ uij;k ðx; yÞqðyÞuj ðyÞuk ðyÞdXðyÞ ZX Z þ uij;j ðx; yÞpðyÞdXðyÞ þ uij ðx; yÞqðyÞbj ðyÞdXðyÞ X ZX @ quj uij ðx; yÞ dXðyÞ @t X
ð8Þ
In order to analyze the fluid flow and heat transfer in fluid mechanics, one needs to consider the conservation laws of the mass, momentum and energy. After some operations on these conservation laws, the governing differential equations in viscous fluid mechanics can be expressed as:
where X denotes the domain of the problem and C the boundary of the domain X. x is referred to as the source point and y as the field point, and ð Þ;k ¼ @ð Þ=yk . The fundamental solutions appearing in integral Eq. (8) can be expressed as:
@ q @ qui þ ¼0 @t @xi
ð1Þ
uij ðx; yÞ
@ qui @ qui uj @ rij þ ¼ þ qbi @t @xj @xj
ð2Þ
tij ðx; yÞ ¼
@ qE @ @T @wi þ ¼ k þ qbi ui @t @xi @xi @xi
ð3Þ
In the above equations, q is the fluid density, t the time, ui the ith velocity component, xi the ith component of spatial coordinates at point x, bi the body force per unit mass (e.g. the gravity force), rij the stress tensor, E the total energy per unit mass, T the temperature, k the thermal conductivity, and the repeated subscripts stand for summation. Other quantities included in Eq. (3) can be determined using
E ¼ e þ ui ui =2
ð4aÞ
( ¼
1 f7dij 16apl
lnð1r Þ þ r;i r;j g for 2D
1 f7dij 16aplr
þ r;i r;j g
for 3D
ð9Þ
1 3ðni r;j nj r;i Þ þ ðbr;i r;j þ 3dij Þnk r;k 8apr a
ð10Þ
uij;k ðx; yÞ ¼
1 7dij r ;k dik r ;j djk r;i þ br;i r ;j r ;k 16aplr a
ð11Þ
uij;j ðx; yÞ ¼
3r ;i 8aplr a
ð12Þ
where a = b 1 with b = 2 for 2D and b = 3 for 3D problems, r is the distance from the source point x to the field point y, i.e. r = ky xk, and
r;i ¼
@r y xi ¼ i @yi r
ð13Þ
185
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Differentiating Eq. (8) with respect to the source points x and using the singular integral separation technique, the velocity gradient boundary-domain integral equations can be derived as follows [4]: Z Z @ui ðxÞ ¼ uij;l ðx; yÞtj ðyÞdCðyÞ þ t ij;l ðx; yÞuj ðyÞdCðyÞ @xl C C Z þ uij;l ðx; yÞnk ðyÞqðyÞuj ðyÞuk ðyÞdCðyÞ ZC Z @ quj uij;l ðx; yÞqðyÞbj ðyÞdXðyÞ þ uij;l ðx; yÞ dXðyÞ @t X X Z Z uij;kl ðx; yÞqðyÞuj ðyÞuk ðyÞdXðyÞ uij;jl ðx; yÞpðyÞdXðyÞ X
þ
X
1 3 fðb þ 7cÞdij dkl 2ðdik djl þ djk dil ÞgqðxÞuj ðxÞuk ðxÞ þ dil pðxÞ 8bcl 4bl ð14Þ
in which c = b + 2 and the used fundamental solutions are
1 3ðni djl þ nl dij nj dil Þ þ b½3ðnj r;i ni r;j Þr;l 8par b þnl r;i r;j þ br;m nm ½dil r;j þ djl r;i 3dij r;l cr;i r;j r;l
t ij;l ðx; yÞ ¼
uij;kl ðx;yÞ ¼
ð15Þ
1 dil djk þdik djl 7dij dkl 16aplrb bðdik r;j r;l þdjk r;i r;l 7dij r;k r;l þdjl r;i r;k þdkl r;i r;j þdil r;j r;k Þ þbcr;i r;j r;k r;l ð16Þ 3 fdil br;i r;l g 8aplr b
uij;jl ðx; yÞ ¼
ð17Þ
Incorporating the indices i and l in Eq. (14) and substituting the resulting velocity divergence boundary-domain integral equation into continuity Eq. (1), the integral equation for evaluating the pressure p can written as [4]:
3 pðxÞ ¼ 4l
Z
Z uij;i ðx; yÞtj ðyÞdCðyÞ t ij;i ðx; yÞuj ðyÞdCðyÞ C C Z uij;i ðx; yÞnk ðyÞqðyÞuj ðyÞuk ðyÞdCðyÞ ZC þ uij;ki ðx; yÞqðyÞuj ðyÞuk ðyÞdXðyÞ ZX Z @ quj þ uij;i ðx; yÞqðyÞbj ðyÞdXðyÞ uij;i ðx; yÞ dXðyÞ @t X X 3 1 @q @q qðxÞui ðxÞui ðxÞ u ðxÞ þ 4bl qðxÞ i @xi @t
other integrals included in Eqs.(8), (14), and (18) are strongly or weakly singular and are interpreted in the Cauchy principal value sense. The singularity separation technique [16] is applied to evaluate the domain integrals appearing in Eq. (14) for kernels uij;kl and uij;jl , and in Eq. (18) for kernel uij;ki , since they are strongly singular even for internal collocation points. For the weakly singular integrals, the element/cell sub-division techniques [16] are used to remove the integral singularities. Eqs. (8), (14), and (18) are general boundary-domain integral equations for the analysis of viscous fluid flows. For incompressible flows, the pressures appearing in the discretized algebraic equations from Eq. (8) can be eliminated by employing Eq. (18), such that the final system of equations only contains the velocities/tractions as unknowns. Using the results obtained by Eqs. (8) and (18), the velocity gradients can be solved through the evaluation of integral Eq. (14). They are only applied to analyze isothermal viscous flows. In order to solve convective heat transfer problems, the energy equation needs to be considered. For this purpose, a new boundary-domain integral equation for evaluating the temperature will be derived in the next section. 4. Temperature boundary-domain integral equation In order to make the symbols appearing in this section consistent with those used in Eqs. (8), (14), and (18), let yi represent xi in Eq. (3), and then the energy equation can be expressed as follows:
@ qE @ @T @wi þ ¼ k þ qbi ui @t @yi @yi @yi
Multiplying both sides of the above equation by a weight function T⁄ and integrating over the domain O bounded with the boundary C, it is obtained that
Z
T
X
uij;i ðx; yÞ ¼
ð19Þ
3r ;j 8aplra
ð20Þ
uij;ki ðx; yÞ ¼
3 8aplr
djk br;j r;k b
Z
T
X
Z X
ð21Þ
When the field points y approaches the source points x, the fundamental solutions in Eqs. (9)–(12), (15)–(17), and (19)–(21) are singular. From the kernel functions t ij;l included in Eq. (15) and tij;i in Eq. (19), it can be seen that the second boundary integrals in Eqs. (14) and (18) are hyper-singular and only exist in the sense of the Hadamard finite-part. Therefore, integral Eqs (14) and (18) are valid only for internal collocation points. For boundary points, the traction-recovery technique as done in literatures [3,16] is used to evaluate the velocity gradient and pressure. Except for the above hyper-singular integrals, all the
Z @ @T @wi dX þ T T k dX @y @y @yi X X i i Z þ T qbi ui dX Z
ð23Þ
Integrating the former two domain integrals on the right-hand side of Eq. (23) by parts and using Gauss’s divergence theorem, it follows that
where
3 fdjl br;j r;l gnl 4apr b
@ qE dX ¼ @t
X
ð18Þ
t ij;i ðx; yÞ ¼
ð22Þ
T
Z Z Z @ @T @T @T dX ¼ k T k dC Tk dC þ Tk @yi @yi @n @n C C X @ @T dX @yi @yi @wi dX ¼ @yi
Z
T wi ni dC
Z
C
X
wi
@T dX @yi
ð24aÞ
ð24bÞ
If the weight function T⁄ is the Green’s function T⁄(x, y) which satisfies the following equation:
@ @T ðx; yÞ þ dðy xÞ ¼ 0 @yi @yi
ð25Þ
where d(y x) is the Dirac delta function of y at the point x, then in terms of the integration property of the Dirac delta function [17], the last domain integral of Eq. (24a) can be written as:
Z X
TðyÞk
Z @ @T dXðyÞ ¼ TðyÞkðyÞdðy xÞdXðyÞ ¼ kTðxÞ @yi @yi X ð26Þ
Substituting this equation into Eq. (24a) and the result together with Eq. (24b) into Eq. (23), it follows that
186
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Z
T ;n ðx; yÞkTðyÞdCðyÞ Z Z þ T ðx; yÞni ðyÞwi ðyÞdCðyÞ T ;i ðx; yÞwi ðyÞdXðyÞ X ZC Z @ qE þ T ðx; yÞqðyÞbi ðyÞui ðyÞdXðyÞ T ðx; yÞ dXðyÞ @t X X ð27Þ
kTðxÞ ¼
T ðx; yÞqðyÞdCðyÞ
Z
C
C
where T ;n ðx; yÞ ¼ @T ðx; yÞ=@n, T ;i ðx; yÞ ¼ @T ðx; yÞ=@yi and q(y) is the heat flux, i.e.,
qðyÞ ¼ k
@TðyÞ @n
ð28Þ
The fundamental solution T⁄ and its derivatives in Eq. (25) can be written as:
(
T ðx; yÞ ¼
1 2p
1 ln rðx;yÞ for 2D
1 4prðx;yÞ
for 3D
1 T ;i ðx; yÞ ¼ r;i 2par a T ;n ðx; yÞ
1 ¼ r ;i ni 2par a
ð29Þ
@ui @uj 2 @uk uj l þ ui ðp þ qui ui =2Þui þ qcv ui T 3 xk @xj @xi
5. Numerical implementation for steady incompressible flows Eqs. (8), (14), (18), and (27) are general integral equations for solving convective heat transfer problems. To be considered in this section is the numerical implementation of these equations in steady incompressible flows. In this case, the density q is constant and the time-related domain integrals in Eqs. (8), (14), (18), and (27)) disappear. The details of the numerical implementation of integral Eqs. (8), (14), and (18) can be found in Refs. [2,3]. To numerically implement Eq. (27), the boundary C and the domain X of the problem should be discretized into boundary elements and internal cells, respectively. For each boundary element or internal cell, the temperature, heat flux and wi can be expressed in terms of their nodal values through shape functions as follows:
T¼
X g g N T
ð33aÞ
g
ð30Þ q¼
X g N qg
ð33bÞ
g
ð31Þ wi ¼
The thermal conductivity k is assumed to be constant in the derivation process of integral Eq. (27). If the thermal conductivity to be considered is the spatial function of coordinates y, a domain integral will appear in the resulting integral equation [18]. As a result, the complexity of the solution of convective heat transfer problems will be raised. The varying thermal conductivity problems will be treated with in the future research. Eq. (27) is a general energy boundary-domain integral equation valid for steady, unsteady, compressible, and incompressible viscous fluids. It is bounded only for internal collocation points, since the second integral associated with the kernel function T ;n ðx; yÞ in Eq. (27) is strongly singular when the source point x approaches the field point y. For boundary collocation points, the rigid body motion strategy can be used to treat the strong singularity as in conventional BEM [16]. Except for the second integral, all other integrals in Eq. (27) are weakly singular and therefore can be evaluated accurately by using the element/cell sub-division technique [16]. Eq. (27) together with Eqs. (8) and (18) described in the above section constitute a boundary-domain integral equation system for analyzing convective heat transfer problems using the primitive variables. For compressible viscous flows, since the density is not constant, the number of unknowns being velocity/traction, pressure, temperature, and density is one more than the number of integral Eqs. (8), (18), and (27). Therefore, an equation of state, such as p = qRT for perfect gases (here R is the specific gas constant), is required to close the final system of equations. For incompressible viscous flows, the system of Eqs. (18), (18), and (27) is closed. Since l is assumed constant as done in Section 2, integral Eqs. (8) and (18) are uncoupled from the temperature integral Eq. (27) and thus the velocity and pressure can be solved independently on the temperature, after which the temperature can be solved from the integral Eq. (27). It is noted that the temperature integral equation itself is not independent of the velocity and pressure, embodied in Eq. (27) through the term wi. Substituting Eq. (5) into Eq. (4a), and the result together with Eq. (6) into Eq. (4b) yields:
wi ¼ l
From Eq. (32), it can be seen that the quantity wi appearing in integral Eq. (27) includes the velocity, velocity gradient, pressure, and temperature.
ð32Þ
X g g N wi
ð33cÞ
g
where Ng are the shape functions for boundary elements and internal cells [16], and Tg represents the value of temperature T at node g. Substituting Eqs. (33) into the discretized form of Eq. (27) and collocating x for all boundary nodes (the rigid body motion is used for determination of diagonal terms) yield the following algebraic matrix equation:
½HfTg ¼ ½Gfqg þ fbg þ ½Eb fwg
ð34Þ
where {T} and {q} are vectors consisting of temperatures and heat fluxes at all boundary nodes, respectively, and {b} is a constant vector from body forces. The vector {w} is given through the evaluation of Eq. (32) or Eqs. (4) and (5) at all boundary and internal nodes. [Eb] is the coefficient matrix determined by both the third boundary integral and the fourth domain integral on the right-hand side of Eq. (27). At each boundary node, either temperature or heat flux is specified as a boundary condition. So after applying boundary conditions to Eq. (34) and rearranging the equation, it follows that
½Ab fXg ¼ fY b g þ ½Eb fwg
ð35Þ
where {X} is a vector consisting of unknown temperatures and unknown heat fluxes, and {Yb} is a known vector. The square matrix [Ab] is invertible and thereby
fXg ¼ fY b g þ ½Eb fwg
ð36Þ
where
fY b g ¼ ½Ab 1 fY b g ½Eb ¼ ½Ab 1 ½Eb
ð37Þ
Similarly, for internal nodes, Eq. (27) gives
fT I g ¼ ½AI fXg þ fY I g þ ½EI fwg
ð38Þ
where {TI} is the vector consisting of all temperatures at internal nodes. Substituting Eq. (36) into Eq. (38) yields:
fT I g ¼ fY I g þ ½EI fwg where
ð39Þ
187
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Table 1 The number of boundary elements, internal cells, boundary nodes, internal nodes and total nodes in four BEM models.
Cell240 Cell320 Cell480 Cell640
Boundary element
Internal cell
Boundary node
Internal node
Total node
76 96 92 112
240 320 480 640
76 96 92 112
203 273 435 585
279 369 527 697
Fig. 1. Geometry and boundary conditions for Poiseuille flow.
" I
b
fY g ¼ fY I g þ ½AI fY g
ð40Þ
½EI ¼ ½EI þ ½AI ½Eb
Upon close inspection of Eq. (32), it can be found that the vector {w} appearing in Eqs. (36) and (39) consists of the constant values and unknown temperatures over all nodes after the velocities, velocity gradients, and pressures are computed using integral equations (8), (14), and (18). Therefore, substituting the computational results into Eq. (32), the vector fwg can be written as:
fT b g fwg ¼ fcg þ sB; ½Dt fT I g
ð41Þ
where {c} is a constant vector consisting of both the term qcvuiT with known boundary temperatures and all other terms in Eq. (32), and {Tb} is the vector consisting of unknown temperatures at boundary nodes. The coefficient matrices [B] and [D] corresponding to the vectors {Tb} and {TI}, respectively, are determined by the term – qcvui. It is noted that after using boundary conditions, the columns of matrix [B] corresponding to the known boundary temperature nodes should be zero. Substituting Eq. (41) into Eqs. (36) and (39), it follows that
e bg ½Ab fXg þ ½ e E b fT I g ¼ f Y
ð42Þ
e Ig E I fT I g ¼ f Y ½AI fXg þ ½ e
ð43Þ
where
8 b b > > < ½A ¼ ½I ½E ½B ½e E b ¼ ½Eb ½D > > : eb f Y g ¼ fY b g þ ½Eb fcg
ð44Þ
8 I I > < ½A ¼ ½E ½B ½e E I ¼ ½I ½EI ½D > : I e g ¼ fY I g þ ½EI fcg fY
ð45Þ
Eb ½Ab ½ e I ½A ½ e EI
fXg fT I g
(
¼
e bg fY e Ig fY
) ð46Þ
The above equation is the final system of equations with boundary unknowns {X} and internal temperatures {TI} as unknowns. Since Eq. (46) is a linear equation about {X} and {TI}, it can be solved without iterative computation. 6. Numerical examples To verify the correctness of the proposed method in this paper, two numerical examples for two-dimensional steady incompressible flows are presented in this section. The first one is the wellknown Poiseuille flow which has analytical solutions to verify, and the second one is a 2D pipe flow with curvature and the comparison is made to FLUENT CFD results. 6.1. Poiseuille flow It is assumed that the fluid flows between two parallel plates without body forces (Fig. 1). The ‘no-slip’ condition is applied to the upper and lower sides. The left and right ends are insulated, while the upper and lower sides are imposed with the temperatures of 40 and 10, respectively. Other boundary conditions applied are t x ¼ p ¼ 50; uy ¼ 0 on left end and tx ¼ p ¼ 0; uy ¼ 0 on right end. Therefore the viscous fluid moves with a uniform x-component of velocity ux, and the velocity and temperature vary only in the y-direction. Referring to literature [19], the analytic solutions for the horizontal velocity and temperature can be expressed as:
ux ¼
in above equations, [I] is the identity matrix. Combining Eqs. (42) and (43) yields
#
T¼
p0 yðH yÞ 2l
p’2 ð2y HÞ4 15 p’2 H4 þ ð2y HÞ þ 25 þ H 192lk 192lk
ð47aÞ
ð47bÞ
where H = 1 is the distance between the upper and lower sides, p0 = dp/dx = 10 the gradient of pressure, k = 1 is the thermal conductivity and l = 1.
Fig. 2. BEM models with the different number of internal cells: (a) 240; (b) 320; (c) 480; (d) 640.
188
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
Table 2 Computed temperatures along y-direction at x = 2.5. y
0.125
0.25
0.375
0.5
0.625
0.75
0.875
Cell240 Cell320 Cell480 Cell640 Analytical
14.118778 14.119476 14.108114 14.108603 14.106038
18.010649 18.011593 17.992470 17.993120 17.988281
21.797131 21.798141 21.774462 21.775145 21.768799
25.551184 25.552202 25.527060 25.527737 25.520833
29.297131 29.298143 29.274578 29.275251 29.268799
33.010651 33.011597 32.992710 32.993340 32.988281
36.618771 36.619469 36.608498 36.608940 36.606038
Table 3 Relative errors (%) between the numerical results and analytical solutions along ydirection at x = 2.5. y
0.125
0.25
0.375
0.5
0.625
0.75
0.875
Cell240 Cell320 Cell480 Cell640
0.0903 0.0953 0.0147 0.0182
0.1243 0.1296 0.0233 0.0269
0.1301 0.1348 0.0260 0.0292
0.1189 0.1229 0.0244 0.0271
0.0968 0.1003 0.0197 0.0220
0.0678 0.0707 0.0134 0.0153
0.0347 0.0367 0.0067 0.0079
Fig. 5. BEM model for 2D pipe with curved part.
35
Temperature
Fig. 3. Computed temperatures along y-direction at x = 2.5.
Current
30
Fluent
25
20
15
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
Length of centreline Fig. 6. Computational temperature along the center line of the pipe from the inlet to the outlet.
Fig. 4. Geometry and boundary conditions for 2D pipe flow.
To use the formulations presented in the paper, the computational domain and its boundary are required to be discretized into the internal cells and boundary elements, respectively. As shown in Fig. 2, four BEM models with the different number of internal cells (240, 320, 480 and 640) are considered to examine the influences of the number of internal cells on the computational results. The number of boundary elements, internal cells, boundary nodes, internal nodes, and total nodes for each BEM model is given in Table 1. Table 2 lists the computed temperatures and analytical solutions along y-direction of x = 2.5, the relative errors of which are shown in Table 3. Fig. 3 is a plot of the temperature profile. The
189
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
0.3
0.4
0.5
Y-coordinate
2426 22
Fig. 7. Computational temperature on the outlet of 2D pipe.
30
20
Temperature
Current 25
12
Fluent
34 32 30
14
28 26 24 22 16
20
18
16
14 12
0.2
20
0.1
18
10 0.0
2 8 30 3 2
15
12
28 30
20
2426 22 18
34
25
36
Temperature
Fluent
20
32
Current 30
14 1 6
35
16 14
12
20
Fig. 9. Contour plot of temperature computed by the current method.
15
The exact solutions for the heat flux q on the upper and lower sides of the fluid are 25.8333 and 34.1667, respectively, and the computed results using Cell640 BEM model are 25.8657 and 34.1359. The relative errors are 0.13% and 0.09%, respectively.
10 0.0
0.1
0.2
0.3
0.4
0.5
Length of line L Fig. 8. Computational temperature along the line L.
6.2. 2D angled pipe flow
computational results for velocity are not given in this paper as detailed description can be found in Ref. [2]. It is noted that the problem is independent of the density q and the specific heat at constant volume cv (see the analytical solution (47)). From Tables 2 and 3 and Fig. 3, it can be seen that the computational results with the use of 240, 320, 480, and 640 internal cells are in good agreement with the analytical solutions. The comparison of the results denoted by Cell240 and Cell480 to the analytical solutions shows that the results of 480 internal cells is closer to the analytical solutions than that of 240 internal cells. However, from the results denoted by Cell240 and Cell320, we can find that the relative errors between the results of 240 cells and analytical solutions are smaller than that of Cell320 results. These mean that the fining of the internal cells along y-direction (Fig. 2(a) and (c)) can make the numerical results more accurate, while the mesh refinement along x-direction improves computational accuracy very little (see Fig. 2(a) and (b)). The phenomenon is mainly due to the fluid velocity and temperature varying along y-direction, not along x-direction.
The second example deals with a convective heat transfer over 2D pipe with a curved part (Fig. 4). The curvature of the curved part is determined by a radius of R = 0.5. The fluid with the property of l = 1, q = 1 and cv = 0.4930966 is subjected to a pressure of p = 50 on the inlet of the pipe and is pressure free on the outlet. The pipe is approximated using 140 linear boundary elements with 140 boundary nodes and 600 internal cells with 531 internal nodes (Fig. 5). The boundary conditions applied are tx = 0, ty = 50, T = 30 on the inlet of the pipe, tx = ty = 0, q = 0 on the outlet, and ux = uy = 0, q = 50 on the upper wall and ux = uy = 0, T = 10 on the lower wall. Fig. 6 shows the computed temperatures along the center line of the pipe from the inlet to the outlet. Figs. 7 and 8 are the distributions of temperature along the y-axis on the outlet of the pipe and line L (see Fig. 4), respectively, and specific results are listed in Table 4, along with the relative errors. For the purpose of verification, this problem is also computed using the commercial software FLUENT, where a fine mesh with 2000 quadrilateral
Table 4 Temperatures on the outlet of pipe and the line L (see Fig. 4). Distance
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Outlet
Line L
Current
FLUENT
Error (%)
Current
FLUENT
Error (%)
12.4586 14.9283 17.4042 19.8831 22.3658 24.8533 27.3456 29.8417 32.3413 34.7998
12.5502 15.0785 1735955 20.1103 22.6288 25.1536 27.6841 30.1266 32.7439 34.9444
0.73 0.80 1.09 1.13 1.16 1.19 1.22 1.24 1.20 0.41
11.5245 13.1173 14.7776 16.5162 18.3448 20.2757 22.3215 24.4935 26.7990 29.2373
11.5275 13.1123 14.7654 16.4977 18.3206 20.2463 22.2871 24.4546 26.7568 29.1906
0.026 0.038 0.083 0.112 0.132 0.145 0.154 0.159 0.158 0.160
30 26 24 22 18 32
30 28
14
16
26 24 22 20
12
20
28
2426 22
18
14 16
20
28
34
12
14 16 2 8
X.-W. Gao et al. / International Journal of Heat and Mass Transfer 63 (2013) 183–190
32
190
18
14
16
12
Fig. 10. Contour plot of temperature computed by FLUENT.
elements and 2121 nodes is used and pressure boundary conditions are applied on the inlet and outlet. The computational results using FLUENT are also shown in Figs. 6–8 and Table 4. Figs. 9 and 10 are the contour plots of the temperature obtained by BEM and FLUENT, respectively. From Figs. 6–8 and Table 4, it can be seen that the BEM results are close to the results by FLUENT with a good accuracy. However, there is a small discrepancy between these two results near the outlet of the pipe. Table 4 shows the maximum error on the outlet of 2D pipe being 1.24%. 7. Conclusions A novel boundary-domain integral equation for computing the temperature in viscous fluids is presented based on the conservation form of the energy equation. The derived formulation is general and applicable to steady, unsteady, compressible, and incompressible flows. Together with the formulations for the velocity, pressure, and velocity gradient proposed in Refs. [2–4], a boundary-domain integral equation system for solving convective heat transfer problems is established. For incompressible flows, the system of equations is closed and the velocity and pressure integral equations are uncoupled from the temperature integral equation. However, the temperature integral equation involves the velocity, pressure, and velocity gradient, which need
to be computed in advance. Two numerical examples for convective heat transfer in steady incompressible fluids have been provided to demonstrate the correctness of the derived formulations. Acknowledgements The authors gratefully acknowledge the National Natural Science Foundation of China under Grant NSFC Nos. 10872050 and 11172055 for financial supports to this work.
References [1] A. Bourchtein, Exact solution of the generalized Navier–Stokes equations for benchmarking, Int. J. Numer. Methods Fluids 39 (2002) 1053–1071. [2] X.W. Gao, A boundary-domain integral equation method in viscous fluid flow, Int. J. Numer. Methods Fluids 45 (2004) 463–484. [3] X.W. Gao, A promising boundary element formulation for three-dimensional viscous flow, Int. J. Numer. Methods Fluids 47 (2005) 19–43. [4] X.W. Gao, Explicit formulations for evaluation of velocity gradients using boundary-domain integral equations in 2D and 3D viscous flows, Int. J. Numer. Methods Fluids 54 (2007) 1351–1368. [5] M. Ikeuchi, K. Onishi, Boundary element solutions to steady convective diffusion equations, Appl. Math. Model. 7 (1983) 115–118. [6] Y. Shi, P.K. Banerjee, Boundary element methods for convective heat transfer, Comput. Methods Appl. Mech. Eng. 105 (1993) 261–284. [7] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part I: General formulations and 1D implementation, Comput. Methods Appl. Mech. Eng. 192 (2003) 4281–4298. [8] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part II. 2D implementation, Comput. Methods Appl. Mech. Eng. 192 (2003) 4299–4312. [9] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part III: Numerical examples, Comput. Methods Appl. Mech. Eng. 192 (2003) 4313–4335. [10] G.F. Dargush, P.K. Banerjee, A boundary element method for steady incompressible thermoviscous flow, Int. J. Numer. Methods Eng. 31 (1991) 1605–1626. [11] P.K. Banerjee, K.A. Honkala, Development of BEM for thermoviscous flow, Comput. Methods Appl. Mech. Eng. 151 (1998) 43–62. [12] K. Onishi, T. Kuroki, M. Tanaka, An application of a boundary element method to natural convection, Appl. Math. Model. 8 (1984) 383–390. [13] N. Tosaka, N. Fukushima, Integral Equation Analysis of Laminar Natural Convection Problems, BEM VIII, Computation Mechanics Publication, SpringerVerlag, Southampton, Berlin, 1986. [14] C.P. Rahaim, A.J. Kassab, Pressure correction DRBEM solution for heat transfer and fluid flow in incompressible viscous fluids, Eng. Anal. Bound. Elem. 18 (1996) 265–272. [15] K. Vajravelu, A. Kassab, A. Hadjinicolaou, BEM solution to transient free convection heat transfer in a viscous, electrically conducting, and heat generating fluid, Numer. Heat Transfer Part A: Appl. 30 (1996) 555–574. [16] X.W. Gao, T.G. Davies, Boundary Element Programming in Mechanics, Cambridge University Press, Cambridge, 2002. [17] M.D. Greenberg, Application of Green’s Functions in Science and Engineering, Prentice-Hall, Englewood Cliffs, New Jersey, 1971. [18] X.W. Gao, A meshless BEM for isotropic heat conduction problems with heat generation and spatially varying conductivity, Int. J. Numer. Methods Eng. 66 (2006) 1411–1431. [19] R.A. Granger, Fluid Mechanics, Dover, New York, 1995.