### 5.13.4. The Thermal Phase Change Model

This model describes phase change induced by interphase heat transfer; it may be used to simulate boiling and condensation, or melting and solidification. For example, it may be used to model condensation of saturated vapor bubbles in sub-cooled liquid, or evaporation of saturated bubbles in superheated liquid.

This section discusses the theory. For modeling information, see Thermal Phase Change Model.

It is essential to consider the heat transfer processes on each side of the phase interface. Hence, the Two Resistance model for interphase heat transfer must be used in conjunction with the Thermal Phase Change model. For details, see The Two Resistance Model.

Recall that, in this case, the sensible heat flux to phase from the interface is:

| (5–248)

and the sensible heat flux to phase from the interface is:

| (5–249)

- and are the phase and phase heat transfer coefficients respectively.

- The interfacial temperature is determined from considerations of thermodynamic equilibrium. Ignoring effects of surface tension on pressure, assume , the saturation temperature.

Note that this is in contrast to the case where there is no mass transfer. In that case, the interfacial temperature is determined from the sensible heat balance .

In the case of interphase mass transfer, the interphase mass transfer is determined from the total heat balance, as follows.

Total heat flux to phase from the interface:

| (5–250)

Total heat flux to phase from the interface:

| (5–251)

- denotes mass flux into phase from phase .

- and represent interfacial values of enthalpy carried into and out of the phases due to phase change, see below for details.

The total heat balance now determines the interphase mass flux:

| (5–252)

The secondary heat flux term must be modified in order to take account of the discontinuity in static enthalpy due to latent heat between the two phases. This is achieved using a modification of the upwind formulation Equation 5–227, due to Prakash [53]. In this formulation, the bulk fluid enthalpy is carried out of the outgoing phase, as in the default upwind formulation. However, the saturation enthalpy is carried into the incoming phase. Thus:

| (5–253)

#### 5.13.4.1. Wall Boiling Model

Wall boiling starts when the wall temperature achieves a temperature that is sufficiently large to initiate the activation of wall nucleation sites. This activation temperature is typically a few degrees above the saturation temperature. However, at this stage, the average temperature of the liquid in the vicinity of the heated wall is still well below the saturation temperature, hence in the sub-cooled boiling regime.

Evaporation starts in the microscopic cavities and crevices, which are always present on the solid surface. Liquid becomes supersaturated locally in these nucleation sites, leading to the growth of vapor bubbles at the sites. The bubbles become detached from the sites when they are sufficiently large that external forces (inertial, gravitational, or turbulent) exceed the surface tension forces that keep them attached to the wall. As the bubbles depart from the wall, they are displaced by superheated liquid in the vicinity of the nucleation sites, after which the nucleation site is free to create another bubble. In regions of the wall not affected by bubble growth, the wall heat transfer to the liquid may be described by single phase convective heat transfer.

It is clear that detailed physics of bubble growth is very complex, and occurs on very small length scales in the vicinity of the wall. It is unrealistic to model the detailed physics in a phase-averaged Eulerian multiphase model. It is also unrealistic to resolve the small length scales with ultra-fine meshes. The model described here is a so-called mechanistic model, which aims to model the important physical sub-processes using engineering correlations. The model is a sub-grid scale model, in the sense that the complex physics is assumed to take place in a vicinity very close to the wall that is smaller than the mesh resolution at the wall.

The first and most well-known model of this kind was formulated by Kurul and Podowski (1991) [159], from the Rensselaer Polytechnic Institute. It is known as the RPI model. In this model, a number of the sub-models of the overall mechanistic model were taken from correlations originally developed for exploitation in one-dimensional thermo-hydraulic simulation methods. This model was implemented in CFX-4, with the near-wall distance taking the place of the centerline wall distance in the one-dimensional models. Unfortunately, this lead to results that were strongly mesh-dependent.

In the implementation in CFX-5.7.1, Egorov et al. formulated modifications of the one-dimensional correlations with the aim of restoring mesh independence to the results. This modification of the original RPI model is described here.

The following topics are discussed in this section:

- Partitioning of the Wall Heat Flux

- Sub-models for the Wall Boiling Model

For more information about the wall boiling model and its usage, see Wall Boiling Model in the CFX-Solver Modeling Guide.

##### 5.13.4.1.1. Partitioning of the Wall Heat Flux

A fundamental feature of the mechanistic model of wall nucleation is the algorithm for deciding how the wall heat flux is to be partitioned among the separated physical processes of evaporation and sensible heating of the liquid phase. In regions of the wall not influenced by nucleation sites, it is sufficient to consider the wall heat flux as contributing solely to single-phase liquid convective heat transfer. However, in the vicinity of the nucleation sites, some of the heat contributes to vapor production, and the remainder to super heating of the liquid phase as it displaces the rising bubbles. This latter process is known as quenching.

