Removal of tetracycline hydrochloride from aqueous solution by three 3D uranyl-organic frameworks

Applied Mathematical Modelling 73 (2019) 514–525

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Modeling transient fully-developed natural convection in vertical ducts – Method of eigenfunction superposition C.Y. Wang Departments of Mathematics and Mechanical Engineering, Michigan State University, East Lansing, MI 48824, USA

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 27 August 2018 Revised 25 March 2019 Accepted 3 April 2019 Available online 16 April 2019

This paper models the transient natural convection problem for the hydrodynamical and thermal fully-developed open vertical ducts. The powerful method of eigenfunction superposition yields analytic solutions for general duct cross sections. It is found that the induced flow rate approaches a constant in the steady state, but the transient depends on the Prandtl number (Pr). The rise time to reach 95% of steady state is 3/(Pr λ1 ) if Pr < 1, and is 3/λ1 if Pr > 1, where λ1 is the lowest eigenvalue of the Helmholtz equation. The theory is then applied to the circular, semi-circular, rectangular and equilateral triangular ducts. For low Pr at small times, the transient induced velocity is larger in the velocity boundary layer near the walls but the maximum velocity is near the corners of the duct. © 2019 Elsevier Inc. All rights reserved.

Keywords: Free convection Transient Eigenfunction superposition

1. Introduction Natural convection in open ducts is important for devices which utilize buoyancy as a method for heat removal or fluid transport. Examples include electronic equipment, nuclear reactor, insulating panel, solar heating etc. We shall not discuss the literature on the steady-state natural convection, but consider the less-studied transient natural convection. Previous theoretical works mostly considered transient natural convection between two vertical plates, with various unsteady imposed thermal boundary conditions. Numerical finite difference method has been used for plates with a finite height, including entrance effects [1–4]. Bejan [5] gave an estimation of the thermal steady-state boundary layer growth on a vertical plate. Let Y be the vertical distance from the lower edge and l be the boundary layer thickness. Then



Y ∼

RaY 1/4 l, (RaY Pr )1/4 l,

Pr > 1 Pr < 1

(1)

Here RaY is the Rayleigh number based on Y and Pr is the Prandtl number defined as

RaY =

gβ Y 3 T 

αν

,

Pr =

ν α

(2) is the thermal expansion coefficient, T is the difference between the wall

where g is the gravitational acceleration, β temperature and the ambient temperature, α is the thermal diffusivity and ν is the kinematic viscosity. Eq. (1) can be rewritten as

Y l 3 gβ T  = max{1, Pr−1 } l α2 E-mail address: [email protected] https://doi.org/10.1016/j.apm.2019.04.008 0307-904X/© 2019 Elsevier Inc. All rights reserved.

(3)

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

515

Nomenclature A b bj cj dj Fj g j Jn kj l L m n nˆ Nu Pr Q r RaY S t T Ta Tc Tw w x y Y

area aspect ratio coefficients Eq. (22) coefficients Eq. (20) coefficients Eq. (30) sum Eq. (48) gravitational acceleration (m/s2 ) integer Bessel function of first kind coefficients Eq. (45) half width (m) characteristic length (m) integer integer normal coordinate Nusselt number Prandtl number normalized flow rate Eq. (11) radial coordinate Rayleigh number based on vertical distance Eq. (2) non-dimensional energy integral Eq. (12) normalized time non-dimensional temperature Eq. (5) ambient temperature (K) cup temperature (K) wall temperature (K) normalized axial velocity Cartesian coordinate Cartesian coordinate vertical distance (m) α thermal diffusivity (m2 /s) β thermal expansion coefficient (K−1 )  difference λj jth eigenvalue φ eigenfunction θ angle coordinate ν kinematic viscosity (m2 /s) prime dimensional quantity overbar steady state tilde transient state

