Revisiting the Schrage Equation for Kinetically Limited Evaporation and Condensation

The Schrage equation is widely used by engineers to estimate heat/mass transfer rates of evaporation and condensation in the kinetically limited regime. This high-rate regime is defined by the dominance of the kinetic resistance at the liquid–vapor interface, which can occur when diffusional mass transfer resistances are absent. For example, the diffusion of vapor through the air is rate-limiting, but this resistance can be minimized in a hermetically sealed device. Thus, the Schrage equation holds great importance for modeling and design of systems that operate at high rates of evaporation/condensation, as found in applications such as power generation, refrigeration, electronics cooling, and distillation.

Schrage published his dissertation in 1953 [1], formulating what we refer to as the “Schrage equation” to correct the description of the half-space problem of evaporation and condensation previously discussed by Hertz [2] and Knudsen [3]. Specifically, Schrage accounted for the nonzero macroscopic velocity of the vapor phase resulting from the net evaporation/condensation flux. A similar treatment was given by Kucherov and Rikenglas [4], but here we focus on Schrage's formula. In 1991, Barrett and Clement [5] noted that Schrage's formulation violates the conservation of momentum and energy between the liquid–vapor interface and the far-field vapor. In fact, the seldom-mentioned fifth chapter of Schrage's thesis acknowledges this failure and seeks to rectify the issue of conservation. Schrage developed a solution based on conservation of mass, momentum, and energy very similar to the “moment solutions” later applied by Anisimov [6] and Ytrehus [7], although their choice of ansatz was better-informed than Schrage's.

A few recent works have sought to validate the Schrage equation via molecular dynamics (MD) simulations of the classic two-surface problem of evaporation/condensation between two parallel liquid surfaces, typically spaced ∼10 mean free paths apart [8–10]. Yu and Wang [11] and Liang et al. [12] observed good agreement between the fluxes predicted by the Schrage equation and their MD simulations for evaporation/condensation of argon, as did Bird et al. [13] for n-dodecane. However, Chandra and Keblinski's [14] MD study of water found that the Schrage equation overpredicted heat and mass fluxes by ∼20%. These studies present an interesting quantification of the Schrage equation's applicability to the two-surface problem; however, the two-surface problem [8–10,15–20] is physically distinct from the half-space problem [8,15,21] and is thus treated differently in the kinetic theory literature. Graur et al. [22] recently used the S-model kinetic equation to quantify the inaccuracy of the Schrage equation and demonstrated that kinetic equations can provide similar accuracy to molecular dynamics and direct simulation Monte Carlo methods with significantly less computation time. However, this study was limited to the evaporation of monatomic fluids.

A large body of literature treats the half-space problem via various approximate solutions to the Boltzmann transport equation (BTE), thus explicitly conserving mass, momentum, and energy. These methods are more tedious than the classically cited treatment by Schrage, and the engineering community has largely preferred the Schrage equation, especially the approximate formulation derived by Nabavian and Bromley [23] and given in Carey's 1992 textbook [24]. This approximate formulation is attractive since it expresses the heat/mass flux of liquid–vapor phase change as an explicit function of temperatures and pressures of the condensed and vapor phases. Further, it was derived for the low Mach number flow conditions typically seen in thermofluid engineering devices. Here, we revisit the Schrage equation through the lens of decades of kinetic theory literature and more recent experimental data to provide: (1) quantification of the accuracy of the Schrage equation for modeling evaporation and condensation; (2) an analysis of fundamentally accurate approaches for calculating heat/mass fluxes; and (3) recommendations for higher-accuracy estimations of evaporation/condensation fluxes.

where the subscript “$v$” denotes properties of the equilibrium vapor flow in the far-field and $uv$ is the net velocity of the vapor stream. Commonly, the vapor phase is modeled as a semi-infinite half-space without flow tangential to the liquid surface, such that $uv=uvx̂$⁠, where $x̂$ is the unit vector perpendicular to the liquid surface and pointing into the vapor phase (see Fig. 2).

where $ρe$ is the equilibrium density of the saturated vapor at the liquid temperature $TL$ (see Fig. 1 for illustration). One consequence of assigning Eq. (3) as a boundary condition is nonequilibrium (i.e., a discontinuity) at the liquid–vapor interface when net evaporation/condensation occurs. This loss of equilibrium is evident in that $fL$ can only be an equilibrium distribution function if $fL−(ξ)=fL+(ξ)$ for $ξx$ < 0, which results in zero net velocity. The half-range Maxwellian form of $fL+$ has been corroborated via molecular dynamics simulations of monatomic [25,26] and polyatomic [27] species when $TL$ is well below the critical temperature. Frezzotti and Barbante also observed the half-range Maxwellian from their Enskog–Vlasov model calculations [28]. We note that other MD studies appear to contradict these findings [29–31], but the discrepancy appears to arise only from the choice of interphase boundary position. References [26] and [27] evaluated the distribution function at the center of the transition layer (the angstrom-scale region where the density transitions from liquid to vapor) while Refs. [30] and [31] explicitly chose a “vapor boundary” in the vapor region outside of the transition layer. In fact, Refs. [26] and [27] also evaluated the distribution function at points similar to those chosen by Refs. [30] and [31] and the results match well. Furthermore, nonequilibrium thermodynamic analysis has shown that the temperature of the liquid–vapor interface is very close to the temperature of the bulk liquid adjacent to the interface [32].

