## Abstract

In two-phase flows of steam, when the velocity is between the equilibrium and frozen speeds of sound, the system is fundamentally unstable. Because any disturbance of the system, e.g. imposition of a small supercooling on the fluid, will cause condensation, the resulting heat release will accelerate the flow and increase the supercooling and thus move the system further from thermodynamic equilibrium. But in high-speed flows of a two-phase mixture, dynamic changes affect the thermodynamic equilibrium within the fluid, leading to phase change, and the heat release resulting from condensation disturbs the flow further and can also cause the disturbances to be amplified at other Mach numbers. To investigate the existence of instabilities in such flows, the behaviour of small perturbations of the system has been examined using stability theory. It is found that, although the amplification rate is highest between the equilibrium and frozen speeds of sound, such flows are temporally unstable at all Mach numbers.

## 1. Introduction

Addition of heat to a flowing compressible fluid is essentially an unsteady process as the changes of the flow properties brought about by the addition of heat will propagate into the stream by wave action. Flow of a two-phase fluid, e.g. wet steam, is a special case of the above problem as dynamic changes experienced by the fluid will affect the thermodynamic equilibrium within it. To regain thermodynamic equilibrium, there will be phase change, and the heat release resulting from condensation will affect the flow properties, e.g. pressure and velocity, and the changes will propagate in the medium. Generally, thermodynamic changes lag the dynamic ones and, under some circumstances, the resulting heat release can cause the imposed disturbances to be amplified and the flow to become unstable. The nomenclature is given as table 1.

An extreme example of the problem occurs when the flow velocity is between the equilibrium and frozen speeds of sound. It is well known that, when a perfect gas flows in a frictionless constant area duct, addition of heat to the gas raises its temperature except in the Mach number range when it will cause it to cool. In the case of a two-phase fluid, phase change is accompanied by heat release and the relevant parameter is the supercooling of the flow. This is because addition of heat to the fluid will also affect its pressure and hence the local saturation temperature which has to be considered. It can be shown that addition of heat to a two-phase fluid flowing in a frictionless constant area duct will reduce its supercooling except in the Mach number range when it will cause it to increase and the system is fundamentally unstable. Here, *T*_{s}(*P*) denotes the saturation temperature at pressure *P*. This is because any disturbance of the flow, e.g. imposition of a small supercooling on the fluid, will result in condensation and the resulting heat release will accelerate the flow, increase the supercooling and thus move the system further from equilibrium (Bakhtar & Young 1978). The lower Mach number limit corresponds with the equilibrium speed of sound. This process provides a positive feedback mechanism for disturbances to be amplified and the flow to become unstable. Such oscillations were experienced in the course of investigations on the behaviour of wet steam in the blowdown equipment at Birmingham (Bakhtar *et al.* 1997*a*,*b*). In steam turbines, the flow in the gaps between stators and rotors is generally in this range and any possibility of resonance between the flow and the blading is of considerable practical importance. Steam turbines play a dominant role in the generation of the mains electrical power supply. Although instability associated with the above conditions was the main motivation for the present study, release of heat following a departure from equilibrium can cause instability of wet steam flows at other Mach numbers. The problems are not specific to steam and could arise in any two-phase mixture. Their examination in relation to steam is because of the background of the authors.

Some unsteady conditions are experienced in nucleating flows of steam in convergent–divergent nozzles, when the zone of rapid condensation occurs at a short distance downstream of the throat and the pressure rise associated with heat release is of a magnitude to propagate upstream (e.g. Barschdorff 1971; Yousif *et al.* 1972; Saltanov 1976; Skillings & Jackson 1987; White & Young 1993; Adam & Schnerr 1997; Schnerr 2005), but these lead to steady oscillations. The unstable conditions considered in the present paper differ from these in character and form and occur over a much wider range of conditions.

Many unstable flow conditions have been studied in single-phase fluids in the literature and a comprehensive review of the subject is given by Drazin & Reid (1993). The present problem is more akin to these but differs from them in that the mechanism driving the instability is heat release associated with phase change. To investigate the existence of instabilities, the behaviour of small perturbations of the system will now be examined using stability theory.

