4.1 Estimation of the phase slips in the wave of compression

Let us show that, in compression waves, the phase slip can reach appreciable values. The characteristic time for a wavy processes, $\unicode[STIX]{x1D70F}^{\ast }$ , can be evaluated as

where $L^{\ast }\sim 1~\text{m}$ is the characteristic length of the compression wave; $a_{l}^{\ast }\sim 1000~\text{m}~\text{s}^{-1}$ is the speed of sound in liquid; $u_{l}^{\ast }$ is the liquid velocity behind the wave front.

Under condition (4.1), the convective components of substantial derivatives with respect to time in the equation for single bubble motion (1.9) can be neglected. Neglecting vapour density compared with that of the liquid $(\unicode[STIX]{x1D70C}_{g}\ll \unicode[STIX]{x1D70C}_{l})$ , too, we obtain the equation for the phase slip, $\unicode[STIX]{x0394}u=u_{g}-u_{l}$ , in the form

In compression waves, the two first terms on the right-hand side of (4.2) lead to an increase in the absolute value of the difference in the phase velocities. The first term, which is the rate of the velocity decreasing in the compression wave, results in an increase of the absolute value of the negative phase slip. The second term represents how, in a compression wave, bubbles collapse ( $\unicode[STIX]{x2202}a/\unicode[STIX]{x2202}t<0$ ), lose their virtual mass and the absolute value of the difference in phase velocities increases. The last term, the interfacial drag force, always diminishes the phase slip. The difference in the phase velocities $\unicode[STIX]{x0394}u$ will stop increasing as soon as the last term reaches a value equal to the sum of absolute values of the first two terms.

When the difference in phase velocities has the same order of magnitude as the velocity,

all three terms on the right-hand side of (4.2) are of the same order:

Estimation (4.4) is obtained taking into account supposition (4.3) and evaluations (4.1) for $\unicode[STIX]{x1D70F}^{\ast }$ , $u_{l}^{\ast }\approx u^{\ast }$ . The bubble radius is estimated to be $a\sim 10^{-4}~\text{m}$ from experimental liquid overheats and $\unicode[STIX]{x1D70F}^{\ast }$ . The bubble Reynolds number is estimated as $Re\sim 10^{4}$ from the bubble radius and $\unicode[STIX]{x0394}u$ ; $C_{\unicode[STIX]{x1D707}}(Re)\sim 4\times 10^{-2}$ is taken from Schlichting (Reference Schlichting1964).

Therefore, the interfacial drag forces are rather small, and the difference in the phase velocities in the compression wave can reach values comparable with the flow velocity. The situation $|\unicode[STIX]{x0394}u|\sim u$ did take place in calculations with the model of Ivashnyov et al. (Reference Ivashnyov, Ivashneva and Smirnov2000) shown in figure 5(c). Thus, as it is based on the assumption that $|u_{g}-u|\ll u$ , the model goes beyond its validity limits. In order to calculate compression waves in bubbly media one needs to abandon the assumption.

Unfortunately, the known models for two-velocity two-phase flow with a compressible carrying phase do not permit one to simulate correctly the flow with phase slip comparable to the flow velocity. These models are named by different authors as ‘non-hyperbolic’ (Nigmatulin Reference Nigmatulin1991), ‘non-evolutionary’ (Kulikovskii & Shikina Reference Kulikovskii and Shikina1996), ‘incorrect’ (Crayko Reference Crayko1982), ‘non-stable’ (Gear Reference Gear1971; Lahey et al. Reference Lahey, Cheng, Drew and Flaherty1980), etc.

4.2 The requirements for the model of a compressible medium in mathematical terms

For any compressible medium, the velocity of the small perturbations spreading, the sound velocity, must be finite. In mathematical terms, this finiteness means that the model is hyperbolic. A hyperbolic model is one that has only real mathematical characteristics. In essence, the characteristics are the lines in $(t,x)$ coordinates along which small high-frequency perturbations are spread. Any hyperbolic system can be rewritten in an equivalent characteristic form, i.e. as a system of the equations for characteristic lines and invariants preserved along them.

For the simplest model of ideal gas adiabatic flow,

the characteristic form is the system of Riemann invariants,

where $\unicode[STIX]{x1D70C}$ , $P$ and $u$ are the density, pressure and velocity of an ideal gas; $\unicode[STIX]{x1D6FE}$ is the adiabatic curve exponent; $a_{g}$ is the speed of sound in gas.

It follows from (4.6) that, along the first characteristic, the line $\text{d}x/\text{d}t=u+a_{g}$ , the first Riemann invariant is preserved, the complex of parameters $2a_{g}/(\unicode[STIX]{x1D6FE}-1)+u=\text{const.}$ Alternatively,

The second invariant, $2a_{g}/(\unicode[STIX]{x1D6FE}-1)-u=\text{const.}$ , is preserved along the second characteristic $\text{d}x/\text{d}t=u-a_{g}$ ,

4.3 Characteristics of a ‘classic’ model for two-phase two-velocity flow

Following Lahey et al. (Reference Lahey, Cheng, Drew and Flaherty1980), let us consider the generally adopted ‘dirty water’ model of a fluid with incompressible particles. Let the fluid be compressible with a linear equation of state. The model consists of two equations for mass and two ‘symmetrical’ equations for momentum:

where $a$ is the particle radius; $\unicode[STIX]{x1D70C}_{l}$ and $\unicode[STIX]{x1D70C}_{g}$ are the densities of the fluid and particles; $u_{l}$ and $u_{g}$ are the phase velocities; $\unicode[STIX]{x1D6FD}$ is the volumetric vapour content; $P$ is the pressure in the fluid; $M_{l}$ and $M_{g}$ are the interfacial area forces for the liquid and particles; $F_{w}$ is the wall shear force per unit mixture volume; $n$ is the number of particles per unit mixture volume; $S$ is the drag force per unit velocity given by (1.10); $f_{vm}$ is the virtual mass force; $a_{vm}$ is the virtual mass acceleration; $\unicode[STIX]{x1D6FD}_{l}$ is the coefficient of the liquid compressibility; $\text{D}_{l}/\text{D}t$ and $\text{D}_{g}/\text{D}t$ are the material derivatives following the liquid and the particles.

A standard procedure of characteristic searching gives an algebraic equation of the fourth order:

Here, $\text{d}x/\text{d}t$ is the angle of the characteristic inclination in the plane $(t,x)$ .

For a bubbly flow $\unicode[STIX]{x1D70C}_{g}\ll \unicode[STIX]{x1D70C}_{l}$ with a ‘small’ difference in the phase velocities $|\unicode[STIX]{x0394}u|\ll a_{l}$ , where $a_{l}=\unicode[STIX]{x1D6FD}_{l}/\unicode[STIX]{x1D70C}_{l}$ is the speed of sound in the carrying fluid, we search for the roots of (4.16) by expending in a power series in a small parameter $\unicode[STIX]{x1D700}=\unicode[STIX]{x0394}u/a_{l}\ll 1$ . The mathematical characteristics have the form within accuracy to the first non-zero terms

Thus, the system (4.10), (4.11) has two real and two complex characteristics and is attributed to a hybrid hyperbolic–elliptical type.

An analysis of two-speed models of a compressible fluid with disperse particles analogous to (4.10)–(4.15) given by Ramshaw & Trapp (Reference Ramshaw and Trapp1978), Crayko (Reference Crayko1982) and Nigmatulin (Reference Nigmatulin1991) also showed the presence of two complex mathematical characteristics.

4.4 The link of the solution stability and the model type

Let us consider the evolution of a harmonic perturbation in a medium described by the simplest elliptical equation:

Equation (4.20) has two complex characteristics:

Solution of (4.20) is usually sought inside the field at which the boundary function $\unicode[STIX]{x1D711}$ or its gradient is specified. However, in our case, time is one of the variables, and determination of the function $\unicode[STIX]{x1D711}$ at the boundary means, from a physical point of view, specification of the result of the solution in advance, i.e. determination of the state that the system will come to before solving the problem. From this follows the term ‘incorrectness’ of the model (Crayko Reference Crayko1982).

Let us pose a typical problem about the evolution of a small harmonic perturbation with wavelength $\unicode[STIX]{x1D706}$ in an infinite area $-\infty <x<\infty$ , $t>0$ . Since (4.20) is an equation of second order, besides the function $\unicode[STIX]{x1D719}_{0}(x)$ , it is necessary to specify its derivative with respect to time at $t=0$ . Let the initial conditions be

The solution of (4.20) that satisfies the initial conditions (4.22) is given by

It follows from solution (4.23) that the initial perturbation separates into two parts: the first one is damping and the second one is increasing exponentially. In total, the perturbation infinitely increases. From this, there arise the other names for the deficiency of the system, e.g. its ‘non-evolutionary character’ (Kulikovskii & Shikina Reference Kulikovskii and Shikina1996).

For a medium described by the simplest hyperbolic equation,

the solution of the same problem about the evolution of the perturbation with the same initial conditions (4.11) describes the propagation of a perturbation having a constant amplitude,

i.e. the solution is stable.

For the model of bubbly flow, a propagation analysis for small-amplitude waves was carried out in Ramshaw & Trapp (Reference Ramshaw and Trapp1978). The system was rewritten in a characteristic form and it has been shown that complex characteristics lead to instable solutions.

4.5 Known approaches to improvement of the non-hyperbolic model

Stability of non-hyperbolic model solutions is usually attained by introducing a pseudo-viscosity (Lahey et al. Reference Lahey, Cheng, Drew and Flaherty1980). The coefficients of pseudo-viscosity are specified to ensure computational stability and are limited by nothing. From the theory of the boundary layer in a viscous fluid, it is known that even small coefficients at high-order derivatives can change the solution essentially (Schlichting Reference Schlichting1964). Therefore, the numerical solution does not correspond to the model equations: any solution of a set of equations with complex characteristics is always instable, and the very solutions are always stable according to numerical calculations.

Developing the idea of Gear (Reference Gear1971) that an appropriate modelling of a virtual mass can improve the model stability, Lahey et al. (Reference Lahey, Cheng, Drew and Flaherty1980) made an attempt to eliminate the defect by choosing an arbitrary parameter $\unicode[STIX]{x1D706}$ in the general formula for the virtual mass acceleration proposed by Drew, Lahey & Cheng (Reference Drew, Cheng and Lahey1979):

However, it occurred that no coefficients permitted them to get rid of the instability in the field of low volumetric content of a dispersed phase $\unicode[STIX]{x1D6FD}<0.3$ , i.e. in the field of a bubbly flow.

Considering that the model of two-phase flow did not take into account particle interactions, Nigmatulin (Reference Nigmatulin1984), as earlier Wallis, Richter & Kuo (Reference Wallis, Richter and Kuo1976), suggested addition of the momentum equation for a disperse medium with the gradient of the ‘pressure of particles’ proportional to their volumetric content,  $\unicode[STIX]{x1D6FD}$ :

In the work cited, calculations were not presented. However, our tests have shown that even when the model is specified so that the pressure of the particle interactions is equal to the pressure in the liquid, it is not successful in getting rid of the defect for the volumetric contents of the dispersed phase $\unicode[STIX]{x1D6FD}<0.1$ (for a bubbly flow).

4.6 A correct hyperbolic model for two-velocity flow

The forces of interphase interaction are on the right-hand sides of the momentum equation (4.11). Their expressions contain derivatives with respect to time. To search characteristics, these derivatives with respect to time should be transported to the left-hand sides of the equations. Thus, the terms characterizing the interphase interaction are the only ones that can influence the characteristics and change the type of the system (Bourè Reference Bourè1975; Ramshaw & Trapp Reference Ramshaw and Trapp1978; Lahey et al. Reference Lahey, Cheng, Drew and Flaherty1980).