The net emitted flux may be smaller than predicted by Eq. (3) as a result of nonidealities at the liquid–vapor interface. To account for this effect, the evaporation coefficient $σ̂e$ is introduced such that the actual emitted distribution is $σ̂efL+$ [21,26,33]. A similar allowance for diffuse reflection of a fraction of the condensing stream is made via the condensation coefficient, $σ̂c$⁠. At equilibrium, Eq. (1) holds everywhere in the vapor, and therefore $σ̂e=σ̂c=σ̂$ [21,24], with the latter often referred to as the “accommodation coefficient.” The common practice extends this equality of evaporation and condensation coefficients to nonequilibrium analyses [21,34]; we will follow this precedent here under the caveat that this assumption is only certain to be accurate near equilibrium, i.e., at low fluxes. The accommodation coefficient is a temperature-dependent property of a given fluid [5], bounded by 0 and 1. We note that the accommodation coefficient has a firm physical meaning, as described above, despite the fact that it is often treated as a fitting parameter for model equations. Much debate exists around the values of the accommodation coefficient for particular fluids (especially water) as it is very difficult to measure and likely highly sensitive to trace amounts of noncondensable gases in the vapor and impurities at the liquid–vapor interface. Accommodation coefficients for various fluids were compiled by Paul [35], and many studies on the accommodation coefficient of water were reviewed by Eames et al. [33]. More recent experimental measurements of the accommodation coefficient of water were briefly discussed by the authors [36].

For polyatomic species, $fL$ must also be integrated over the internal energy space and the energy conservation expression (Eq. (4c)) changes to account for the transport of internal energy; examples can be found in Cercignani's textbook [38] and Barrett and Clement's article [5].

where $pe$ is the equilibrium saturated vapor pressure at $TL$ and $J$ is defined as positive for net evaporation. $Γ$ is often called a “correction factor” that accounts for the net velocity in the vapor phase and is a function of the speed ratio $Sv=uv/2RTv$⁠, where $2RTv$ is the most probable thermal speed. Note that we have assumed ideal gas behavior to replace mass density with pressure in Eq. (6). More detailed derivations of Eq. (6) can be found in Schrage's thesis [1] and Carey's textbook [24].

For complete accommodation (⁠$σ̂$ = 1), Eq. (8) predicts only half the flux of Eq. (7). Ytrehus [21] noted that Eq. (7) is the more accurate expression for evaporation, except very close to the sonic condition.

Despite the widespread use of the Schrage equation in parts of the thermofluid engineering community [24,39–41], the kinetic theory community has pointed out flaws in Schrage's derivation that lead to an incorrect physical picture [5,21]. The root of the problem is the assumption that the distribution function of the incident stream is exactly the $ξx$ < 0 portion of the shifted Maxwellian $fv$ (Eq. (2)), which is equivalent to assuming $fvξx<0$ describes the incident vapor stream as an infinitesimal distance away from the liquid surface. This condition entirely overlooks the transition from nonequilibrium at the liquid–vapor interface to the state of thermal equilibrium achieved in the far-field vapor. The transition region is known as the Knudsen layer (see Fig. 2), and only beyond this region is the gas described by the equilibrium distribution $fv$ [5,21].

While Schrage's choice of $fL−$ is a reasonable first approximation, it neglects collisions between molecules of the incoming stream from the equilibrium region with molecules from the stream emitted from the liquid–vapor interface. This collision process enables relaxation of the nonequilibrium condition at the interface to the equilibrium state of the far-field vapor. The paradox of Schrage's physical picture is that he prescribes a nonequilibrium condition at the interface using a piece of the far-field equilibrium distribution, yet allows for no mechanism by which relaxation to the equilibrium state can be achieved (in fact, his kinetic models do not incorporate the collision integral from BTE).

As a result of this conceptual shortcoming, Schrage's derivation satisfies mass conservation (Eq. (4a)) but violates momentum and energy conservation at the liquid–vapor interface (Eqs. (4b) and (4c)). Barrett and Clement [5] mathematically demonstrated that the explicit Schrage equation violates energy conservation; however, the lack of explicit enforcement of momentum and energy conservation in the derivation of the Schrage expressions is sufficient to expect this conclusion. In fact, Schrage himself noted that his choice of $fL−$ “cannot possibly be a rigorous description of those molecules with $ξx$ < 0 just at the phase interface” (adjusted to the notation in this work, Chapter V of Ref. [1]). Schrage proceeded to develop a mass flux expression by assuming instead $fL−=(1−Bξx)fL+,ξx<0$⁠, where $fL+,ξx<0$ indicates the same Maxwellian distribution as $fL+$⁠, but for $ξx<0$ instead of $ξx>0$⁠, and $B$ is an arbitrary constant. This approach does allow for the conservation of mass, momentum, and energy, and it resembles the moment methods then in use for rarefied gas dynamics analysis and later applied to evaporation and condensation (see Secs. 2.3 and 4). However, we will not discuss this method further as we have not seen it used elsewhere in the literature and because a better choice of ansatz for $fL−$ was introduced later.

In Secs. 2.3, 3.1, and 4, we will briefly introduce more rigorous analyses of the half-space problem developed in the kinetic theory community and discuss the ramifications of ignoring the Knudsen layer.

where $I(f,f)$ is the binary collision integral which is described in detail in the literature [21,37,42]. The boundary conditions for the half-space problem are $fL+$ and $fv$ (as specified in Sec. 2.1) at the interface and at the equilibrium edge of the Knudsen layer, respectively. The thickness of the Knudsen layer is not known a priori, but analysis of the BTE shows that its thickness scales with the mean free path of the gas phase (⁠$λ$⁠, average distance gas particles travel between collisions). In contrast to Schrage's analysis, $fL−$ is a result of solving the BTE subject to the half-space boundary conditions, rather than a boundary condition to be enforced. The BTE is difficult to solve for practical problems due to the complicated nature of the collision integral [38]; however, decades of theoretical and experimental work have led to a variety of approximate solution methods.

As previously discussed in the literature, the half-space problem for a given vapor species is fully characterized by three parameters: the pressure ratio $pK*$⁠, temperature ratio $TK*$⁠, and Mach number $MK=uK/γRTK$ [21,42,43], where $γ$ is the ratio of specific heats. Focusing our analysis on the temperature range where vibrational excitation can be neglected, $γ$ relates to the number of rotational degrees-of-freedom (⁠$j$⁠) as $γ=(5+j)/(3+j)$⁠. As illustrated in Fig. 2, we use the subscript “$K$” instead of “$v$” to denote parameters evaluated at the equilibrium edge of the Knudsen layer to stress that these parameters are boundary conditions for the continuum region. The standard approach is to treat the Knudsen layer problem in the half-space framework and then couple the Knudsen layer solution to the continuum region. We briefly discuss this coupling in sec. 5. As heat and mass fluxes are of more interest than the Mach number for most engineering applications of liquid/vapor phase change, we will refer to the dimensionless flux as the third parameter, rather than the Mach number (the two are related by the expression $J*=2πγMKpK*/TK*$⁠).