Now for an open channel, an estimate of the order of the entrance length Y is when the boundary layers on both sides covers the whole channel, i.e. when l is half of the distance between the plates. For example, consider an open channel of 1 cm height and 2 mm distance (e.g. in electronic devices). For a temperature difference of 10 K, Eq. (3) shows Y/l = O(1) for air at 20 °C. Thus, the entrance length is about 1/10 of the channel width, and most of the air in the channel is thermally fully-developed, √ i.e. the temperature becomes independent of the axial distance. Bejan [5] also showed the velocity boundary thickness is Pr = 0.84 times the thermal boundary layer thickness, or the velocity and thermal entrance lengths are of the same order. The entrance lengths for the developing transient problem or for more confined open cylinders are even smaller than that of the steady-state for parallel plates. Similarly for water in a tube of 0.2 mm diameter and 1 mm height (e.g. in microfluidics), one can also show Y/l = O(1). Experimentally the fully-developed state can be detected when the longitudinal velocity or temperature no longer varies with distance from the entrance. The fully-developed transient natural convection between vertical plates was solved analytically using Green’s function [6] or Laplace transform [7–9]. The transient in a vertical annulus was studied numerically [10,11] and analytically [12]. There are few reports for fully-developed duct flow of other cross-sectional shapes. Of importance are the works of Al-Nimr and El-Shaarawi [6] who used Green’s function to solve the suddenly heated circular cylinder and Morini and Spiga [13] who studied the rectangular cylinder using double Fourier series.

516

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

In this paper we present the eigenfunction superposition method applied to uniform wall temperature (UWT) transient natural convection in general open vertical ducts. The efficient method is then illustrated with the circular and rectangular duct, and also applied to the semicircular duct and the equilateral triangular duct. The last two cases have never been studied before, and existing methods would prove difficult to implement. 2. Mathematical formulation of the transient natural convection inside a duct of arbitrary cross section The assumptions are as follows. (1) (2) (3) (4)

The vertical duct of constant cross section is open-ended, and there is no applied axial pressure gradient. The temperature variation is small, such that the Boussinesq approximation is valid, and the flow is laminar. The fluid and thermal properties are constant. The system is quiescent and in thermal equilibrium and until the walls of the duct are suddenly raised to a constant higher temperature. (5) The ensuing natural convection is fully-developed, i.e., the flow and temperature are independent of axial position. We assume originally the fluid and the duct are at ambient temperature Ta then suddenly the walls are raised to a constant temperature Tw . In the fully-developed state where the variables are all independent of the axial distance and lateral velocities are zero, the convection terms are also zero. The energy equation is then reduced to the simpler diffusion equation

∂T  = α∇ 2 T  ∂t Here Tʹ is the temperature, tʹ is the time and ∇ dimensional temperature

(4) 2

is the Laplacian operator over the cross section. Define a non-

T  − Ta Tw − Ta

T =

(5)

and a non-dimensional time

αt 

t=

(6)

L2

where L is a characteristic length which normalizes all other lengths and ∇  2 = L−2 ∇ 2 . Then Eq. (4) becomes

∂T = ∇ 2T ∂t

(7)

In the beginning T = 0, then suddenly T = 1 is imposed on the walls. This problem is the same as the unsteady heat conduction problem of a solid cylinder. The momentum equation is

  ∂ w  = ν∇ 2 w + gβ T  − Ta  ∂t

(8)

Here wʹ is the axial velocity, ν is the kinematic fluid viscosity and similarly due to the fully-developed assumption, the convection terms are zero. The last term in Eq. (8) is the forcing term due to buoyancy. Define a normalized velocity

w=

w L2 gβ (Tw − Ta )/ν

(9)

Eq. (8) becomes

1 ∂w = ∇ 2w + T Pr ∂ t

(10)

The velocity w is zero at t = 0 and also on the walls. The induced flow rate, normalized by L4 gβ (Tw − Ta )/ν , is



Q=

wdA

(11)

where the integral is over the cross sectional area A. Also of interest is the normalized energy integral, representing the amount of heat carried off per time



S=

wT dA

(12)

The energy integral is also the normalized “mixing cup” temperature Tc

S=

Tc

(Tw − Ta )/Q

(13)

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

517

The local Nusselt number, evaluated on the wall, is

   ∂T   Nu =  ∂ nˆ 

(14)

Here nˆ is the normal to the boundary. However, for fully developed natural convection, the temperature is decoupled from the induced velocity, and the Nusselt number is only related to the transient heat conduction, and has less significance. 3. General solution using the eigenfunction superposition method The eigenfunction superposition method [14,15] has long been known as a means to solve problems with a Laplacian operator. It requires the availability of the eigenfunctions and eigenvalues of the Helmholtz equation. The method was used to solve unsteady heat conduction problems [16]. This paper would be the first to apply this method to transient natural convection. Briefly, the eigenfunction superposition method is sketched as follows. Consider the Helmholtz equation

∇ 2 φ j + λ j φ j = 0, φ j = 0 on boundary

(15)