## 2. The governing equations

The equations are essentially the standard unsteady gas dynamic conservation equations of mass, momenta and energy as applied to two-dimensional inviscid flow of a compressible fluid with heat addition in which the source of heat is the energy resulting from phase change. To simplify the solution, it is assumed that the vapour behaves as a perfect gas, the liquid is at saturation temperature corresponding to the local pressure, inter-phase slip is negligible and the contribution of the liquid to the density of the mixture is small. These assumptions have been used widely in the literature in many investigations of nucleating and two-phase flows of high-quality low-pressure steam and, with a small sacrifice of accuracy, yield meaningful solutions. With these assumptions, the conservation equations may be written as

Continuity:(2.1a)Momentum in the *x*-direction:(2.1b)Momentum in the *y*-direction:(2.1c)Energy equation:(2.1d)where *e* is the specific internal energy of the fluid; *P* is the pressure; *ρ* is the density; and *u* and *v* are the components of velocity in *x*- and *y*-directions, respectively.

Equation of state:(2.1e)where *R* is the gas constant and *T* is the temperature. The main consequence of phase change is the release of latent heat which directly affects the energy equation. To deal with this effect, *e* is written as(2.2)where *w* is the wetness fraction; *E* is the specific internal energy of phase change; and suffixes *L* and *G* refer to the liquid and vapour phases, respectively. Considering *E* as a constant and writing(2.3)substitution into equation (2.1*d*) after some rearrangement yields(2.4)The last term in equation (2.4) is the rate of condensation which, in wet steam after the nucleation stage, may be written as(2.5)where ; *N* is the number of droplets per unit mass; *r* is the droplet radius; *L* is the enthalpy of the phase change (assumed constant); and *α*_{r} is the heat transfer coefficient. Substitution in the energy equation yields(2.6)The above equations are sufficient to describe two-dimensional flows of two-phase fluids and in their more complete form without the introduction of the simplifying assumptions have been used for the solution of steady flows of nucleating and wet steam by Bakhtar & Mahpeykar (1997), among others. The treatment was originally based on Denton's numerical scheme, but it has subsequently been refined and the algorithm based on a Runge–Kutta scheme (Zamri 1997; Bakhtar *et al.* 2007). In the present study to investigate the stability conditions considered, it has been necessary to use, as the starting point, the steady-state solution of the basic flow configuration. For this purpose, the available treatment has been used, but for consistency the same simplifications have been introduced into the codes (Zamri 1997).

## 3. Linearization of the equations

To investigate the most amplified mode of the disturbances, all the fluid and flow properties in the above equations can be taken as the sum of a mean value and a small perturbation. Substitution into equations (2.1)–(2.6), cancelling terms describing the mean flow properties which satisfy the steady form of the equations and neglecting second-order perturbation terms, yields the linearized form of the governing equations.

Because the droplets are assumed to be at the saturation temperature, perturbations in droplet temperature can be expressed in terms of the perturbations in pressure *P*′ via the Clausius–Clapeyron relation, i.e.(3.1)Further simplifications may be obtained by selecting a flow configuration in which the mean flow is approximately parallel to, and has small gradients in, the *x-*direction. Consequently, the *y* component of the velocity *v* can be regarded as small in comparison with *u* and is neglected. Furthermore, the axial gradient of all mean variables can be neglected. With these assumptions, the linearized equations reduce to(3.2a)(3.2b)(3.2c)and(3.2d)A bar over a symbol denotes its mean value and a dash indicates its fluctuating component.

### (a) Stability equations