In order to find the reason for the defect of the two-phase two-velocity flow model (4.10)–(4.15), let us ‘reconstruct’ its derivation. At first, we limit ourselves to the consideration of the mixture of an incompressible fluid ( $\unicode[STIX]{x1D70C}_{l}=\text{const.}$ ) with dispersed particles of a constant radius, $a=\text{const.}$

The scheme of the classic statement of the problem of a single particle overflowing by a forward flow is shown in figure 7(a). The force acting on a representative particle in an inviscid incompressible flow is (Lamb Reference Lamb1957)

Here, $u_{l\infty }(t)$ is the velocity of the liquid at distance; $u_{g}(t)$ is the velocity of the test particle; $\text{d}/\text{d}t$ are ordinary derivatives with respect to time. The terms on the right-hand side are usually grouped into expressions for the forces of buoyancy and virtual mass, $f_{b}$ , $f_{vm}$ .

In a two-phase medium, the particle velocities vary with space and time: $u_{g}(t,x)$ . To apply formula (4.28) obtained for a single particle for a particle in ensemble, the derivatives are to be taken along the trajectory of motion of just this particle:

It is shown in appendix B that the liquid velocity at distance $u_{l\infty }(t)$ is equal to the average velocity of the liquid in the space surrounding this particle, $u_{l}$ . Therefore, the derivatives of the liquid velocity could be taken along the trajectory of particles as well:

Thus, the force acting on a particle in a two-phase medium, $f_{bl}$ , takes the form

To generalize this expression for the case of particle movement in a viscid fluid, the Stokes force is added to the sum of forces acting on the particle (on the right-hand side of (4.31)). The expression for this force does not contain derivatives and does not affect the characteristics of the model.

where $S$ is the drag force per unit velocity given by (1.10).

Figure 7.

The statement of the problem on a particle overflowing by a forward (a) or gradient (b) flow.

From the mass conservation equation for the liquid phase one from (4.10a ) for the flow of an incompressible fluid ( $\unicode[STIX]{x1D70C}_{l}=\text{const.}$ ),

it follows that the divergence of the average velocity arises due to change in the particle distribution in time and space. On the other hand, the formula for the force of the interphase interaction (4.33) is obtained under the assumption that the divergence is zero since the field of the velocity is symmetrical relative to the particle midsection according to the problem statement (figure 7 a). Therefore, closing expression (4.33) is only valid for the model of a rarefied mixture where the particles do not effect the fluid and formula (4.34) is reduced to

For the model with the mass conservation equation (4.34), it needs the expression for the interphase force that considers the presence of neighbouring particles. To deduce such an expression, one must state the problem so that the gradient of the average micro-flow velocity around a representative particle is equal to $\unicode[STIX]{x2202}u_{l}/\unicode[STIX]{x2202}x=0$ in macro-flow.

The scheme of the problem statement of a particle overflowing by a gradient flow is shown in figure 7(b). As a flow of an incompressible fluid can be gradient only due to a change in the particle distribution, this statement considers the interaction of the neighbouring particles (with the same degree of accuracy as the mass equation (4.34) does). Voyinov, Voyinov & Petrov (Reference Voyinov, Voyinov and Petrov1973) solved the problem on a sphere overflowing by a gradient flow. (Since the equations of a multiphase flow model were obtained earlier than the solution of Voyinov et al. (Reference Voyinov, Voyinov and Petrov1973), an analogous solution had been obviously obtained by the authors of the model independently.) In appendix C, we give the solution for a particle of variable radius $a(t)$ (which will be required later on). If the radius is constant ( $a=\text{const.}$ ), the expression for the interfacial force has the form

where the expressions for the buoyancy force $f_{b}$ and virtual mass force $f_{vm}$ are given by (4.32). Comparing (4.33) and (4.36), we see that the presence of neighbouring particles results in the appearance of an additional term in the expression for the interphase force. This term can be named as the force of bubble interaction, $f_{bb}$ ,

The field of stresses in the ‘gas of disperse particles’, $\unicode[STIX]{x1D6F1}(t,x)$ , can be introduced as

where $\unicode[STIX]{x1D6F7}(t)$ is an indefinite time function.

The expression for the force of interphase interaction (4.36) accounting for the mutual influence of dispersed particles can be rewritten in the form

where $f_{b}$ and $f_{vm}$ (4.32) have been obtained from the solution of the problem of a particle overflowing by a forward flow (figure 7 a).

In the form (4.39), the force $f_{bb}$ has been ‘divided’ into two components. One forms a part of $f_{b}^{\prime }$ and the other is a part of $f_{vm}^{\prime }$ . It should be noted that there is an insignificant, at first sight, distinction between the forces $f_{b}$ and $f_{b}^{\prime }$ , and between $f_{vm}$ and $f_{vm}^{\prime }$ as well: all of the derivatives are taken along the trajectory of particle movement in the first ones; in the second ones the derivatives of the liquid velocity are taken along the trajectory of the liquid phase movement.

On substituting the expression (4.39) for the interphase interaction force $f_{bl}^{\prime }$ into the one-dimensional version of the conservation equations for mass and momentum averaged throughout the volume (Ishii & Hibiki Reference Ishii and Hibiki2006), we obtain

where the expressions for $f_{b}^{\prime }$ and $f_{vm}^{\prime }$ are given by formulae (4.41); $n$ is the number of particles per unit mixture volume.

After substituting expressions (4.44), (4.41) for the interphase interaction force $f_{bl}^{\prime }$ into the right-hand sides of the momentum equation (4.43), transferring the buoyancy force $f_{b}^{\prime }$ into the left-hand side, we are coming to a ‘classic’ system of equations for a bubbly flow with a ‘symmetrical’ momentum equation (4.10)–(4.15).

Thus, the ‘classic’ model (4.10)–(4.15) considers a mutual influence of dispersed particles (4.18), but for the particles located in an incompressible liquid, $\unicode[STIX]{x1D70C}_{l}=\text{const.}$ The model (4.5)–(4.7) is correct for the case of an incompressible carrying phase but is not correct on being generalized to the case of a compressible liquid since the system of equations for a compressible medium is closed by the ratio obtained for an incompressible liquid which supposes the velocity of perturbations spread between particles to be instantaneous. Therein lies the self-contradiction of the model, its defect.

Let us correct the system by replacing the expression for the force of the particle interaction in an incompressible medium by the force in a compressible one, $f_{bb}^{\prime \prime }$ . The expression for $f_{bb}^{\prime \prime }$ is taken in a relaxation approximation,

Here, $\unicode[STIX]{x1D6F1}_{x}=\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}/\unicode[STIX]{x2202}x$ is the pressure gradient in the dispersed particle medium; $\unicode[STIX]{x1D70F}^{\ast }$ is the time within which the gradient flow is settled around the test particle, i.e. the time necessary for the wave to run the distance between the particles; $a/\unicode[STIX]{x1D6FD}^{1/3}$ is the distance between the bubbles; $a_{l}$ is the speed of sound in liquid.

On replacing $f_{bl}^{\prime }$ by $f_{bl}^{\prime \prime }$ in the momentum equations (4.43), we obtain the following system of equations:

The model (4.47)–(4.51) has five real mathematical characteristics in the whole range of volumetric contents. For the first couple, exact analytical expressions are found:

and the other three characteristics are determined from the solution of the cubic equation:

where the parameters $y$ and $k$ are determined by (4.17).

Searching for the roots of (4.53) in the form of the expansion in a small parameter $\unicode[STIX]{x1D700}=\unicode[STIX]{x0394}u/a_{l}\ll 1$ , we obtain

Thus, the model corrected for a compressible fluid with particles of constant radius (4.47)–(4.51) is hyperbolic in the whole range of parameters.

In spite of the fact that the time of installation of a stationary gradient around a disperse particle $\unicode[STIX]{x1D70F}^{\ast }$ is negligibly small and much less compared with the characteristic times of any problem, these small corrections inserted into the classic model can change the solutions qualitatively. Critical regimes can be realized in stationary flows under a third characteristic becoming equal to zero. Moreover, waves spreading with the velocity of the third characteristic can be formed during transient processes. Let us give some examples of such flows.

4.7 Simulation with a hyperbolic model of a two-velocity steady-state flow in a diffuser

Let us consider a flow in which calculations with a ‘classic’ model were followed by ‘numerical instabilities’ in solutions (Lahey et al. Reference Lahey, Cheng, Drew and Flaherty1980). For a channel with a variable cross-section, the stationary version of the improved model (4.47)–(4.51) taking into account gravity and neglecting the drag force between the channel walls and flow is

where $F$ is the cross-sectional channel area; $\unicode[STIX]{x1D714}$ is the volume-average velocity of the mixture; $g$ is the acceleration due to gravity; $a_{l}$ is the speed of sound in the liquid; $S$ is the drag force per unit velocity given by (1.10); $C_{\unicode[STIX]{x1D707}}$ is the drag coefficient between a solid particle and fluid taken from Lahey et al. (Reference Lahey, Cheng, Drew and Flaherty1980). The liquid is compressible with the linear equation of state (4.13).

Equations (4.57)–(4.59) in a form suitable for numerical integrating are

where $\unicode[STIX]{x0394}u=u_{g}-u_{l}$ is the difference in phase velocities.

Let lightweight particles, $\unicode[STIX]{x1D70C}_{g}=1~\text{kg}~\text{m}^{-3}$ , of a constant radius $a=0.15~\text{cm}$ be injected with velocity $u_{g0}=0.048~\text{m}~\text{s}^{-1}$ into the narrow bottom end of the diffuser (figure 8 a). The volumetric rate of the particles is $Q_{g}=\unicode[STIX]{x1D6FD}_{0}u_{g0}F_{0}=0.547\times 10^{-6}~\text{m}^{3}~\text{s}^{-1}$ ( $\unicode[STIX]{x1D6FD}_{0}=0.01$ ). Water is flowing in the opposite direction (downward). The pressure at a bottom nozzle section is equal to atmospheric pressure.

Figure 8(b–d) shows the solutions obtained under different values of the fluid velocity at the bottom cross-section, $u_{l0}$ . Among the solutions obtained, there are continuous, $\text{OA}_{1}$ , $\text{OA}_{2}$ , $\text{OA}_{3}$ , and broken-off, OSB, $\text{OC}_{1}$ , $\text{OC}_{2}$ , $\text{OC}_{3}$ , ones. In continuous regimes, which are realized under small liquid velocity, $-0.25\,845~\text{m}~\text{s}^{-1}<u_{l0}<0$ , the particles rising in a downward liquid flow ( $u_{g}>0$ , $u_{l}<0$ ) reach the top of the channel (figure 8 d). In these regimes, the liquid velocity decreases in absolute value and the velocity of particles increases, tending to the velocity of a single particle floating up in an infinite volume of a viscous motionless liquid, $u_{max}$ , along the nozzle height.

Figure 8.

The flow in the diffuser in the frame of a hyperbolic model. (a) The channel geometry. (b–d) The profiles of (b) the volumetric particle content, (c) the third (mass) characteristic, and (d) the liquid and particle velocities. Here, $u_{max}$ is the velocity of a single particle floating up in an infinite volume of liquid.

Under liquid velocities $u_{l0}\leqslant -0.25\,846~\text{m}~\text{s}^{-1}$ , the parameter gradients exceed the limiting value and the solutions break off (points C $_{1}$ –C $_{3}$ ). Let us write the equation for $\text{d}\unicode[STIX]{x1D6FD}/\text{d}z$ (4.62a ) substituting coefficients B $_{1}$ , B $_{2}$ , C $_{1}$ , C $_{2}$ (4.64)–(4.66) in it:

From (4.67), it is seen that $\text{d}\unicode[STIX]{x1D6FD}/\text{d}z\rightarrow \infty$ if the denominator of the equation on the right-hand side is tending to zero.

Let us note that the expression in square brackets of the denominator in (4.67) is the numerator of (4.54), which is the equation for the third (mass) characteristic. The bound between the gradients tending to infinity and the characteristic $(\text{d}z/\text{d}t)_{3}$ approaching zero is clearly demonstrated by comparison of panels (b) and (c) of figure 8.