We assume the eigenvalues λj and the eigenfunctions φ j are known for the duct cross section. It can be shown that the eigenfunctions are complete, and different eigenfunctions are orthogonal [14]. Normalize the eigenfunctions such that



φ 2j dA = 1

(16)

We also order the eigenvalues such that λ1 < λ2 < . Now first consider the steady-state (t → ∞) fully-developed natural convection problem. From Eq. (7) and the boundary conditions the steady-state temperature is uniform

T¯ = 1

(17)

From Eq. (10) the steady-state velocity is governed by

∇ 2 w¯ = −1

(18)

This equation is the same as that which governs the steady flow due to a constant pressure gradient. Expand unity in a sum of eigenfunctions

1=

∞ 

c jφ j

(19)

j=1

Using orthonormal properties, the coefficients are



cj =

φ j dA

(20)

Expand also the steady velocity in eigenfunctions

w¯ =

∞ 

b jφ j

(21)

j=1

Eqs. (15) and (18) yield the coefficients

bj =

cj

(22)

λj

Thus the steady state induced flow rate is

 Q¯ =

w¯ dA =

∞ c2  j j=1

λj

(23)

Eq. (23) also represents the steady flow rate due to a constant pressure gradient. For the unsteady natural convection, let

˜ T = 1 − T˜ , w = w¯ − w

(24)

where the tilde denotes the transient part which is zero on the boundary and decays to zero as t → ∞. Such a decomposition seems to be well known, see the early editions of Carslaw and Jaeger [17]. Eq. (7) gives

∂ T˜ = ∇ 2 T˜ ∂t

(25)

518

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

with the condition T˜ = 0 on the boundary and the initial condition

T˜ = 1 at t = 0

(26)

The problem is analogous to the transient conduction problem [17], but we shall solve it using eigenfunction superposition. The solution is simply

T˜ =

∞ 

c j φ j e −λ j t

(27)

j=1

For the transient velocity, Eqs. (10), (18) and (24) give

˜ 1 ∂w ˜ = T˜ − ∇ 2w Pr ∂ t

(28)

Using the Helmholtz equation, the solution is found to be composed of a homogeneous solution plus a particular solution ∞ 

˜ = w

a j φ j e−λ j Pr t +

∞ 

j=1

d j φ j e−λ j t (Pr = 1 )

(29)

j=1

Eqs. (26)–(28) yield

dj =

Pr c j

(30)

(Pr −1 )λ j

˜ = w¯ at t = 0, Eqs. (21), (22), (29) and (30) give Since w

aj = bj − dj =

−1

cj

(31)

(Pr −1 ) λ j

So the induced velocity is

˜ = w = w¯ − w

∞ 

φ j (b j − a j e−λ j Pr t − d j e−λ j t )

j=1

 ∞  c jφ j = 1+ λj

j=1

Pr 1 e−λ j Pr t − e −λ j t (Pr −1 ) (Pr −1 )

(32)

Thus the general transient natural convection problem is solved. The induced flow rate from Eqs. (11), (20) and (32) is

 Q=

wdA =

∞ c2   j 1+ j=1

Pr 1 e−λ j Pr t − e −λ j t (Pr −1 ) (Pr −1 )

λj

(33)

The energy integral from Eqs. (12), (20), (24), (32) and (33) is



w(1 − T˜ )dA = Q −

S= =Q− =

  ∞

∞ c2   j j=1

λj

  ∞ 

b j φ j − a j φ j e−λ j Pr t − d j φ j e−λ j t

ci φi e−λi t dA

i=1

j=1



c j b j − a j e−λ j Pr t − d j e

 −λ t j

e−λ j t φ 2j dA

j=1

1+

∞ 

Pr 1 e−λ j Pr t − e −λ j t (Pr −1 ) (Pr −1 )



1 − e −λ j t



(34)

Notice we used the orthonormal properties to reduce the product of two sums into one sum, and avoided the area integration. Some asymptotic observations are as follows. For large Prandtl numbers

Q∼

∞ c2   j j=1

λj



1 − e −λ j t , S ∼

∞ c2   j j=1

λj

1 − e −λ j t

2

(35)

Note that the transient flow rate Q for infinite Prandtl numbers is analogous to the transient flow rate due to a suddenly applied pressure gradient. For small Prandtl numbers

Q∼

∞ c2   j j=1

λj



1 − e−λ j Pr t , S ∼

∞ c2   j j=1

λj

1 − e−λ j Pr t



1 − e −λ j t



(36)

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