The stability of the flows described in §3 can now be considered. It is noted that the flow quantities are functions of the axial position *x* and the cross-channel coordinate *y.* In order to make progress, it is assumed that the instabilities develop over very short scales and thus it is pertinent to study the local stability of the flow. The flow properties are taken to be of the formwhere the superscript T denotes the transpose. Intrinsic in the assumption that the flow is parallel is the postulate that flows in the cross-channel direction can be neglected. Now infinitesimal perturbations to this system are considered and these are assumed to take the formIn this expression, *α* is the stream-wise wavenumber and *c* is the phase speed. The real part of the product *αc* gives the frequency of the modes, *Ω*. Since the general case is being considered, the analysis for both the cases of temporal and spatial stabilities will be presented. Perturbation of the two-phase Euler equations (3.2*a*)–(3.2*d*) yields the following equations governing the normal modes:(3.3a)(3.3b)(3.3c)and(3.3d)These equations can be manipulated to obtain a single second-order differential equation, namely(3.4)where(3.5a)(3.5b)and(3.5c)Thus provided the steady flow variables across the channel are available, equation (3.4) can be solved using a shooting technique. This is preferred to the alternative of using an analytical solution which would have to employ hypergeometric functions and would still need to be evaluated numerically. In fact, there would be no additional insight available from the expression of the solutions in this form and the direct numerical solution is less complex. For the present analysis, equation (3.4) may be written in a finite-difference form, i.e.(3.6)and(3.7)where *n* is the station identifier in the transverse direction. Substitution into equation (3.4) yields(3.8)where(3.9a)(3.9b)and(3.9c)When the basic mean flow variables are specified, equation (3.4) contains two variables, namely *α* and *c*. If the temporal growth is sought, then *α* is assumed to be real and known. With *Y* the vertical extreme of the channel, using the normalization condition at *y*=0 and the boundary conditions of zero flow into the solid wallsEquation (3.4) can be solved for *c*. If spatial growth is sought, *ω* is assumed to be real and known and the equations furnish *α*.

When solving for temporal growth, for a known value of *α*, a value for the complex variable *c* is guessed. Thus, at each point along the profile the values of the terms on the r.h.s. of equation (3.9*a*)–(3.9*c*) and hence values of *A*_{n}, *B*_{n} and *C*_{n} are known. Using normalization at *y*=0 and at *y*=0 and *y*=*Y*, the equations can be solved for values of for all *y*. Once has been determined, the disparity of at *y*=0 from zero can be used to iterate on the value of *c*, using a secant method. The calculations are repeated for a range of values of *α*. The procedure is similar when solving for spatial growth, in which the case *ω* is assumed to be known and the equations are solved for *α*.

## 4. Application to wet steam flows

To investigate the existence of instabilities in wet steam flows, the above treatment was applied to a number of specific wet steam flow conditions. These were generated theoretically by applying an existing two-dimensional two-phase treatment to the flow in two convergent–divergent nozzles. Treatments for such two-phase flows have been under study by the group at Birmingham for many years and the quality of agreement obtained with a wide range of experimental measurements has been very good. More recently, the numerical scheme has been refined further and the algorithm modified to use a Runge–Kutta scheme (Bakhtar *et al.* 2007). Using the treatment, steady-state solutions obtained have been used as the starting point for the present analysis. But, as already indicated to cater for consistency with the instability calculations, the same simplifying assumptions have been introduced into the two-dimensional code. Wet steam flow conditions that have been analysed are taken from specific axial positions in the two convergent–divergent nozzles. The shape of the converging section was the same in both nozzles and in both sets of solutions; the inlet stagnation conditions were set at 1.008 bar, 374 K, wetness fraction of 0.5% and droplet radius of 0.05 μm. In both nozzles, the flow was choked close to the equilibrium speed of sound at the throat. To ensure that the gradients of the flow properties in the direction of flow were small, the angle of divergence of the first nozzle was 0.1°. Thus, the flow velocity increased slowly from the equilibrium to the frozen speed downstream of the throat and the outlet Mach number was slightly above unity.

The droplet radius adopted is a little smaller than the measured droplet sizes in the operating turbines reported. This size was selected to avoid the possibility of secondary nucleation experienced in the equipment with larger droplet diameters (Bakhtar *et al.* 1999). On the other hand, the wetness fraction adopted is lower than that experienced in the wet stages of steam turbines. As the dominating factor is the surface area presented by the droplets, the combination of droplet size and wetness fraction adopted yields droplet surfaces which are larger than those experienced in operating turbines by a small factor. The effect of this disparity on the extremely high calculated amplification rates is expected to be small.