Hence, the mechanistic heat partitioning model has the following general structure:

where is consumed for heating the sub-cooled liquid, and is the evaporation heat flux. is further assumed to be partitioned into two sub-processes:

where is the heat flux corresponding to convective heat transfer, and is the heat flux corresponding to quenching. Hence, in total, there is a three-way partition:

In this model, vapor is assumed to be saturated everywhere, and no part of the wall heat flux is arranged for superheating of the vapor phase.

The heat partitioning model considers the whole wall surface as being separated into two fractions:

- Fraction is influenced by the vapor bubbles, formed on the wall.

- Fraction is the rest of the wall surface, with .

The wall area fraction represents the part of the wall surface that does not 'feel' the presence of the vapor phase. The wall heat flux for this part of the surface is modeled in a similar way as for the single-phase flow of pure liquid, by using the turbulent wall function procedure in the case of the turbulent flow. This part of , called the convective heat flux , can be correlated for turbulent flow as:

where is the temperature of the solid wall, is the temperature of the liquid at the wall, and is the turbulent heat transfer coefficient, which depends on the velocity field and on the near-wall grid cell size. Kurul and Podowski (1991) [159] modeled using a one-dimensional Stanton number correlation. Egorov and Menter (2004) [163] modeled using the turbulent wall function.

The wall area fraction represents the remaining part of the wall surface, which exchanges heat with both phases. The already mentioned heat flux comes from this part of the surface and is consumed for evaporation of the initially sub-cooled liquid:

| (5–254)

where is the evaporation mass transfer rate per unit wall area, and and are the specific enthalpies of the saturated vapor and sub-cooled liquid, respectively.

A part of the heat flux to the liquid, coming from the wall area fraction , is transported between the bubble departure and the next bubble formation at the same nucleation site. This additional mechanism of heating the liquid phase is called quenching, and is modeled as:

The area fraction values and play an important role in the heat-partitioning model. They are related to the nucleation site density per unit wall area and to the influence area of a single bubble forming at the wall nucleation site. The latter value is modeled by introducing the bubble departure diameter value , which can generally serve as a length scale of the wall boiling mechanism. The RPI model assumes that the diameter of the bubble influence zone is twice as large as , therefore the non-dimensional area fraction of the bubble influence is:

| (5–255)

with defined later in Equation 5–258.

Here the implicitly assumed factor of 2 for the diameter of the bubble influence can be user-specified.

It should be noted that, for the evaporation rate and corresponding heat flux (Equation 5–254), the upper limit of Equation 5–255 for bubble influence area fraction does not apply. Instead of using a factor of as defined by Equation 5–255, the value is directly modeled as being proportional to the nucleation site density .

##### 5.13.4.1.2. Sub-models for the Wall Boiling Model

As mentioned in the previous section, the two most important parameters governing the heat partitioning model are the nucleation site density and the bubble departure diameter . In the RPI wall-boiling model they are correlated to the wall superheat and to the near-wall liquid sub-cooling , respectively.

The following topics will be discussed:

- Wall Nucleation Site Density

- Bubble Departure Diameter

- Bubble Detachment Frequency

- Bubble Waiting Time

- Area Influence Factors

- Convective Heat Transfer

- Quenching Heat Transfer

- Evaporation Rate

###### 5.13.4.1.2.1. Wall Nucleation Site Density

The model for wall nucleation site density adopted in the RPI model is that of Lemmert and Chawla (1977) [164]:

Note that the wall superheat in the above equation cannot be negative. Its negative value means that the wall temperature drops below the saturation temperature, where there is no boiling and the heat partitioning model is not used.

This model was implemented in CFX-4 with parameter , giving an overall correction factor of . Egorov and Menter (2004) [163], conjecturing that this was a deliberate alteration by the RPI group, related to the corresponding factor in the bubble waiting time model. Egorov and Menter also reformulated the correlation as follows:

where and . This formulation avoids fractional powers of physical dimensions, hence is more amenable to the use of different unit systems.

Alternative correlations have been proposed for nucleation site density by Kocamustafaogullari and Ishii (1983) [160].

###### 5.13.4.1.2.2. Bubble Departure Diameter

Kurul and Podowski (1991) [159] adopted the correlation for bubble departure diameter due to Tolubinski and Kostanchuk (1970) [165]:

| (5–256)

The parameters of the model are dimensional (, , ) and are chosen to fit pressurized water data. Hence the model is clearly not universal. Negative liquid sub-cooling is possible here, where it means onset of the bulk boiling. Limiting the bubble departure diameter applies to this situation and prevents from growing too high.

