Hydrodynamic instability of solar thermosyphon water heaters

S. C. Du

Graduate Assistant.

B. J. Huang

Professor.

R. H. Yen

Professor.

Department of Mechanical Engineering, National Taiwan University, Taipei, Taiwan 10764, Republic of China

Hydrodynamic Instability of Solar

Thermosyphon Water Heaters

The flow instability of a solar thermosyphon water heater is studied analytically. A system dynamics model is derived by means of a one-dimensional approach and a linear perturbation method. The characteristic equation is obtained and the Nyquist criterion is used to examine the flow stability. The parameter M is a dimensionless parameter of system stability. The stability maps are plotted in terms of 14

param-eters. The occurrence of hydrodynamic instability is determined by comparing the stability curves and the designed values ofM. Flow instability is shown not to occur in most of solar water heaters commercially available, because the loop friction is relatively high in the design and because solar irradiation in field operation is still not high enough to cause flow instability.

1 Introduction

The phenomenon of reverse flow has been observed in solar thermosyphon water heaters. Reverse flow can reduce the over-all efficiency and thus should be avoided. Two kinds of reverse-flow mechanisms are noted: namely, thermosyphon saturation and hydrodynamic instability. The reverse flow due to mosyphon saturation results mainly from reverse of the ther-mosyphon head during periods of low irradiation or at a relatively high rate of cooling of the collector side. Reverse flow of this kind can be suppressed by increasing the vertical distance between the tank bottom and the collector top (Vax-man and Sokolov, 1986; Morrison, 1986).

Reverse flow caused by hydrodynamic instability for ther-mosyphon loops of simple geometry has been studied (for example, Greif, 1988; Welander, 1967; Keller, 1966; Creveling etal., 1975; Bau and Torrance, 1981; Huang and Zelaya, 1987). Flow instability may occur and finally results in reverse flow in the thermosyphon loops under certain operating conditions. It had thus been suspected that hydrodynamic instability may also exist and causes flow reverse in solar systems as the solar thermosyphon water heater is a kind of thermosyphon loop. To study analytically the stability of the solar thermosyphon water heaters, Zvirin et al. (1978) used linear system theory and assumed linear distribution functions for the perturbed (transient) components of the temperatures in the collector and the tank (or the heat exchanger situated in the tank) which leads to a simple characteristic equation of the system. This assumption is not exactly true since it implies that the response of the collector temperature to the variation of solar radiation is instantaneous. Zvirin and Greif (1979) used the same as-sumption to study a simple thermosyphon loop that was in-vestigated by Welander (1967), but they failed to verify the hydrodynamic instability predicted by Welander and attributed this error to the assumption of linear distribution function for the perturbed temperatures.

The present study was to investigate the existence of

hydro-dynamic instability in the thermosyphon water heater. Instead of using the approximation of Zvirin et al., a system dynamics model is derived analytically by using a one-dimensional ap-proach and linear perturbation method. A system transfer function, which treats the solar irradiation as the system input and the flow rate as the system output, is derived to represent the system dynamics behavior of solar thermosyphon water heaters. The characteristic equation is then obtained and the Nyquist criterion is used to examine the flow instability. The analysis is then tested by the numerical solutions in the time domain. Dimensionless parameters related to flow instability are derived and stability maps are constructed and used to determine the occurrence ,of instability.

2 Governing Equations

The solar thermosyphon water heater studied is a closed loop, similar to that of Mertol et al. (1981) (Fig. 1). A heat exchanger is designed in the tank to extract indirectly the solar energy absorbed in the collector so that the freezing problem can be avoided by using a working fluid of low frozen point inside the natural circulation loop. In the derivation of the governing equations we have made the following assumptions: (1) A two-node dynamic model is used for the collector. That is, the collector is assumed to be composed of a solid phase (absorber plate) and a fluid phase.

(2) The heat capacity effects of the connecting pipe and exchangers walls are ignored.

(3) The convective heat-transfer coefficients and the phys-ical properties of the circulating fluid (water), except the den-sity in the derivation of buoyancy term in the momentum equation, are all constant.

(4) The axial conduction of the plate and of the fluid along the flow direction is negligible.

From assumption (3), the conservation of mass yields m = m(t). The energy equations derived for the loop are Contributed by the Solar Energy Division of THE AMERICAN SOCIETY OF

MECHANICAL ENGINEERS for publication in the ASME JOURNAL OF SOLAR ENERGY

ENGINEERING. Manuscript received by the ASME Solar Energy Division, Apr. 29, 1992; final revision Aug. 3, 1993. Associate Technical Editor: Z. Lavan.

PpApCPl dTp

' dt

= Q-U

a

(T

p

-T

a

)

UW(TP-T„), for collector plate (1)

Journal of Solar Energy Engineering FEBRUARY 1994, Vol. 116/53

rU„(Tp-Tw), where A -Aw for fluid in the collector -UP(TW-Ta), where A=AU for fluid in the riser

•Ue(Tw-T,), where A = Ae for fluid in the exchanger UP(TW-Ta), where A =Ad for fluid in the downcome

(2) The momentum equation of the loop is

F

»»~I7

=

®

T

"

ev •^dy-Fpsm' where

F

lm