519

Fig. 1. Cross sections of the ducts. Table 1 Convergence of the steady state flow rate Eq. (42). N

2

4

6

8

10

12

Q

0.3893

0.3922

0.3925

0.3926

0.3927

0.3927

In the rare case that Pr = 1, taking the limit on Eqs. (33) and (34) yields

Q=

∞ c2

 j j=1

λj



1 − ( 1 + λ j t ) e −λ j t ,

S=

∞ c2

 j j=1

λj



1 − ( 1 + λ j t ) e −λ j t ( 1 − e −λ j t )

(37)

In what follows we shall apply the general solution to some basic duct shapes, i.e. the circular duct, the semi-circular duct, the rectangular duct and the equilateral triangular duct (Fig. 1). 4. The circular duct and the semi-circular duct Transient natural convection in the circular duct (Fig. 1(a)) was solved by Al-Nimr and El-Shaarawi [6] using Green’s functions. The flow rate was given in double infinite series. Since the circular duct has some importance, we shall use our method to investigate it in more detail. The solution to the Helmholtz equation for a circular boundary is well known, e.g. in membrane vibrations. The axisymmetric eigenvalues are the roots of

J0

 λj = 0

(38)

where Jn is the Bessel function of the first kind. The ordered eigenvalues are

{λ j } = {5.7832, 30.4713, 74.887, 134.04, 222.93, 326.56, . . .}

(39)

The orthonormal eigenfunctions are

φj

 λ jr   π J1 λ j

J0 = √



(40)

where r is the radial coordinate. We also derive

c j = 2π



1

0



φ j rdr = 2

π λj

(41)

From Eq. (23) the steady state flow rate is truncated

Q¯ =

N c2  j j=1

λj

(42)

where N is number of terms in the truncated infinite series. Table 1 shows the convergence of the series Eq. (42) compared to the exact closed form steady-state Poiseuille flow rate π /8 = 0.39270. We see that about 10 terms would guarantee a four-figure accuracy. Eqs. (41) and (42) also suggest Q¯ ∼ λ−2 since the subsequent λj in Eq. (39) are much larger than λ1 . 1 The energy integral S has similar properties. In the transient solutions Eqs. (33) and (34) the decaying exponentials in time are also dominated by the first eigenvalue, and do not affect the convergence rate. Fig. 2 shows the transient growth of Q and S for typical Prandtl numbers 0.7 (air) and 7 (water). Both start from zero and approach the steady state value of π /8. To estimate the starting time (to reach 95% of steady state) set the smaller of λ1 Pr t or λ1 t to 3. Thus, for Pr = 0.7, the starting time is 0.74, and for Pr = 7 the starting time is 0.52. These results are indeed reflected in Fig. 2.

520

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

Fig. 2. Transient increase of flow rate Q (solid curves) and energy integral S (dashed curves) for the circular duct. Top: Pr = 7, bottom: Pr = 0.7. Table 2 Comparison of transient flow rates. t

Q (Pr = 0.7) Eq. (33)

Q (Pr = 7) Eq. (33)

Q (pressure) Eq. (43)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.75 1 2 ∞

0 0.04501 0.1105 0.1755 0.2314 0.2759 0.3096 0.3440 0.3735 0.3923 0.3927

0 0.1472 0.2548 0.3154 0.3493 0.3684 0.3791 0.3870 0.3913 0.3927 0.3927

0 0.1813 0.2745 0.3264 0.3555 0.3719 0.3810 0.3878 0.3915 0.3927 0.3927

On the other hand, let us compare starting flow due to natural convection and starting flow due to a constant pressure gradient. In the steady state both cases are governed by Eq. (18) but the transients are different. For the pressure-driven flow in a circular tube Szymanski [18] used separation of variables to obtain (in our variables)

Q=

π

− 4π

8

∞  e −λ j t j=1

(43)

λ2j

Table 2 shows the starting flow rate due to natural convection caused by a step change of boundary temperature Eq. (33) and that due to a step change in the pressure gradient Eq. (43). Notice that although all three transient flow rates start from zero and approach π /8, the rise is fastest for the flow due to constant pressure, and slowest for natural convection at low Prandtl numbers. Al-Nimr and El-Shaarawi [6] did not publish any numerical values for comparison. The semi-circular duct has not been studied before. Using cylindrical coordinates (r, θ ) as in Fig. 1(b), separating variables and applying the boundary conditions, it can be shown that the orthonormal eigenfunctions are