The series of solutions turning from a hidden line to the right ( $\text{OA}_{1}$ – $\text{OA}_{3}$ ) and left ( $\text{OC}_{1}$ – $\text{OC}_{3}$ ) point out that the line is a separatrix with a saddle singular point ‘S’ in it. There is a unique solution OSB with flow transition from a ‘supercritical’ to a ‘subcritical’ one. (Using the terms of gas dynamics, it can be said that the continuous regimes $\text{OA}_{1}$ are ‘supercritical’ ones, since the characteristic $(\text{d}z/\text{d}t)_{3}$ is always positive (figure 8 c).) The separatrix OSB cannot be obtained by a fourth integration of the system of ordinary differential equations (4.62) since it is not possible to guess an exact, with an infinite number of characters after the point, value of the liquid velocity $u_{l0}$ . At the slightest deviation from this exact value, the solution turns either left, as $\text{OA}_{1}$ do, or right, like $\text{OC}_{1}$ do (figure 8 b). To obtain the separatrix OSB, a special numerical method was required.

In a subcritical sector of the separatrix, SB, quick compaction of particles occurs. At that point, the absolute value of the negative stress gradient in the dispersed particle medium $\unicode[STIX]{x1D6F1}_{z}$ increases sharply; $\text{d}\unicode[STIX]{x1D6FD}/\text{d}z\rightarrow \infty$ (4.67), and the separatrix ‘breaks off’ (point B). The effect of the increase in the absolute value of $\unicode[STIX]{x1D6F1}_{z}$ with particle content $\unicode[STIX]{x1D6FD}$ increasing follows from (4.59), (4.67). The increase in $\unicode[STIX]{x1D6FD}$ results in the liquid velocity decreasing (4.56a ), which, according to (4.59), causes a subsequent increase in the absolute value of $\unicode[STIX]{x1D6F1}_{z}$ , which, since the regime is ‘subcritical’ and the denominator in (4.67) is negative, leads to the increase in the particle content $\unicode[STIX]{x1D6FD}$ . The rate of changing $\unicode[STIX]{x1D6FD}$ and the other parameters in the subcritical sector of the separatrix SB is determined by the value of $\unicode[STIX]{x1D70F}^{\ast }$ (the time within which a gradient flow is settled around a particle).

Thus, the calculations of liquid–particle countercurrent flow with a hyperbolic version of the ‘dirty water’ model (4.47)–(4.51) have shown the existence of a limiting regime $\text{OA}_{3}$ with a maximum by absolute value of liquid velocity, $u_{l0}\approx -0.25\,845~\text{m}~\text{s}^{-1}$ , for which a continuous solution with the particles rising up to the top nozzle section can be obtained.

4.8 The comparison of calculations with a hyperbolic and ‘classic’ model

A stationary version of the ‘classic’ model (22)–(27) of Lahey et al. (Reference Lahey, Cheng, Drew and Flaherty1980), taking into account gravity and neglecting the liquid compressibility ( $\unicode[STIX]{x1D70C}_{l}=\text{const.}$ ) and the drag force between the channel walls and the fluid $F_{w}$ , for a channel of variable cross-section is

After transformations, we obtain equations of the same form as (4.62), but with other expressions for B $_{1}$ , B $_{2}$ , C $_{1}$ and C $_{2}$ which do not contain the pressure due to neglect of liquid compressibility. Therefore, the volumetric particle contents and phase velocities can be calculated with the help of only (4.62a ), which, after substituting B $_{1}$ , B $_{2}$ , C $_{1}$ and C $_{2}$ , takes the form

The denominator of (4.71) vanishes only at a stagnation point, $u_{g}=u_{l}=0$ , and, therefore, discontinuous solutions are not possible.

Figure 9 compares calculations of the flow in the diffuser with the improved and classic models. For the continuous regimes $\text{OA}_{1}$ – $\text{OA}_{3}$ , the calculations do not differ. The discontinuous solutions $\text{OC}_{2}$ and $\text{OC}_{3}$ coincide up to the points of break off (C $_{2}$ , C $_{3}$ ). However, the classic model builds the solutions further.

Thus, the classic model rejects the existence of the limit regime $\text{OA}_{3}$ predicted by the hyperbolic model. ‘Numerical instabilities’ of solutions of the classic model can be due to solution bifurcation in the vicinity of a saddle singular point indefinable by this model.

Figure 9.

The flow in a diffuser in the frame of the ‘classic’ and hyperbolic models. (a) The channel geometry. (b,c) The profiles of (b) the volumetric particle content, (c) the liquid and particle velocities. Here, $u_{max}$ is the velocity of a single particle floating up in an infinite volume of liquid.

4.9 Superslow ‘creeping’ shock waves in a two-phase two-velocity medium

Consider a transient upward flow of solid particles. A cylindrical vessel of 20  cm height is occupied with water (figure 10 a). At $t=0$ , solid particles of radius $a=0.15~\text{cm}$ start being injected from the bottom of the tube with velocity $u_{g0}$ . The liquid velocity at the bottom is zero, $u_{l0}=0$ . The pressure at the top tube cross-section is equal to atmospheric pressure.

Let us consider the flow in the frame of the hyperbolic ‘dirty water’ model (4.47)–(4.51), in which the force of gravitation is additionally considered and the flow wall friction is neglected. Then, equations (4.48)–(4.50) take the form

For calculations, the ‘double triangle’ numerical method was used, described briefly in § 6, and in appendix B of Ivashnyov et al. (Reference Ivashnyov, Ivashneva and Smirnov2000) in detail.

Figure 10.

High-frequency vibrations in an upward flow of solid particles. (a) The flow scheme. (b–e) The oscillograms of (b,c) the pressure, (d) the volumetric particle content and (e) the volumetric mixture flux.

First, we will analyse the case when the particles are injected with a small velocity $u_{g0}=1~\text{cm}~\text{s}^{-1}$ . The volumetric rate of the injected particles is $\unicode[STIX]{x1D714}_{0}=\unicode[STIX]{x1D6FD}_{0}u_{g0}=10^{-3}~\text{m}^{3}~\text{s}^{-1}$ ( $\unicode[STIX]{x1D6FD}_{0}=0.1$ ).

The particle injection causes waves of two types: fast and slow ones. Figure 10(b–e) shows the evolution of fast waves. The particle injection leads to the formation of a compression wave which propagates up to the vessel top. On reaching it, the wave is reflected from the free surface as a wave of rarefaction and, by the instant $t\approx 2H/a_{l}=0.26~\text{ms}$ , returns to the bottom cross-section. On reaching the bottom end, the wave is reflected and goes upward. After reflecting from the open end, it returns into the vessel as a compression wave. As a result, there arise pressure vibrations with period $T\approx 4H/a_{l}=0.52~\text{ms}$ (figure 10 b).

As the particles fill up the vessel (figure 10 d), the ‘fork’ ideality is broken and the vibration amplitude decreases gradually (figure 10 b,c). To characterize the intensity of the acoustic vibrations, let us use a relative deviation of the volumetric mixture rate $\unicode[STIX]{x1D714}$ . Summing equations (4.47) of the phase continuity rewritten for constant density of the liquid, we obtain the integral

Thus, in the case of fluid component incompressibility, the volumetric mixture rate $\unicode[STIX]{x1D714}$ is constant and equal to its value at the entrance,

and therefore the amplitude of vibrations $\unicode[STIX]{x1D714}$ characterizes the mixture compressibility. As is seen from figure 10(e), by the time moment $t=0.5~\text{s}$ , the change in $\unicode[STIX]{x1D714}$ is comparable to its absolute value, $|\unicode[STIX]{x0394}\unicode[STIX]{x1D714}|\sim \unicode[STIX]{x1D714}_{0}$ .

Figure 11.

The evolution of slow waves in an upward flow of solid particles. (a,b) For seven instants (shown in the plot), the profiles of (a) the particle volumetric content and (b) the particle mass flow rate, $\unicode[STIX]{x1D714}_{g}=\unicode[STIX]{x1D6FD}u_{g}$ . (c,d) For $t=0.4~\text{s}$ , (c) the velocity profiles of the liquid, particles and third characteristic, and (d) the pressure profile. Here, $u_{max}\approx 0.28~\text{m}~\text{s}^{-1}$ is the velocity of a single particle floating up in an infinite volume of liquid.

Simultaneously with the acoustic waves, two slow waves whose velocities are determined by the phase velocities, $u_{g}$ and $(\text{d}z/\text{d}t)_{3}$ , travel up the vessel. In figure 11(a), a stepwise increase in the volumetric particle content rising to the vessel top is seen. Figure 11(c) shows the profiles of the velocities of the two mass characteristics. These velocities, practically coinciding, $u_{g}\approx (\text{d}z/\text{d}t)_{3}$ , decrease along the wave. Since the characteristic velocity decreases, backward perturbations lag behind forward ones and the waves expand (like a rarefaction wave in an ideal gas). Thus, the wave shown in figure 11(a) represents a combination of two expanding waves.

Since both characteristic velocities are positive, $u_{g}>0$ , $(\text{d}z/\text{d}t)_{3}>0$ , the wave moves only upward and, after its passage, the particle volumetric rate equalizes within the tube height, $\unicode[STIX]{x1D714}_{g}=\unicode[STIX]{x1D6FD}u_{g}$ (figure 11 b).

High-frequency oscillations do not have a visible effect on slow-wave evolution, obviously due to the fact that rarefaction and compression quick waves compensate each other’s influence. By the time moment $t=0.7~\text{s}$ , when slow waves have reached the upper cross-section, the vibrations caused by the quick-wave propagation attenuate, and a stationary flow satisfying the integral ratio

is settled.

Figure 12.

The upward flow of solid particles in the case when the velocity of particle injection exceeds the velocity of a single particle floating up in an infinite volume of liquid, $u_{max}$ (shown in c). (a–b) The profiles of the particle (a) volumetric content and (b) mass flow rate, $\unicode[STIX]{x1D714}_{g}=\unicode[STIX]{x1D6FD}u_{g}$ , at different time moments. (c) The velocity profiles of the liquid, particles and the third characteristic at $t=0.4~\text{s}$ . (d) The pressure profiles. Here, $D$ is the velocity of the shock (second) wave determined as the travelled distance/time ratio.

Let us now analyse the second case, when the particles are injected with a velocity exceeding $u_{max}$ , the velocity of a single particle floating up in an infinite volume of liquid,

The volumetric rate of particle injection is $\unicode[STIX]{x1D714}_{0}=\unicode[STIX]{x1D6FD}_{0}u_{g0}=0.005~\text{m}^{3}~\text{s}^{-1}$ ( $\unicode[STIX]{x1D6FD}_{0}=0.01$ ).

The profile of volumetric particle content (figure 12 a) shows a peak formed by two slow waves, a ‘plume of particles’ rising up to the vessel top. The profiles of the velocities of the mass invariant spreading are given in figure 12(c). The velocities coincide up to the moment of slow-wave arrival and then separate, decreasing along the first wave and increasing along the second one. Therefore, the first wave expands and the second one steepens, and turns into a shock wave.

Any perturbation spreads with a characteristic velocity. As the velocities of the first two characteristics exceed that of the third one (figure 12 c),

the velocity of the forward going stretching wave is determined by the first characteristic, $(\text{d}z/\text{d}t)_{1,2}=u_{g}$ , and the movement of the lagging shock one is defined by the third characteristic, $(\text{d}z/\text{d}t)_{3}$ .

The bound between the third characteristic and the shock wave is determined by the condition of jump evolution. Passing on to a moving frame linked with the jump and neglecting the change in the liquid density, we obtain five ratios preserved across the jump:

where $v_{l}$ and $v_{g}$ are the liquid and particle velocities in the moving frame; $\unicode[STIX]{x1D70C}_{l}^{\ast }=\unicode[STIX]{x1D70C}_{l}(P_{atm})=\text{const.}$ represents the liquid density.

There are 11 unknown parameters in five equations (4.79)–(4.83): $\unicode[STIX]{x1D6FD}$ , $P$ , $v_{l}$ , $v_{g}$ , $\unicode[STIX]{x1D6F1}_{x}$ are unknown on both sides of the shock wave, and the wave velocity $D$ is not known either. Therefore, six additional ratios should be added to the five equations to resolve them. Invariants preserved along five characteristics are such ratios, but since there are only five characteristics, one of them must come to the jump twice: from the fields ahead of and behind it (figure 13). It is this characteristic of coming twice with which the jump is bounded. Then, the jump velocity $D$ must satisfy the condition

where $[(\text{d}x/\text{d}t)_{3}]_{1}$ and $[(\text{d}x/\text{d}t)_{3}]_{2}$ are characteristic velocities before and behind the front of the shock wave.

Figure 13.

The bound between characteristics and a shock wave.

Figure 12(c) clearly shows that condition (4.86) is satisfied.

5.1 Basic model assumptions

A model of a three-phase medium is suggested. The phases are fluid, immovable $\unicode[STIX]{x1D6FC}$ -bubbles growing at walls and $\unicode[STIX]{x1D6FD}$ -bubbles moving in the flow (figure 6).

1. (i) Boiling is considered to start at the channel wall surface. The number of centres of nucleation per unit wall surface area $N$ is a free model parameter. In the flow volume, boiling centres are initially lacking and appear only after the fragmentation of wall $\unicode[STIX]{x1D6FC}$ -bubbles.

2. (ii) The fragmentation of both $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles is induced by Kelvin–Helmholtz instability. The critical value of the Weber number, $We^{\ast }$ , is the second of two free model parameters.

3. (iii) The size and velocity distribution of $\unicode[STIX]{x1D6FD}$ -bubbles across a channel section is neglected.

4. (iv) The $\unicode[STIX]{x1D6FC}$ -bubbles are assumed to be in the form of a hemisphere and the $\unicode[STIX]{x1D6FD}$ -bubbles are ball-shaped. (The supposition about regular geometry of the bubbles is used for the interphase drag and heat transfer models.)

5. (v) The pressures in the phases are considered to be equal, $P_{l}=P_{\unicode[STIX]{x1D6FC}}=P_{\unicode[STIX]{x1D6FD}}=P$ .

   This assumption is accepted on the grounds of the evaluation of the time of pressure equalization between a motionless vapour bubble and the evaporating fluid surrounding it. This process studied using the model for a spherical bubble in an inertial liquid has shown that the pressures in the bubble and liquid equalize within only 2–3 bubble oscillations. (The evaporation is an essential damping factor.) The oscillation period can be evaluated by the formula

   (5.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{b}\sim a\sqrt{\frac{\unicode[STIX]{x1D70C}_{l}}{P}}.\end{eqnarray}$$

   The characteristic dimensions of $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles are determined by the fragmentation process and can be evaluated from the expression for the Weber number (1.8). Assuming that the critical Weber number $We^{\ast }=1$ , using the evaluation for the flow velocity $u_{l}^{\ast }\sim 10~\text{m}~\text{s}^{-1}$ (4.1) and accepting the difference in the phase velocities $\unicode[STIX]{x0394}u_{l}^{\ast }=1~\text{m}~\text{s}^{-1}$ , we can evaluate $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubble dimensions as

   (5.2a,b ) $$\begin{eqnarray}a_{\unicode[STIX]{x1D6FC}}\leqslant \frac{\unicode[STIX]{x1D70E}We^{\ast }}{2\unicode[STIX]{x1D70C}_{g}(u_{l}^{\ast })^{2}}\sim 10^{-6}~\text{m},\quad a_{\unicode[STIX]{x1D6FD}}\leqslant \frac{\unicode[STIX]{x1D70E}We^{\ast }}{2\unicode[STIX]{x1D70C}_{g}(\unicode[STIX]{x0394}u_{l}^{\ast })^{2}}\sim 10^{-4}~\text{m}.\end{eqnarray}$$

   Thus, the oscillation period is $\unicode[STIX]{x1D70F}_{b}\sim 10^{-6}~\text{s}$ even for large $\unicode[STIX]{x1D6FD}$ -bubbles, that is several orders of magnitude less than the characteristic time of the pressure change in the compression wave, which does not exceed $\unicode[STIX]{x1D70F}_{P}\sim 10^{-3}~\text{s}$ (figure 1). Since $\unicode[STIX]{x1D70F}_{b}\ll \unicode[STIX]{x1D70F}_{P}$ , we can neglect the inertia of the liquid surrounding a bubble and adopt the assumption about the equality of the pressures in the phases (v).

6. (vi) The parameters in the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles are uniform over their radii and equal to the parameters on the line of saturation, $T_{g}=T_{s}(P)$ , $\unicode[STIX]{x1D70C}_{g}=\unicode[STIX]{x1D70C}_{gs}(P)$ .

7. (vii) The vapour density is much less than the liquid density, $\unicode[STIX]{x1D70C}_{g}\ll \unicode[STIX]{x1D70C}_{l}$ .

   From this assumption, the equality of the mean-mass mixture velocity $u$ and the liquid velocity $u_{l}$ follows:

   (5.3) $$\begin{eqnarray}u=(1-\unicode[STIX]{x1D712})u_{l}+\unicode[STIX]{x1D712}u_{g}\approx u_{l},\quad \unicode[STIX]{x1D712}=\frac{\unicode[STIX]{x1D70C}_{g}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D70C}}.\end{eqnarray}$$

   Here, $\unicode[STIX]{x1D712}$ is the mass fraction of $\unicode[STIX]{x1D6FD}$ -bubbles in the flow; $u_{g}$ is the velocity of $\unicode[STIX]{x1D6FD}$ -bubbles.

8. (viii) The velocity of fragments separating from $\unicode[STIX]{x1D6FC}$ -bubbles is assumed to be zero, $u_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}=0$ . (After the separation, the fragments instantly acquire the $\unicode[STIX]{x1D6FD}$ -vapour velocity (according to assumption (iii)).) The velocity of the liquid vaporizing into $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles is considered to be equal to the velocity of $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -vapour, $u_{l\unicode[STIX]{x1D6FC}}=0$ , $u_{l\unicode[STIX]{x1D6FD}}=u_{g}$ .

5.2 Equations of mass balance

The one-dimensional mass equations for the liquid, $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -vapour have the form

where $\unicode[STIX]{x1D6FC}$ is the volumetric fraction of $\unicode[STIX]{x1D6FC}$ -vapour in the whole volume of the three-phase mixture; $N=\text{const.}$ is the number of $\unicode[STIX]{x1D6FC}$ -bubbles per unit area of the wall surface; $D_{ch}$ is the inner diameter of the channel; $n$ and $c$ are the numbers of $\unicode[STIX]{x1D6FD}$ -bubbles per unit volume and mass of the moving mixture (not including $\unicode[STIX]{x1D6FC}$ -bubbles); $\unicode[STIX]{x1D70C}$ is the density of the moving mixture; $J_{l\unicode[STIX]{x1D6FC}}$ and $J_{l\unicode[STIX]{x1D6FD}}$ are the liquid masses evaporating into the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles per unit volume of all the mixture in unit time; $J_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ is the mass of $\unicode[STIX]{x1D6FC}$ -vapour converted into $\unicode[STIX]{x1D6FD}$ -vapour per unit mixture volume in unit time; $a_{\unicode[STIX]{x1D6FC}}$ and $a_{\unicode[STIX]{x1D6FD}}$ are the radii of $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles.

The intensities of liquid evaporation into $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles are determined by

where $j_{l\unicode[STIX]{x1D6FC}}$ and $j_{l\unicode[STIX]{x1D6FD}}$ are the rates of liquid evaporation into each bubble of the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -phases; $\unicode[STIX]{x1D706}_{ls}$ and $D_{ls}$ are the coefficients of thermal conductivity and thermal diffusivity of the fluid at the curve of saturation.

The dependence of the dimensionless heat flux into a bubble, the Nusselt number, on the dimensionless liquid overheat, $Nu(Ja)$ , is determined from a self-similar solution of the problem of heat bubble growth in a superheated liquid (Scriven Reference Scriven1959) with formula (1.5).

The intensity of the phase transition between the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -vapour bubbly phases is

where $\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FC}}$ is the intensity of $\unicode[STIX]{x1D6FC}$ -bubble breakup: the number of $\unicode[STIX]{x1D6FC}$ -bubble fragments passing into the flow from a unit of the wall surface per unit time.

5.3 Equations of momentum balance

The momentum equations for the fluid and the $\unicode[STIX]{x1D6FD}$ -bubbles moving in it are

Here, $\text{D}/\text{D}t$ and $\text{D}_{g}/\text{D}t$ are the material derivatives following the liquid and $\unicode[STIX]{x1D6FD}$ -bubbles; $R_{\unicode[STIX]{x1D6FC}l}$ and $R_{\unicode[STIX]{x1D6FD}l}$ are the specific forces of the vapour–liquid interaction for the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -vapour phases per unit mixture volume; $F_{w}$ is the drag force between the walls and the flow per unit mixture volume.

On generalizing the solution for a single bubble with a changeable radius in a gradient flow of an incompressible non-viscous fluid (appendix C) to a compressible viscid one, we obtain the expression for the interphase force the same way as was done in § 3.6:

Here, $\unicode[STIX]{x1D6F1}_{x\unicode[STIX]{x1D6FD}}$ is the gradient of the stresses in the $\unicode[STIX]{x1D6FD}$ -bubble medium; $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}^{\ast }$ is the characteristic time of the $\unicode[STIX]{x1D6FD}$ -bubble interaction set-up, the characteristic time within which a perturbation runs the distance between neighbouring bubbles; $S$ is the drag force per unit velocity for $\unicode[STIX]{x1D6FD}$ -bubbles; $c_{\unicode[STIX]{x1D707}}$ is the interfacial drag coefficient; $Re$ is the bubble Reynolds number; $\unicode[STIX]{x1D707}_{l}$ is the coefficient for the liquid dynamic viscosity.

The first term in square brackets in (5.12) is the buoyancy force; the sum of the following two terms is the virtual mass force. The first component of the virtual mass force is determined by the bubble acceleration; the second one is determined by the change in the bubble radius: the virtual mass of a growing bubble increases and it slows down in contrast to a collapsing bubble which loses part of its virtual mass and accelerates (appendix A in Ivashnyov et al. Reference Ivashnyov, Ivashneva and Smirnov2000).

The term $R_{\unicode[STIX]{x1D6FC}l}$ , the specific interfacial force between the liquid and the $\unicode[STIX]{x1D6FC}$ -bubbles, is obtained from expression (5.12) for $R_{\unicode[STIX]{x1D6FD}l}$ by replacing the $\unicode[STIX]{x1D6FD}$ -bubble velocity by the velocity of $\unicode[STIX]{x1D6FC}$ -bubbles, i.e. zero, and the $\unicode[STIX]{x1D6FD}$ -bubble volume by that for $\unicode[STIX]{x1D6FC}$ -bubbles which is the volume of a hemisphere:

Here, $\unicode[STIX]{x1D6F1}_{x\unicode[STIX]{x1D6FC}}$ is the stress gradient in the $\unicode[STIX]{x1D6FC}$ -bubble medium; $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FC}}^{\ast }$ is the characteristic time of the $\unicode[STIX]{x1D6FC}$ -bubble interaction set-up; $f_{\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FC}}$ is the drag force for an $\unicode[STIX]{x1D6FC}$ -bubble; $c_{\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FC}}=c_{\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FC}}(Re_{u})$ is the drag coefficient for an $\unicode[STIX]{x1D6FC}$ -bubble; $Re_{u}=2a_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D70C}_{l}|u|\unicode[STIX]{x1D707}_{l}^{-1}$ is the flow Reynolds number.