= —J^

}• ( ^C i-'U A wg0\Aw + Au + Ae + Ad ' fm 0 Ffm = —. (3) (4)

The derivation of Eq. (3) uses relation for the frictional head that is expressed by a semi-empirical formula Hj= cmd, where c and d are coefficients determined by a loop friction test suggested by Huang and Hsieh (1985).

The above equations can be separated into the steady-state and perturbed parts according to the linear perturbation method and expressed in the normalized form:

(a) for the steady-state:

0 = q- ua6p- (dp- 6W), collector plate

3dx , - 8m collector - updw, riser -ue(6w-8i), exchanger • UpvWi downcomer 1 d„ex*ezdx=Wj-mfd, momentum N o m e n c l a t u r e (5) (6) (7) LU= l-uv+ Luh Ld = Ldv+ Ldh Fig. 1 Closed-loop thermosyphon solar hot water system

(b) for the perturbed state:

ddB i i i i

rp-f- = q - ufip -(6P-6W), collector plate ddw , 86w -d6„

f \.j \.j —

or dx dx

dp-d'w, where r = rw, collector -up6w, where r = ru, riser

ue6m where/• = /•„, exchanger

Wtm df _ 1 wfm dr v/f„

(8)

(9) -up6„, where r = rd, downcomer 'wex*ezdx-d*fd~lf, momentum. (10) Ac = A = CP = c = d = e = F = FR = F' =

f =

g = H = Hf = h =

/ =

L = M = Ls = m =

K =

ne = Q = area of collector, m

cross-section area of pipe, m2 specific heat, kJ/kg °C coefficients of friction loss head

coefficients of friction loss head

unit vector

coefficient in momentum equation

plate heat removal factor plate efficiency factor flow rate, dimensionless acceleration of gravity, m/s2 height of thermosyphon system, m

loop friction head, m thermosyphon head of each component, dimensionless solar irradiation on collector slope, W/m2

length, m

stability parameter, dimenr sionless

length of whole loop, m mass flow rate, kg/s Nusselt number

number of tubes in exchanger solar energy input per unit

length of collector, W/m; Q = IAAra)e/Lc q = heat input per unit length,

dimensionless

r = time constant, dimensionless T = temperature, °C

t = time, s

U = heat-transfer coefficient per unit length, W/m °C UL = heat loss coefficient of

collec-tor, W/m2 °C

u = heat-transfer coefficient per unit length, dimensionless W = width of the thermosyphon

system, m

w = coefficient in momentum equation, dimensionless x = coordinate along the loop,

dimensionless

y = coordinate along the loop, m p = density, kg/m5

7 = collector tilt angle, deg (3 = thermal expansion coefficient

of water, ° C_ 1 4> = dimensionless variable \j/ = time constant, s T = time, dimensionless 6 = temperature, dimensionless (rot) = transmittance-absorptance product of collector AZ = A?

relative height between ex-changer and collector, m relative height between exchanger and collector, dimensionless steady-state value perturbed value

Subscripts

a

c

cr

d

dv

e

fm

n

P

ref

t

tm

u

uv

w

X

y

z

= ambient

= collector

= critical state

= downcomer

= vertical part of downcomer

= heat exchanger, effective

= friction head in momentum

equation

= solar normal incidence

= collector plate, pipe of riser,

and downcomer

= reference state

= tank

= inertia in momentum

equa-tion

= riser

= vertical part of riser

= water

= loop tangent direction, for

dimensionless coordinate

= loop tangent direction

= gravity direction

Equations (5) to (10) are the governing equations. The di-mensionless variables are defined in Appendix A. The water in the tank (different from the working fluid inside the ex-changer) is assumed to be well mixed with a uniform temper-ature since both the mass of water contained in the tank and the thermal mixing effect is large.

3 Analytical Solution

3.1 Steady-State Solutions. An implicit solution for the

flow-rate of the working fluid / is obtained by solving the steady-state equations

WfnF=KG)+h«(f)-he(f)-h

d

(f). (11)

The variables h's represent the thermosyphon heads in the

various components which are presented in Appendix B. Equa-tion (11) is solved for / numerically by iteraEqua-tion.

By assuming a linear steady-state temperature distribution for the working fluid in the collector and exchanger (Ong, 1974; Zvirin et al., 1977), we obtain an explicit solution for the flow rate:

f=M

i/(rf+i). _where M-lc(q-uad,)

(1 + ua) + lctia/leue' Az =

-hiv + Id

(12) Here M is a parameter that determines the flow rate of the working fluid in the steady-state; 4> is related to heat transfer in both the collector and the heat exchanger; Az represents the height of the collector relative to the exchanger or tank; vtym represents the system friction in the steady-state.

3.2 System Dynamics Model. Taking Laplace transform

of the perturbed Eqs. (5) to (10) with zero initial conditions and solving analytically, we obtain the transfer function of the system

/ ' ( s ) Gn(s)

q'(s) (S> Gd(s)'