φ j = k j cos[(2n − 1 )θ ]J2n−1 ( λ j r ) 

k j = (−π /4 )J2n−2 (

(44)

−1/2 λ j )J2n ( λ j )

(45)

The eigenvalues λj are obtained from

J2n−1 (

λj) = 0

(46)

For each n, a sequence of roots of Eq. (46) is computed and the procedure is repeated for many n values and, finally, all roots are listed in ascending order (as indicated in Table 3). The value of n is determined by ordering λj in increasing magnitude. By ordering in terms of λj (instead of n, say) the series sums would have more consistent convergence properties. We also need the coefficients

 cj =

φ j dA =

 π /2  −π / 2

1

0

φ j r drdθ =

n+1

2(−1 ) k j Fj 2n − 1

(47)

Here the integrals Fj can be evaluated numerically or by exact series [19]

 Fj =

1 0

r J2n−1

λ j dr =

∞  4(2n − 1 )(n + l )J2n+2l (

l=0

λj)

λ j (2n + 1 + 2l )(2n − 1 + 2l )

(48)

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

521

Table 3 Eigenvalues λj and the corresponding cj of the semi-circular duct. j

λj

n

cj

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

14.682 40.707 49.219 79.939 95.278 103.50 122.91 152.24 169.40 177.52 178.34 21.967 243.04 246.50 263.20

1 2 1 3 2 1 4 3 2 1 5 4 6 3 2

1.0058 −0.3327 −0.1751 0.1927 −0.01026 0.3112 −0.1327 0.02472 −0.1220 −0.1121 0.09979 −0.02500 −0.07914 0.07860 0.004334

Table 4 The steady asymptote for Q (and S) for the semi-circular duct. N

10

20

30

40

50

60

70

80

90

100 Q¯

7.397

7.422

7.428

7.433

7.434

7.435

7.436

7.437

7.437

Fig. 3. Transient increase of flow rate Q (solid curves) and energy integral S (dashed curves) for the semi-circular duct. Top: Pr = 7, bottom: Pr = 0.7.

Fig. 4. 3D Velocity profiles for the semi-circular duct. Pr=0.7, (a) t = 0.02. (b) t = 0.1 (c) t = 0.5.

Table 3 shows the first fifteen eigenvalues. Fig. 3 shows the transient flow rate and the energy integral. The convergence of the steady state asymptote Eq. (42) for all curves is given in Table 4. In comparison to the transient of the circular duct, the induced flow rate for the semi-circular duct is lower, and the convergence of the series slower. We used N = 80 for a four-digit accuracy. Fig. 4 shows some transient velocity profiles for Pr = 0.7. Notice that for small times the velocity near the walls is higher than that of the interior. In contrast, the maximum velocity is always at the center for the starting flow due to a pressure gradient. Thus this phenomenon is greatly diminished at large Prandtl numbers since the flow approaches that due to a suddenly applied pressure gradient. 5. The rectangular duct Transient natural convection in the rectangular duct was solved by Morini and Spiga [13] using somewhat tedious double Fourier series and Laplace transform. We shall complement their work. Fig. 1(c) shows Cartesian axes are placed at the center of a rectangle of dimensions 2L and 2bL, where b is the aspect ratio. Again by separating variables and applying the

522

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525 Table 5 The first fifteen eigenvalues λj the corresponding (n,m), and the cj values with asterisks. j \b 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

0.25

0.5 ∗

41.946 (1,1)0.8106 61.685 (2,1)−0.2702∗ 101.16 (3,1)0.1621∗ 160.38 (4,1)−0.1158∗ 239.34 (5,1)0.09006∗ 338.03 (6,1)−0.07369∗ 357.78 (1,2)−0.2702∗ 377.52 (2,2)0.09006∗ 416.99 (3,2)−0.05404∗ 456.47 (7,1)0.06235 476.21 (4,2)0.03860∗ 555.17 (5,2)−0.03002∗ 594.64 (8.1)−0.05404∗ 653.86 (6,2)0.02456∗ 752.56 (9,1)0.04768∗

0.75 ∗

12.337 (1,1)1.1463 32.077 (2,1)−0.3821∗ 71.555 (3,1)0.2293∗ 91.297 (1,3)−0.3821∗ 111.04 (2,2)0.1274∗ 130.77 (4,1)−0.1638∗ 150.52 (3,2)−0.07642∗ 209.73 (5,1)0.1274∗ 209.73 (4,2)0.05459∗ 249.22 (1,3)0.2293∗ 268.96 (2,3)−0.07642∗ 288.69 (5,2)−0.04246∗ 308.43 (6,1)−0.1042∗ 308.44 (3,3)0.04585∗ 367.65 (4,3)−0.03275∗