The problem of weak (or low Mach number) evaporation/condensation has been solved using linear perturbation theory [8,44,45]. To describe strong evaporation, Anisimov [6,46] applied a moment method based on the collision-invariant moments of the BTE (i.e., mass, momentum, and energy), adapting methods developed for the study of the shock wave structure. Ytrehus [7] developed this moment method further and with Alvestad [47] extended it to condensation. Shankar and Marble [48] also treated the half-space problem using a similar moment analysis, but their approach only provided a steady solution for condensation. Moment methods require certain assumptions to be made about the incident stream, yet these assumptions are much looser than Schrage's and, importantly, allow the simultaneous conservation of mass, momentum, and energy. The BGK method (developed by Bhatnagar, Gross, and Krook [49] and Welander [50]) is widely used as a starting point to achieve analytical/numerical solutions to the BTE by approximating the collision integral as $I(f,f)→νc(fM−f)$⁠, where $νc$ is the characteristic collision frequency and $fM$ is the local Maxwellian described by the local temperature, density, and macroscopic velocity of the gas. This model is sometimes referred to as the “relaxation time approximation” as 1/$νc$ is a timescale over which $f$ relaxes to the local Maxwellian via collisions.

Several numerical methods have been developed to solve the BGK approximation of the BTE for the half-space problem [51,52]. These numerical models require no additional assumptions about the form of $fL−$⁠, and have thus been used to validate the moment methods, which have the advantage of offering closed-form analytical solutions. Numerical solutions to the BGK equation have been developed to the point of describing monatomic and polyatomic gas species, but they are not capable of accounting for any further properties specific to the species in question. Direct simulation Monte Carlo (DSMC) solutions provide the most precise description of the half-space problem owing to the development of collision models that do not require the relaxation time approximation. Further, these collision models have parameters that can be tuned to replicate the experimentally determined viscosity particular to unique gas species (see Bird [53,54] for a detailed description of the DSMC method and its collision models). DSMC calculations for hard-sphere gases, numerical methods using the relaxation time approximation, and the moment method all provide consistent mass flux calculations for the half-space problem for evaporation [43] and condensation [55].

Ytrehus [21] provided simple guidelines for selecting the appropriate model for the half-space problem according to the magnitude of the dimensionless driving pressure, $Δp*=(pe−pK)/pe$⁠. The weak mass transfer regime is described by $|Δp*|≪$ 1 and is appropriately described by the linear perturbation theory solutions. Strong mass transfer occurs when $|Δp*|$ ∼ 1 and should be modeled using the moment methods, high-order numerical approximations (e.g., numerical solutions to the BGK equation), or DSMC. Problems of moderately strong mass transfer (⁠$|Δp*|$ < 0.1) can alternatively be modeled using a linearized solution to the moment methods that will be discussed in Sec. 4.

The existing kinetic theory literature describes some key differences between the descriptions of evaporation/condensation from BTE analysis and Schrage's equation. Recalling Eq. (9), Schrage's explicit dimensionless flux expression for both evaporation and condensation, we see that Schrage's analysis states that both the dimensionless temperature and pressure must be specified to determine the nondimensional flux. In contrast, the BTE-based methods all agree that only one of these can be a free parameter for evaporation: e.g., $TK*=TK*(pK*)$⁠, $J*=J*(pK*)$⁠, etc., as already discussed. On the other hand, the BTE methods do show that condensation requires two parameters to be specified; this fundamental difference between evaporation and condensation is not captured by the Schrage equation. Furthermore, Ytrehus showed that the explicit Schrage equation (Eq. (7)) overestimates the dimensionless flux in the case of $σ̂$ = 1 [21].

We believe that it would be valuable to the community to provide quantitative examples of the deviations in fluxes calculated using the Schrage equation (both in full and explicit form) from theoretically accurate values for common engineering problems. In this section, we describe a numerical method from the literature for approximately solving the BTE for polyatomic gases and compare fluxes calculated from this method to calculations using the Schrage equation. Much of this analysis is performed for evaporation and condensation of water at low temperatures (0–100 °C), but we also show that these insights apply to other fluids and different operating conditions.

We found that our numerical solver could reproduce the same relationships among pressure ratio, temperature ratio, and Mach number as the original references for evaporation [43] and condensation [55] of polyatomic gases (see Appendix  B4). Further, our results for monatomic gases matched those computed by Sone [15]. However, we note that the formulation is incapable of differentiating between species that have the same number of rotational degrees-of-freedom [43]. Thus, to validate the application of this analysis to the common engineering problem of evaporation/condensation of water, we compared our numerical results to the experimental results and DSMC calculations of Lu et al. [58]. They used the variable soft sphere collision model with parameters that replicate the viscosity of water vapor, meaning their results are specifically accurate for water. Our present model matches their DSMC and experimental results very well (see Appendix  B4), suggesting the hard-sphere model and relaxation time approximation offer sufficient accuracy for real gases.

The simplifications to the collision term of the Holway model result in a kinetic description that distinguishes between species only by their number of rotational degrees-of-freedom. Recalling that evaporation of a pure fluid is fully described by one dimensionless parameter (⁠$pK*$⁠, $TK*$⁠, or $J*$⁠), the Holway model provides a single curve per specified $j$ relating any two parameters to each other. That is, the relation $Jj=0*=J*(pK*)j=0$ is universal for monatomic species and separate relations exist for diatomic (⁠$j$ = 2) and nonlinear polyatomic (⁠$j$ = 3) species. In the case of condensation, there is a two-dimensional parameter space per value of j (i.e., $Jj=0*=J*(pK*,TK*)j=0$⁠). To reduce the parameter space, it is useful to consider the case where $pK*$ and $TK*$ are not independent parameters but rather are constrained by the saturation relation of the fluid, i.e., the vapor is in a saturated state, a widely relevant condition in engineering applications [24,59–61]. When this special case is considered, each fluid has a unique relation (⁠$J*(pK*)$ or $J*(TK*)$⁠) corresponding to its saturation curve. Schrage's analysis of the case of a saturated vapor state provides fluid-specific solutions for both evaporation and condensation, with the former being an erroneous outcome as evaporation should be described by one parameter only.