5.4 Interfacial force generalization to the case of fragmenting bubbles

The breaking was not taken into account in the expression for the force of the interfacial interaction between the $\unicode[STIX]{x1D6FD}$ -bubbles and the liquid, $R_{\unicode[STIX]{x1D6FD}l}$ (5.12). Formula (5.12) predicts the acceleration of breaking bubbles since their radii diminish. For our flow, the change in bubble radii occurs not only due to their growth but due to their breakup as well. However, this is in contradiction with physical reality as the total volume of bubbles and, consequently, the total virtual mass must not change under the breaking process. Thus, the breakup should not lead to bubble acceleration. Formula (5.12) should be modified for these types of flows.

In order to exclude the non-physical effect, we will pass on from the parameters characterizing a single bubble to mean parameters of the vapour phase. From the mass balance equation for a single (not fragmenting) bubble,

we obtain

Taken separately, either of (5.12) or (5.19) is not valid for breaking bubbles, but their composition becomes valid after passing on to the parameters that do not change under breaking. Multiplying (5.19) by $\unicode[STIX]{x1D70C}_{l}(u_{g}-u)/2$ and summing up the result with (5.12), we obtain

Equation (5.20) excludes the abovementioned defect of a non-physical bubble acceleration.

Expression (5.16) for the force of the interfacial interaction between the liquid and the $\unicode[STIX]{x1D6FC}$ -bubbles is generalized similarly:

The momentum equation for the $\unicode[STIX]{x1D6FD}$ -vapour (5.11) can be simplified. The order of magnitude of the term on the left-hand side of (5.11) can be evaluated as

where $W$ is the characteristic acceleration. The order of magnitude of the first component of $R_{\unicode[STIX]{x1D6FD}l}$ is

Since $\unicode[STIX]{x1D70C}_{g}\ll \unicode[STIX]{x1D70C}_{l}$ , the term on the left-hand side of (5.11) can be neglected and the equation can be rewritten as

On substituting the expressions for $R_{\unicode[STIX]{x1D6FC}l}$ and $R_{\unicode[STIX]{x1D6FD}l}$ from (5.21) and (5.20) into (5.10) and (5.24), we obtain the equations for the mixture momentum and phase slip in the forms

Here,

is the sum drag force between the wall with $\unicode[STIX]{x1D6FC}$ -bubbles located on it and the fluid per unit mixture volume.

5.5 The dependence of the drag force between the flow and the walls

Let us evaluate the friction losses in the flow under consideration. At the first stage, the wave of the pressure increase with amplitude $\unicode[STIX]{x0394}P_{1}^{\ast }\approx 1~\text{MPa}$ and length $L\approx 1~\text{m}$ forms behind the pressure undershoot.

The force of the wall resistance to the flow motion is evaluated as

where $\unicode[STIX]{x1D709}$ is the dimensionless friction coefficient for the tube approximated as a dependence on the channel Reynolds number,

For $u\sim 10~\text{m}~\text{s}^{-1}$ , $\unicode[STIX]{x1D70C}_{l}=741~\text{kg}~\text{m}^{-3}$ , $\unicode[STIX]{x1D707}_{l}=10^{-4}~\text{kg}~\text{m}^{-1}~\text{s}^{-1}$ , $D_{ch}=0.075~\text{m}$ , we obtain $Re_{D}\approx 5\times 10^{6}$ .

Under such $Re$ values, a universal (logarithmic) Prandtl dependence is valid for smooth tubes (Schlichting Reference Schlichting1964). From the dependence, we find $\unicode[STIX]{x1D709}(Re)\approx 5\times 10^{-2}$ . From (5.28), we have

Estimation (5.30) shows that, for the processes at the stage of fast waves, which we are going to study, the friction forces are negligibly small compared with inertia ones and can be neglected.

To verify our model, we are going to calculate the second main stage as well. This stage is characterized by a ‘slow wave of boiling’, the velocity of which is determined from the condition of Weber number equality to its critical value ahead of the wave (Ivashnyov & Ivashneva Reference Ivashnyov and Ivashneva2012). The pressure before the wave front is practically uniform and the inertia forces there are not large and can be evaluated as

That is, friction forces can be comparable to inertia ones ahead of the wave of boiling.

A characteristic time of the main stage is $\unicode[STIX]{x1D70F}_{2}^{\ast }\sim 0.1~\text{s}$ . In the work of Kucienska, Seynhaeve & Giot (Reference Kucienska, Seynhaeve and Giot2008), an oscillatory process with period ${\sim}0.01~\text{s}$ was calculated with stationary and non-stationary friction models. The discrepancy was ${\approx}10\,\%$ . Therefore, for modelling our, order of magnitude slower, flow, we can use the expressions for the drag force coefficient $\unicode[STIX]{x1D709}$ for stationary flows in tubes (Schlichting Reference Schlichting1964):

5.6 Equation for the energy balance

As distinct from the previous model (Ivashnyov et al. Reference Ivashnyov, Ivashneva and Smirnov2000) for fluxes with a negligibly small phase slip, the present one is intended to be applied for flows where the difference in the phase velocities is comparable to the mean-mass velocity. In this case, the kinetic energy of the small-scale motion of liquid drawn into the movement by bubbles can be essential and should be taken into account.

The components of the interfacial force $R_{\unicode[STIX]{x1D6FD}l}$ may be divided into two groups. In the first one, there are two first components in the curly braces of formula (5.20) obtained from the solution of the problem of a particle overflowing by an ideal liquid:

The third and fourth components added to account for the liquid non-ideality are attributed to the second group:

Since the work of the forces of the first group must not cause an increase in the liquid temperature, we will suppose the work of $R_{\unicode[STIX]{x1D6FD}l}^{\prime }$ to transform into the energy of a small-scale liquid motion $k_{l}$ , $[k_{l}]=\text{m}^{2}~\text{s}^{-2}$ , in order to avoid a non-physical variation in temperature. Assuming that the work of only $R_{\unicode[STIX]{x1D6FD}l}^{\prime \prime }$ forces changes the liquid heat content, $i_{l}$ , the caloric equation for the liquid is written as

where $i_{l}$ is the liquid enthalpy; $i_{ls}$ is the liquid enthalpy on the line of saturation; $Q_{l\unicode[STIX]{x1D70E}}$ is the heat flux from the liquid (‘ $l$ ’) to the phase interface (‘ $\unicode[STIX]{x1D70E}$ ’); $A$ is the sum work of the second group components $R_{\unicode[STIX]{x1D6FC}l}^{\prime \prime }$ and $R_{\unicode[STIX]{x1D6FD}l}^{\prime \prime }$ of the interfacial force, which is always positive, $A>0$ .

The first two terms on the right-hand side of (5.35) are the heat flux to the liquid. The second one, $(J_{l\unicode[STIX]{x1D6FC}}+J_{l\unicode[STIX]{x1D6FD}})(i_{l}-i_{ls})$ , is the ‘retained heat’: since the enthalpy $i_{ls}$ of the liquid layer adjacent to the interfacial surface differs from the average enthalpy $i_{l}$ , part of the heat flux is wasted to heat the liquid itself. Only a part, $Q_{l\unicode[STIX]{x1D70E}}$ , of the heat flux reaches the interfacial surface and is wasted on evaporation.

The equation for the balance of heat fluxes on the interfacial surface is

Here, $Q_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70E}}$ and $Q_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70E}}$ are the heat fluxes from the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles to the phase interface per unit volume and unit time; $h$ is the specific heat of vaporization; $i_{g}$ is the vapour enthalpy.

The equations of conservation for $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -vapour are

The balance equation for the kinetic energy of a small-scale liquid motion has the form

Then, the law of total energy conservation can be written in the form

Applying the second law of thermodynamics to multiphase flow models, Ishii & Hibiki (Reference Ishii and Hibiki2006) have shown the necessity of ensuring positiveness of the interfacial entropy generation, $\unicode[STIX]{x0394}$ . Using Gibbs ratios,

in the equations of heat energy for the phases (5.35), (5.38), we obtain the balance equation for the mixture entropy:

Here $S_{l}$ and $S_{g}$ are the liquid and vapour entropies; $T_{s}$ is the temperature on the line of saturation, which is equal, according to the adopted assumption, to the vapour temperature $T_{s}=T_{g}$ (for vaporizing liquid, $T_{s}\leqslant T_{l}$ ).

5.7 Interfacial area transport equation

As the number of wall surface bubbles does not change, the IAT equation of the model has the form of only the $\unicode[STIX]{x1D6FD}$ -bubble number balance:

where $\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FC}}$ and $\unicode[STIX]{x1D713}_{\unicode[STIX]{x1D6FD}}$ are the fragmentation rates of the $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles, which are defined with relaxation ratios

where $N^{\ast }$ and $n^{\ast }$ are the numbers of bubbles that would be formed if the fragmentation were instantaneous; $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FC}}^{\ast }$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}^{\ast }$ are the characteristic fragmentation times for $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles.

The parameters $N^{\ast }$ , $n^{\ast }$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FC}}^{\ast }$ , $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}^{\ast }$ are specified using the solution of the problem on the rise in the amplitude of a small harmonic perturbation with wavelength $\unicode[STIX]{x1D706}$ arising on a plane interphase boundary (Lamb Reference Lamb1957; Birkhoff Reference Birkhoff1960; Ishii & Hibiki Reference Ishii and Hibiki2006):

where $\unicode[STIX]{x1D70E}$ is the surface tension coefficient.

From the condition that a bubble is divided by the wave with the fastest growing amplitude,

and the condition that the length of the wave cannot exceed the bubble diameter, $\unicode[STIX]{x1D706}\ast \leqslant 2a$ , one obtains the criterion of surface stability of the bubble overflowing by fluid, the Weber number,

and the evaluation of its critical value, $We^{\ast }=3\unicode[STIX]{x03C0}$ .

The number of fragments appearing as a result of a bubble breakup ( $N^{\ast }/N$ , $n^{\ast }/n$ ) is determined as the ratio of the bubble diameter to the length of the wave causing breakup:

The characteristic time of breakup, $\unicode[STIX]{x1D70F}^{\ast }$ , is adopted as the ‘ $e$ ’-time ( $e=2.7\ldots$ ) rise in the amplitude of the perturbation with wavelength $\unicode[STIX]{x1D706}^{\ast }$ :

5.8 Equations for phase states

Neglecting the dependence of the liquid density on the temperature and considering the isotherms in pressure-specific volume (PV) coordinates to be straight lines near the curve of saturation (Scripov Reference Scripov1972), one can write the equation of state for the liquid as

where $P_{s}(T_{0})$ , $\unicode[STIX]{x1D70C}_{ls}(T_{0})$ and $\unicode[STIX]{x1D6FD}_{l}(T_{0})$ are the pressure, the liquid density and the coefficient of isothermal liquid compressibility on the line of saturation at initial liquid temperature $T_{0}$ .

The caloric equation for the liquid follows from the equation of its state (5.51):

where $i_{l}$ is the liquid enthalpy; $i_{ls}$ is the liquid enthalpy on the line of saturation.

The isobaric specific heat capacity of the liquid, $c_{l}$ , is determined by differentiating (5.53):

The liquid entropy is a function of its temperature and is calculated as follows:

where $S_{l0}$ is the liquid entropy at the temperature $T_{0}$ .

Within a wide range of temperatures, $450~\text{K}\leqslant T_{s}\leqslant 590~\text{K}$ , simple relationships

describe the vapour density and enthalpy on the line of saturation with a maximal relative error of 2 %. Here, $a_{g}^{2}=2\times 10^{5}~\text{m}^{2}~\text{s}^{-2}$ is the approximation constant.

5.9 General form and characteristics of the model