1 ∗

6.3228 (1,1)1.450 26.062 (2,1)−0.4833∗ 37.166 (1,2)−0.4833∗ 56.965 (2,2)0.1611∗ 65.540 (3,1)0.2900∗ 96.384 (3,2)−0.09667∗ 98.853 (1,3)0.2900∗ 118.59 (2,3)−0.09667∗ 124.76 (4,1)−0.2071∗ 155.60 (4,2)0.06905∗ 158.07 (3,3)0.05800∗ 191.38 (1,4)−0.2071∗ 203.72 (5,1)0.1611∗ 211.12 (2,4)0.06905∗ 217.29 (4,3)−0.4143∗

4.9349 (1,1)1.6211∗ 24.674 (2,1)−0.5404∗ 24.675 (1,2)−0.5404∗ 44.414 (2,2)0.1801∗ 64.153 (3,1)0.3242∗ 64.154 (1,3)0.3242∗ 83.892 (3,2)−0.1081∗ 83.893 (2,3)−0.1081∗ 123.37 (4,1)−0.2316∗ 123.37 (1,4)−0.2316∗ 123.37 (3,3)0.06485∗ 143.11 (4,2)0.07720∗ 143.11 (2,4)0.07720∗ 182.59 (4,3)−0.0463∗ 182.59 (3,4)−0.0463∗

Table 6 Convergence of Q¯ from Eq. (23). Asterisked values from Shah and London [20] and hatted values from Morini and Spiga [13]. b \N

10

20

30

40

50

60

Refs [20,13]

0.25 0.5 0.75 1

0.01748 0.1141 0.3038 0.5610

0.01753 0.1142 0.3043 0.5619

0.01754 0.1143 0.3044 0.5621

0.01754 0.1143 0.3044 0.5621

0.01755 0.1143 0.3045 0.5622

0.01755 0.1143 0.3045 0.5622

0.01755∗ 0.01755^ 0.1143∗ 0.1144^ 0.3045∗ – 0.5623∗ 0.5624^

Fig. 5. Transient increase of flow rate Q (solid curves) and energy integral S (dashed curves) for the rectangular duct. For each set, top curves are for Pr = 7, and bottom curves for Pr = 0.7 (a) top set b = 0.5, bottom set b = 0.25 (b) top set b = 1, bottom set b = 0.75.

boundary conditions the orthonormal eigenfunctions are

φ j = cos[(n − 0.5 )π x] cos[(m − 0.5 )π y/b]/ b

Here n, m are integers, determined by the {[(n − 0.5)2 + (m − 0.5)2 /b2 ]π 2 }. The coefficients cj are

 cj =

√ n+m 4 b(−1 ) φ j dA = (n − 0.5 )(m − 0.5 )π 2

(49) ordering

of

the

eigenvalues

λj

chosen

from

the

set

(50)

Table 5 shows the first 15 eigenvalues and related information. Notice, due to symmetry, some eigenvalues are repeated but that does not affect orthogonality [14]. The frequencies are only a subset of those of the vibrating membrane, since only the even modes are included. As for the mode shape, there are 2(n − 1) internal nodal lines in the x direction, and 2(m − 1) internal nodal lines in the y direction. The convergence of the steady state flow rate Eq. (23) is shown in Table 6. We used N = 50 for a four-digit accuracy. Our values are compared to Shah and London [20] and Morini and Spiga [13] after converting their results using our normalization. Fig. 5 shows the transient flow rates and energy integrals for various aspect ratios. Our curves agree with those of Morini and Spiga [13]. In order to investigate the transient development of velocity in more detail, some level curves are plotted in Fig. 6 for small aspect ratio b = 0.25. It can be seen that for short times the maximum velocity occurs near the short ends,

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

523

Fig. 6. Constant velocity curves for rectangular duct b = 0.25, Pr = 0.7 (a) t = 0.025, from outside, w = 0.0 02, 0.0 04, 0.0 06 (b) t = 0.1, from outside, w = 0.005, 0.01, 0.015, 0.02, 0.025.