This analysis of the Schrage equation's accuracy first considers water with $TL$ = 25 °C and a saturated vapor state as a case study and then briefly discusses how these results can be applied to other fluids and operating conditions, including nonsaturated states. Here we show the solutions as functions of the dimensionless driving pressure (⁠$Δp*$⁠) as has been done by others [21,55,58].

where $Q$ is a calculated quantity the subscripts “mod” and H denote calculations from the model in question and the Holway model, respectively. This error definition is used consistently throughout this work and has the property of being positive when $Qmod$ overpredicts the absolute value of the quantity, e.g., the flux error is positive when $Jmod*$ overpredicts the net evaporation or net condensation rate. Consistent with previous studies, the explicit Schrage equation overpredicts the evaporative flux [21]. The full expression (Eq. (6)) overpredicts the evaporative flux to a lesser degree and consistently overpredicts the flux by ∼15% for both evaporation and condensation across the range considered here (Fig. 3(b)). Interestingly, Eq. (6) is less accurate than the explicit expression for condensation. This trend of the explicit expression performing better for condensation and worse for evaporation is a coincidence resulting from the fact that $Γ′$ is always less than $Γ$ (see Appendix  C), artificially shrinking the contribution of the condensing stream and thus resulting in $JEq.(6)*$ < $JEq.(7)*$⁠. Observing that Eq. (6) always overpredicts the absolute value of the flux, this inequality means that the approximation of $Γ′$ can help to offset the inaccuracy of Eq. (6) under net condensation but exaggerates this inaccuracy under net evaporation.

where $(pK*)σ̂=1$ is the pressure ratio corresponding to the model results calculated with complete accommodation. We note that Eq. (13) was derived assuming fully diffuse reflections. Conceptually, Eq. (13) maps the pressure ratio of a real fluid to that of a pseudo-fluid with the property $σ̂=1$ such that $SK(pK*)=SK((pK*)σ̂=1)$ and $TK*(pK*)=TK*((pK*)σ̂=1)$⁠. Therefore, the relations among these parameters need only be calculated for complete accommodation, then the dataset can be transformed to describe an arbitrary value of $σ̂$ by replacing the pressure ratio with Eq. (13) at virtually no computational expense compared to the numerical Holway solver. While this method incurs no error for evaporation, transforming condensation calculations causes $TK*(pK*)$ to deviate from the saturation relation initially enforced in the complete accommodation calculation. However, the condensation flux is only weakly dependent on $TK*$ when holding $pK*$ constant. We calculated the relative error incurred by assigning constant temperature ratios $TK*=Tsat*$(⁠$Δp*$ = 0$) =1$ or $TK*=Tsat*(Δp*$= −0.5$)=$1.023 instead of using the saturation relation to calculate the condensation flux of water at $TL$ = 25 °C from −0.5 ≤$Δp*$ < 0 for complete accommodation. We found the error caused by slightly deviating from the saturation curve to be ≈ 1% at most (see Appendix  D).

Figure 4 compares the explicit Schrage equation predictions to the Holway model calculations for the same conditions as Fig. 3 but with varying values of $σ̂$⁠. The Schrage expression properly captures the trend that decreasing $σ̂$ suppresses the flux of net evaporation and condensation. We observe that the explicit Schrage equation is much more accurate for $σ̂$ < 1 than for full accommodation; this finding agrees well with the observation of Lu et al. [58] that their experimental data for evaporation of water into its pure vapor could be described well by DSMC calculations or the explicit Schrage equation with the same $σ̂$ = 0.31. However, our data show that the relative error does not converge to zero as $σ̂$ → 0, but rather that Eq. (7) can underpredict the flux for sufficiently small values of $σ̂$ (Fig. 4(b)). We found that Eq. (6) shows the same trend with $σ̂$ as Eq. (7).

Despite speculation that the full, nonlinear Schrage equation (Eq. (6)) would not suffer from the same inaccuracy widely attributed to the more popular explicit form (Eq. (7)) [12,14], our results show that Eq. (6) still fails to accurately describe the kinetics of evaporation and condensation. This failure can be attributed to the lack of enforcement of momentum and energy conservation in Schrage's derivation, allowed by neglecting molecular collisions and setting the incident stream to the $ξx$ < 0 portion of the far-field vapor distribution function. With quantification and cause of its inaccuracy, the question remains: under what conditions has Schrage's approximation given reasonable accuracy despite its flaws? The answer ought to lie in the accommodation coefficient. As shown in Fig. 4(b) the error of the explicit Schrage expression (and Eq. (6) as well) is relatively low for fluids with small $σ̂$⁠. Polar fluids like water and alcohols typically fall into this range, and these fluids are perhaps the most well-studied [24,35,62]. Furthermore, $σ̂$ is a difficult property to measure accurately, and the literature reflects much disagreement about its value, particularly for polar fluids like water [33–35]. Thus, $σ̂$ is often treated as a fitting parameter in evaporation/condensation experiments [35,58,63–65]. Because the Schrage expressions overestimate the flux when the true $σ̂$ of the fluid is large (⁠$σ̂$ ≳ 0.5, from Fig. 4), a reasonable value of $σ̂$ (i.e., bounded by 0 and 1) for Eq. (6) or Eq. (7) can always be fitted to experimental data. For instance, we performed nonlinear least squares regressions to fit the Schrage expressions to the same $σ̂$ = 1 curve for the Holway model as shown in Fig. 3(a). Using $σ̂$ = 0.94, the Schrage equation predicts the $σ̂$ = 1 Holway curve with less than 10% error and even better accuracy for low flux conditions (Fig. 5). In practice, $σ̂$ is unlikely to be fit to data spanning such a large range of condensation and evaporation conditions; rather, the experimental data is likely to cover either condensation or evaporation for a limited set of flux conditions [35,58,63–65]. In this case, it is likely that fitting $σ̂$ will reasonably reproduce the data within the bounds of uncertainty. For example, considering the $σ̂$ = 1 curve for the Holway model (Fig. 3(a)) restricted to evaporative fluxes from 0–100 W/cm² (for water at $TL$ = 25 °C) results in less than 0.5% and 2% error from Eqs. (6) and (7), respectively.