Thus, the boiling flow model accounting for wall surface nucleation and bubble breakup includes 11 differential equations in partial derivatives: the mass balance equations for the three phases (5.4), (5.5); the equations for the mixture momentum (5.25) and the phase slip (5.26); the equations for the stresses in the disperse particle media (5.13), (5.17); heat energy balance equations for each of the three phases (5.35), (5.38); and the IAT equation (5.44). All 11 mathematical characteristics are real and always remain real under any vapour content. Eight of them are found exactly:

The other three, two ‘sonic’ and one mass one, are determined from the solution of the cubic equation

Here, $c_{f}$ is the parameter characterizing the compressibility of the three-phase mixture under a ‘frozen’ heat exchange:

where $\unicode[STIX]{x1D6FE}$ is the summary volumetric content of $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -phases in the channel; $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6F4}}$ and $V$ are the density and specific volume of the three-phase mixture:

Here, $\unicode[STIX]{x1D712}_{l}$ , $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D6FC}}$ and $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D6FD}}$ are the mass contents of the liquid, $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles in the three-phase mixture:

$\text{A}_{0}$ , $\text{A}_{1}$ , $\text{B}_{0}$ , $\text{B}_{1}$ and $\text{C}_{1}$ are dimensionless coefficients:

where $\unicode[STIX]{x1D6EC}$ is a dimensionless complex of parameters.

Searching the roots of (5.58) in the form of an expansion in the small parameter $\unicode[STIX]{x1D700}=\unicode[STIX]{x0394}u/c_{f}\ll 1$ , we have

6.1 Idea of the method

Any hyperbolic system of equations can be rewritten in a characteristic form. For example, for the equations of an ideal gas (4.5), it is the form of Riemannian invariants (4.6), (4.7):

where $I_{i}$ is the invariant; $a_{i}$ is the velocity of $I_{i}$ -invariant propagation.

Depending on the sign of $a_{i}$ , a ‘left triangle’ ( $a_{i}>0$ ) or ‘right triangle’ ( $a_{i}<0$ ) scheme is applied to (6.1). The left triangle scheme has the form

where $j$ is the number of the circuit node; $k$ is the number of the time layer; $\unicode[STIX]{x0394}t$ and $\unicode[STIX]{x0394}x$ are the integration steps in time and coordinate.

The only condition for stability of this scheme is to satisfy the Courant criterion,

where $C$ is the Courant number, which is the ratio of the physical and scheme velocities of information spreading. It is verified that the scheme ensures numerical solution convergence to the exact solution under diminution of the mesh size $\unicode[STIX]{x0394}x$ . A more detailed description of the method with calculation examples is given in appendix B of Ivashnyov et al. (Reference Ivashnyov, Ivashneva and Smirnov2000).

7.1 Boundary conditions

The law of the pressure change at the exit cross-section has been chosen so as to achieve agreement between the experimental and theoretical rates of the pressure drop in the first wave of rarefaction:

where $P_{atm}$ is the atmospheric pressure.

Near the closed vessel end, the condition of equality of the flow velocity to the middle velocity of the layer of wall bubbles, $v_{\mathit{cs}}$ , is applied, which is expressed in terms of the radial velocity of the wall bubbles as

where $N$ is the number of bubbles per unit wall area; ${\dot{a}}_{\unicode[STIX]{x1D6FC}}$ is the radial velocity of an $\unicode[STIX]{x1D6FC}$ -bubble.

Substituting ${\dot{a}}_{\unicode[STIX]{x1D6FC}}$ deduced from the balance equation for the $\unicode[STIX]{x1D6FC}$ -bubble mass,

we obtain

7.2 Choice of the free model parameters

The model has two free parameters: the number of boiling centres per unit wall area $N$ and the critical value of the Weber number $We^{\ast }$ . First, we ‘froze’ the bubble fragmentation by specifying $We^{\ast }\rightarrow \infty$ , and chose the number of centres $N$ to achieve agreement between the minimum pressure levels in the calculated and experimental oscillograms. We obtained $N=4\times 10^{6}~\text{m}^{-2}$ (figure 14 a).

This number recalculated in terms of the mixture mass is

Then, having fixed $N$ , we chose the critical Weber number giving the best agreement between the experimental and theoretical widths of the compression wave in the calculations and experiment (figure 14 b). We obtained $We^{\ast }=2$ .

Figure 14.

The adjustment of the free model parameters.

7.3 Model validation and verification

With no change in the free parameters chosen to simulate an early stage of decompression, the model is tested for its ability to simulate ‘boiling shocks’ at the main stage of blowdown. These ‘shocks’ were simulated by our previous model accounting for an in-volume nucleation (Ivashnyov et al. Reference Ivashnyov, Ivashneva and Smirnov2000). Figure 15 compares the calculated and experimental large-scale pressure oscillograms for the two experiments of Edwards & O’Brien (Reference Edwards and O’Brien1970). The depressurization process is calculated with four models: equilibrium, boiling at constant number of bubbles ( $c=10^{6}~\text{kg}^{-1}$ ) and two models accounting for bubble breakup, assuming boiling up in the flow volume ( $c_{0}=4\times 10^{5}~\text{kg}^{-1}$ , $We^{\ast }=1$ ) and at the channel walls ( $N=4\times 10^{6}~\text{m}^{-2}$ , $We^{\ast }=2$ ). For the second experiment, the calculations with our previous model accounting for an in-volume nucleation coincide with those obtained using the equilibrium model (dotted and dash–dot lines in figure 15 b), since the previous model predicts conversion of the mixture into an equilibrium state as early as at a primary stage of the decompression (figures 3 d, 5 a,b).

Thus, the model accounting for a wall surface nucleation (thin continuous line in figure 15) is ‘able’ to model boiling shocks at the main stage of efflux for both experiments of Edwards & O’Brien (Reference Edwards and O’Brien1970).

Figure 15.

Large-scale oscillograms of pressure and volumetric vapour content under the conditions of the two experiments of Edwards & O’Brien (Reference Edwards and O’Brien1970): (a,b) $P_{0}=7~\text{MPa}$ , $T_{0}=515~\text{K}$ ; (c,d) $P_{0}=10.3~\text{MPa}$ , $T_{0}=558~\text{K}$ . All of the oscillograms are given for a cross-section of $x=1.5~\text{m}$ , except that the pressure oscillograms for the second experiment are shown for a cross-section of $x=2~\text{m}$ (there are no experimental data for $x=1.5~\text{m}$ ). For $\unicode[STIX]{x1D6FE}$ , see formula (5.60).

With the same parameters, specified for one-dimensional flow under a long tube depressurization, the model supposed was applied for an essentially two-dimensional flow through a short tube with a sharp right-angled entrance. A comparison of the calculation results with the experiment and the predictions of the three models referred to above is given in appendix A.

Table 1 gives the relative errors of the four model predictions for three experiments: the two experiments of Edwards & O’Brien (Reference Edwards and O’Brien1970) on a high-pressure tube depressurization and the experiment of Tikhonenko et al. (Reference Tikhonenko, Kevorkov and Lutovinov1979) on saturated water efflux through a short tube with a right-angled entrance. The free parameters of all of the models were specified to simulate early stages of vessel depressurization (figures 3 a,c and 14). Therefore, in the table, the errors are only presented for the main stages of decompression (figure 15). The errors are defined from the mean deviations of the theoretical large-scale pressure oscillograms from the experimental ones. For the experiments of Tikhonenko et al. (Reference Tikhonenko, Kevorkov and Lutovinov1979), the errors are determined from the accuracy of the mass flow rate predictions.

It is seen that the model supposed ensures the closest correspondence to experimental data. Besides, it simulates with practically the same error long-tube decompression and a flow through a short nozzle, i.e. the predicative accuracy does not change under an essential change in the type of flow simulated.

7.4 Estimation of the number of wall nucleation centres from an experiment of another type

In the paper of Tokita et al. (1997), copper plates were electrically heated in a vessel filled with distilled degassed water. Parallel gashes with a rectangular cross-section of 1 mm depth were dented in the plates. The gashes were spaced by 1mm. The width of the gashes was varied from 0.2 to 1.0 mm. Photographs showed that bubbles were formed only in the gashes of the plates. A dependence of the specific heat flux  $q$ on the plate overheating $\unicode[STIX]{x0394}T=T_{w}-T_{s}(P)$ was measured, where $T_{w}$ is the wall temperature and $T_{s}$ is the saturation temperature of the liquid. The value of $\unicode[STIX]{x0394}T$ permits judgements to be made about the intensity of heat exchange between the plate and the adjoining mixture of fluid and vapour. The less $\unicode[STIX]{x0394}T$ is, the more intensive the heat exchange is.

The value of $\unicode[STIX]{x0394}T$ differed depending on the width of the gashes under a fixed  $q$ . Under a minimal width, 0.2 mm, the wall overheating was minimal. The overheating increased with widening of the gashes, attaining a maximum value coinciding with $\unicode[STIX]{x0394}T$ for a smooth plate when the gash width was increased to 1.0 mm. This permits estimation of the distance between wall nucleation centres on a smooth plate as ${\sim}1~\text{mm}$ . In terms of the surface density of centres  $N$ , it is ${\sim}10^{6}~\text{m}^{-2}$ , which corresponds in order of magnitude to the value chosen in § 7.2 ( $N=4\times 10^{6}~\text{m}^{-2}$ ).

Kolev (Reference Kolev2005) collected experimental dependences of the number of nucleation centres on wall superheating for different heated surfaces, polished copper, brass and steel, in motionless or slowly flowing water (under pool boiling). According to data from 14 papers, the range of active nucleation centre densities lies between $10^{2}$ and $10^{7}~\text{m}^{-2}$ . Our number specified for wall nucleation centres is close to the upper value, obviously due to the fact that the liquid overheats after the depressurization, 10–30 K, corresponding to an upper bound on wall overheating of ${\sim}100~\text{K}$ .

8.1 The structure in calculations with the full model

Figure 16(a) compares experimental and calculated pressure profiles at 1, 2, 3 and 4 ms. For the instant 3 ms, the calculated channel length distributions of different parameters are shown in figure 16(b–d). Thin vertical lines show the channel cross-sections of the beginning of the initial decompression wave (A–A), its end coinciding with the beginning of the compression wave (B–B) and the compression wave end (C–C) for the instant $t=3$  ms.

Figure 16.

The parameter distribution along the channel under the conditions of the second experiment of Edwards & O’Brien (Reference Edwards and O’Brien1970):  $P_{0}=10.3~\text{MPa}$ , $T_{0}=558$  K. (a) The experimental (points) and theoretical (solid lines) pressure profiles at moments 1, 2, 3 and 4 ms. (b–d) For $t=3$  ms, the profiles of (b) the volumetric contents of $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -vapour, number of boiling centres, radii of $\unicode[STIX]{x1D6FC}$ - and $\unicode[STIX]{x1D6FD}$ -bubbles; (c) the mixture density and temperatures of the liquid and vapour phases; (d) the liquid velocity, the difference in the phase velocities, the velocity of the 11th invariant propagation and the Weber numbers. The cross-sections of the wave origin, and the minimum and maximum of the pressure profile at the instant 3 ms are shown by thin vertical lines.

After vessel depressurizing, the initial wave of decompression is formed and moves to the closed channel end with the speed of sound in pure water. It converts the fluid into a metastable state. On the pressure decreasing below saturation pressure, boiling starts at the wall centres of nucleation. The growth of the $\unicode[STIX]{x1D6FC}$ -bubbles compensates the pressure decrease in the initial decompression wave. As a result, the amplitude of the wave steadily decreases while it propagates towards the closed tube end (figure 16 a).