Fig. 7. Constant velocity curves for rectangular duct b = 0.5, Pr = 0.7 (a) t = 0.03, from a corner towards the interior, w = 0.0 015, 0.0 03, 0.0 045, 0.0 055, 0.0 065, 0.0 065, 0.0 055, 0.0 045, 0.0 03, 0.0 015 (b) t = 0.3, from outside, w = 0.01,0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08. Table 7 Convergence of Eq. (54) when truncated to N terms. N=2

4

6

8

10

12

14

16

18

20

0.7651

0.7769

0.7786

0.7790

0.7792

0.7793

0.7793

0.7794

0.7794

0.7794

while for larger times the maximum velocity moves to the center. For moderate aspect ratios, Fig. 7(a) shows the maximum velocity for short times actually occurs near the corners, while the interior has low velocity. This phenomenon was not reported before. 6. The equilateral triangular duct Fig. 1(d) shows the cross section of the equilateral triangular duct. This geometry is not describable in separable coordinates, and existing methods using separation of variables fail. The eigenvalues were found by Lame and presented by Schelkunoff [21]. Since the equilateral triangle has three-fold rotational symmetry, Wang [16] reduced the eigenvalues to a simpler form. Let Cartesian coordinates be placed at the centroid. The eigenvalues are

λj =

4π 2 j 2 9

(51)

The orthonormal eigenfunctions are

φj

      √ 2 2 j π (2 + x ) jπ y j π (2 + x ) = 5/4 sin cos √ − cos 3 3 3 3

(52)

Eq. (20) yields

cj =

√ 3 5/4 2 jπ

(53)

The steady state flow rate can be summed exactly

Q¯ =

∞ c2  j j=1

λj

=

√ ∞ √ 81 3  1 9 3 = = 0.7794 20 2π 4 j4

(54)

j=1

where we have used a summation formula [22]. The result of Eq. (54) also agrees with the exact solution for the pressure induced flow in an equilateral triangular duct [20]. Table 7 shows the convergence rate is much faster. We used N = 16 for our results. Fig. 8 shows the transient Q,S integrals and Fig. 9 shows the constant velocity lines. From the caption, we see at small times (t = 0.03) the velocity rises from the wall until a maximum is reached near the corners, then the velocity increases to a very low value in the interior. For intermediate times (t = 0.3) the velocity profile gradually becomes fuller, but the maximum is still near the corners, and the local minimum is at the center. For larger times (t = 0.7) the velocities tend to the steady state, with maximum at the center.

524

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

Fig. 8. Transient increase of flow rate Q (solid curves) and energy integral S (dashed curves) for the equilateral triangular duct. Top: Pr = 7, bottom: Pr =0.7.

Fig. 9. Constant velocity curves for the equilateral triangular duct Pr = 0.7 (a) t = 0.03, from a corner towards the interior, w = 0.0 015, 0.0 03, 0.0 045, 0.0 06, 0.0 075, 0.0 06, 0.0 045, 0.0 03, 0.0 015 (b) t = 0.3, For a corner toward the interior, w = 0.015, 0.03, 0.045, 0.06, 0.075, 0.075 (c) t = 0.7, For a corner toward the interior, w = 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 2.

7. Discussions and conclusions This paper models the transient natural convection in fully-developed vertical ducts. This difficult thermo-fluids problem is formulated in general under the afore-mentioned assumptions. If the eigenvalues and eigenfunctions of the Helmholtz equation are available, which is indeed the case for the cross-sectional geometries considered, the method of eigenfunction superposition is simple and efficient. For other cross sections, eigenfunctions may not be available in closed form, but may be found by semi-numerical means such as the variational method. The eigenfunction superposition method is similar to the integral transform method (and its generalizations [23,24]) in the sense that eigenfunctions of the Helmholtz equation is used. However, the eigenfunction superposition method is quite different and easier to implement than the integral transform method, since integral transform pairs, integration of the differential equation, and inversion of results are not needed. We derived analytic solutions for arbitrary open-ended ducts. The results show the following physics (some may be already known but most are new revelations) ∗ The transient temperature is same as that due to heat conduction in a solid by a suddenly applied surface temperature. ∗ The induced transient velocity consists of a sum of two different decaying exponentials. ∗ The induced flow rate (and the energy integral) approaches the steady state flow rate caused by a constant pressure gradient. Due to an extra exponential factor, the energy integral rises slower than the induced flow rate. ∗ The transient rise to steady state also differs. For large Prandtl numbers the induced flow rate matches the flow rate due to a suddenly applied pressure gradient. The rise is slower at low Prandtl numbers. ∗ An estimate of the rise time (to reach 95% of steady-state value) depends on the first eigenvalue. It is 3/(λ1 Pr )if Pr > 1, and is 3/λ1 if Pr < 1. Therefore the rise time is same for different gases, but not for different liquids. ∗ The solutions are expressed in a single infinite series (instead of double series). The convergence of the series depends on the geometry of the cross section. ∗ For the flow caused by a suddenly applied pressure gradient, the velocity maximum is always at or near the center. For transient natural convection, the velocity is larger near the boundaries, especially at low Prandtl numbers and small times. The maximum velocities occur near the corners of the cross section. These ideal theoretical results would be useful for the practicing engineer.