where Gd(s) and Gn(s) are the denominator and numerator of G{s) defined as Gd(s)=cAwfl Wfm S+df l)SCiC4C5C6C1(l-piP2PiP4) {hc + hu -he -hd\ (13) G „ ( 5 ) = c3. ( %c + ^ „ - "ge-"gd)- ( 1 4 ) •ha All variables in the above equations are presented in Appendix

C. The flow rate / ' (s) is the system output which is induced by solar radiation q' (the input).

System stability is then determined according to the unstable roots of Gd(s). That is, the characteristic equation of the system is

Gd(s)=0. (15) 3.3 Stability Analysis. According to the Nuquist

crite-rion, if Gd(s) has roots in a finite region bounded by a closed contour Cs along a clockwise direction in the s-domain, the mapping contour of Gd(s) from Cs will encircle the origin in a clockwise direction. As Gd(s) is a transcendental function, the number of roots are infinite. However, the flow response in solar thermosyphon water heaters is generally slow; that is, a solar thermosyphon system behaves as a low-pass filter. The responses at high frequencies are therefore beyond the range of practical applications. It is thus unnecessary to examine the unstable roots in the high-frequency regions. A finite region bounded by a Cs contour that covers a frequency up to 0.5 rad/s is used to examine the existence of unstable roots. It will

be shown later that the frequencies of all unstable roots are typically smaller than 0.1 rad/s.

3.4 Wave-Equation Form of the Governing Equations. The

term ddw/dx in the perturbed Eq. (9) is obtained from the steady-state solutions. Equation (9) is further written as in a wave equation form

90vv dx <adJ J<i Xlvexp -UaX + l\uexp(-upx/f) ^j- = '\eexp(-uex/f) ^Kdexp(-upx/f) ri„, collector 0, riser 0, exchanger 0, downcomer; (16) where

<xw=f/rw, a„ =///•„, ote=f/re, ad=f/re £w= l / >w. L = up/ru, £e = ue/re, £d=up/rd

\w=—k\/rvl, \,= —k2/ru,

Ke = — kj/fej Ad = — K^/Vd - e x p [ - ( l + ua)T/rp]

(17)

x j exp[(l+ ua)v/rp][q (v)+6w(x,v)]dv\. (18)

Relative to the standard one-dimensional wave equation with wave speed a,

— + a — = 0,

3T dx

(19) Eq. (16) has additional terms for heat loss and variation of flow rate. Here aw, au, ac, and ad represent the speeds of temperature wave fronts in each component. According to Eq. (17) variations of the thermal parameters r„, ru, re, rd result in the variation of wave speeds a„, au, ae, ad. This important property is used later to explain the results.

3.5 Time Solutions of the Governing Equations. Except

the first equation for the collector, the equations in Eq. (16) for other components are solved by the method of Laplace transform. The time solutions for the temperature distributions are 6w(x, r) = exp(-ux/f) where a) J0 ( T - X / C < ) f'(v)dv •e(T-x/ct)-\\ f (v)dv (20)

u = up, Bi =62, X = X„, a = a„, riser

u = ue, dj = 03, X = Xe, a. = ae, exchanger (21) u = up, 6, =04, X = Xd, cx = ad, downcomer Here e(j) is the unit step function.

According to the time solution, the temperature at the inlet of each component 0/ propagates downstream with individual wave speeds. The flow r a t e / ' ( r ) appears in the integration parts. The exponential term in Eq. (20) represents a damping factor. Thus the amplitude of the temperature waves is damped as x proceeds. This damping effect increases with increasing heat loss coefficient u and decreasing flow r a t e / . In this case,

the thermosyphon system tends to become more stable. In

contrast, if/is large or u is small, the damping effect is small

and the system tends to become more unstable.

Equation (10) is also solved by Laplace transform to yield

the time solution for flow rate:

f

( T ) =

— exp

$;

exp

l

df

0 0'(x, v)

w

t

x

't

z

dx dv\.

(22)

Here/' (T) is determined by integrating the instantaneous

ther-mosyphon head from T = 0 to r.

An explicit solution of/' ( 0 is evaluated numerically because

of the complexity of the above equations. A finite difference

method is used here to solve Eq. (16). According to the von

Neumann numerical stability analysis for the hyperbolic-type

wave equation, we used the "time-centered implicit" scheme

with second-order accuracy, i.e., the Crank-Nicolson scheme

(Anderson et al., 1984). Several numerical tests using the

schemes of unwind, center difference, leap frog, and

Crank-Nicolson, etc., were performed. Crank-Nicolson scheme was

proved to be the best.

C N _ d I

-_

o Q ) CO \ £ -CT>0~ ^ d : o tu 5 QC -as s Flo w 0.00 5

s :

s

: 0. 0 Case Case Case

t L

E i i i i 11 1 : W/H 2 : W/H 3 : W/H A / +

/^P*

/& 3 - E = 1.5 / i ^ 3 - L = 1.0 ^ = 0.5 ^ ^ y ^ 2-E

^** ^—r=°

i

-

L lc=0.2 le=0.1 7 = 27° uo= 0 . 3 1 6 Linear Approximation ue = 1.660 Exact i i i I i i i i i i 500 Solution C=1.17, d = 1.0 I ! 1 | 1 1 i 1 I I i I 1 | I I I i 1 I 1 1 I | I 1 1 1 1 I 1 1 1 | 1000 1500 2000 2500 Fig. 2 exact Q, ( W / m )

Solution of the flow rate in the steady-state: comparison between solution and linear approximation solution

4 Analytical Results and Discussions

If not stated explicitly, the data presented in Appendix D,

which represents a conventional design of solar thermosyphon

water heaters commercially available, will be used in the

fol-lowing analysis.

4.1 Steady-State Solution. Equation (11) is an implicit

solution for flow rate / in the steady state; Eq. (12) is an

approximate solution. Both solutions are similar for various

Q and W/H (aspect ratio of the solar system) (Fig. 2), especially

for smaller W/H (i.e., the system is shorter). The approximate

solution, Eq. (12), is thus approximately valid. The parameter

i defined in Eq. (12) represents a dimensionless variable to

orrelate the steady-state flow rate/and becomes an important

arameter in the following stability analysis.

4.2 Stability Parameters. The system parameters that

af-act the stability of a solar thermosyphon system are divided

Ho two categories according to the governing Eqs. (5)-(10)

nd Eq. (22):

a) the design parameters, including the geometric parameters,

4, l

c

, W/H, and y; the heat-transfer and friction coefficients

u

e>

u

a

, u

p

, and W/

m

; and the time constants r

p

, r

w

, r

e

, r

u

,

and r

d

. Here, W/H and y are used for convenience instead

of W and H. w,

m

is included since it is function of the

geometric parameters l

e

, l

c

, W/H, etc.

(b) the operation parameters including the input q and the

tank temperature 0, in the steady state.

All the above dimensionless parameters are related to the

system stability since they appear in the system dynamics model,