When the $\unicode[STIX]{x1D6FC}$ -bubbles reach a critical size (see the $\mathit{We}_{\unicode[STIX]{x1D6FC}}$ profiles in figure 16 d), they begin to break up. The fragments of wall $\unicode[STIX]{x1D6FC}$ -bubbles, after their separation and passing into the flow, compose a $\unicode[STIX]{x1D6FD}$ -vapour (figure 16 b). Thus, at the cross-section of the initial decompression wave termination (B–B), there are a significant number of $\unicode[STIX]{x1D6FD}$ -bubbles, $c\approx 10^{5}~\text{kg}^{-1}$ , but since the bubbles are of ${\approx}20$ micron radius (figure 16 b), the volumetric vapour content is negligibly small, $\unicode[STIX]{x1D6FD}\approx 10^{-6}$ .

Along the compression wave, $\mathit{We}_{\unicode[STIX]{x1D6FC}}$ stays higher than the critical value, $\mathit{We}^{\ast }$ , $\mathit{We}_{\unicode[STIX]{x1D6FD}}$ does not increase to $\mathit{We}^{\ast }$ , and therefore $\unicode[STIX]{x1D6FD}$ -bubbles do not fragment. However, due to $\unicode[STIX]{x1D6FC}$ -bubble breakup, the number of $\unicode[STIX]{x1D6FD}$ -bubbles increases from $c=10^{5}$ to $10^{7}~\text{kg}^{-1}$ , the volumetric $\unicode[STIX]{x1D6FD}$ -vapour content increases to $\unicode[STIX]{x1D6FD}\approx 2\times 10^{-4}$ (figure 16 b) and the mixture density, $\unicode[STIX]{x1D70C}$ , decreases (figure 16 c). The liquid overheating decreases in the compression wave to ${\approx}15$  K.

Behind the cross-section C–C, the pressure decrease is negligibly small. This corresponds to the experimental results showing that the compression wave is followed by a ‘pressure plateau’.

The pressure plateau is terminated right next to the discharge vessel end where the pressure drops by $\unicode[STIX]{x0394}P\approx 3$  MPa (figure 16 a). This essential pressure drop causes the flow acceleration and the increase in the phase slip. As a result, the Weber number reaches its critical value for in-volume $\unicode[STIX]{x1D6FD}$ -bubbles (figure 16 d) and their chain breakup starts. The number of boiling centres increases to $10^{14}~\text{kg}^{-1}$ (figure 16 b), the liquid boils up in a shock mode (the vapour fraction reaches maximal values $\unicode[STIX]{x1D6FD}\sim 1$ ) and the phase temperatures equalize (figure 16 c). This structure with 3 MPa pressure drop is a ‘slow’ wave of boiling which, later on, starts to move upstream at a speed of ${\approx}10~\text{m}~\text{s}^{-1}$ (see § 7.3 and figure 15 c,d).

8.2 The compression wave structure in calculations with stationary equations

The configuration consisting of two waves, an initial decompression wave and a compression one, moves deep into the channel with the velocity close to the speed of sound in a pure liquid, $a_{l0}$ (figure 16 a). Let us go over to the coordinate system ( $t$ , $y$ ) moving with the speed $a_{l0}$ ( $y=x+a_{l0}t$ ). In § 5.9, it was shown that the system (5.4), (5.5), (5.13), (5.17), (5.25), (5.26), (5.35), (5.38), (5.44) had only real characteristics and could be rewritten in a characteristic form. For the moving frame, the form is

where $I_{i}=I_{i}$ ( $P,\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\ldots$ ) are the invariants, the expressions that are preserved along characteristic curves if $Q_{i}=0$ ; $W_{i}$ is the velocity of $i$ -invariant propagation (the tangent of inclination of $i$ -characteristic in the ( $t,y$ ) plane); $Q_{i}$ are the source terms.

Let us estimate $W_{i}$ . The equations for characteristics in a fixed coordinate system ( $t,x$ ) were found in § 5.9. Nine of them are mass ones and two are sonic:

The characteristics in the ( $t,y$ ) plane are shown schematically in figure 17(a). The mass characteristics 1–9 pass through the field between cross-sections A–A and C–C with the velocity

The 10th one crosses the field with the speed

Figure 17.

(a) Characteristics in the fields of an initial wave of decompression and a compression wave. (b) Pressure profiles along a 100 m tube.

In the field of the initial decompression wave between cross-sections A–A and B–B (figure 17 a), the 11th invariant propagation velocity is close to zero, $W_{11}\approx 0$ , as the wave velocity is equal to the speed of sound in a pure liquid before the boiling up, $|(\text{d}x/\text{d}t)_{11}|=-a_{l0}$ , and remains close to $a_{l0}$ after the boiling begins.

For the field of the compression wave between cross-sections B–B and C–C, the propagation velocity of the 11th invariant ( $W_{11}$ ) is (figure 16 d)

Thus, the 11th characteristic is the ‘slowest’ one (in the moving frame), and the time interval within which it travels through the compression wave field determines the time of a stationary regime onset, $\unicode[STIX]{x1D70F}^{\ast }\sim l^{\ast }/W_{11}=10^{-2}~\text{s}$ , where $l^{\ast }\approx 1$  m is the characteristic length of the compression wave.

The time of the parameter variation can be estimated at the compression wave boundaries as the time within which the pressures at cross-sections B–B and C–C change by the value of the wave amplitude, $\unicode[STIX]{x0394}P_{w}\sim 1$  MPa. In figure 17(b), the pressure profiles in a 100 m long tube are shown. It is seen that the pressure change takes place within the time interval $\unicode[STIX]{x0394}T\sim 0.1$  s. Since $\unicode[STIX]{x1D70F}^{\ast }\ll \unicode[STIX]{x0394}T$ , we can use the stationary analogue of the model to clear up the wave mechanism.

We are looking for a self-similar solution of a step type moving with the velocity $D$ being constant relative to a fixed coordinate system:

Replacing the coordinate $x$ by $z$ in a stationary version of the model proposed in (5.4), (5.5), (5.13), (5.17), (5.25), (5.26), (5.35), (5.38) and (5.44), and neglecting insignificant flow friction on the channel walls, we obtain

where $v=u+D$ is the flow velocity in a moving frame; $\unicode[STIX]{x0394}v=u_{g}-u$ is the phase slip. (The classic analysis of shock waves in media with energy release is based on the conservation of mass, impulse and energy across the wave front. By eliminating the velocity from the conservation laws, the Hugoniot curve with energy release is obtained. This analysis for a one-velocity non-equilibrium model of boiling liquid was carried out in our previous works (Ivashnyov et al. Reference Ivashnyov, Ivashneva and Smirnov2000; Ivashnyov & Ivashneva Reference Ivashnyov and Ivashneva2012). Because of the peculiarities of the given model in accounting for wall surface nucleation, it is not possible to carry out the classic analysis of shock waves and use its standard conclusions.)

At point $z=0$ , all of the parameters are taken from the numerical calculation at cross-section B–B at instant $t=3$  ms (figure 16): $P^{\ast }=4.171$  MPa, $T^{\ast }=T_{0}=558$  K, $u^{\ast }=(P_{0}-P^{\ast })/\unicode[STIX]{x1D6FD}_{l}=8.5~\text{m}~\text{s}^{-1}$ , $\unicode[STIX]{x1D6FC}^{\ast }=3.20\times 10^{-6}$ , $\unicode[STIX]{x1D6FD}^{\ast }=5.64\times 10^{-6}$ , $c^{\ast }=2.40\times 10^{5}~\text{kg}^{-1}$ , $\unicode[STIX]{x0394}v^{\ast }=-0.824~\text{m}~\text{s}^{-1}$ . (The phase slip is negative since $\unicode[STIX]{x1D6FC}$ -bubble fragments do not have time to speed up to the flow velocity.)

The wave velocity $D$ is varied to find step-like solutions satisfying the condition behind the wave,

Figure 18 shows the solutions obtained under different $D$ values. The wave velocity $D$ is given in terms of the initial sound velocity in liquid, $a_{l0}=985.5~\text{m}~\text{s}^{-1}$ . Looking for the waves of a compression type, we have limited ourselves to analysis of the regimes with a positive pressure gradient at point $z=0$ , $\text{d}P/\text{d}z>0$ . Among the solutions obtained there are continuous, $\text{OA}_{1}$ , $\text{OA}_{2}$ , $\text{OA}_{3}$ , OBD, and broken-off ones, $\text{OC}_{1}$ , $\text{OC}_{2}$ , $\text{OC}_{3}$ (figure 18 a).

The broken-off solutions $\text{OC}_{1}$ , $\text{OC}_{2}$ , $\text{OC}_{3}$ are realized under wave velocities $D<0.9996848112\cdot a_{l0}$ . The reason for the break is zeroing of the denominator (8.12) on the right-hand side of (8.10) at points $\text{C}_{1}$ , $\text{C}_{2}$ , $\text{C}_{3}$ . The denominator is proportional to the product of the velocities of the 10th and 11th invariant propagation in the moving frame. Writing the denominator in terms of characteristic velocities in fixed coordinates, we obtain

Since $(\text{d}x/\text{d}t)_{10}\approx a_{l0}$ , $(\text{d}x/\text{d}t)_{11}\approx -a_{l0}$ and $D\sim a_{l0}>0$ , the first multiplier always differs from zero, $|(\text{d}x/\text{d}t)_{10}+D|\approx 2a_{l0}$ . However, the second one is close to zero,

From figure 18(a,d) it is seen that, at point $\text{C}_{3}$ , where $(\text{d}x/\text{d}t)_{11}=-D$ (figure 18 d), the pressure gradient tends to infinity and the solution breaks off.

Under wave velocities $D\geqslant 0.9996848112\cdot a_{l0}$ , continuous regimes OBD, $\text{OA}_{1}$ , $\text{OA}_{2}$ , $\text{OA}_{3}$ are realized. All of the regimes satisfy condition (8.20). Under $D=0.9996848112\cdot a_{l0}$ , a limiting regime OBD is realized. In this regime, there is also the point at which $\text{Den}$ (8.12) becomes equal to zero (point B), but since $\text{Num}$ (8.11) is also equal to zero at this point, the solution does not break off.

Figure 18.

The structure of the compression wave calculated with the stationary equations (4.67)–(4.75) under different values of the wave propagation velocity $D$ (shown in a). (a) Pressure profiles. (d,e) The velocity of the 11th invariant propagation in the regimes (d)  $\text{OC}_{3}$ , OBD, $\text{OA}_{3}$ and (e)  $\text{OA}_{2}$ . (b,c,f,g) Distributions of the other parameters in the regimes (b,c) OBD and (f,g)  $\text{OA}_{2}$ .

Let us analyse $\text{Num}$ (8.11). Its first two items are determined by the intensity of liquid evaporation into the bubbles, $J_{l\unicode[STIX]{x1D6FC}}$ and $J_{l\unicode[STIX]{x1D6FD}}$ , and are always positive:

The third one, $\text{Num}_{3}$ , is proportional to the product of $\unicode[STIX]{x1D6FC}$ -bubble fragmentation intensity and the velocity of the $\unicode[STIX]{x1D6FD}$ -phase, $J_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}u_{g}$ , which is the flux of impulse wasted to speed up initially motionless fragments of $\unicode[STIX]{x1D6FC}$ -bubbles to the velocity $u_{g}$ . It is always negative:

All of the solutions are waves of compression and, along the wavelength, the change of pressure is positive, $\text{d}P/\text{d}z>0$ , the flow slows down, $\text{d}v/\text{d}z<0$ (8.14), and the Weber number for $\unicode[STIX]{x1D6FC}$ -bubbles, $\mathit{We}_{\unicode[STIX]{x1D6FC}}$ (5.48), decreases. At the OBD, $\mathit{We}_{\unicode[STIX]{x1D6FC}}$ drops below a critical value (figure 18 b). The fragmentation of $\unicode[STIX]{x1D6FC}$ -bubbles comes to an end, and $J_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ and $\text{Num}_{3}$ become equal to zero. This leads to the increase in the pressure gradient and the flow deceleration.

