Linear stability analysis of double-diffusive solar pond

Energy ('onvers. Mgmt Vol. 23. No. 2. pp. 67 71. 1983

0196-8904/83 $3.00+0.00 Copyright ~" 1983 Pergamon Press t.td

Printed in Great Britain. All rights reserved

LINEAR STABILITY ANALYSIS OF DOUBLE-DIFFUSIVE SOLAR POND M. S. S O D H A and A J I T K U M A R Department of Physics, Indian Institute of Technology, Hauz Khas New Delhi 110 016, India (Received

18 M a y 1982)

Abstract In this communication, the stability of the double-diffusive solar ponds has been investigated in the linear approximation. The corresponding linearized system of equations of motion is reduced to a single integro-difl'erential equation using the Green-function technique. In contrast to the conclusions of Veronis that, initially, the instability occurs as an oscillatory mode and at a value of R 7 (Rayleigh number for temperature) greater than Rs the motion becomes steady, the present analysis shows that, initially, as R r increases from zero but remains considerably less than R s, exponentially growing and decaying modes (steady motion) occur first; for a value of R r more than a critical value R'r, the motion becomes oscillatory. This oscillatory motion may, due to the basic non-linear dynamics of the system, even become aperiodic. Further, for R.~ , zc, the minimum value of R r for which steady motions can occur tends to .~f ~,2. Rs ' where ,~ = K s / K r where K s and K r are diffusivity coefficients for salt and temperature, respectively; as a contrast, according to Veronis, R r ~ a - ~ R s ; a = v / K r, v being the kinematic viscosity.

Stability

Solar pond

INTRODUCTION

certain observations of solar pond instabilities show that the gradient tends to break into small sublayers over the entire gradient region, although the instability occurs only in a small section of the gradient. It follows then that only local properties and local gradients play a key role in the occurrence of instabilities, and that the linear constant gradient solution should be adequate in predicting stability if applied to the appropriate sublayers instead of the entire non-convective zone. Therefore, in what follows, we shall consider only the local gradients and properties of the double-diffusive system and follow the procedure for deducing the basic equations adopted by Zangrando [2].

One of the main problems connected with effective functioning of double-diffusive solar ponds is its instability with respect to vertical motion of the fluid. Disturbances that cause vertical transfer of the fluid will. as well, transfer heat and thus destroy the insulating effect of the solution. The vertical motion of the solution depends upon its thermodynamical and hydrodynamical parameters. The convective motion may be oscillatory or steady. Veronis [1] made a normal mode analysis of this problem and came to the conclusion that the instability in a double-diffusive system initially starts as oscillatory modes and then, at a value of Rr (Rayleigh number for temperature) greater than R s (Rayleigh number for salinity), the motion becomes steady. He also conjectured that, for R s ~ ac, the minimum value of RT for which steady convection can occur tends to R s. Our analysis of stability leads to the conclusion that the instability in a doublediffusive system occurs first in the form of steady convective modes. Afterwards, a transition into oscillatory modes takes place. These oscillatory modes, in turn. due to the nonlinearity of the system, might even become aperiodic. Double-diffusive fluid systems with linear gradients and constant properties have been extensively studied. The conclusion, which one draws while going through the literature, is that neither linear gradients nor constant properties can be assumed over the entire gradient region of a solar pond [2-4] as the temperature variation is fairly large and transport coefficients are fairly temperature sensitive. Further, H u 23"2

MODEL STABILITY ANALYSIS

Consider a layer of fluid of depth D and infinite horizontal dimensions. The layer is subjected to non-constant gradients of temperature and salinity. which are established in the layer at time t = 0: z is positive downwards. Following the standard method of linear theory [5, 6] one may express temperature and salinity as T(x,y,z,t)=

T,(z)+

T*(x,y,z,t).

S ( x , y , z, t ) = S , ( z ) + S * ( x ,

y , :, t )

(1)

(2)

where T c and S,. correspond to the conductive distribution and T* and S* are the convective perturbations. The density distribution in the layer may 67

68

SODHA and KUMAR: LINEAR STABILITY ANALYSIS OF DOUBLE-DIFFUSIVE POND

also be expressed as p(x,y, z, t) = pc(z) + p*(x,y, z, t)

(3)

p,(z)=po{l-~(T,.- To)+fl(Sc-So)} p*(x,y,z,t) =po{-~T* +flS*}

(4)

with

(5)

where

Here, v is the kinematic viscosity, Kr and Ks the kinematic diffusivities for temperature and salinity, respectively, and 1~ is the unit vector along the positive direction of the z-axis, g is the gravitational acceleration. To eliminate the pressure in equation (9), one may express q. Vq = grad { q[------~:} I - curl (curl q)

# aT s,r, # = =1 a~S r.e, # O is the mean density and the suffix zero refers to values at some reference level, e and # may be constant throughout the layer. The fluid may be assumed to be confined to a small layer, so that p * Ap, --~<--=E P0 Po

< 1

--

-curl~ =gl~ - 1 gradp + vWq. P

By taking the curl of this expression, we get

t x gradp + v W q .

Now, as we assume the convective perturbation to be small

1 _l p p,.+p,

1 [lP*+ofp*'~2 ~ p,~ l o p\I ~,, / j .

grad {~} x gradp = grad {~} x gradpc + grad {~} x gradp*

(i) Continuity equation

6p 6-7 + pV.q = 0

(7)

where q is the convective velocity in the fluid with components (u, V, W). As 6p/6t = 0, equation (7) reduces to

aV

aW

divq = ~xx + fffy + - ~ - = 0.

aq + q V ' q = g k - p 1 v P + vV2q. &-

8~ ÷ q" VT = KrWT.

+ grad

-grad

(9)

(iii) Equation for energy conservation aT

= grad {~} x gpck x grad p * ~

xgradp*.

(8)

(ii) Equation for momentum conservation

In our approximation, the second and third terms of the above approximation can be neglected. So, the linearized form of equation (9) reduces to a _ v v2 t

grad p *

(10) Taking the curl of this equation, we get

(iv) Equation for salt conservation aS

(12)

Therefore,

-po ~-

is the Lagrangian derivative of the relative density variation (P*/Po) of a parcel of fluid. The variation of this term, in our approximation, can be at most of order e and is, therefore, neglected in the governing equations. The basic equations are

au

~-+grad

(6)

where Ap, is the maximum variation of the conductive distribution over the whole layer. Further, following Veronis [1], it is assumed that velocity q at the onset of the convective state is sufficiently small so that it squares and products can be neglected. Similarly, we neglect the second order terms in the temperature and salinity perturbations. The term

po[dt + g r a d p . q

and introduce vorticity t / = curl q. Then from equation (9), we get

a~ ÷ q" VS = KsV2S.

(11)

{6 _ vV2}

(-V:q) =g--'V pc

× (Vp* × 1~).

(13)

SOl)HA and KUMAR:

LINEAR STABILITY ANALYSIS OF DOUBLE-DIFFUSIVE POND

velocity must vanish. Further, due to the continuity condition

The z-component of this equation is (14)

72 W - ?z - = 02

where

p2 {72 V -I = - -?...2 + 03,-

z=0,1

at

(22)

must also hold. For simplicity, one may put

is the Laplacian in the two horizontal dimensions. Similarly, the equations for temperature and salt concentration are written as

?T*

69

(w, 0, z ) = [ w(z, r), O(z, r), Z(z, ~)} •e x p [ i ( K , . x + k , . y ) ] .

(23)

Then from equations (18), (19) and (20), we obtain the following set of equations

dT,, + Ww-~ = KrV2T * dZ

(15)

?S* d& ~ - + W ~-~ = KsV2S *.

(16)

=Rra20-Rsa2"£

(24)

Let us now non-dimensionalize the equations (14), (15) and (16). For this purpose, one may put D2 t=--r. Kr

Kj W=--W D

-- X

T - : 'tT,.(D)- 7;.(0)}.0 = A T . 0 (17)

The depth is scaled by D. Dropping the star specifying the perturbed quantities and defining two Rayleigh numbers Rr and Rs for temperature and salinity, respectively, by

Rs -

V2} V2W

D dT,. " W = AT~''d~

~ -£f'V 2}E=

&

-

(18)

D dS, AS,E'dz'W--GzW

where G~ and G2 are the dimensionless gradients of the conductive state, a = v/Kr is the Prandtl n u m b e r and ~ = Ks/K r is the Schmidt number. The boundary conditions for this problem are at

z=0,1

(26)

(21),

as the temperature and salinity are kept constant on the boundaries, the normal component of convective

(27)

dro FI(z, r; Zo%)

ro

(28)

where F~ is the Green-function satisfying the equation

~,r-

-a:

FI = 6 ( z - z o ) 6 ( ~ - r o )

F~=0

~ G[ W

(20)

W=O=E=O

-a2w

x { - W(zo, to)}

(19)

{

dz o

v. K s

-RTV~O+RsV~E

-4 t i"'fr-V2 0=-)

;f }

equations (14), (15) and (16) can be written in dimensionless tbrm

{l(j. ~T ~

=

Further, let us consider the case of constant gradients (G~ = G2 = 1). The solution for equation (25) can be written as 0(z,r) =

V • K I

Z

a 2 = K~. + K~.

fl "g "AS, "D3

Rr --

-- a 2

where

S -= {S,(D) -- S,(O).X =AS,..E.

7 "g "AT,.'D 3

~

at

z=0,1

(29)

(30)

Solving (29) by a standard method [7], satisfying the boundary conditions (30),and using the method of reflection, one obtains F I = (1;'2 w.~/~(~ - %) x [exp { - ( z - zo)214(z - ri0} - exp [ - ( z + zo)2/4(z - z0)}].

(31)

Similarly for equation (26), the solution is given as

£(z, r) =

i' dz o f ndr o F2(z, r; zo, %)

,do

,d~o

× { - w(:0, ~0)} (32)

70

SODHA and KUMAR:

LINEAR STABILITY ANALYSIS OF DOUBLE-DIFFUSIVE POND

with F2 = (J{" '/2/2 x/rr(T - To) x [exp { - ( z

-

2o)2/4j{(r

-

z0) I

- exp {--(z + zo)'-/4,~(z - r0)}].

(33)

Finally, substituting for 0 and Y~from equations (28) and (32) into equation (18) and finding the appropriate Green-function for the operator

Let J{, R s and a be fixed. Then from equations (37) and (38), it is clear that, initially, as Rr is increased from zero but remains considerably less than Rs, A will be positive and equation (37) will have the solution which is the superposition of exponentially growing and decaying modes. It means that, in the beginning, steady convective modes are generated in the system. If Rr is increased further for values of Rr greater than the critical, viz. (39)

R~r = a-,/2 + ~ - , / 2 R s

we arrive at the following integro-differential equation which is equivalent to the system of equations (18), (19) and (20) V2W=

dz

dr0 F3(z, T; z;, z0)

)

d ro

x {aRra20(zo, "c;,) - a R s a 2 S (z(l, %)}

(34)

with F3=(a

'n/2~-%)) x {exp [-(z

- Z o ) 2 / 4 a ( r - %)]

- exp [ - ( z + Zo)2/4ff (r

-- T0)]} .

(35)

Now, we proceed to study the integro-differential equation (34) for the onset of convection. For this purpose, consider r --, %. In this case, the full Greenfunction reduces to the retarded Green-function (Green-function for outgoing waves only). Further

!im°:0dz°

0

the steady convection motion will change into an oscillatory motion. If R r is further increased, the amplitude and the period of these oscillations will increase because of the increasing influence of the temperature field. Expressed in terms of the motion of a typical fluid particle, the explanation is that, during its oscillatory displacement, the particle experiences a restoring force which decreases as R~. increases and hence the period increases. Later on, due to the nonlinearity of the system, these oscillatory modes may even become aperiodic. Further, it is clear from (39) that, for Rs ~ oo, the m i n i m u m value of Rr for which steady monotonic convective modes can occur tends to R s" 9ff -v2. According to the normal mode analysis made by Veronis, it tends to ~r J. Rs. Graphically, these results can be represented in the R s - R r plane as in Fig. 1. The numerical calculations correspond to the values a = 1, g ( = 0.5. The curve 1 corresponds to the initial steady monotonic modes while the curve 2 belongs to the oscillatory regime of motion which occurs when Rr increases beyond the critical value R~-~ 201. The curves 3 and 4 represent the results obtained by Veronis [1] and correspond to the oscillatory and steady monotonic regimes of

d~0 ~ - -

x exp

x / 4 ~ ( r - r0)

I

{

4(z - r0)J W(zo, c 0 ) - ' W ( z )

6

equation (34) reduces to the following differential equation for the z-component of the velocity W

Rr:

Rs/

4

Qt"

d2W dz 2

Rr * R e

5OOO

a2{1 -- a v2" R r + a l / 2 ~ -'/2" Rs} W = 0

(36)

or d2W ~dz -

a2AW= 0

(37)

where A =l-a

I/2. Rr+a'/2Off '/2"R s.

(38)

As A depends upon the parameters 24(-, or, R r and R s, equation (37) can have the following solutions (i) exponential solution if A is positive; (ii) oscillatory solutions if A is negative.

R rc= 201

~1~ iooo

e,

50o0

Fig. 1. The numerical calculations correspond to the values (r = 1, • = 0 . 5 . The curve I corresponds to the initial steady monotonic modes while the curve 2 belongs to the oscillatory regime of motion which occurs when R r increases beyond the critical value R'r ~ 201. The curves 3 and 4 represent the results obtained by Veronis [1] and correspond to the oscillatory and steady monotonic regimes of motion respectively. Finally the curves 5 and 6 represent the results obtained by Huppert and Moore [8].

SODHA and KUMAR: LINEAR STABILITY ANALYSIS OF DOUBLE-DIFFUSIVE POND motion, respectively. Finally, curves 5 and 6 represent the results obtained by Huppert and Moore [8]. According to them, for 0 < R s < R~, 6, no oscillations take place because the salt gradient is too small. Curve 5 belongs to the oscillatory regime of motion. Finally, the conduction state is destabilized by monotonic modes as Rr exceeds R~n6.

CONCLUSION The analysis presented in this communication is based on the exact reduction of the basic system of linearized equations of motion to a single integrodifferential equation for the z-component of the velocity. The analysis of this basic equation, for the onset of convection, leads to the conclusion that, in the beginning, as Rr is increased from zero but remains considerably less than Rs, the instability occurs as exponentially growing or decaying modes. If Rr is increased beyond some critical value R'r, these exponential modes change into oscillatory ones. Further, for infinitely increasing Rv, the amplitude and the period of these oscillations increase because of the growing influence of the temperature field. Finally,

71

due to the nonlinearity of the system, these oscillatory modes might become aperiodic. All these results are in qualitative agreement with the results obtained by Huppert and Moore [8] via nonlinear computational analysis of the problem. In fact, this agreement can be understood to some extent because our treatment takes into account the Brownian character of the motion of the fluid molecules and hence indirectly takes into account the weak nonlinearity of the system. REFERENCES

I. G. Veronis, J. Fluid Mech. 34, 315 (1968). 2. F. Zangrando, Ph.D. Thesis, Univ. of New Mexico (1979). 3. H. E. Huppert and P. C. Manins, Deep-Sea Res. 20, 315 (1973). 4. M. S. Sodha, N. D. Kaushik and S. K. Rao, Int. J. Energy Res. 5, 32I (1981). 5. S. Chandrasekhar, 14ydrodynamic and ttydromagnetic Stability. Clarendon Press, Oxford (1961). 6. J. S. Turner, Buoyancy Effects in Fluids. Cambridge Univ. Press (1973). 7. P. M. Morse and H. Feshbach, Methods qf Theoretical Physics, Part I. McGraw-Hill, New York (1953). 8. H. E. Huppert and D. R. Moore, J. Fluid Mech. 78, 821 (1976).