Eqs. (13) and (14). However, the three parameters q, 6„ and

W/

m

can be grouped together to yield a new parameter M which

is defined in Eq. (12). To prove this, the neutrally stable value

of Mare evaluated for three design cases with (0, = 0) or without

(6,9^0) heat losses in the tank and varied loop friction. The

heat losses from the tank to the surroundings are assumed to

be zero (0, = 0) for Cases 1 and 3, the frictional parameters F

tm

and Fj

m

for Case 3 are twice those in Case 1. The heat loss

from the tank is not zero (0, = 0.5) for Case 2, but F,

m

and F

fm

for both Cases 1 and 2 are identical.

The neutrally stable curves for three design cases are shown

in Fig. 3 to be identical. Other calculations for various designs

show the same trend. Hence, M has a similarity property and

can be used as a parameter to characterize the stability of a

solar thermosyphon system. Therefore, 14 parameters

char-acterize the system stability and the neutrally stable value of

M, denoted by M

cn

for a given solar system, is written in the

form

S r l - 8- •8-Unstable Region r p= 7 5 0 . 8 6 , re = 3 0 1 . 9 5 ru= rd = 4 5 . 6 7 * » » • » * Case 1 •j + Case 2: - A CaBG 3: Cose B * Stable Region » Case C 8t- 0 B,-.5

3, = 0 Doubled friction of Cass 1 and Case 2

„ wl m/ wf m= 1 0 Cass D Case E j r = 0 2 le=o!l , W / H = 0 . 4 7 = 2 7 ' ua = 0.316 Ue = 1.660 * • » « « Real Operation^ Region i i i i \ \ \ \ i i > i I i \ \ i i M \ i \ i i i i i I i i i i i i 50 100 150 .200 250 300 350 Fig. 3 Effect of r„

(

W M

h, 4. 7J» y> ue> ua, up; r„, rw, re, rU! rd,—)= 0. (23)

n wf„ A neutrally stable curve is drawn in terms of the function

parameters in Eq. (23) that divide the system performance into stable and unstable operation according to the stability map (Fig. 3). Mcr represents the maximum permissible flow rate in the steady-state for the solar system to become hydrodyn-amically stable, s i n c e / i s function of M, Eq. (12). Basically, M represents the operating conditions of the solar system be-cause the functional relation includes terms for solar irradia-tion and heat loss from the tank. According to Fig. 3, for a system operated at a larger value of M, larger / is expected with smaller damping effect (Section 3.5). Hydrodynamic in-stability occurs when the system parameter M exceeds its crit-ical value Mcr ( M > Mcr), i.e., when Mis located in the unstable region of the stability map. The hydrodynamic stability of a solar thermosyphon water heater is determined by verifying the location of M in a stability map such as Fig. 3.