To investigate the stability of the flow at higher Mach numbers and to verify that, at the range of conditions considered, the effect on flow stability of the gradients of flow properties in the direction of the flow were small, the second nozzle had an angle of divergence of 4° and exit Mach number of approximately 1.44.

Temporal stability calculations were carried out for the data from both nozzles. At the same Mach numbers, the stability characteristics of the flow conditions from both nozzles were identical. The angle of divergence of the second nozzle was considerably larger than that of the first nozzle. Consequently, the gradients of the flow properties in both the flow direction and across the channel were higher for the second nozzle. The identity of the results demonstrates that the effect of the gradients of the flow properties in both directions on the stability of the flow is small and their neglect in the analysis is justified. Consequently, the spatial stability calculations were carried out for the data from the first nozzle only.

### (a) Temporal stability of the flow

#### (i) Flow in first nozzle

The variations of Mach number along the centre line of the first nozzle are shown in figure 1. Six axial stations, indicated in figure 1, were chosen for the stability analysis. The midstream Mach numbers at each station are given in the figure caption.

The results of the temporal growth calculations for the first nozzle are plotted as amplification factor and disturbance frequencies as functions of wavelength in figure 2*a*,*b*, respectively, and the variations of amplification rate with frequency are plotted in figure 2*c*. The numbers on the graphs refer to positions in the nozzles and the corresponding Mach numbers are given in the figure caption.

It will be seen from figure 2*a* that, at any specific Mach number, there exists a critical wavelength below which the flow is stable, but, once the critical wavelength is exceeded, the flow amplifies disturbances very highly at all Mach numbers. The amplification rate increases with increases in Mach number in the subsonic range. It is highest in the Mach number range between the equilibrium and frozen speeds of sound but it then decreases with increases in Mach number in the supersonic range; even where the Mach number is 0.424 the amplification rate is 39 000 s^{−1}.

The variation of frequency of oscillations with wavelength is shown in figure 2*b*. It will be noted that, for each wavelength, there is only one frequency of oscillations which causes the flow to become unstable. Conversely, for each frequency of disturbances there is only one wavelength which leads to instability.

There is need to investigate these conditions experimentally, but disturbances present in the flow would be spread over a range of frequencies; however, for a component wave of a disturbance to be amplified, its wavelength should match the value given by the solutions.

For the given dimensions of individual machines, the oscillations may develop characteristic frequencies. Should this be the case, the blade natural frequencies can be designed to avoid them.

It is also seen that as would be expected the frequency decreases with increases in wavelength. For each wavelength, the frequency of the oscillations increases with increases in the flow Mach number.

The variations of the amplification rate with frequency are shown in figure 2*c*. The lower frequency waves have a higher amplification rate. At each Mach number there exists a frequency above which the flow is stable. These limiting frequencies correspond with wavelengths below which the flow is stable.

#### (ii) Flow in second nozzle

The midstream Mach number distribution for this nozzle is shown in figure 3. Four stations, indicated on the figure, were chosen for stability analysis. The mid-passage Mach numbers at these stations are given in the figure caption. The flow is subsonic at station 1, within the range of equilibrium and frozen speeds at stations 2 and 3 and at a high supersonic Mach number at station 4.

Only temporal growth calculations were performed for this nozzle. The results are shown in figure 4. The amplification factor is plotted against wavelength in figure 4*a*. The variations of the frequency of oscillations with wavelength are shown in figure 4*b*, while variations of amplification rate with frequency are given in figure 4*c*. The results are very similar to those obtained for the first nozzle. At station 4, which is at a higher supersonic Mach number, the temporal amplification rate is lower than those at stations 2 and 3. These show that a two-phase flow is unstable at all Mach numbers but instability is greatest between the equilibrium and frozen speeds of sound. The stability characteristic is independent of the shape of the nozzle. The determining factor of flow stability is the mean flow properties, particularly the Mach number.

#### (iii) Marginal stability curve