C.Y. Wang / Applied Mathematical Modelling 73 (2019) 514–525

525

References [1] C.F. Kettleborough, Transient laminar free convection between heated vertical plates including entrance effects, Int. J. Heat Mass Trans. 15 (5) (1972) 883–896. [2] H.M. Joshi, Transient effects in natural convection cooling of vertical parallel plates, Int. Comm. Heat Mass Trans. 15 (2) (1988) 227–238. [3] A. Andreozzi, B. Buonomo, O. Manca, Transient natural convection in vertical channels symmetrically heated at uniform heat flux, Num. Heat Trans. A Appl. 55 (5) (2009) 409–431. [4] B. Buonomo, O. Manca, Transient natural convection in vertical microchannel heated at uniform heat flux, Int. J. Therm. Sci. 56 (2012) 35–47. [5] A. Bejan, Convection Heat Transfer, Forth ed., Wiley, New Jersey, 2013. [6] M.A. Al-Nimr, M. El-Shaarawi, Analytical solution for transient laminar fully-developed free convection in vertical channels, Heat Mass Trans. 30 (4) (1995) 241–248. [7] A.K. Singh, H.R. Gholami, V.M. Soundalgekar, Transient free convection flow between two vertical parallel plates, Heat Mass Trans. 31 (5) (1996) 329–331. [8] B.K. Jha, A.O. Ajibade, Transient natural convection flow between vertical parallel plates: one plate isothermally heated and the other thermally insulated, Proc. Inst. Mech. Eng. 224 (E4) (2010) 247–252. [9] M. Spiga, P. Vocale, Step response for free convection between parallel walls, Heat Mass Trans. 51 (2015) 1761–1768. [10] M.A.I. El-Shaarawi, M.A. Al-Attas, Unsteady natural convection in open-ended vertical concentric annuli, Int. J. Num. Meth. Heat Fluid Flow 2 (1992) 503–516. [11] S. Wang, A. Faghri, T.L. Bergman, Transient natural convection in vertical annuli: numerical modeling and heat transfer correlation, Num. Heat Trans. A Appl. 61 (11) (2012) 823–836. [12] M.A. Al-Nimr, Analytical solution for transient laminar fully-developed free convection in vertical concentric annuli, Int. J. Heat Mass Trans. 36 (9) (1993) 2385–2395. [13] G.L. Morini, M. Spiga, Transient laminar natural convection along rectangular ducts, Int. J. Heat Mass Trans. 44 (2001) 4703–4710. [14] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, McGraw-Hill, New York, 1953. [15] K.F. Riley, M.P. Hobson, S.J. Bence, Mathematical Methods of Physics and Engineering, third ed., Cambridge Univ Press, Cambridge, 2006. [16] C.Y. Wang, Transient diffusion in triangular cylinders, Adv. Appl. Math. Mech. 7 (2015) 818–826. [17] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, 2nd ed., Clarendon, Oxford, 1959. [18] P. Szymanski, Quelques solutions exactes des equations de l’hydrodynamique du fluide visqueux dans le cas d’un tube cylindrique, J. Math. Pures Appl. 11 (1932) 67–107. [19] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1965. [20] R.K. Shah, A.L. London, Laminar Flow Forced Convection in Ducts, Acad Press, New York, 1978. [21] S.A. Schelkunoff, Electromagnetic Waves, Van Nostrand, New York, 1943. [22] L.B.W. Jolley, Summation of Series, second ed., Dover, New York, 1961. [23] M.D. Mikhailov, M.N. Ozisik, Unified Analysis and Solutions of Heat and Mass Diffusion, Wiley, New York, 1984. [24] R.M. Cotta, Integral Transforms in Computational Heat and Fluid Flow, CRC, Boca Raton, FL, 1993.