_For fixed M, the heat input Q increases with increasing (T,-Ta) and friction factor v/fm. Hence, the solar system tends to be unstable for larger solar irradiation (# high), smaller tank loss and loop friction (ua, 9„ wfm small). This result coincides with that of Huang and Zelaya (1987) for a simple rectangular loop.

Regarding the effect of rw on stability, according to Fig. 3,

Mcr=1.132, q'=5.0E-06 Case A ; M= 1.576 Case B : M=1.132 Case C : M=0.727 lc=0.2, le=0.1, W/H=0.4 . p =750.86, rw =94.36, wtm /w,m =10, 7 = 2 7 ° =301.95, ru = rd =45.67, u„=0.31 6, u, =1 ,i 600 I I I I I I | 1200 1600 2400 Time ( sec ) 3600

Fig. 4 Numerical solutions in the time domain

values of Mcr first increase with increasing rw, then decline after reaching a maximum value. This phenomenon is ex-plained from the numerical solution in the time domain. Flow responses / ' (T) for five design cases shown in Fig. 3 were carried out to verify the stability for a step input q'. Figure 4 shows that Case C is stable as the oscillation is damped out; Case A is unstable as the amplitude oscillates and increases; Case B is a neutrally stable as oscillation persists. The frequency of the oscillation for unstable cases is about 0.1 rad/s identical with the unstable roots found in the Nyquist analysis.

The time solutions for Cases A, D, and E are presented in Fig. 5. The solutions for Case A and E reveal unstable oscil-lations predicted by the stability map in Fig. 3 since the design points are located in the unstable region. The time solution for Case D is stable. The slight oscillation is due to numerical error as the M value for Case D is near Mcr. The stability map shown in Fig. 3 is thus verified by the time solutions of Fig. 5.

A maximum value of Mcr in the neutrally stable curves of Fig. 3 needs explanation. Welander (1967) found that insta-bility of a thermosyphon loop resulted from the phase differ-ence between the thermosyphon buoyancy and friction forces; both are related to the wave speeds. Variations of r„ alters the wave speed a„ (=f/rw). The phase difference between the thermosyphon buoyancy and friction forces is then altered. For given M, continuously increasing r„ results in continuous phase shift from out-of-phase to in-phase, then to out-of-phase mode again ( + 2ir from the first out of phase). The operating point thus moves along a line of constant M line from an unstable to stable region, then again to an unstable region (Fig. 3).

From these maps, we can determine the stability of a solar thermosyphon system by comparing the value of M in the operation with Mcr. The values of M for most solar water heaters commercially available are calculated here. In practice, the maximum solar irradiation / is about 1200 W/m2, Q is thus < 2400 W / m (for a two-panel collector with standard size 2 m x 1 m each). For loop friction coefficients chosen from the value determined by Huang et al. (1990) and the design parameters presented in Appendix D, M = 0.073, which is far below the neutrally stable curve in Fig. 3. The system designed with parameters given in Appendix D is thus in the stable region. Although the operating parameter Mmay be increased by decreasing the loop friction, the loop friction in practical solar thermosyphon hot water, heaters is always finite and cannot be decreased greatly. The practical design of solar

ther-3.OE-0O6 - , O.0E+0OO M = 1.576, q'=5.OE-06 Cose A : Mcr=1.132, rw =94.36 Case D : Mcr=1.772, r „ =226.47 Case E : Mcr=0.813, r , =283.03 lc=0.2, le=0.1, W/H=0.4 wt m/ w( m= 1 0 , 7 = 2 7 ° Case E l M | l l l | l l l | l l l | l l i | l l l | l l l | l l i [ I l l | l l l | l I l | l l l | 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 2400 Time ( sec )

Fig. 5 Numerical solutions in the time domain

mosyphon water heaters cannot diminish more than half the

friction presented in Appendix D. This behavior indicates that

the values of M of practical solar systems are invariably far

below the neutrally stable curve. A hydrodynamic instability

is thus unexpected for most solar thermosyphon water heaters.

4.3 Effect of Parameter Variations on Stability. The

par-ametric study of stability for solar thermosyphon systems is

complicated according to the characteristic equation, Eq. (23).

Here, only two of 14 parameters are varied in a stability map

and the others are kept constant. The effects of the latter 12

parameters are discussed separately.

(1) Effect ofr

e

, r

u

, and r

d

. The magnitude of parameters,

r

e

, r

u

and r

d

affect the wave speeds transferred in the loop. As

stated above, the wave speeds vary with the phase differences

between the thermosyphon buoyancy and friction forces and

also affect the stability. The effect of r

e

is shown in Fig. 6.

Three design cases were studied. The wavelike stability curve

for Case 3 is noted, explained previously, but rapid variation

of phase difference results in two maxima for Case 3.

Although r

e

in Case 3 is larger than in Cases 1 and 2, the

stability curve of Case 3 is not necessarily higher than that of

Case 2 because of the variation of wave speeds and phase

differences. Although M

cr

varies rapidly with r

e

and r

w

, the

operating value of M for practical solar water heaters is still

far below the neutrally stable curves as shown in Fig. 6. The

solar water heaters commercially available are thus invariably

hydrodynamically stable.

The effects of r

u

and r

d

are similar to those of r

w

and r

e

. As

the range of variation are narrow, their stability curves are

similar to that in Fig. 3 and are not presented here.

(2) Effect of w,

m

/w/

m

. In general, w„„ is larger than w

/m

.

Five ratios w

tm

/Wf,„ were studied. According to Fig. 7, the

larger is the ratio w

lm

/w/

m

, the lower are the stability curves.

Hence, the system tends to become more unstable (the stable

region is reduced). This property is indicated by the exponential

damping term in Eq. (22). As large ratio w„„/Wf

m

has a small

damping effect, the system tends to become more unstable.

For Case 5 in Fig. 7, the stability curve is low and the system

tends to become more unstable. In practical designs, the ratio

w

tm/w/m is invariably finite and the values of M for practical

operations of thermosyphon water heaters are thus always far

below the neutrally stable curves, therefore solar systems are

hydrodynamically stable.

(3) Effect of W/H, l

c

and l

e

. Figure 8 shows the stability

curves for four ratios W/H. The larger is the ratio W/H, the

higher is the stability curve. Hence, the system tends to be

rp=750.86 A A A a i D a • n D D D ^ * A A A + + + + + + + + Unstable Region * 4 . o a o ° ¥ * + + + Stable Region A A A A ft a D D a + + + ru= rd= 4 5 . 6 7 A* wt m/ wf m= 1 0 A ' ' ' l c = 0 . 2 le=0.1 n° oa W / H = 0 . 4 D 7 = 27° ua= 0 . 3 1 6 u e = 1 . 6 6 0 + Case 1 : re= 3 0 1 . 9 5 A Case 2 : re= 6 7 9 . 4 0 a Case 3 : re = 1 2 0 7 . 8 1

Real Operation Region 5 0 i I i i i 150 200 1 I I 250 350 Fig. 6 Effect of re Case Case Case Case Case wtm/wfm = | w tm / w fm = 8 W tm / W fm = I 0 W tm / W fm = 1 5 W tm / W fm = 2 0 Unstable Region Case 2 na aa a aD Cose 3 +4- + -I- + + + 4 Cose 4 A A A A A A A A Cose 5 " o o o o o o a Real Operotioin Region

rp= 7 5 0 . 8 6 , re = 3 0 1 . 9 5 r = rH = 4 5 . 6 7 l c = 0 . 2 le=0.1 W / H = 0 . 4 7 = 27° ua= 0 . 3 1 6 ue=1.660 • t t t t Stable Region * • * . ' I ' M 50 100

PF

l | l l i i | l i i i \ l i i i ( i 150 200 250 300 T - r - 1 350 Fig. 7 Effect of wlm/m,m

Case 1 Case 2 Case 3: Case 4 W/H=0.6 W/H=0.5 W/H=0.4 W/H=0.3 Unstable Region + * 3 Case 1 + + + A - A T- Cos. 2 . - * " " , ^ - a A A A A a K' 2 : « " * I Cose 3 „ „ x * . Stable Region ° u ,

Real Operation Region

r „ = 7 5 0 . 8 6 , r. = 3 0 1 . 9 5 ru= rd = 4 5 . 6 7 wt r„ / w ,n, = 1 0 lc=0.2 le=0.1 7 = 27° ua=0.316 ue=1.660 • + + I A A I I I I | I I I I | I I I I | I I I I | I I I I | I I I I | I I I I | 0 50 100 150 200 250 300 350 Fig. 8 Effect of W/H

Bode plot of Solar T h e r m o s y p h o n System Dynamics G(s)

io-« 10-5 10< 103 i o -Frequency, r a d / s

10-1 100

10-6 10-5 10-* 10-3 10-2 Frequency, r a d / s

Fig. 9 Bode plot of transfer function G{s)

10-1 10°

more stable (the stable region is broadened). This behavior coincides with the conclusion of Huang and Zelaya (1987), who found that a larger distance between the heater (collector) and the cooler (exchanger) (larger Az, smaller W/H) results in a system becoming more unstable. This property is explained by the damping effect in Eq. (20). As a larger Az results in larger M and larger / , the damping effect is smaller and the system is more unstable.

The effects of lc and /„ were examined with a stability map similar to Fig. 6 but not presented here. The solar system is shown to be more stable for larger /<, or le because the longer collector or heat exchanger creates a larger damping effect.

The neutrally stable curves (Fig. 8) for various W/H, lc, and le are far above the values of M for practical conditions of most commercial solar water heaters. This property indicates that the solar system is hydrodynamically stable.

(4) Effect ofue, ua, rp andy. Stability curves for various ue are similar to Fig. 3 and are not presented here. The system tends to be more stable for larger ue as the damping effect shown by Eq. (20) is greater for large ue. Again, most com-mercial solar water heaters operate at values of M far below the stability curves, hence no flow instability occurs.

The effects of 7, rp, and u„ on the hydrodynamically stability are minimal since the variation of these parameters is small in

practical design. Commercial solar water heaters operate at values of M far below Mcr, thus no flow instability occurs.

According to the above discussion, commercial solar water heaters invariably operate at values of M far below the neutrally stable value Mcr and are thus invariably hydrodynamically stable. This property is mainly due to larger loop friction and small irradiation in practice.

4.4 System Transfer Model of a Solar System. The dy-namic behavior of a solar thermosyphon system is clearly shown by the system dynamics model, Eqs. (13) and (14), using a Bode plot of G(s). The design parameters presented in Ap-pendix D are used in the calculation and the Bode plot presented in Fig. 6 shows that the solar water heater is an overdamped system with low-pass behavior. The gain drops at 20 db/decade approximately for frequencies up to 0.01 rad/s. Furthermore, a more rapidly decreasing rate is seen at frequencies greater than 0.01 rad/s. The phase lag approaches - 100 deg at to = 0.01 rad/s, and - 2 7 0 deg at w>0.1 rad/s. Hence, the dynamic behavior of the solar system is approximately second order with a delay at frequency co<0.01 rad/s, but alters to a system of higher order at frequency w>0.1 rad/s. Nevertheless, the dynamic system is stable.

According to this typical Bode plot, the upper limit fre-quency (<0.5 rad/s) used in Section 3.3 to map s contour to

Journal of Solar Energy Engineering

FEBRUARY 1994, Vol. 116/59

Gd{s) plane is adequate. Although unstable roots may exist in a region of higher frequency, this possibility can be ignored since the system possesses low-pass properties and the gain drops to a value much smaller than - 4 0 db. Morrison et al. (1980) showed that a step variation of solar radiation to a thermosyphon hot water system resulted in a low-frequency oscillation of flow that was rapidly damped. The system was shown to be stable similar to the present case shown in Fig. 9.

5 Conclusion

Reverse flow commonly observed in solar thermosyphon water heaters was suspected to be due to hydrodynamic insta-bility. A theoretical analysis was carried out with a one-di-mensional approach and a linear perturbation method to derive a system dynamics model. A transfer function that treats solar irradiation as the system input and flow rate as the system output was derived to describe the dynamic behavior of solar thermosyphon water heaters. The characteristic equation was then obtained and the Nyquist criterion was used to examine the flow instability. According to the governing equations and the steady-state solutions, the parameter M can be a dimen-sionless parameter for system stability. The stability maps are then plotted in terms of 14 parameters. The occurrence of hydrodynamic stability is determined by comparing the sta-bility curves with the designed values of M. Flow instasta-bility cannot occur in most solar water heaters commercially avail-able, because the loop friction is high in the design and solar irradiation in field operation is still not great enough to cause flow instability. The reverse flow observed in some solar ther-mosyphon water heaters is shown not to be due to hydrody-namic instability, it may be caused by thermosyphon saturation.

Acknowledgment

The present study was supported by Energy Commission, Ministry of Economic Affairs, Taiwan, R.O.C. through Grant No. 782J1 and by National Science Council of the Republic of China, through Grant No. NSC78-0401-E002-18.

References

Anderson, D. A., Tannehill, J. C , and Pletcher, R. H., 1984, Computational

Fluid Mechanics and Heat Transfer, McGraw-Hill, New York.

Bau, H. H., and Torrance, K. E., 1981, "On the Stability and Flow Reversal of an Asymmetrically-Heated Open Convection Loop," J. Fluid Mech., Vol. 106, pp. 417-433.

Creveling, H. F., De Paz, J. F., Baladi, J. Y., and Schoenhals, R. J., 1975, "Stability Characteristics of a Single-Phase Free Convection Loop," J. Fluid

Mech., Vol. 67, Part 1, pp. 65-84.

Greif, R., 1988, "Natural Circulation Loops," ASME Journal of Heat

Trans-fer, Vol. 110, pp. 1243-1258.

Huang, B. J., and Hsieh, C. T., 1985, " A Simulation Method for Solar Thermosyphon Collector," Solar Energy, Vol. 35, pp. 31-43.

Huang, B. J., and Zelaya, R., 1987, "Stability Analysis of a Thermosyphon Loop," ISES Solar World Congress, Hamburg.

Keller, J. B., 1966, "Periodic Oscillations inaModel of Thermal Convection,"

J. Fluid Mech., Vol. 26, Part 3, pp. 599-606.

Mertol, A., Place, W., and Webster, T., 1981, "Detailed Loop Model (DLM) Analysis of Liquid Solar Thermosyphons with Heat Exchangers," Solar Energy, Vol. 27, No. 5, pp. 367-386.

Morrison, G. L., and Ranatunga, D. B. J., 1980, "Transient Response of Thermosyphon Solar Collectors," Solar Energy, Vol. 24, pp. 55-61.

Morrison, G. L., 1986, 'Reverse Circulation in Thermosyphon Solar Water Heaters," Solar Energy, Vol. 36, No. 4, pp. 373-379.

Ong, K. S., 1974, " A Finite-Difference Method to Evaluate the Thermal Performance of Solar Water Heater," Solar Energy, Vol. 16, pp. 137-147.

Vaxman, M., and Sokolov, M., 1986, "Effects of Connecting Pipes in Ther-mosyphonic Solar Systems," Solar Energy, Vol. 37, No. 5, pp. 323-330.

Welander, Pierre, 1967, "On the Oscillatory Instability of a Differentially Heated Loop," J. Fluid Mech., Vol. 29, pp. 17-30.

Zvirin, Y., Shitzer, A., and Grossman, G., 1977, "The Natural Circulation Solar Heaters—Models with Linear and Nonlinear Temperature Distributions,"

Int. J. Heat Mass Transfer, Vol. 20, pp. 997-999.

Zvirin, Y., Shitzer, A., and Bartal-Bornstein, A., 1978, "On the Stability of the Natural Circulation Solar Heater," Proc. 6th Int. Heat Transfer Conference, Toronto, Canada.

Zvirin, Y., and Greif, R., 1979, "Transient Behavior of Natural Circulation Loops: Two Vertical Branches with Point Heat Source and Sink," Int. J. Heat

Mass Transfer, Vol. 22, pp. 499-504.

A P P E N D I X A

Definitions of Parameters and Dimensionless Variables

pp, i P\v^w^—pw ^ . Pw^e^pw, , y/w — y y/e —

u

w

u

w

'

'/'» = ; Wd=—-— (A 1)

u

w

u

w T

-Q-.

a

2_

\JW i refL-/iv T ' T ' / ' 1 ref J-'s ^s Ua_.

u'

Ue ~UU

u

t • up = -f (A2) Ut. t_. _^P_. ___^V. _JU_. _}KL. _TU_ , . , v 7 — . i "p — . ) r«i — . > re — f > r"~ f ' d ~ f \A i)

'ref 'ref ^ref ^ref ^ref ^ref

f_mCpw V>tm = Fin, Uw

F

f

J-

fuX\

d

*-7?w**ref* •* re W/m = " Trrt-L, -. (A4) A P P E N D I X B

Relations of Parameters and Coefficients in the Steady-State Solution (11) hc(f)= ( I - 6 1 ) _{* . - *}

M

(sin 7); *a J hu(f)=e2-(l-b5) (Bl) up

^(7)=-(i-&

3

)(0

3

-0,)+0,4; M ? ) = 0 4 - ( i - * 6 ) (B2)