The variations of the critical wavelength with Mach number from both nozzles are plotted together, as a graph of marginal stability in figure 5*a*. It will be noted that both sets of points fall on the same curve. For wavelengths below the curve the flow is stable. But, even at low Mach numbers, the critical wavelength is relatively short in comparison with relevant dimensions of full-scale steam turbines. Thus, if wet steam flows at a relatively constant Mach number for any distance then it will be unstable and will amplify disturbances imposed on it. The maximum amplification rate is plotted against Mach number in figure 5*b*.

It will also be noted that amplification rates indicated by the present results are extremely high. The character of the disturbances will become nonlinear very quickly.

### (b) Spatial stability of the flow

As already discussed, spatial growth calculations were carried out for the first nozzle only and the results are shown in figure 6. The spatial amplification rate is plotted against frequency in figure 6*a*. For conditions at stations 1 and 2, where the flow is subsonic, the flow remains stable at all angular speeds. At other Mach numbers the flow is unstable. The behaviour at conditions of stations 5 and 6 was very similar, so the graph labelled 6 represents both sets of solutions. These amplification rates are lower than those associated with the temporal growth rates but, as plotted, they are measured per unit length. To measure the rates per unit time they should be multiplied by the velocity, but they are lower than those of the temporal growth rates. Nevertheless they are substantial. It is also noted that surprisingly the amplification rates are higher for the supersonic flows than for those between the equilibrium and frozen speeds. This is attributed to the fact that in the two supersonic flows considered the Mach number is only slightly above unity and very sensitive to disturbances which could make them flip from supersonic to subsonic and back.

The variations of the frequency of oscillations with the wavelength are shown in figure 6*b*. It will be noted that these frequencies are lower than those associated with the temporal growth of the disturbances and for flows in the Mach number range between the equilibrium and frozen speed of sound; the rise in frequency with increasing wavelength is very sharp. Thus, it would appear that these occur at a different frequency from those of the temporally unstable oscillations. The reflection of pressure waves in a turbine is complicated owing to the complexity of the geometry, but if any of the wave paths coincide with multiples of the wavelengths there will be potential for resonance which would be undesirable.

#### (i) Dependence on variation of properties across channel

In the previous analysis, it was assumed that the stability characteristics were only weakly dependent on the variations of the mean fluid and flow properties across the channel. The validity of this assumption is examined further in appendix A. It will be seen from equation (3.4) that the cross-channel variations of properties are represented by the term . In the absence of these variations, this term can be put to zero and the resulting equation can be solved analytically. A comparison between the resulting solutions using this approach and the numerical procedure for a typical case is given in appendix A. It will be seen that the effect of neglecting the cross-channel variations on the results is small.

## 5. Stability of dry flows

Returning to the main flow equations, the influence of two-phase effects on the flow is through the term *ξ* on the energy balance in equation (2.6). In dry flows there will be no condensation and *ξ*=0. If this is substituted in the stability equation, then neglecting the cross-channel variations it will reduce to(5.1)which gives(5.2a)and(5.2b)Hence, it follows that regardless of the flow conditions a dry flow will always be stable. This is not surprising; for example, in parallel shear flows, Raleigh's criterion requires that, for instability to occur, the flow must contain at least one inflection point in the distribution of variables in the *y*-direction. Thus, in the case of dry flows the profile of the variables determines the stability of the flow. In contrast, in wet steam flow, this component is swamped by the effect of *ξ*, which is of the order of 10^{7}.

## 6. Discussion and conclusion

It has been demonstrated that, in the flow of wet steam, in the Mach number range between the equilibrium and frozen speeds of sound, due to the effect of condensation on supercooling, there exists a positive feedback mechanism which will amplify any disturbance within the flow and cause it to become unstable. The results also show, interestingly, that flowing wet steam is temporally unstable at all Mach numbers. Instability is greatest between the equilibrium and frozen speeds of sound but the flow is unstable at all other Mach numbers. This is because, when the frequency is high, the rate of internal phase change within the system will not be sufficient to maintain thermodynamic equilibrium and a disturbance will move the system from it. The heat release associated with the consequential phase change will amplify the disturbance. Should this type of instability occur in real turbines, it would have serious repercussions on blade performance and durability. It has been reported that the majority of blade failures in steam turbines occur in the low-pressure stages where the steam is wet (Dewey *et al.* 1983). A major proportion of the flow in the stator–rotor gaps is between the equilibrium and frozen speeds of sound. In the rotor–stator gaps the flow is in moderate subsonic range and still potentially highly unstable. Furthermore, the flow in the throat regions of control valves is also in this range.