While the dimensionless relation between flux and driving pressure for evaporation of a fluid with $j$ = 3 applies to any absolute $TL$⁠, the relation does not strictly hold for condensation for different $TL$ values or different working fluids as the saturation relation varies. We calculated the condensation flux relations for water at $TL$ = 1 °C and $TL$ = 100 °C and found that the relation $J*(Δp*,TK*sat;TL$ = 25 °C) from Fig. 3(a) deviates by less than 0.2% and less than 0.5% from $J*(Δp*,TK*sat;TL$ = 1 °C) and $J*(Δp*,TK*sat;TL$ = 100 °C), respectively (see Appendix  D). Furthermore, fixing either $TK*$ = 1 or $TK*=$ 1.023 for $TL$ = 25 °C, rather than using the saturation curve, results in a maximum error of ∼1%, as discussed in Sec. 3.3. Thus, we infer that the condensation flux of water is only weakly dependent on $TK*$ in this range and saturated conditions for 1 °C ≤ $TL$ ≤ 100 °C and that condensation flux can be accurately modeled using any single calculated flux relation, e.g., the data shown in Fig. 3(a). Recalling the dimensionless form of the explicit Schrage equation (Eq. (9)), we see that the Schrage equations are also only weakly dependent on $TL$ for the same range since $Tv*$ (using the Knudsen-layer ignorant notation) is very close to unity. It follows that the error in flux prediction at other $TL$ in this range should be like that of Fig. 3(b).

Other fluids with $j$ = 3 have condensation flux relations that differ from water only by their saturation relations. We also calculated the flux relation with $TK*$ = 1.1 fixed and observed less than 5% deviation from the curve for saturated water at $TL$ = 25 °C (Appendix  D). This represents a relatively extreme case: for saturated water at $TL$ = 25 °C this corresponds to $Δp*$ = −3.9. As for monatomic (⁠$j$ = 0) and diatomic (⁠$j$ = 2) species, we found that the evaporation and condensation fluxes deviate little from the $j$ = 3 case (Appendix  D). This result agrees well with Frezzotti and Ytrehus [55], who found that $J*$ as a function of $Δp*$ calculations for different values of $j$ collapse to a single curve for condensation at the same $TK*$⁠. We note here that the choice of normalizing by the mass flux emitted from the liquid (⁠$JL+$⁠) is justified by this result: $J*$ is not independent of $j$ if its normalization involves the sonic speed (which depends on $j$ via the ratio of specific heats). Synthesizing these observations about the dependence of the dimensionless flux on $TK*$ and $j$⁠, the results of Fig. 3 provide a reasonable understanding of the general accuracy of the Schrage equations. For example, we calculated the flux relation for argon using $TL$ = 100 K for the condensation curve and compared this result with the Schrage equations (Fig. 6). The error of the Schrage equations for argon is nearly identical to their error for water.

The same trimodal ansatz cannot be applied to the problem of half-space condensation two of the macroscopic flow variables must be specified in that case. Thus, Ytrehus and Alvestad [47] introduced a fourth mode, $a4−(x)f4−(ξx)$⁠, representing the molecules reflected back toward the interface as a result of the first collisions of the evaporated molecules, where $f4−(ξx)$ is a known distribution calculated from the BGK model. $f4−(ξx)$ and the three conservation equations for condensation using four modes are shown in Appendix  E (Eqs. (E8)–(E10)). The fourth mode adds additional unknown to the control volume analysis, $a4−(0)$⁠, thus requiring that two of the macroscopic flow variables be specified to describe the problem. As we are unaware of any extension of the moment method to condensation of polyatomic species [55], our moment method calculations for net condensation do not distinguish species by internal degrees-of-freedom. These same analyses for evaporation and condensation can be applied to fluids with arbitrary $σ̂$ via Ytrehus' transform [21] by substituting all $pK*$ terms in the conservation equations with the $(pK*)σ̂=1$ expression that results from rearranging Eq. (13). Therefore, the flux relation for evaporation/condensation of any fluid can be approximated by numerically solving the corresponding system of equations (at a very small computational cost compared to the Holway model).

where $ω=32π/(32+9π)$ and $ω′=(23π−32)/4π$⁠.

We compare the flux prediction performance of the moment method, Eq. (16) (L. moment method), and the explicit Schrage equation (Eq. (7)) to the Holway model in Fig. 7. Here we considered a monatomic species with $TK*$ = 1 fixed in the net condensation regime. This moment method reproduces the Holway model flux calculations with high accuracy: within 0.4% and −4% for evaporation and condensation, respectively. The discontinuity of the relative error of the moment method at $Δp*$ = 0 results from the additional approximation of the fourth mode required for condensation. The linearized moment expression outperforms the explicit Schrage equation in predicting the evaporative flux and the lower condensation flux regime. Equation (16) is an especially good choice of model for low flux evaporation/condensation, providing less than ±5% error in the range −0.1 < $Δp*$ < 0.1. This behavior is expected as Eq. (16) is an expansion of the highly accurate moment method around the equilibrium (zero flux) condition. On the other hand, the Schrage expressions do not have the guarantee of being especially accurate at low fluxes, as can be seen in Fig. 3. Along with this contrast, Eq. (16) is comparable to the explicit Schrage equation in complexity (both are explicit, linear in $Δp*$⁠, and facilitate computation for any value of $σ̂$⁠), making it a particularly attractive substitute for Eq. (7) in computational fluid dynamics (CFD). We also observed that, like the Schrage expressions, the moment method and its linearized formulation tend to increase in accuracy as $σ̂$ decreases (Appendix  E2). We also analyzed the accuracy of using the moment method and Eq. (16) for a polyatomic fluid with $j$ = 3. As one would expect following our observation in 3.5 that the flux is impacted little by internal degrees-of-freedom, we found that the relative errors of these two models increase marginally compared to the data in Fig. 7(b). For example, the maximum error of the moment method for evaporation increases to 1% (see Appendix  E2).