Ue Up

0 , = [ ? / Uo( l - & , ) . Z ?2. V & 4 + 0 / ( l - 6 3 ) - 6 4 1 / ( 1 - V & 2 ' & 3 « & » )

(B3) e2={q/ua(l-bi) + e,>bl>(l-bi)-b4}/(l-brb1'bi.b4) (B4) e3={q/ua(l-bi)-b2 + et-bi-b2-(l -b3)'b4]/(l-brb2'b3'b4) (B5) e4=lq/ua(l-bl).b2.bi + el(l-b3)}/(.l-brb2-b3-b4) (B6) bx = exp{ - u jc/ \ f { \ + ua)})\ b2 = e x p ( - uplu/ f ) ; &3 = e x p ( -M e/e/ 7 ) . (B7) b4 = exp(-upld/f); b5 = exp(-upluv/f); b6 = exp(-upldv/f)

(B8)

A P P E N D I X C

Definitions of Parameters in Eqs. (13) and (14)

ki = - "• [ei-q/ua); k2 = 62up/f;

J\l + Ua)

ki = ut{fii-B,)/f\ k4 = 64Up/f (Cl)

Ci = V p S2 + [rw (ua + 1) + /-Js + M„; 02 = ^ 5 + (ua+ 1); c ^ o , - . ? (C2) c4 = r„rp(l + ua)s + [/•„,( 1 + ua)2 + r„]; c5 = /vs + up; c6 = res+ue; c1 = rds + up (C3) d'i (s) = [ ( c , z3) • / ' + ( c3u3) ^ ' ] / 2 i ; «2 (^) = [ ( c i z4) ' / + ( c3«4) ' « r ' ] / « i (C4) eUs) = l(c1z5)'f + (,c3v5)'g']/zi; Oi(sy=Uc1zs)-f + {c3v6)'q']/zi (C5) Vi = {\-P\)PiPiPi\ v4 = (l-p1); V5 = (l-Pi)Pi; v6 = (l-pl)p2pi (C6) Pi = e x p [ - (cilc)/(Jc2)]; p2 = exp(-cslu/f); p3 = e x p ( - c6/e/ 7 ) ; (C7) p4 = exp(-c1ld/f); p5 = e x p ( - c5/ „!, / 7 ) ; p6 = exp(-cyldv/f) (C8) *1 = (1 -PlP2PlP4)CiC3 (C9) " * 2 ,