Note that the model is strongly dependent on a liquid temperature scale. In the original experimental data, this was taken as the pipe center-line temperature. CFX-4 used the cell-center value in near wall cells, but this proved to give mesh-dependent results. Egorov and Menter restored mesh independence by using the logarithmic form of the wall function to estimate the liquid temperature, , at a fixed value of . ANSYS CFX uses this method, with a user-modifiable value of .

###### 5.13.4.1.2.3. Bubble Detachment Frequency

The computation of the evaporation rate requires an additional model parameter, namely the frequency of the bubble detachment from the nucleation site. The model adopted by Kurul and Podowski is that due to Cole (1960) [166]:

Note that, due to its dependence on gravity, this correlation is taken from pool boiling. It is simply estimated as the bubble rise velocity divided by the bubble departure diameter. The drag coefficient factor was taken to be unity by Ceumern-Lindenstjerna (1977) [169].

###### 5.13.4.1.2.4. Bubble Waiting Time

Kurul and Podowski (1991) [159] employed the model of Tolubinski and Kostanchuk (1970) [165]. This fixes the waiting time between departures of consecutive bubbles at 80% of the bubble detachment period:

| (5–257)

The numerator in this equation is adjustable in ANSYS CFX.

###### 5.13.4.1.2.5. Area Influence Factors

Recall from Partitioning of the Wall Heat Flux that Kurul and Podowski assumed a diameter of influence of a nucleating bubble equal to twice the bubble departure diameter . Encoding this as a user-modifiable parameter (default value = 2), the area fraction of the bubble influence is given by:

with defined later in Equation 5–258.

The area fraction subjected to single phase liquid convective heat transfer is limited from below by a small value, so its actual form is:

###### 5.13.4.1.2.6. Convective Heat Transfer

As discussed in Partitioning of the Wall Heat Flux, single phase convective heat transfer to the liquid phase is modeled using the turbulent wall function (Egorov and Menter 2004 [163]).

This replaces the mesh dependent Stanton Number correlation originally employed by Kurul and Podowski (1991) [159].

###### 5.13.4.1.2.7. Quenching Heat Transfer

As discussed in Partitioning of the Wall Heat Flux, quenching heat transfer to the liquid phase in the area of influence of the vapor phase is modeled using a quenching heat transfer coefficient:

In order to close the model, the quenching heat transfer coefficient , participating in the quenching heat flux to liquid, must be defined. This value depends on the waiting time between the bubble departure and the next bubble formation.

With this value, the quenching heat transfer coefficient is correlated as:

where is the liquid temperature conductivity coefficient. (Mikic and Rohsenow 1969 [167], Del Valle and Kenning 1985 [168]).

As for the case of bubble departure diameter, Egorov and Menter (2004) [163] used the logarithmic form of the wall function to estimate the liquid temperature at a fixed value of 250.

###### 5.13.4.1.2.8. Evaporation Rate

Knowing the bubble departure frequency, as well as the bubble size and the nucleation site density, one can obtain the evaporation rate as a product of the bubble mass, the detachment frequency and the site density:

This was the form adopted by Kurul and Podowski (1991) [159].

Egorov and Menter expressed the evaporation rate in terms the non-limited area fraction :

(unlike , can exceed 1). In this case, the evaporation rate obtains the form:

which is an extension of the Kurul and Podowski expression and coincides with it when the Parameter is equal to 2 (default value).

In the final form of the evaporation rate, the area fraction factor is limited by:

| (5–258)

defaults to 0.5.

##### 5.13.4.1.3. Determination of the Wall Heat Flux Partition

The wall heat flux partition discussed in Partitioning of the Wall Heat Flux, together with the user’s choice of submodels described in Sub-models for the Wall Boiling Model, gives the total heat flux as a complex nonlinear function of the wall superheat

In the case of an isothermal wall, the three components of the heat flux are computed explicitly as functions of the given wall superheat. However, in the case of a wall with specified total heat flux, or specified wall heat transfer coefficient and outside temperature, the above equation determines implicitly as a function of the total heat supplied to the wall. In this case, is computed by solving the nonlinear equation using the binary chop method. The resulting value of is subsequently used to determine the three individual components of the wall heat partition.

Important: If you want to use the User Specified option for any of the sub-models described in Sub-models for the Wall Boiling Model, and if you want any of these to be functions of the wall superheat, then it is important to use the correct name for this variable. As the wall boiling model is described as a phase-pair model, this variable is a phase-pair based variable named Tsuperheat. For example:

= water | steam.Tsuperheat

It is not appropriate to use a phase-specific temperature field at the wall, such as water.Temperature, to determine as this variable does not participate in the binary chop algorithm to determine the wall superheat.

Release 18.2 - © ANSYS, Inc. All rights reserved.