We note that applying Eq. (16) to condensation is physically flawed as it describes the flux as a function of $Δp*$ only when $TK*$ should be an additional parameter. However, Ytrehus' discussion [21] of linearization shows that the temperature ratio drops out of the mass flux expression, even for condensation. This finding is in agreement with our analysis in Sec. 3.5: our data shows that the condensation flux has a weak dependence on $TK*$⁠. Further, the use of Eq. (16) to model net condensation seems no more inappropriate than the use of the Schrage expressions for evaporation, where the specification of both $Δp*$ and $TK*$ is physically flawed. Equation (16) outperforms the explicit Schrage equation for $Δp*≳−0.3$⁠.

As noted in Sec. 2.3, solving the BTE for the half-space problem provides boundary conditions for the continuum region, in which the Navier–Stokes equations aptly describe the flow. In engineering applications, complicated modeling (such as two-phase CFD) may be needed for the liquid and the continuum vapor region. It is useful to regard $pK*$⁠, $TK*$⁠, and $J*$ (or $MK$⁠) as jump conditions that form the boundary conditions at the interface for the Navier–Stokes equations. Here we assume that the system has a continuum region after the Knudsen layer, requiring that the characteristic length scale (⁠$d$⁠) is much larger than the mean free path (i.e., that the Knudsen number Kn $≡λ/d$ is small). We note that $d$ should be taken as a geometric length scale or macroscopic boundary layer thickness, whichever is smaller [21,69]. For small Knudsen numbers, the solutions to the half-space problem discussed in Secs. 2.3, 3.1, and 4 give accurate jump conditions for practical evaporation/condensation scenarios in the moderately strong and strong regimes, and often in the weak mass transfer regime as well.

where $k$ is the thermal conductivity of the vapor and $T(0)$ is the external (Navier–Stokes) solution evaluated at the Knudsen layer boundary. When heat supply/removal occurs in the liquid or via radiation, the mass transfer rates can be stronger, and it has been shown that the half-space assumption for the Knudsen layer is valid since the external gradients become negligible [21,69]. Equations (18) and (19) converge to the half-space solutions under the condition $|Δp*|≪$ 1, provided the external temperature gradient is small. Here, we focus on connecting the half-space approaches to the vapor flow.

The Schrage expressions provide a mass flux boundary condition at the liquid–vapor interface that depends on the pressure ratio and the temperature ratio, that is, $pK/pe$ and $TK/TL$⁠, consistent with the Knudsen layer-aware kinetic theory analysis of net condensation (Sec. 2.3). However, the implementation of the boundary condition must be different for net evaporation, as the kinetic theory dictates that only the pressure ratio or temperature ratio can be prescribed to give the mass flux. Choosing to prescribe the pressure ratio, , for instance, a boundary condition must be imposed on both the mass flux and $TK$⁠. In this way, the temperature of the vapor at the boundary is still related to the other variables, despite not appearing in the mass flux expression. Alternatively, the temperature ratio could be prescribed in the boundary condition and $pK$ would then be output alongside the mass flux. We note that because evaporation is a one-parameter problem, the vapor stream at the Knudsen layer boundary is not required to be in a saturated state. That is to, say, given an input $Δp*$⁠, the temperature ratio could correspond to a saturated, superheated, or supersaturated state, depending on the saturation curve of the fluid (Ytrehus [21] discusses this fact at greater length).

We remark here that when dealing with the classic case of condensation from a saturated vapor reservoir, $pK*$ is often considered to be the pressure of the saturated vapor, but Eq. (17) does not guarantee that $TK=Tsat(pK)$⁠. Therefore, the saturated condition occurs further upstream in the continuum region. This feature is not unique to the linear approximation of condensation; the same result is seen when coupling the fully nonlinear moment method to the Navier–Stokes description of the continuum region [21]. Depending on the fluid, it is even possible for an inverted temperature gradient to exist in the continuum region (with the vapor hotter closer to the liquid); this phenomenon is discussed in greater detail by Ytrehus [21] and Cercignani [38].

The final point we will comment on regarding the external continuum flow is that engineered systems commonly have a flow with a nonzero velocity component tangential to the liquid–vapor interface. The tangential velocity inside the Knudsen layer must be zero in the case of evaporation but not in the case of condensation [21,70,71]. However, the shear flow has little effect on the normal component of the condensing flow in the Knudsen layer unless the tangential Mach number (⁠$MτK$⁠) is greater than 3 [21,71]. Thus, the effect of tangential flow on the Knudsen layer solution for condensation can safely be ignored when $MτK≪$ 1.

Many more interesting aspects of evaporating and condensing flows have been studied by the kinetic theory community. We recommend Ytrehus' very thorough treatment of the half-space problem [21] (which includes a worked example for implementing Eqs. (16) and (17) in CFD), Kogan's review [69], Cercignani's textbook [38], Sone's textbook [15], and Frezzotti and Barbante's recent review [28]. These texts offer a plethora of information on the subject and point to many other enlightening publications.