Loop Friction Test : In H f = In c + d In m System A : c = 8 . 2 7 , d = 1.35 System ~ ' " System System C D m = 0 . 0 2 m = 0 . 0 4 1 1 1 i 1 I 1 1 1 -3.8 - r - r - p r -3.3 = 0.08 ~>~r -2.1 In m, m : Kg/sec Fig. 1A Friction loss head

z3 = kxc2{\ + uB) (61 -pi )p2p3p4 + _'(b2-P2)P3P4 k3 k4 + — (b3-p3)p4 + — (b4-p4) Q (CIO) /i? c = sin 7 (1 - / ? , ) p2p3p4fc2 lcCi-fc2(l-pi) z4 = ^ic2(l + M„) (61 - . p , ) +pi (b2-p2)p3p4 + ~ (b3-p3)p4 + — (b4-p4) ~e rd Zs=piP2)—(b3-pi)p4 + —{b4-p4)\c4 ire I'd C4 ( C l l ) + klc2(l + ua)(bl-pi)p2 + - (b2-p2)c4 (C12) z6=PiP2P3—(b4-p4)c4 + kiC2(l + ua)(bl-pl)p2p3 I'd + -(.b2-p2)p3c4 + -(b3-p3)c4 (C13) (1 -pip2p3p4) c5c6c7 (C18) / V = ( 1 - A ) ( 1 - / ? 5 ) 7 < W7 (C19) * C = (1 - P i ) ( l -p3)P2fc\CSc1 (C20) / V = ( l - A ) ( l -Pe)P2P3fc]c5c6 (C21) A P P E N D I X D S y s t e m D e s i g n D a t a