According to (8.15), under the deceleration ( $\text{d}v/\text{d}z<0$ ), the difference in the phase velocities decreases, $\text{d}(\unicode[STIX]{x0394}v)/\text{d}z<0$ . Initially negative $\unicode[STIX]{x0394}v$ decreases further, increasing in absolute value, which causes the increase in the Weber number for $\unicode[STIX]{x1D6FD}$ -bubbles (figure 18 b). Conditions for $\unicode[STIX]{x1D6FD}$ -bubble fragmenting arise. Due to the increase in the interfacial area caused by the fragmentation, boiling intensifies. Since the heat flux into a $\unicode[STIX]{x1D6FD}$ -phase is proportional to the $\unicode[STIX]{x1D6FD}$ -bubble number, $J_{l\unicode[STIX]{x1D6FD}}\sim n$ , the second positive component of the numerator ( $\text{Num}_{2}$ ) increases along the wave and causes the increase in the pressure gradient (8.10). The increase in $\text{d}P/\text{d}z$ causes the decrease in $\text{d}v/\text{d}z$ and $\text{d}(\unicode[STIX]{x0394}v)/\text{d}z$ (8.14), (8.15). The absolute value of the negative difference in the phase velocities $|\unicode[STIX]{x0394}v|$ increases. The increase in $|\unicode[STIX]{x0394}v|$ is why the Weber number for $\unicode[STIX]{x1D6FD}$ -bubbles (5.48) does not decrease after breakup in spite of the diminishing of bubble radius, and so breakup repeats. Thus, one breakup creates the conditions for the next one, i.e. a chain fragmentation proceeds. The bubble fragmentation comes to an end at point ‘B’ where $J_{l\unicode[STIX]{x1D6FD}}$ and $\text{Num}$ (8.11) on the right-hand side of (8.10) tend to zero due to the mixture approaching an equilibrium state ( $T_{g}\rightarrow T_{l}$ , figure 18 c), and so $\text{d}P/\text{d}z\rightarrow 0$ (8.10). Along the segment BD, under maintenance of heat equilibrium, equilibrium in the phase velocities is established (figure 18 b).

Thus, the regime OBD is characterized by the greatest, among the continuous solutions, pressure gradient and the smallest width of the shock wave, ${\sim}1$  cm. This is the only regime in which $\unicode[STIX]{x1D6FD}$ -bubble fragmentation takes place.

Let us consider one of the continuous regimes of the OA series, regime $\text{OA}_{2}$ , in which the amplitude of the compression wave is equal to that obtained in calculations with a full model for the moment $t=3$  ms (figure 16 a). The regime is obtained under the wave velocity $D=1.1\cdot a_{l0}$ , which is greater than the velocity of the 11th invariant propagation along the wavelength; that is, the regime $\text{OA}_{2}$ is supersonic ( $\text{Den}>0$ (8.22)). Compared with the OBD, the pressure and velocity gradients in the regime $\text{OA}_{2}$ are smaller (figure 18 a,b,f). The flow velocity remains great enough for the Weber number of $\unicode[STIX]{x1D6FC}$ -bubbles to be higher than its critical value along the wavelength, $\mathit{We}_{\unicode[STIX]{x1D6FC}}>We^{\ast }$ (figure 18 f). As a result, the third, negative, term of the numerator, $\text{Num}_{3}$ , is ‘preserved’ proportional to the flux of impulse wasted to accelerate the $\unicode[STIX]{x1D6FC}$ -bubble fragments up to the flow velocity, $J_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}u_{g}$ . The presence of the third, negative, term in $\text{Num}$ (8.25) is the reason why the pressure reaches a plateau level, $\text{d}P/\text{d}z=0$ , earlier than the mixture attains equilibrium, when the first two terms of the numerator (8.23), (8.24) become equal to zero, $\text{Num}_{1}=\text{Num}_{2}=0$ .

A small pressure gradient results in a small difference in the phase velocities $\unicode[STIX]{x0394}u$ , which does not permit $\unicode[STIX]{x1D6FD}$ -bubbles to fragment ( $\mathit{We}_{\unicode[STIX]{x1D6FD}}<We^{\ast }$ ). Their number increases only due to $\unicode[STIX]{x1D6FC}$ -bubble fragments entering into the flow (figure 18 f).

Thus, the regime $\text{OA}_{2}$ , compared with the OBD, is characterized by lower pressure gradients and a greater wavelength, ${\sim}1$  m. In this regime, only $\unicode[STIX]{x1D6FD}$ -bubble breakup proceeds.

The regime $\text{OA}_{3}$ , with the velocity of the compression wave equal to the velocity of the initial wave of decompression, $D=a_{l0}$ , is a limiting one for series OA. The regime $\text{OA}_{3}$ is analogous to the OBD: the breakup of $\unicode[STIX]{x1D6FC}$ -bubbles stops at some point and the fragmenting of $\unicode[STIX]{x1D6FD}$ -bubbles starts, which proceeds in a chain mode with the pressure going over on a plateau level close to $P_{s}$ ( $T_{0}$ ).

The comparison of the experiment and calculations with a stationary model shows that the experimental compression wave profile corresponds to the $\text{OA}_{2}$ regime: the wavelength is approximately 1 m and the plateau level behind the wave is ${\approx}1$  MPa lower that $P_{s}$ ( $T_{0}$ ). This can evidence that, in the experiment, the pressure rise behind the initial decompression wave is caused only by wall bubbles fragmenting.

8.3 Compression wave regimes in calculations with a full model

The pressure profile is close to the regime $\text{OA}_{1}$ at instant $t=1$  ms, and to $\text{OA}_{2}$ at $t=3$  ms. As is seen from the wave structure at $t=3$  ms, the wavelength is ${\approx}1$  m (figures 16 a, 18 a) and the fragmentation of only $\unicode[STIX]{x1D6FC}$ -bubbles takes place (see the $\mathit{We}_{\unicode[STIX]{x1D6FC}}$ and $\mathit{We}_{\unicode[STIX]{x1D6FD}}$ profiles in figures 16 c, 18 f).

Each of the regimes OA corresponds to its own velocity of the compression wave $D$ which is greater than the velocity of the initial wave of decompression $a_{l0}$ . Under the transfer from $\text{OA}_{1}$ to $\text{OA}_{2}$ , etc., the wave velocity decreases and the pressure level behind its front increases. Since the compression wave cannot overrun the initial decompression wave, the increase in $D$ is limited to its rise up to $a_{l0}$ . In the limiting regime $\text{OA}_{3}$ with $D=a_{l0}$ , the pressure in the compression wave increases up to a level close to $P_{s}(T_{0})$ (figure 18). The regime $\text{OA}_{3}$ was not obtained in calculations of a 4 m length tube depressurization with a full model. To obtain it we had to lengthen the tube up to 100 m. In figure 17(b), it is seen that the pressure behind the compression wave actually verges towards $P_{s}(T_{0})$ in the course of time. (The exact equality to $P_{s}(T_{0})$ cannot be attained due to temperature wastage on evaporation.)

8.4 Modelling of pressure oscillations

The pressure rise practically up to the $P_{s}(T_{0})$ level by the time $t=4$  ms, corresponding to the regime $\text{OA}_{3}$ , was fixed in both the experiments of Edwards & O’Brien (Reference Edwards and O’Brien1970) (see the experimental oscillograms in figure 2(a,b) and the pressure profile at $t=4$  ms in figure 16 shown by crosses). It is possible that an accelerated realization of the $\text{OA}_{3}$ regime was provoked by interaction of the compression wave with an obstacle. Aluminium discs for void fraction measurements fitted at the pressure sensor at a cross-section of $x=1.5$  m (from the closed end) might be such an obstacle.

In the experiments of Lienhard et al. (Reference Lienhard, Alamgir and Trela1978), the pressure oscillation with the pressure rise up to an initial saturation value (figure 2 c) had been fixed in only one from the experiments presented. Since it was recorded under the use of a heavy plug weighing 18 g (instead of a 3 g plug), it could be caused by the conditions at the exit cross-section.

In all of the regimes of OA type, including $\text{OA}_{3}$ with the pressure rise up to $P_{s}(T_{0})$ , the wave velocity $D$ exceeds the velocity of the perturbations going upstream, $D>|(\text{d}x/\text{d}t)_{11}|$ ; that is, all the regimes are supersonic ones (figure 18 d,e). Therefore, the parameters behind the compression wave have no influence on the wave, and its evolution is only determined by the flow parameters ahead of the wave. Since we could influence the flow only via boundary conditions, to model the $\text{OA}_{3}$ regime, we specified a short pressure perturbation at the discharge tube end before the moment of compression wave formation, at $t=0.55$  ms. The perturbation amplitude was taken to be equal to the amplitude of the compression wave in the regime $\text{OA}_{3}$ :

The results of the problem solution with boundary conditions (8.26) are presented in figure 4(b). The short impulse provokes the formation of a compression wave corresponding to regime $\text{OA}_{3}$ (figure 18). The pressure in it increases to a pressure close to $P_{s}(T_{0})$ and the bubble number increases to $10^{6}~\text{kg}^{-1}$ .

Reflection. On approaching the closed tube end, the initial wave of depressurization is reflected and a return wave arises. The incident wave accelerates the flow from zero velocity to

where $\unicode[STIX]{x1D6FD}_{l}=\unicode[STIX]{x1D70C}_{l}a_{l}$ is the coefficient of isothermal liquid compressibility and $\unicode[STIX]{x0394}P^{\ast }$ is the amplitude of the incident wave.

The reflected wave retards the flow from velocity $u^{\ast }$ to zero,

where $\unicode[STIX]{x1D6FD}_{mix}=\unicode[STIX]{x1D70C}_{mix}a_{mix}$ is the effective compressibility of boiling liquid; $a_{mix}$ is the effective sound velocity in the two-phase mixture; $\unicode[STIX]{x0394}P_{ref}$ is the amplitude of the reflected wave.

Equating the right-hand sides of (8.27) and (8.28) yields

where $k_{mix}$ is the coefficient of the boiling liquid compressibility.

Therefore, the boiling liquid compressibility can be judged by the amplitude of the reflected wave. According to calculations using the model of boiling at a constant number of centres nucleated in the flow bulk, $k_{mix}\approx 0.2$ (figure 4 a). (For the cross-section at 1 m distance from the closed end, this evaluation is fairly rough since it does not consider the damping of the reflected wave while it propagates through boiling liquid. This effect is studied in Ivashnyov et al. (Reference Ivashnyov, Ivashneva and Smirnov2000).) In calculations with the model accounting for nucleation in the bulk of the liquid and bubble breakup, the amplitude of the reflected wave, together with $k_{mix}$ , is close to zero. This is explained by the fact that, in the mixture behind the compression wave, the bubble number is $10^{12}~\text{kg}^{-1}$ (figure 5 b), which is sufficient to maintain the state of equilibrium boiling. An equilibrium mixture is characterized by a low speed of sound, $a_{e}\sim 10{-}100~\text{m}~\text{s}^{-1}$ , and the coefficient of its compressibility is negligibly small, $k_{mix}\sim a_{e}/a_{l}\sim 0.01{-}0.1$ . In calculations with the model accounting for wall nucleation and the bubble breakup, $k_{mix}\approx 0.3$ , this corresponds to the experiment. This indicates that the model proposed simulates the interfacial area formation in waves adequately.

On the arrival of the rarefaction wave reflected from the closed tube end (the moment $t=5$  ms in figure 4 b), the liquid returns to a metastable state and the mixture becomes non-equilibrium again. This reversion is bound with the fact that, after the pressure decrease, the interfacial area is insufficient to maintain equilibrium boiling. It is the composition of the wave of compression and the reflected rarefaction wave that forms ‘a pressure oscillation’, i.e. the pressure rise to $P_{s}(T_{0})$ and its subsequent drop (figure 2 a,b).