The solutions predict the existence of a critical wavelength, below which the flow is stable at all Mach numbers. In the Mach number range between the equilibrium and frozen speeds of sound, the predicted critical wavelength is approximately 47 mm. The critical wavelengths are a little longer at Mach numbers both above and below this range, but these are still smaller than typical dimensions of turbine stator–rotor gaps. For unstable flows, the frequency of the oscillations decreases with increases of the wavelength. The frequency is highest for the shortest wavelengths. It is probable that the frequencies associated with the shorter wavelengths are too high to have serious practical consequences, but very high amplification rates are also predicted for the relatively lower frequencies which could lead to problems. Oscillations, once developed, can propagate in the rest of the flow field. For given machine dimensions, the oscillations may develop characteristic frequencies. If so, they can be avoided in the design of the natural frequencies of the blades.

The dissipation associated with the internal heat transfer will enhance the damping due to viscous shear and heat conduction within the flow. These effects have been neglected in the present treatment. The inclusion of the additional terms in the equations will increase the mathematical complexities considerably. The predicted frequencies will not be affected by the neglect of the viscous and other dissipation terms but the amplification factors will be lower than predicted. There is need to investigate these conditions experimentally.

The present results also have implications for the stability of boundary layers. Stability characteristics of boundary layers are based on the spatial amplification of disturbances. The present spatial stability calculations predict wet steam flows to be stable at velocities below the equilibrium speed of sound. But these flows are shown to be temporally very unstable at all Mach numbers. It may be argued that the instability mechanism associated with wet vapour will also be the dominating factor for boundary layers and cause earlier transition to a turbulent regime with consequences for the associated profile losses. On the other hand, flows over turbine blading are generally accelerating and the distances over which the Mach number may be regarded as constant are very short and consequently the flow will be stable. The problem warrants a separate experimental investigation.

## Acknowledgments

The investigations reported in this paper were carried out in the Schools of Mechanical and Manufacturing Engineering and of Mathematics and Statistics of the University of Birmingham. M.Y.Z., one of the authors, is grateful for financial support from Tenaga Nasional Berhad which enabled him to take part in the investigation.

## One-dimensional stability properties

If the cross-channel variation of the basic quantities can be neglected, equation (3.4) becomes(A1)This equation can be solved very simply to obtain(A2)where the boundary condition that at *y*=0 has been exploited and *Ψ* is given by(A3)

Using the boundary condition at *y*=*Y*, it is found that(A4)where . Now, it can be inferred that and for *n*∈*Z*.

Thus,(A5)The dependency of the solutions on cross-channel variations can be examined by comparing the values of the eigenfunction , obtained numerically by solving equation (3.4) and using equation (A 1). The eigenfunction was found to be almost identical at all stations. Selecting, as an example, the conditions at station 1 of the first nozzle, which had the largest variation in properties across the channel, the comparison between the real and imaginary parts of the resulting values of , using the two methods of solution, are given in figure 7*a*,*b*, respectively. In this the value of is normalized to unity. It can be seen that the distributions of the real parts are almost identical. Equation (A 2) contains only the real component. Some discrepancy is seen in the distribution of the imaginary part. However, the absolute value of the difference is very small, the maximum being 4×10^{−4}, in comparison with the value of the real part. It can thus be concluded that the stability characteristics of wet steam flow equations is virtually independent of the cross-channel variations of the properties.

## Footnotes

- Received June 18, 2007.
- Accepted November 14, 2007.

- © 2007 The Royal Society

## References

## Notice of correction

Equation (2.2), the penultimate paragraph on page 546 and figure 5*a* are now present in their correct forms.

A detailed erratum will appear at the end of volume 464.

31 January 2008