T h e working fluid in t h e loop is water. T h e basic d a t a for the design p a r a m e t e r s are given here.

(1) Ua a n d U„.

T h e heat-transfer coefficients, Ua a n d U„, a r e evaluated

from the collector test results according t o A N S I / A S H R A E 93-1986 S t a n d a r d . As quasi-steady Ua a n d Uw is a s s u m e d , b y

c o m p a r i n g t h e steady-state terms of the energy e q u a t i o n with the Hottel-Whillier-Bliss e q u a t i o n , we o b t a i n hc = z3fc2(\ -p{) + &ic2(l + ua)

7(1 +"J

( l - 6 i ) c i „ AeUL F' (ACUL Ua~ Lr ' Uw~l-F'\ L, ( D l ) -fc2(l-Pi) k2 (1 -pip2p3p4) 'C^c-, sin 7 (C14) hu=]z4f(l-p5) + -0-bi)cs u„

where F' a n d UL are related t o t h e p a r a m e t e r , FR(roi)„ and

FRUL, defined in t h e Hottel-Whillier-Bliss equation t h r o u g h

the following relation:

F'U, = - GnCnw In 1 -Al-Ps) C4(l ~P\P2PT,Pi) CtCtCj (C15) he=)zif(l-p3 )+-f(l-b3)

-70

-to)

-c6 c4(l -Pip2p3p4) C1C5C7 (C16) hd= \z6f(l-p6) + r / [ I - 6 e rd

-7(1 -ft)

-c7 C4(l ~PlP2PlP4) CiC5C6 (C17) F ' ( r a ) „ = ^ ^ F ' [ / L (D2)

where G0 is the mass flow rate per unit area of collector. (TOC)„

is t h e t r a n s m i t t a n c e - a b s o r p t a n c e p r o d u c t of t h e collector at solar n o o n . F o r t h e solar collector commercially available that

we u s e d , FR(T<X)„ = 0.642, FRUL = 4.22 W / m2 ° C , ^c= 3 . 7 8 5

m2,Lc= 1.96 m , a n d (ra)„ = 0.85. H e n c e , U„= 11.0 W / m °C

and U„ = 34.8 W / m ° C .

(2) I/,.

F o r simplification, we neglect the heat loss from riser and d o w n c o m e r , i.e., UP = Q.

(3) Ue.

Assume that there are n

e

tubes in the exchanger, that the

flow is a fully developed laminar in the tubes, and that the

tank temperature is constant. We have Nusselt number

N

u

= 3.66 from convection theory; hence £4 = 7.2 • /?

e

W/m ° C.

Here the thermal conductivity of water is taken as 0.628 W/

m °C (at 40°C). With n

e

=&, U

e

= 51.6 W/m °C.

(4) c and d in frictional loss head H

f

.

The frictional loss tests for four typical solar thermosyphon

water heaters were measured (Fig. Al). For System C, c= 1.17,

d = 1.07. For simplicity, let d= 1.

(5) Other design parameters.

L

s

=Wm,L

c

/L

s

= 0.2, L

e

/L

s

= 0.l, W/H=0A,

7 = 27 deg, T

a

= 20°C, T, = 20°C,

Pp

A

!

, = 30kg/m,

C

pp

= 870 J/kg °C. A

W

= 0.00151 m

2

,

A

u

= A

d

= 0.0003S m

2

, A

e

= 0.00251 m

2

.

I w J The American Society of

LHy® Mechanical Engineers

GET ASME COUPON BOOKS

Easy to Use & Economical

Use coupons to purchase all ASME publications ~ including special

publica-tions, codes and standards, and technical papers (preprints). Use coupons to

save money. Technical papers cost less when you purchase them with coupons!

One coupon may be redeemed for one technical paper. That's a savings o f$.50

for members, $1.00 for non-members (off the regular price for preprints).

Technical Papers (preprints) Coupon Book

Contains 10 coupons

Order No. CB0001

$40 (ASME Members) / $80 (Non-Members)

Publications Coupon Book

Contains 10 coupons ($10 each)

Order No. CI30002

$100 (Member & Non-Member)

To order: write ASME Information Central, 22 Law Drive. Box 2300, Fairfield, NJ 07007-2300 or call 800-THE-ASME (843-2763) or fax 201-882-1717.