Regarding the implementations we have discussed, the moment method correctly requires both $pK*$ and $TK*$ to compute $J*$ for condensation, while needing only one input to determine the other two quantities for evaporation. It is important to consider that small errors in calculating $TK*$ for evaporation can often result in large discrepancies in the temperature difference $TL$ – $TK$ because $TK*$ is a ratio of absolute temperatures. For example, a 1% error in $TK*$ at low fluxes and $TL$ = 300 K produces an error of ∼3 K in $TL$ – $TK$⁠. The effect of internal degrees-of-freedom must also be considered: although we demonstrated in Sec. 3.5 that the flux is mostly insensitive to $j$⁠, Frezzotti [43] pointed out that the same is not true of the temperature ratio. In comparison to the Holway model (with complete accommodation), the moment method predicts $TK*$ for evaporation with at most 0.6% relative error (Fig. 8(a), the error would be significantly worse if not using Cercignani's method to account for internal degrees-of-freedom). The linearized moment method should be used cautiously when an accurate calculation of the temperature difference is required: its temperature ratio error exceeds 1% for $Δp*$ > 0.25 for $j$ = 0 and $j$ = 2, and for $j$ = 3 it exceeds 1% when $Δp*$ > 0.14. As only the value of $j$ affects the $TK*(pK*)$ relation for evaporation, Fig. 8(a) can be used to determine the $TK*$ accuracy of these two methods for complete accommodation of any fluid in evaporation conditions up to $Δp*$ = 0.5. Both models see a reduced error for $σ̂$ = 0.25 (Fig. 8(b)), suggesting that the temperature ratio prediction error decreases with decreasing $σ̂$ in the same manner as the flux prediction error. The temperature ratio error (and flux error) of these models for arbitrary values of $σ̂$ can be directly calculated from the evaporation datasets we have provided (see Appendix  F) by transforming the pressure ratio using Eq. (13) and comparing the transformed data to Eq. (16) and the moment method.

To quantify the temperature ratio prediction accuracy of the moment method for net condensation, we computed the temperature ratio predicted by the Holway model as a function of dimensionless driving pressure at the dimensionless flux given by the linearized moment method, i.e., $TK*(Δp*,JEq. (16)*(Δp*))$⁠. We assessed the accuracy of the moment method and linearized moment method by comparing $TK*(Δp*,JEq. (16)*(Δp*))$ to $TK,Moment*(Δp*,JEq. (16)*(Δp*))$ and $TK,Eq. (17)*(Δp*)$⁠, respectively, where the former is the temperature ratio predicted by the moment method at the dimensionless flux calculated by Eq. (16). These evaluations are shown in Fig. 9(a), where we find that the moment methods are much less accurate at temperature ratio prediction for condensation than evaporation. Further, we observe that these monatomic-specific models become more inaccurate when used to model species with two and three internal degrees-of-freedom.

However, one would expect from our findings in Sec. 3.5 that since the flux is weakly dependent on the temperature ratio, the temperature ratio ought to be strongly dependent on the flux if $Δp*$ and $J*$ are taken as the free parameters for condensation. Thus, the results of Fig. 9(a) should be highly specific to the condition $J*=JEq.(16)*(Δp*)$⁠. As shown in Fig. 9(b), this condition leads to large temperature ratios when the driving force for condensation is large. Comparing the temperature ratio calculated by the Holway model (solid black line in Fig. 9(b)) to the saturation curve of water at $TL$ = 25 °C (dashed-dotted black line in Fig. 9(b)), it is apparent that these conditions typically correspond to heavily superheated vapor. Further, we caution that the moment method can be more unreliable for temperature ratio prediction at high condensation fluxes than indicated in Fig. 9(a), see Appendix  E3.

Another calculation of interest is to find the vapor saturation state needed to produce a certain condensing flux on a liquid at $TL$⁠. In contrast to the case demonstrated in Fig. 9, both $Δp*$ and $TK*$ are adjustable parameters to match the flux, but are subject to the constraint $TK*=Tsat*(Δp*;TL)$⁠. We quantified the accuracy of the moment method and Eq. (16) for determining these saturated vapor states by comparing the predicted temperature ratios to that of the Holway model. We performed this analysis for saturated water vapor with $TL$ = 25 °C and $σ̂$ = 1 (Fig. 10). For the linearized moment method, we solved Eq. (16) for $Δp*$ and selected the temperature ratio from the saturation curve (using Eq. (17) results in a higher error and a nonsaturated state for the water vapor). In this scenario, the moment method and linearized moment method achieve less than 0.1% and 0.5% temperature prediction error, respectively.

The Schrage equations for evaporation or condensation rates of pure fluids have long been understood to satisfy neither momentum nor energy conservation at the liquid–vapor interface. We have quantified the inaccuracy of the Schrage equations via comparison to numerical solutions of the Holway model, an extension of the BGK model to polyatomic gases. Our conclusions summarized below apply to the range of conditions −0.5 ≤ $Δp*$ ≤ 0.5, covering the low to moderately high flux regimes that are typically encountered in thermofluid engineering. Our key findings are:

1. The complete Schrage expression (Eq. (6)) overpredicts evaporation/condensation fluxes by ∼15% in the case of perfect accommodation.

2. The explicit Schrage equation (Eq. (7)) overpredicts evaporation fluxes by more than Eq. (6), yet is more accurate than Eq. (6) in predicting condensation fluxes.

3. For polar fluids such as water where the accommodation coefficient is likely less than 0.5, the explicit Schrage equation is typically accurate within 10%. Its accuracy tends to be higher for smaller values of $σ̂$⁠, although its error does not asymptotically tend to zero as $σ̂$ → 0.

   As our numerical solver for the Holway model is rather computationally expensive, we also analyzed the accuracy of the moment method (Eqs. (E1)–(E3), (E8)–(E10)) and the linearization of this moment method (Eqs. (16) and (17)) described by Ytrehus [21] and Cercignani [38]. Our findings and recommendations are:

4. The moment method is a highly accurate model for mass flux, providing less than 0.4% flux prediction error for evaporation and less than 4% for condensation of monatomic fluids with complete accommodation. The temperature prediction error of the moment method for evaporation is less than 0.6% for any number of rotational degrees of freedom.

5. The linearized moment method (Eqs. (16) and (17)) is quite accurate in the low flux regime, achieving less than 5% flux prediction error in the range −0.1 < $Δp*$ < 0.1 for monatomic fluids with complete accommodation. These equations are comparable in complexity to the explicit Schrage equation, which overpredicts fluxes by 16%–25% in the same range. For low flux evaporation (0 ≤ $Δp*$ < 0.1), Eq. (17) provides a temperature ratio prediction error of less than 0.75% for any number of rotational degrees of freedom. Further, Eq. (16) is more accurate than the explicit Schrage equation for any $Δp*$≳−0.3. We recommend using Eqs. (16) and (17) instead of the explicit Schrage equation.

6. While the flux prediction accuracies mentioned in conclusions (4) and (5) apply to monatomic fluids, little additional error arises in applying these methods to polyatomic gases: the moment method relative errors increase to 1% and 4.1% for evaporation and condensation, respectively. For water—and other polyatomic species with low accommodation coefficients—these methods are more accurate than the values stated in (4) and (5) because their accuracy improves as $σ̂$ decreases.

7. Neither the moment method nor its linearization provides high-accuracy predictions of the temperature ratio for net condensation. Both models are only accurate for monatomic species at low fluxes. More theoretical work is needed to address the lack of computationally inexpensive, high-accuracy models for condensation. However, both models accurately predict saturated vapor states when the saturation constraint is enforced.

As liquid–vapor phase change devices reach increasingly high heat fluxes, high-accuracy models of the interfacial kinetic resistance become essential to the calculation of transport rates. Here, we have evaluated a selection of models based on a comparison to theoretical results. However, the moderate to strong mass transfer regimes of kinetically limited evaporation and condensation are difficult to access experimentally; thus, there is little data available for experimental validation. More experiments pushing further into the kinetically limited regime would be useful for such validation and clarification of the accommodation coefficients of various fluids. To facilitate the adoption of higher-accuracy alternatives to the Schrage equation, we have: (a) provided select sets of numerical data from our Holway model calculations (see Appendix  F and Supplemental Materials on the ASME Digital Collection); (b) fit our data to engineering equations (see Appendix  G); and (c) provided code for moment method calculations of evaporation and condensation for arbitrary values of $σ̂$⁠.² These approaches, as well as Eqs. (16) and (17), are valid for modeling evaporation/condensation in systems with vanishingly small Knudsen numbers and a single species in the vapor phase.

We are grateful to Professor Gang Chen for his insights and stimulating discussions.

• Air Force Office of Scientific Research (Grant No. FA9550-19-1-0392; Funder ID: 10.13039/100000181).

$a$ =

hard sphere diameter, m

$a(x)$ =

scalar functions in Eq. (14)

$B$ =

arbitrary constant

$CJ$ =

convergence criterion for flux

$CR$ =

convergence criterion for residuals

$C1−C4,D1−D4,$$K1−K4$ =

constants for engineering equations

$d$ =

characteristic length scale, m

$d1$ =

constant used in spatial grid definition

$d2$ =

constant used in spatial grid definition

$f$ =

velocity distribution function, kgs³/m⁶

$f′$ =

velocity distribution function for a polyatomic gas, s³/m⁶ J

$F$ =

reduced distribution function, s/m⁴

$F̂$ =

reduced local equilibrium distribution function, s/m⁴

$F−$ =

integration factor

$G$ =

reduced distribution function, 1/m²/s

$Ĝ$ =

reduced local equilibrium distribution function, 1/m²/s

$G−$ =

integration factor

$H$ =

reduced distribution function, J s/m⁴

$Ĥ$ =

reduced local equilibrium distribution function, J s/m⁴

$H−$ =

integration factor

$Hj−$ =

modified integration factor

$I(f,f)$ =

binary collision integral

$j$ =

number of internal degrees of freedom

$J$ =

mass flux, kg/m²/s

$k$ =

thermal conductivity, W/m/K

$kB$ =

Boltzmann constant, J/K

Kn =

Knudsen number

$L$ =

normalized length of domain

$m$ =

particle mass, kg

M =

Mach number

$n$ =

number density of molecules having rotational energy $ϵ$⁠, 1/m³/J

$N$ =

number density, 1/m³

$Ne′$ =

modified number density defined by Eq. (B17), 1/m³

$p$ =

pressure, Pa

$Q$ =

calculated quantity

$R$ =

specific gas constant, J/kg/K

$S$ =

speed ratio

$t$ =

time, s

$T$ =

temperature, K

$u$ =

net velocity vector, m/s

$x$ =

cartesian coordinate, m

$x̂$ =

unit vector normal to the liquid/vapor interface directed toward the vapor phase

$z¯$ =

fraction of inelastic collisions

Greek Symbols

$α−$ =

free parameter of the moment method for net condensation

$β−$ =

free parameter of the moment method

$γ$ =

ratio of specific heats

$Γ$ =

Schrage's correction factor

$Γ′$ =

Schrage's correction factor, linearized

$Γ¯$ =

Gamma function

$Δp*$ =

dimensionless driving pressure

$ϵ$ =

rotational energy, J

$η$ =

proportionality constant between dimensionless temperature and pressure jumps

$λ$ =

mean free path, m

$μ$ =

kinematic viscosity, Pa s

$ν$ =

collision frequency, 1/s

$ξ$ =

molecular velocity space coordinate

$ρ$ =

density, kg/m³

$σ̂$ =

accommodation coefficient

$σ̂c$ =

condensation coefficient

$σ̂e$ =

evaporation coefficient

$Φ$ =

local equilibrium distribution function

$χ$ =

accommodation prefactor

$ω$ =

constant used in Eqs. (16) and (17)

$ω′$ =

constant used in Eq. (17)

Other Symbols

$∞$ =

infinity, far field

Subscripts

$c$ =

characteristic

$e$ =

equilibrium

$el$ =

elastic

$H$ =

Holway

$hs$ =

hard sphere

$i$ =

index

$in$ =

inelastic

$J$ =

flux

$k$ =

index

$K$ =

Knudsen

$L$ =

liquid

$M$ =

Maxwell–Boltzmann distribution

$mod$ =

model

$r$ =

rotational

$R$ =

relative

$sat$ =

saturation

$t$ =

translational

$tot$ =

total

$v$ =

vapor

$x$ =

vector component

$τ$ =

tangential

$4$ =

fourth mode

Superscripts

$(i)$ =

Cartesian grid point

$(k)$ =

velocity grid point

$*$ =

normalized variable

$+$ =

$ξx>0$ portion of a distribution function or flux

$−$ =

$ξx<0$ portion of a distribution function or flux

Abbreviations

BGK =

Bhatnagar, Gross, and Krook

BTE =

Boltzmann transport equation

DSMC =

direct simulation Monte Carlo