1 Introduction

The attention of many researchers during the last decades focused on the study of heat and mass transfer processes in microchannels and microporous media due to promising perspectives of their use for the thermal control in microelectronics, heat pipes, radiators etc. Microchannel systems can provide effective heat removal at essentially small dimensions of the cooling systems of aerospace electronic equipment (Konovalov et al. 2019). The use of microradiators for the thermal control enables reducing absolute maximum temperatures and temperature gradients on the surfaces of devices subject to high heat fluxes (Konovalov et al. 2019).

Models of porous media are used in studies of processes in shale gas formations (taking into account gas desorption and slippage effects), as well as in the study of fractal porous media. In the paper of Foroozesh et al. (2019), a pore network model was developed to model the gas desorption and the slip flow effects in recovery of shale gas reservoirs. Quantitative evaluation of the effective thermal conductivity of porous media has received wide attention in scientific investigations and engineering. The influence of geometrical factors, porosity and relative surface roughness on the effective thermal conductivity in porous media were discussed and analyzed in Qin et al. (2019).

Studies of heat transfer processes during the flow of two-phase coolants through porous media are relevant for various engineering applications in geothermal systems, nuclear waste disposal systems, storage systems, as well as at transport of thermal energy, see Singh and Myong (2018). Four applications of porous media in the energy sector were considered in the work (Singh and Myong 2018), such as electrokinetic energy converting devices, membrane-based water desalination through reverse osmosis, shale reservoirs and hydrogen storage systems. In addition, this paper provides recommendations on the further development of models of convective and diffusion flows in porous micro- and nanoscale media.

Simulation of heat transfer processes in porous media during boiling and condensation is a rather complicated task. Both mathematical models of two-phase mixtures, and models with a separate description of phases are used for the modelling purpose.

The study of Wang (1997) focused on the modeling two-phase heat transfer processes in porous media characterized by the presence of two-phase and single-phase zones with irregular and moving boundaries between them. The mathematical model included the equations of conservation of mass, momentum, and energy for two-phase mixtures, as well as a system of closing equations characterizing the interaction of phases at the phase boundaries.

The mathematical model developed in Wang (1997) was further used to solve the problem of complete evaporation of a liquid inside porous evaporators made of circular pipes (Alomar et al. 2014a, 2015) and diffusers (Alomar et al. 2014b). The authors used an algorithm for smoothing the effective diffusion coefficient in the regions of transition from convective heat transfer to two-phase flow and further to superheated vapor. The results of calculations demonstrated an influence of the smoothing algorithm of the effective diffusion coefficient on the axial temperature distribution. The influence of smoothing was especially significant in the region of transition from two-phase flow to superheated vapor.

In the works of Wang and Beckermann (1993a, b), a mathematical model of two-phase heat and mass transfer in capillary porous media was presented, which considered two-phase flow as a binary mixture. The model was completed with the mass conservation equation for the fluid. The temperatures of the liquid, vapor, and porous medium were assumed to be the same. The model of a two-phase mixture of Wang and Beckermann (1993a) was used to study boiling flow along a heated surface embedded in a porous medium by Wang and Beckermann (1993b). The problem was solved in the boundary layer approximation. In addition, it was assumed that there was a thin capillary liquid layer with constant saturation on a solid surface. The results of calculations of the saturation profile, phase velocities, and the dependence of the wall heat flux, including the maximum (critical) wall heat flux, on the saturation value are documented in the work of Wang and Beckermann (1993b).

The problem of unsteady heat and mass transfer in a cell of a porous body was studied by Duval et al. (2004). It was assumed that the problem of two-phase fluid flow could be solved independently of the heat transfer problem. The mathematical model of heat and mass transfer included the continuity equations for vapor and liquid phases and the energy equations for vapor, liquid and porous phases. As an example, a stratified unit cell was considered, in which a solid was bordered by a layer of vapor surrounded in its turn by a layer of liquid. The problem was solved numerically. The results of calculations were compared with experimental data for the heating mode of a porous medium, whereas the initial temperature of three phases was equal to the saturation temperature and a volumetric heat flux in a porous body was equal to 107 W/m3.

The model of non-equilibrium two-phase flow in porous media developed by Duval et al. (2004) was used by Bachrata et al. (2012) to simulate high-temperature processes in a porous layer formed during the fragmentation of nuclear reactor fuel in the event of a severe accident with core destruction. The peculiarities of the model are additional hypotheses and a set of the closure equations. The process of flooding the porous layer of fuel fragments with water is essentially non-equilibrium, therefore, the transition region of nucleate boiling was also considered by Bachrata et al. (2012), where the heat transfer rate depends on the local void fraction, mass flow rate and temperature. Experimental data were correlated as an empirical equation that considers these effects. The heat flux in the transition region was determined by the combination of heat transfer at nucleate and film boiling. The computations were performed using the ICARE CATHARE code. In addition, the authors presented a dependence for the speed of movement of the wetting front on the parameters of injection of the coolant.

The results of the study of boiling and condensation in a porous wick of a heat pipe were presented in the work of Hari et al. (2015). Here a mathematical model of a two-phase mixture was used. The problem was solved in the boundary layer approximation. The saturation distributions for boiling and condensation regimes were presented in Hari et al. (2015) and compared with the results of other authors.

A 3D mathematical model of a two-phase flow in a heat pipe with a porous wick was developed by Brahim and Jemni (2012). The authors simulated performance of a heat pipe under critical conditions of high thermal loads. The wick of the heat pipe was modeled as a homogeneous and isotropic porous medium. The problem was solved numerically by the finite volume method. The computations demonstrated that at high heat flux the tangential velocity increases, especially at the liquid–vapor interface, which has a significant effect on phase changes. The 3D analysis yielded more accurate results than the 2D approach and improved visualization, especially for various external conditions.

The work of Kiseev (2001) presents the results of a theoretical and experimental study of heat and mass transfer in capillary porous wicks of heat pipes capable of operating under weightless conditions. Two-phase flows were studied in capillary structures with an average pore diameter of 0.5–15 μm. The authors considered closed loop heat pipes with separated flow of the vapor and liquid phases. It was revealed by Kiseev (2001) that a significant increase in the temperature load causes a pulsation mode of operation, whereas droplets of liquid were ejected from capillaries. The modeling of transport processes in a heat pipe by Kiseev (2001) was reduced to solving a conjugate problem of heat transfer with initial and boundary conditions in the heat source and heat sink regions describing heat transfer with the environment. The system of Navier–Stokes equations was solved numerically. The results of computations were validated against the experimental data obtained in the study. As a result, the authors revealed the main peculiarities of heat transfer in the region of vapor removal channels. In case of high thermal loads, the vaporization front retreated deep into the wick thus causing unfavorable conditions for vaporization.

Results of theoretical and experimental studies of heat transfer in a porous wick of a heat pipe are presented in the work of Hanlon and Ma (2003). Water was supplied at both inlets to the area of the porous wick. The vapor was released through the top surface of the wick. The bottom surface of the wick was heated part. A 2D steady-state problem was solved. Computations indicated that heat transfer during evaporation of a thin film from the surface of the wick depended on the particle size and the porous layer thickness. As an average particle radius of the porous layer decreased, the heat transfer coefficient increased. Computed values of the heat transfer coefficient on the upper surface of the wick were compared with experimental data. Even though the predicted heat transfer coefficient exceeded the experimental data, a qualitative analysis suggested that evaporation from the top surface of the wick played an important role in heat transfer.

Effects of local thermal non-equilibrium (LTNE) on the cooling process in porous media with a phase change of a coolant were considered by Shi and Wang (2011), Yuki et al. (2008), Xin et al. (2014), Alomar et al. (2017, 2018). Presented there are mathematical models of heat and mass transfer of a two-phase mixture in a porous body, which consider temperature difference between the coolant and the porous matrix. Numerical modeling of heat and mass transfer in porous media at high heat fluxes was performed in the work of Yuki et al. (2008), who pointed out that thermal non-equilibrium must be taken into account in modeling of two-phase flows in a porous body at high heat fluxes (above 1 MW/m2) and high flow velocities. Xin et al. (2014) presented results of numerical simulations of heat and mass transfer in a porous medium based on a mathematical model with separate phase flow accounting for thermal non-equilibrium. Computations indicated that the temperature difference between a solid and a liquid was less than 1 K in entire region except for the interface between two-phase flow and vapor.

Two-phase flow in a porous medium was studied by Alomar et al. (2017), (2018) based on a two-phase mixture model with and without account for local thermal non-equilibrium between flow and the solid. The authors analyzed effects of the process parameters on the temperature distribution and saturation temperature.

Results of numerical modeling of complete evaporation of a liquid in an anisotropic porous medium obtained on basis of a two-phase mixture model accounting for the LTNE are presented in the works of Alomar et al. (2019) and Alazmi and Vafai (2004). It was shown by Alomar et al. (2019) that porous medium anisotropy and thermal conductivity of the solid phase significantly affected the beginning and the end of the phase transition process. Alazmi and Vafai (2004) demonstrated that the effect of variable porosity was very significant in the vicinity of the solid boundaries. Effects of local thermal non-equilibrium were enhanced due to thermal dispersion mechanism. Alazmi and Vafai (2004) indicated the range of variation of the process parameters, where effects of thermal dispersion and local thermal non-equilibrium must be taken into account.

A review of these studies demonstrates that they mostly deal with complex unsteady non-equilibrium heat transfer processes at two-phase flow in a porous medium. At modeling of these processes arise difficulties in determining of the heat transfer parameters in the region of transition from a subcooled liquid to boiling and afterwards at the transition to single-phase flow of superheated vapor. At modeling of heat and mass transfer in two-phase flow in porous media, it is necessary to take into account the temperature difference between the porous body and the flow in the area of transition from two-phase flow to vapor at high heat fluxes and high flow velocities.

Similar problems of convective heat and mass transfer at film boiling and condensation of nanofluid on a vertical wall in the approximation of the boundary layer with free outer flow were considered by Avramenko et al. (2014, 2015a, b, c, 2016, 2018). The authors studied the influence of physical properties of the nanofluid and outer flow parameters on the behavior of velocity and temperature profiles, as well on the heat transfer coefficients. In addition, stability of the film was investigated.

As follows from the literature review above, velocities of liquid and vapor phases in porous wicks of heat pipes are rather low so that fluid flow and heat transfer are laminar. Theoretical modelling in the overviewed works was performed using numerical methodology, whereas analytical solutions are practically absent. Analytical solutions, should they be available, provide a powerful tool for the analysis of the peculiarities of fluid flow and heat transfer at boiling in porous media, which clearly indicates physical trends and effects of different parameters and alleviates parametric studies. However, such analytical solutions and their analysis are to the knowledge of the authors not available in the literature so far.

Thus, the objective of the present study is to derive a novel analytical solution of the problem of laminar fluid flow and heat transfer at film boiling of a fluid in a porous medium bordering with a heated solid wall. The porous medium will be modeled in the Darcy–Brinkman and Darcy–Brinkman–Forchheimer approximation. Effects of different thermal boundary conditions on the heated wall (constant temperature and constant heat flux) and different parameters (Forchheimer parameter, Darcy number, permeability of the porous medium etc.) on the Nusselt number will be studied both qualitatively and quantitatively.

Not pretending on modelling of the entire range of non-equilibrium boiling in porous media we restricted ourselves by consideration of steady-state film boiling on a vertical wall under condition of local thermal phase equilibrium due to relatively low flow rates of the coolant.

2 Mathematical Model

The present paper focuses on the problem of heat transfer at boiling of a liquid in a semi-infinite porous medium, one side of which is bordering with a heated vertical solid wall (Fig. 1). The heat flux density on the heated wall is higher than the critical heat flux (CHF), therefore the liquid inside the porous medium evaporates and forms a layer (i.e. a film) of vapor between the liquid phase and the heated wall (Fig. 1).

Fig. 1
figure 1

Schematic of film boiling over a vertical wall

Full size image

In the chosen coordinate system, x coordinate is aligned with the surface of the wall, while the zero point of the y coordinate is located on the wall (Fig. 1). The vapor in the film flows steadily upwards in the direction of the larger x-values. The thickness δ of the vapor film is small compared to the length of the wall, which justifies modelling the film as a boundary layer.

This assumption was confirmed in the study of Kim et al. (2010), where the film thickness was an order of magnitude less than the length of the heated wall (the film thickness was 1 … 1.5 mm, the length was 50 … 100 mm).

Both the vapor and the fluid temperature \(T_{\infty }\) at the outer boundary of the film, i.e. at y = δ, are equal to the saturation temperature at a given pressure. The wall temperature is higher than that of the fluid: \(T_{w} > T_{\infty }\).

The problem will be solved under the following assumptions: porous medium is isotropic; inertia forces in the vapor film can be neglected in comparison with the viscosity and gravity forces; thermal conductivity and convective heat and mass transfer in the vapor film in the streamwise direction (i.e. x-coordinate) are vanishingly small than heat transfer in the direction orthogonal to the wall (i.e. y-coordinate); friction is absent at the boundary between the fluid and the vapor phases; the density of the vapor is much smaller than the density of the liquid; the surface tension effects at the outer boundary of the vapor film are negligible.

Flow regimes in the boiling film in pool boiling are described in detail in the textbook of Çengel (2002). The boundaries of the flow regimes are defined with the help of the heat flux density and the excess temperature (the difference between the wall temperature and the saturation temperature of the liquid at a given pressure). For water, the end of the nucleate boiling and a transition to film boiling occurs for the critical (or peak) heat flux of about 1 MW/m2 and the excess temperature of about 30 °C. According to Çengel (2002), developed film boiling for water is observed for the excess temperatures of above 120 °C. We restrict, however, our analysis with the case, where the wall temperature does not exceed the value of 300 °C, so that the radiation effects remain negligible, see Çengel (2002). Cryogenic applications deal with fluids with very low boiling points (e.g. oxygen, nitrogen or helium), whereas film boiling begins below the melting point of the most of heater materials, so that this regime can be operated without any danger of burnout (Çengel 2002).

Subsequently, fluid flow and heat transfer of the vapor film over a heated wall can be described by the modified Darcy–Brinkman–Forchheimer model. According to Ward (1964), flow regimes in a porous medium can be classified using the Reynolds number based on the average velocity and permeability \(\text{Re}_{K} = u_{\text{av}} \rho \sqrt K /\mu\). According to this classification, for the \(\text{Re}_{K} < 1\) fluid flow can be described exclusively by the Darcy’s law. When \(\text{Re}_{K} > 1\), the Forchheimer quadratic correction must be taken into account. At high permeability values, the Forchheimer correction must also considered, see Nield and Bejan (2006).

According to Nield and Bejan (2006), the full Darcy–Brinkman–Forchheimer model has the following form

$$\rho \frac{1}{\phi }\nabla \left( {\frac{{{\mathbf{V}} \cdot {\mathbf{V}}}}{\phi }} \right) = - \nabla p + \mu \nabla^{2} {\mathbf{V}} - \frac{\mu }{K}{\mathbf{V}} - \frac{{c_{F} \rho }}{\sqrt K }\left| {\mathbf{V}} \right|{\mathbf{V}} + g\Delta \rho {\mathbf{j}},$$
(1)
$$\nabla \cdot {\mathbf{V}} = 0,$$
(2)
$$\left( {\rho c} \right)_{\text{eff}} {\mathbf{V}} \cdot \nabla T = \nabla \cdot \left( {k_{\text{eff}} \nabla T} \right),$$
(3)

where

$$\left( {\rho c} \right)_{\text{eff}} = \left( {1 - \phi } \right)\left( {\rho c} \right)_{s} + \phi \left( {\rho c_{p} } \right),$$
(4)
$$k_{\text{eff}} = \left( {1 - \phi } \right)k_{\text{s}} + \phi k,$$
(5)

ks stands for thermal conductivity of the solid porous medium itself, and k, as before, stands for thermal conductivity of the vapor.

Considering assumptions above, one can simplify system (1)–(3) in the boundary layer approach. The result is

$$\mu \frac{{{\text{d}}^{2} u}}{{{\text{d}}y^{2} }} - \mu \frac{u}{K} - \frac{{c_{\text{F}} \rho }}{\sqrt K }u^{2} = - g\Delta \rho$$
(6)
$$\frac{{{\text{d}}^{2} T}}{{{\text{d}}y^{2} }} = 0,$$
(7)

where cF is the dimensionless Forchheimer form-drag constant, see Nield and Bejan (2006).

The problem is described by the modifed Navier–Stokes equation and the energy equation. The Navier–Stokes equation is written in the Brinkman approximation, so that it considers the component of the Stokes friction according to Nield and Bejan (2006). This ensures fulfillment of the no-slip boundary condition on the walls. Also, this system of equations includes terms that consider the hydraulic resistance caused by the porosity of the medium. Namely, the second term in Eq. (6) describes the linear resistance (Darcy) and the third describes nonlinear resistance (Forchheimer). Since we consider fully developed fluid flow, the advection terms are not considered. In the steady-state process, velocity and temperature derivatives with respect to time are not taken into account.

Two sets of the boundary conditions will be implied while solving Eqs. (6) and (7).

The first set assumes no mechanical interaction between the vapor and the moving liquid. For a pure fluid (no porous medium), the following boundary conditions were considered by Bromley (1950):

$$u = 0 ,\quad T = T_{w} ,\quad {\text{at}} \;y = 0,$$
(8)
$$\left( {\frac{{{\text{d}}u}}{{{\text{d}}y}}} \right)_{y = \delta } = 0,\quad T = T_{\infty } \quad {\text{at}}\;y = \delta$$
(9)

In frames of this set of the boundary conditions, we will additionally study the case where the constant temperature of the heated wall is replaced with the constant heat flux on it

$$u = 0,\quad \left( {\frac{{{\text{d}}T}}{{{\text{d}}y}}} \right)_{y = 0} = - \frac{{q_{w} }}{{k_{\text{eff}} }}\quad {\text{at}}\;y = 0,$$
(10)
$$\left( {\frac{{{\text{d}}u}}{{{\text{d}}y}}} \right)_{y = \delta } = 0 ,\quad T = T_{\infty } \quad {\text{at}}\;y = \delta$$
(11)

which apparently has not been analysed yet in the literature.

For this problem, an increase in the mass flow rate through the vapor film due to the boiling is described by Bromley (1950) using the equation of the mass balance:

$${\text{d}}G = \frac{{q_{w} }}{{L_{\upsilon } }}{\text{d}}x.$$
(12)

The second set of the boundary conditions proposed by Ellion (1954) outlines the condition of a stationary liquid at the outer boundary of the vapor film. Ellion (1954) considered the following boundary conditions

$$u = 0,\quad T = T_{w} \quad {\text{at}}\;y = \, 0$$
(13)
$$u = 0 ,\quad T = T_{\infty } \quad {\text{at}}\;y = \delta$$
(14)

Again, as above, we will additionally study here the case where the constant wall temperature is replaced with the constant heat flux on the wall

$$u = 0,\quad \left( {\frac{{{\text{d}}T}}{{{\text{d}}y}}} \right)_{y = 0} = - \frac{{q_{w} }}{{k_{\text{eff}} }},\quad {\text{at}}\;y = 0,$$
(15)
$$u = 0,\quad T = T_{\infty } \quad {\text{at}}\;y = \delta$$
(16)

For the case of a stationary fluid, Ellion (1954) suggested mass balance equation in a modified form

$$q_{w} = L_{\upsilon } G_{l} .$$
(17)

Here

$$\int\limits_{0}^{x} {G_{l} {\text{d}}x} = \int\limits_{0}^{\delta } {\rho u{\text{d}}y} ,$$
(18)

whereas \(G_{l}\) is the mass flow rate through the film per unit of the surface length.

The advantages of the proposed models consist in that they enable us to analytically analyze a wide range of problems of boiling in porous media under various thermal and hydrodynamic boundary conditions.

3 Pure Fluids

At first let us consider heat transfer subject to two sets of the boundary conditions mentioned above. For the boundary conditions (8)–(11), Bromley (1950) obtained the velocity profile in the following form

$$u = \frac{{g\Delta \rho \delta^{2} }}{\mu }\left( {\frac{y}{\delta } - \frac{{y^{2} }}{{2\delta^{2} }}} \right).$$
(19)

Boundary conditions (8) and (9) (Tw = const) yield the following equations for the vapor film thickness, the temperature distribution and the Nusselt number, see Bromley (1950)

$$\delta = \sqrt[4]{{\frac{{4k_{\text{eff}} \mu \Delta Tx}}{{gL_{\upsilon } \rho \Delta \rho }}}},$$
(20)
$$T = T_{w} - \left( {T_{w} - T_{\infty } } \right)\frac{y}{\delta },$$
(21)
$${\text{Nu}} = \frac{hx}{{k_{\text{eff}} }} = \sqrt[4]{{\frac{1}{4}{\text{RaJa}}}},$$
(22)

where the Rayleigh number Ra is defined as

$${\text{Ra}} = \frac{{gx^{3} \rho \Delta \rho c_{p} }}{{k_{\text{eff}} \mu }},$$
(23)

and the Jakob number is defined as

$${\text{Ja}} = \frac{{L_{\upsilon } }}{{c_{p} \Delta T}}.$$
(24)

The boundary conditions (10) and (11) (qw = const) also yield the velocity profile (19). Then with the help of Eq. (12) one can obtain

$$\frac{{{\text{d}}G}}{{{\text{d}}x}} = \frac{\text{d}}{{{\text{d}}x}}\left( {\int\limits_{0}^{\delta } {\rho u{\text{d}}y} } \right) = \frac{\text{d}}{{{\text{d}}x}}\left( {\frac{{g\rho \Delta \rho \delta^{3} }}{3\mu }} \right) = \frac{{q_{w} }}{{L_{\upsilon } }}.$$
(25)

An integration of this equation under the boundary condition

$$\delta = 0\quad {\text{at}}\; x = \, 0$$
(26)

gives

$$\delta = \sqrt[3]{{\frac{{3q_{w} \mu x}}{{gL_{\upsilon } \rho \Delta \rho }}}}.$$
(27)

For the temperature profile one can obtain

$$T = T_{\infty } + \frac{{q_{w} }}{{k_{\text{eff}} }}\left( {\delta - y} \right).$$
(28)

Based on Eqs. (27) and (28), one can derive a relation for the Nusselt number

$${\text{Nu}} = \sqrt[4]{{\frac{1}{3}{\text{RaJa}}}}.$$
(29)

For the boundary conditions (13)–(16), Avramenko et al. (2015a) obtained the velocity profile in the following form

$$u = \frac{{g\Delta \rho \delta^{2} }}{2\mu }\left( {\frac{y}{\delta } - \frac{{y^{2} }}{{\delta^{2} }}} \right).$$
(30)

Using the boundary conditions (13) and (14) (Tw = const) Ellion (1954) obtained the following equations for the vapor film thickness and the Nusselt number

$$\delta = \sqrt[4]{{\frac{{12k_{\text{eff}} \mu \Delta Tx}}{{gL_{\upsilon } \rho \Delta \rho }}}},$$
(31)
$${\text{Nu}} = \sqrt[4]{{\frac{1}{12}{\text{RaJa}}}}.$$
(32)

The temperature distribution has again the form of Eq. (21).

Boundary conditions (15) and (16) (qw = const) also lead to the velocity profile (30). Using this profile, one can come to the relation for the mass flowrate

$$G = \int\limits_{0}^{\delta } {\rho u{\text{d}}y} = \frac{{g\rho \Delta \rho \delta^{3} }}{12\mu }.$$
(33)

Then, based on Eqs. (17) and (18), one can obtain a relation for the vapor layer thickness

$$\delta = \sqrt[3]{{\frac{{12q_{w} \mu x}}{{gL_{\upsilon } \rho \Delta \rho }}}}.$$
(34)

The use of the temperature profile (28) yields again Eq. (32) for the Nusselt number. So, for the boundary conditions (13), (14) and (15), (16) in view of Eqs. (17) and (18) for the mass flow rate, the Nusselt number remains the same for the boundary conditions Tw = const and qw = const.

All results described above are summarized in Table 1.

Table 1 Characteristics of heat transfer at film boiling of liquid on a vertical heated wall in the absence of a porous medium
Full size table

4 Darcy–Brinkman Approach

According to Avramenko et al. (2015b), in the range of ReK < 1 the quadratic drag (the third term) in the right-hand side of Eq. (6) is much smaller than the linear drag (the second term). Thus, Eq. (6) can be reduced to the following form

$$\mu \frac{{{\text{d}}^{2} u}}{{{\text{d}}y^{2} }} - \mu \frac{u}{K} = - g\Delta \rho ,$$
(35)

which corresponds to the Darcy–Brinkman approach. The general solution of the homogenous Eq. (35) is

$$u = \frac{gK\Delta \rho }{\mu } + C_{1} \cosh \left( {\frac{y}{\sqrt K }} \right) + C_{2} \sinh \left( {\frac{y}{\sqrt K }} \right).$$
(36)

Let us separately consider two cases mentioned above: (1) no mechanical interaction between the vapor and the moving liquid on the outer boundary of the vapor film and (2) a stationary liquid.

4.1 No Mechanical Interaction

In this case described by the boundary conditions (8)–(11), Eq. (36) yields the following velocity profile

$$u = \frac{gK\Delta \rho }{\mu }\left( {1 - \frac{{\cosh \left( {\frac{y - \delta }{\sqrt K }} \right)}}{{\cosh \left( {\frac{\delta }{\sqrt K }} \right)}}} \right) = \frac{gK\Delta \rho }{\mu }\left( {1 - \frac{{\cosh \left( {{\text{Da}}^{ - 1/2} \left( {1 - \eta } \right)} \right)}}{{\cosh \left( {{\text{Da}}^{ - 1/2} } \right)}}} \right).$$
(37)

Here

$${\text{Da}} = \frac{K}{{\delta^{2} }},$$
(38)

is the Darcy number, whereas \(\eta = \frac{y}{{\delta^{{}} }}\) is a dimensionless coordinate.

Profile (37) results in the relation for the mass flow rate through the vapor film

$$G = \int\limits_{0}^{\delta } {\rho u{\text{d}}y} = \frac{{g\rho K\Delta \rho \left( {\delta - \sqrt K \tanh \left( {{\text{Da}}^{ - 1/2} } \right)} \right)}}{\mu }.$$
(39)

We will split the further investigation into two separate subcases: Tw = const and qw = const.

4.1.1 Condition T w = const

For the condition Tw = const, the differential Eq. (12) for the film thickness has the following form

$$\frac{\text{d}}{{{\text{d}}x}}\left( {\frac{g\rho K\Delta \rho }{\mu }\left( {\delta - \sqrt K \tanh \left( {{\text{Da}}^{ - 1/2} } \right)} \right)} \right) = k_{\text{eff}} \frac{\Delta T}{{L_{\upsilon } \delta }},$$
(40)

where the effective thermal conductivity of the porous medium filled with vapor is determined by the relation (5).

A solution of this equation under the boundary condition (26) yields the following transcendental equation for the film thickness

$$K\ln \left( {\cosh \left( {\frac{\delta }{\sqrt K }} \right)} \right) + \frac{{\delta^{2} }}{2} - \delta \sqrt K \tanh \left( {\frac{\delta }{\sqrt K }} \right) = \frac{{k_{\text{eff}} \Delta T\mu x}}{{gL_{\upsilon } K\rho \Delta \rho }},$$
(41)

or

$$\ln \left( {\cosh \left( {\frac{1}{{\sqrt {Da} }}} \right)} \right) + \frac{1}{2Da} - \frac{1}{{\sqrt {Da} }}\tanh \left( {\frac{1}{{\sqrt {Da} }}} \right) = \frac{{k_{\text{eff}} \Delta T\mu x}}{{gL_{\upsilon } K^{2} \rho \Delta \rho }} = \frac{1}{{{\text{Da}}_{x} {\text{RaJa}}}}.$$
(42)

The linear temperature profile gives following equation for the Nusselt number

$${\text{Nu}} = \frac{x}{\delta }.$$
(43)

A combination of Eqs. (41) and (43) gives

$$2{\text{Da}}_{x} \ln \left( {\cosh \left( {\frac{1}{{{\text{Nu}}\sqrt {{\text{Da}}_{x} } }}} \right)} \right) + \frac{1}{{{\text{Nu}}^{2} }} - \frac{{2\sqrt {{\text{Da}}_{x} } }}{\text{Nu}}\tanh \left( {\frac{1}{{{\text{Nu}}\sqrt {{\text{Da}}_{x} } }}} \right) = \frac{2}{{{\text{Da}}_{x} {\text{RaJa}}}} = \frac{2}{{{\text{Ra}}_{\upsilon } {\text{Ja}}}},$$
(44)

or

$$2\ln \left( {\cosh \left( {S^{ - 1} } \right)} \right) + \frac{1}{{S^{2} }} - \frac{2}{S}\tanh \left( {S^{ - 1} } \right) = \frac{1}{{2{\text{Da}}_{x}^{2} {\text{Nu}}_{0}^{4} }},$$
(45)

where\({\text{Ra}}_{\upsilon } = {\text{RaDa}}_{x} ,\)

$$S = {\text{Nu}}\sqrt {{\text{Da}}_{x} } .$$
(46)

Thus, we have a transcendental equation to find the Nusselt number as a function of the local Darcy number \({\text{Da}}_{x}\) and the Nusselt number Nu0 for pure fluids, Eq. (22). Cheng and Verma (1981) used dimensionless Rayleigh and Darcy numbers, \({\text{Ra}}_{\upsilon }\) and \({\text{Ja}}\), to describe heat transfer in porous media. In view of this, we can also recast Eq. (41) as

$$2M^{2} \ln \left( {\cosh \left( {\frac{1}{MN}} \right)} \right) + \frac{1}{{N^{2} }} - \frac{2M}{N}\tanh \left( {\frac{1}{MN}} \right) = \frac{2}{\text{Ja}},$$
(47)

where

$$M = {\text{Da}}_{x} \sqrt {{\text{Ra}}_{\upsilon } } ,\quad N = \frac{\text{Nu}}{{\sqrt {{\text{Ra}}_{\upsilon } } }}.$$
(48)

This is another transcendental equation for the Nusselt number. Unfortunately, this equation cannot be solved analytically in a closed form. It can be however solved analytically in particular limiting cases. At first, we will analyze the limiting case of small permeability of the porous medium: K → 0 and Da → 0.

Case Da → 0

After expanding the left-hand side of Eqs. (41) or (42) in the McLaren series, one can obtain a relation for the film thickness

$$\delta = \sqrt[{}]{{\frac{{2k_{\text{eff}} \Delta T\mu x}}{{gL_{\upsilon } K\rho \Delta \rho }}}}.$$
(49)

It gives an equation for the Nusselt number

$${\text{Nu}} = \sqrt {\frac{1}{2}{\text{Ra}}_{\upsilon } {\text{Ja}}} .$$
(50)

The Nusselt number in porous media can be expressed with the help of the Nusselt number in a pure fluid. Then Eq. (50) can be transformed to

$$\frac{\text{Nu}}{{{\text{Nu}}_{0} }} = \sqrt {2{\text{Da}}_{0} } = \sqrt[4]{{2{\text{Da}}}},$$
(51)

where \({\text{Da}}_{0} = \frac{K}{{\delta_{0}^{2} }}\), and δ0 is the vapor film thickness in fluid flow without a porous medium.

As expected, Eqs. (50) and (51) indicate that heat transfer degenerates, if permeability of the porous medium approaches to zero.

Case Da → ∞

A series expansion of the left-hand side of Eqs. (41) or (42) at Da → ∞ yields

$$\delta^{4} - \frac{{4\delta^{6} }}{9K} = \frac{{4k_{\text{eff}} \mu x}}{{gL_{\upsilon } \rho \Delta \rho }}\Delta T = \delta_{0}^{4} ,$$
(52)

where δ0 is the film thickness for a pure fluid. This equation can be recast as

$$\frac{{\delta^{4} }}{{\delta_{0}^{4} }} - \frac{{4\delta_{0}^{2} }}{9K}\frac{{\delta^{6} }}{{\delta_{0}^{6} }} = \bar{\delta }^{4} - \frac{4}{9}\frac{{\bar{\delta }^{6} }}{{{\text{Da}}_{0} }} = 1.$$
(53)

A substitution

$$\bar{\delta } = \frac{\delta }{{\delta_{0} }} = \sqrt z$$
(54)

enables transforming Eq. (53) to

$$z^{2} - \frac{{4z^{3} }}{{9{\text{Da}}_{0} }} = 1.$$
(55)

This cubic equation has three roots. A physical meaning has the following root

$$\bar{\delta }^{2} = \frac{3}{4}\left( {{\text{Da}}_{0} + \frac{{3^{1/3} \left( { - 1} \right)^{2/3} {\text{Da}}_{0}^{2} }}{F} - 3^{ - 1/3} \left( { - 1} \right)^{1/3} F} \right),$$
(56)

where

$$F = 3{\text{Da}}_{0}^{3} - 8{\text{Da}}_{0} + 4\sqrt {4{\text{Da}}_{0}^{2} - 3{\text{Da}}_{0}^{4} } .$$
(57)

Considering a relation

$$\frac{\text{Nu}}{{{\text{Nu}}_{0} }} = \frac{{\delta_{0} }}{\delta } = \bar{\delta }^{ - 1}$$
(58)

one can further calculate the Nusselt number with the help of Eq. (56).

There exists an easier way to derive a relation for the Nusselt number. For the condition K → ∞, an assumption \({\text{Da}}_{ 0} \approx {\text{Da}}\) is valid. As a result, one can obtain

$$\delta = \frac{{\delta_{0}^{{}} }}{{\sqrt[4]{{1 - \frac{4}{{9{\text{Da}}}}}}}} \approx \frac{{\delta_{0}^{{}} }}{{\sqrt[4]{{1 - \frac{4}{{9{\text{Da}}_{0} }}}}}}.$$
(59)

This leads to the following equation

$$\frac{\text{Nu}}{{{\text{Nu}}_{0} }} = \sqrt[4]{{1 - \frac{4}{{9{\text{Da}}}}}} \approx \sqrt[4]{{1 - \frac{4}{{9{\text{Da}}_{0} }}}}.$$
(60)

The normalized Nusselt numbers \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) calculated by Eqs. (58) and (60) are depicted in Fig. 2.

Fig. 2
figure 2

Dependence of the normalized Nusselt number \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) on the Darcy number Da

Full size image

For Da > 1.5, the relative differences between Eqs. (58) and (60) are less than 3%.

4.1.2 Condition q w = const

For this condition, the differential equation, Eq. (12) for the film thickness looks as

$$\frac{\text{d}}{{{\text{d}}x}}\left( {\frac{{g\rho K\Delta \rho \left( {\delta - \sqrt K \tanh \left( {{\text{Da}}^{ - 1/2} } \right)} \right)}}{\mu }} \right) = \frac{{q_{w} }}{{L_{\upsilon } }}.$$
(61)

A solution of this equation under the boundary condition (26) brings the following transcendental equation for the film thickness

$$\delta - \sqrt K \tanh \left( {\frac{\delta }{\sqrt K }} \right) = \frac{{q_{w} \mu x}}{{gL_{\upsilon } K\rho \Delta \rho }}.$$
(62)

Using (28), one can obtain from Eq. (62)

$$\frac{{\delta^{2} }}{2} - \delta \sqrt K \tanh \left( {\frac{\delta }{\sqrt K }} \right) = \frac{{k_{\text{eff}} \Delta T\mu x}}{{gL_{\upsilon } K\rho \Delta \rho }}.$$
(63)

This leads further to the equation

$$\frac{1}{{S^{2} }} - \frac{2}{S}\tanh \left( {S^{ - 1} } \right) = \frac{2}{{3{\text{Da}}_{x}^{2} {\text{Nu}}_{0}^{4} }}$$
(64)

and finally

$$\frac{1}{{N^{2} }} - \frac{2M}{N}\tanh \left( {\frac{1}{MN}} \right) = \frac{2}{\text{Ja}}.$$
(65)

Case Da → 0

For this case, the film thickness can be expressed as

$$\delta = \frac{{q_{w} \mu x}}{{gKL_{\upsilon } \rho \Delta \rho }},$$
(66)

and the Nusselt number looks as

$${\text{Nu}} = \sqrt {{\text{Ra}}_{\upsilon } {\text{Ja}}}$$
(67)

or

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = 3{\text{Da}}_{ 0} = \sqrt[ 3]{{ 3 {\text{Da}}}}.$$
(68)

Figure 3 depicts the Nusselt numbers calculated by Eqs. (50), (67) and by Cheng and Verma (1981), who considered also nonlinear effects (advection and convection). The differences between the present calculations using a linear approach and the results of the nonlinear approach obtained by Avramenko et al. (2018) increase with the decreasing Jacob numbers.

Fig. 3
figure 3

Dependence of the Nusselt number on the Jacob number \({\text{Ja}}\)

Full size image

Case Da → ∞

For this case, the equation for the film thickness

$$\delta^{3} - \frac{{2\delta^{5} }}{5K} = \frac{{3q_{w} \mu x}}{{gL_{\upsilon } \rho \Delta \rho }} = \delta_{0}^{3}$$
(69)

is transcendental and cannot be solved analytically.

However, we can obtain its solution under the assumption \({\text{Da}}_{ 0} \approx {\text{Da}}\)

$$\delta = \frac{{\delta_{0}^{{}} }}{{\sqrt[3]{{1 - \frac{2}{{5{\text{Da}}}}}}}}.$$
(70)

This yields a solution for the Nusselt number

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt[3]{{1 - \frac{2}{{5{\text{Da}}}}}} \approx \sqrt[3]{{1 - \frac{2}{{5{\text{Da}}_{ 0} }}}}.$$
(71)

4.2 Stationary Fluid

For a stationary fluid under the boundary conditions (13)–(16), Eq. (36) yields the following velocity profile

$$u = \frac{gK\Delta \rho }{\mu }\left( {1 - \frac{{\cosh \left( {\frac{\delta - 2y}{2\sqrt K }} \right)}}{{\cosh \left( {\frac{\delta }{2\sqrt K }} \right)}}} \right) = \frac{gK\Delta \rho }{\mu }\left( {1 - \frac{{\cosh \left( {\frac{{{\text{Da}}^{ - 1/2} }}{2}\left( {1 - 2\eta } \right)} \right)}}{{\cosh \left( {\frac{{{\text{Da}}^{ - 1/2} }}{2}} \right)}}} \right).$$
(72)

In this case, the mass flow rate is

$$G = \int\limits_{0}^{\delta } {\rho udy} = \frac{{g\rho K\Delta \rho \left( {\delta - 2\sqrt K \tanh \left( {\frac{{{\text{Da}}^{ - 1/2} }}{2}} \right)} \right)}}{\mu }.$$
(73)

Again, we will further consider two separate subcases: Tw = const and qw = const.

4.2.1 Condition T w = const

For the condition Tw = const, Eqs. (17) and (18) give

$$\delta^{2} - 2\delta \sqrt K \tanh \left( {\frac{{{\text{Da}}^{ - 1/2} }}{2}} \right) = \frac{{k_{eff} \Delta T\mu x}}{{gL_{\upsilon } K\rho \Delta \rho }}.$$
(74)

Using Eq. (43), one can obtain the following transcendental equation for the Nusselt number

$$\frac{1}{{{\text{Nu}}^{2} }} - \frac{{2\sqrt {{\text{Da}}_{x} } }}{\text{Nu}}\tanh \left( {\frac{1}{{ 2 {\text{Nu}}\sqrt {{\text{Da}}_{x} } }}} \right) = \frac{1}{{{\text{Da}}_{x} {\text{RaJa}}}} = \frac{1}{{{\text{Ra}}_{\upsilon } {\text{Ja}}}},$$
(75)

or

$$\frac{1}{{S^{2} }} - \frac{2}{S}\tanh \left( {\frac{{S^{ - 1} }}{2}} \right) = \frac{1}{{ 1 2 {\text{Da}}_{x}^{ 2} {\text{Nu}}_{ 0}^{ 4} }}.$$
(76)

Another form of this transcendental equation is

$$\frac{1}{{N^{2} }} - - \frac{2M}{N}\tanh \left( {\frac{1}{2MN}} \right) = \frac{1}{\text{Ja}}.$$
(77)

Now let consider limiting cases, where closed form solutions of Eq. (74) can be obtained.

Case Da → 0

A series expansion of the left-hand side of Eq. (74) in the McLaren series yields

$$\delta = \sqrt[{}]{{\frac{{k_{ef} \Delta T\mu x}}{{gL_{\upsilon } K\rho \Delta \rho }}}},$$
(78)

In its turn, Eq. (78) enables obtaining a solution for the Nusselt number

$${\text{Nu}} = \sqrt {{\text{Ra}}_{\upsilon } {\text{Ja}}} .$$
(79)

Finally, one can reduce it to

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt {12{\text{Da}}_{ 0} } = \sqrt[4]{{12{\text{Da}}}}.$$
(80)

Case Da → ∞

A series expansion of the left-hand side of Eq. (74) in the McLaren series at K → ∞ brings

$$\delta^{4} - \frac{{\delta^{6} }}{10K} = \frac{{12k_{\text{eff}} \mu \Delta Tx}}{{gL_{\upsilon } \rho \Delta \rho }} = \delta_{0}^{4} .$$
(81)

Substituting Eq. (54) in Eq. (81), one can obtain

$$z^{2} - \frac{{z^{3} }}{{1 0 {\text{Da}}_{0} }} = 1.$$
(82)

A physically meaningful solution of this equation is

$$\bar{\delta }^{2} = \frac{1}{3}\left( {10{\text{Da}}_{0} + \frac{{20\left( { - 5} \right)^{2/3} {\text{Da}}_{0}^{2} }}{F} - \left( { - 5} \right)^{1/3} F} \right),$$
(83)

where

$$F = 200{\text{Da}}_{0}^{3} - 27{\text{Da}}_{0} + 3\sqrt {81{\text{Da}}_{0}^{2} - 1200{\text{Da}}_{0}^{4} } .$$
(84)

An approximate solution in present case is

$$\frac{\text{Nu}}{{{\text{Nu}}_{0} }} = \sqrt[4]{{1 - \frac{1}{{10{\text{Da}}}}}} \approx \sqrt[4]{{1 - \frac{1}{{10{\text{Da}}_{0} }}}}.$$
(85)

Comparisons of the calculations by Eq. (58), as well as Eqs. (83) and (85) are illustrated in Fig. 4.

Fig. 4
figure 4

Dependence of the normalized Nusselt number \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) on the Darcy number Da

Full size image

In the case of a stationary fluid (Fig. 4), like in the case of a moving fluid (Fig. 2), the normalized Nusselt numbers \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{0} }}} \right. \kern-0pt} {{\text{Nu}}_{0} }}\) decrease with the decreasing Darcy number. However, in the case of a stationary fluid, this decrease is weaker.

4.2.2 Condition q w = const

Based on Eqs. (17) and (18), one can obtain

$$\frac{g\rho K\Delta \rho }{\mu x}\left( {\delta - 2\sqrt K \tanh \left( {\frac{{{\text{Da}}^{ - 1/2} }}{2}} \right)} \right) = \frac{{q_{w} }}{{L_{\upsilon } }},$$
(86)

or

$$\frac{{\delta^{2} }}{2} - 2\delta \sqrt K \tanh \left( {\frac{{{\text{Da}}^{ - 1/2} }}{2}} \right) = \frac{{k_{\text{eff}} \Delta T\mu x}}{{gL_{\upsilon } K\rho \Delta \rho }}.$$
(87)

This equation can be transformed to the following transcendental equation for the Nusselt number

$$\frac{1}{{{\text{Nu}}^{2} }} - \frac{{2\sqrt {{\text{Da}}_{x} } }}{\text{Nu}}\tanh \left( {\frac{1}{{2{\text{Nu}}\sqrt {{\text{Da}}_{x} } }}} \right) = \frac{1}{{{\text{Ra}}_{\upsilon } {\text{Ja}}}}.$$
(88)

This equation also can be rewritten in the following form

$$\frac{1}{{S^{2} }} - \frac{2}{S}\tanh \left( {\frac{{S^{ - 1} }}{2}} \right) = \frac{1}{{12{\text{Da}}_{x}^{2} {\text{Nu}}_{0}^{4} }},$$
(89)

or

$$\frac{1}{{N^{2} }} - \frac{2M}{N}\tanh \left( {\frac{1}{2MN}} \right) = \frac{1}{\text{Ja}}.$$
(90)

Again, let us consider two limiting cases that enable obtaining closed form solutions of Eq. (90).

Case Da → 0

Using Eq. (66) for the film thickness and Eq. (67) for the Nusselt number, one can finally obtain

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt[3]{144}{\text{Da}}_{ 0} = \sqrt[ 3]{{12^{2/3} {\text{Da}}}}.$$
(91)

Case Da → ∞

The equation for the film thickness looks as

$$\delta^{3} - \frac{{\delta^{5} }}{10K} = \delta_{0}^{3} ,$$
(92)

which yields

$$\frac{\text{Nu}}{{{\text{Nu}}_{0} }} = \sqrt[3]{{1 - \frac{1}{{10{\text{Da}}}}}} \approx \sqrt[3]{{1 - \frac{1}{{10{\text{Da}}_{0} }}}}.$$
(93)

All solutions obtained in frames of the Darcy–Brinkman approach are summarized in Table 2.

Table 2 Parameters of heat transfer at film boiling of a liquid on a vertical heated wall for the Darcy-Brinkman approach
Full size table

5 Darcy–Forchheimer–Brinkman Approach

According to Ward (1964), under the condition

$$\text{Re}_{K} \ge 1$$
(94)

the Darcy–Forchheimer regime occurs, i.e. it is necessary to consider both the Darcy linear drag and the Forchheimer quadratic drag. In this case, momentum Eq. (6) holds.

For the subsequent comparisons, it is convenient to present the solution of Eq. (6) in a dimensionless form. For this purpose, we transform Eq. (6) using dimensionless variables, such as

$$\frac{{{\text{d}}^{2} U}}{{{\text{d}}\eta^{2} }} - \frac{U}{\text{Da}} - \frac{\text{Fh}}{{{\text{Da}}^{2} }}U^{2} = - 1,$$
(95)

where

$$\eta = \frac{y}{\delta },$$
(96)
$$U = \frac{u\mu }{{g\Delta \rho \delta^{2} }},$$
(97)
$${\text{Fh}} = \frac{{c_{F} g\rho \Delta \rho K^{3/2} }}{{\mu^{2} }}.$$
(98)

here Fh is a dimensionless parameter.

Equation (95) has a symmetry, which can be described by the following infinitesimal generator

$${\mathbf{v}} = \frac{\partial }{\partial \eta }.$$
(99)

It means that the function U is an invariant of this transformation and, according to Olver (1993), it can be accepted as a new independent argument. Then we can use a substitution

$$\frac{{{\text{d}}U}}{{{\text{d}}\eta }} = s\left( U \right),\quad \frac{{{\text{d}}^{2} U}}{{{\text{d}}\eta^{2} }} = s\left( U \right)\frac{{{\text{d}}s\left( U \right)}}{{{\text{d}}U}}.$$
(100)

Substituting Eq. (100) into Eq. (95) one can obtain

$$s\frac{{{\text{d}}s}}{{{\text{d}}U}} - \frac{U}{\text{Da}} - \frac{\text{Fh}}{{{\text{Da}}^{2} }}U^{2} = - 1.$$
(101)

Having performed a separation of variables and an integration, we have

$$\pm \int\limits_{0}^{U} {\frac{1}{{\sqrt {C_{1} - 2U + \frac{1}{\text{Da}}U^{2} + \frac{2}{3}\frac{\text{Fh}}{{{\text{Da}}^{2} }}U^{3} } }}{\text{d}}U} = \eta + C_{2} .$$
(102)

The integration constant \(C_{2}\) is equal to zero because of the boundary condition

$$u = 0\quad {\text{at}}\; y = 0$$
(103)

The left-hand side of the integral in Eq. (102) is an elliptic integral of the first kind. The inverse function of this solution is

$$U = {\text{root}}_{3} + \left( {{\text{root}}_{2} - {\text{root}}_{3} } \right){\text{sn}}^{2} \left( {Z\left| R \right.} \right),$$
(104)

where sn is the Jacobi elliptic function,

$$Z = \frac{1}{{2{\text{Da}}}}\sqrt { - \frac{{\eta^{2} }}{3}\left( {3{\text{Da}} + 2{\text{Fh}}\left( {{\text{root}}_{2} + 2{\text{root}}_{3} } \right)} \right)} ,$$
(105)
$$R = \frac{{{\text{root}}_{2} - {\text{root}}_{3} }}{{{\text{root}}_{1} - {\text{root}}_{3} }},$$
(106)

\({\text{root}}_{n}\) is the n the ordinal number of a particular root of the equation

$$3{\text{Da}}^{2} C_{1} - 6{\text{Da}}^{2} W + 3{\text{Da}}^{2} W^{2} + 2{\text{Fh}}W^{3} = 0.$$
(107)

It is necessary to say that Eq. (98) is a nonlinear differential equation. According to Kamke (1977), this type of equations can have different solutions. In the present case, Eq. (104) is not acceptable as a solution of the mathematical problem at hand.

Hence, we will try to obtain another approximate solution and then to compare it with the numerical solution. For this purpose, we will perform a linearization of Eq. (95). We will further introduce an unknown mean velocity Um in the last term of the left-hand side of Eq. (95). This yields as a result

$$\frac{{{\text{d}}^{2} U}}{{{\text{d}}\eta^{2} }} - \frac{U}{\text{Da}} - \frac{\text{Fh}}{{{\text{Da}}^{2} }}UU_{m} = - 1.$$
(108)

At first, we will consider case of no mechanical interaction and then stationary fluid.

5.1 No Mechanical Interaction

A dimensionless solution of the linearized Eq. (108) with the boundary conditions (8), (9) or (10), (11) has the following form

$$U = \frac{{{\text{Da}}^{2} }}{{{\text{Da}} + {\text{Fh}}U_{m} }}\left( {1 - \frac{{\cosh \left( {\left( {1 - \eta } \right)\frac{{\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } }}{\text{Da}}} \right)}}{{\cosh \left( {\frac{{\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } }}{\text{Da}}} \right)}}} \right).$$
(109)

The mean velocity Um can be expressed from Eq. (109), such as

$$U_{m} = \int\limits_{0}^{1} {U{\text{d}}\eta = } \frac{{{\text{Da}}^{2} }}{{\left( {{\text{Da}} + {\text{Fh}}U_{m} } \right)^{3/2} }}\left( {\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } - {\text{Da}}\tanh \left( {\frac{{\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } }}{\text{Da}}} \right)} \right).$$
(110)

The condition (94) indicates that it necessary to take into account the Forchheimer drag when \(\text{Re}_{K} \ge 1\). It means that at slow fluid flow in a porous medium we should as usually consider the problem in frames of the model for high Darcy numbers. In this case, Eq. (110) for the mean velocity Um can be presented in the following form

$$U_{m} = \frac{{1 - \frac{2}{{5{\text{Da}}}}}}{{3\left( {1 + \frac{{2{\text{Fh}}}}{{15{\text{Da}}^{2} }}} \right)}}.$$
(111)

A substitution of Eq. (111) into Eq. (109) gives

$$U = \frac{{15{\text{Da}}^{2} + 2{\text{Fh}}}}{{5\left( {3{\text{Da}} + {\text{Fh}}} \right)}}\left( {1 - \frac{{\cosh \left( {\left( {1 - \eta } \right)\sqrt {\frac{{5\left( {3{\text{Da}} + {\text{Fh}}} \right)}}{{15{\text{Da}}^{2} + 2{\text{Fh}}}}} } \right)}}{{\cosh \left( {\sqrt {\frac{{5\left( {3{\text{Da}} + {\text{Fh}}} \right)}}{{15{\text{Da}}^{2} + 2{\text{Fh}}}}} } \right)}}} \right).$$
(112)

A comparison of the numerical integration of Eq. (108) and the linearized solution (108) is presented in Fig. 5. It is evident that the numerical solution of Eq. (108) and the approximate analytical solution (112) agree quite well. The differences between them decrease, if the dimensionless parameters Da and Fh increase.

Fig. 5
figure 5

Comparisons of the numerical integration of Eq. (108) and the analytical solution (112) for the velocity profile in a vapor film depending on the Darcy number and the parameter Fh

Full size image

Of these two parameters, the variation of the Darcy number causes the most noticeable effect. As can be seen from Eq. (108), the square of the Darcy number is included in the denominator of the last term on the left-hand side. With an increase in the Darcy number, the role of this last term declines. Therefore, the difference between the numerical and analytical solutions decreases with the increasing Darcy number.

Expanding Eq. (111) in the McLaren series yields a dimensional form of the mean velocity

$$u_{m} = \frac{{g\Delta \rho \delta^{2} }}{3\mu } - \frac{{2c_{F} g^{2} \rho \Delta \rho^{2} \delta^{6} }}{{45\mu^{3} \sqrt K }} - \frac{{2g\Delta \rho \delta^{4} }}{15\mu K}.$$
(113)

Let us again consider two separate subcases: Tw = const and qw = const.

5.1.1 Condition T w = const

Using Eqs. (12) and (113), one can obtain

$$\delta^{4} - \frac{{4\delta^{6} }}{9K} - \frac{{7c_{F} g\rho \Delta \rho }}{{45\mu^{2} \sqrt K }}\delta^{8} = \delta_{0}^{4} ,$$
(114)

where δ0 is the film thickness for a pure fluid (no porous medium, see Table 1). Using the substitution (54), one can transform Eq. (114) to

$$z^{2} - \frac{4}{{9{\text{Da}}_{0}^{{}} }}z^{3} - \frac{{7{\text{Fh}}}}{{45{\text{Da}}_{0}^{2} }}z^{4} = 1.$$
(115)

A solution of this equation has four roots, which are cumbersome and therefore not presented here.

An approximate simplified approach assumes that \({\text{Da}}_{ 0} \approx {\text{Da}}\). In frames of this approach, one can obtain

$$\delta = \frac{{\delta_{0}^{{}} }}{{\sqrt[4]{{1 - \frac{4}{{9{\text{Da}}}} - \frac{{7{\text{Fh}}}}{{45{\text{Da}}_{{}}^{ 2} }}}}}} \approx \frac{{\delta_{0}^{{}} }}{{\sqrt[4]{{1 - \frac{4}{{9{\text{Da}}_{ 0} }} - \frac{{7{\text{Fh}}}}{{45{\text{Da}}_{ 0}^{ 2} }}}}}}.$$
(116)

Then the approximate solution for the Nusselt number is

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt[4]{{1 - \frac{4}{{9{\text{Da}}}} - \frac{{7{\text{Fh}}}}{{45{\text{Da}}_{{}}^{ 2} }}}} \approx \sqrt[4]{{1 - \frac{4}{{9{\text{Da}}_{ 0} }} - \frac{{7{\text{Fh}}}}{{45{\text{Da}}_{ 0}^{ 2} }}}}.$$
(117)

Figure 6 depicts comparisons of the calculated normalized Nusselt numbers based on Eq. (58) in combination with one of four roots of Eq. (115) against the normalized Nusselt numbers by Eq. (117). It is evident that these curves are in a good agreement with each other.

Fig. 6
figure 6

Comparisons of the numerically and analytically computed normalized Nusselt numbers \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) depending on the Darcy number and the parameter Fh

Full size image

Figure 6 indicates that account for the non-linear Forchheimer resistance is the reason for the Nusselt number to tend to zero at large values of permeability (Da), i.e. when the physical properties of the medium tend closer to those of a pure liquid.

As expected, an increase in the parameter Fh leads to heat transfer deterioration. Obviously, it can be attributed to the additional hydrodynamic drag.

5.1.2 Condition q w = const

For this condition, the differential Eq. (12) yields an equation for the vapor film thickness

$$\delta^{3} - \frac{2}{5K}\delta^{5} - \frac{{2c_{F} g\rho \Delta \rho }}{{15\mu^{2} \sqrt K }}\delta^{7} = \frac{{3q_{w} \mu x}}{{gL_{\upsilon } \rho \Delta \rho }} = \delta_{0}^{3}$$
(118)

Unfortunately, this equation cannot be solved. However, under the assumption \({\text{Da}}_{ 0} \approx {\text{Da}}\) we can obtain an approximate solution in the following form

$$\delta = \frac{{\delta_{0}^{{}} }}{{\sqrt[3]{{1 - \frac{2}{{5{\text{Da}}}} - \frac{{2{\text{Fh}}}}{{15{\text{Da}}_{{}}^{ 2} }}}}}}.$$
(119)

Consequently

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt[3]{{1 - \frac{2}{{5{\text{Da}}}} - \frac{{2{\text{Fh}}}}{{15{\text{Da}}_{{}}^{ 2} }}}} \approx \sqrt[3]{{1 - \frac{2}{{5{\text{Da}}_{ 0} }} - \frac{{2{\text{Fh}}}}{{15{\text{Da}}_{0}^{ 2} }}}}.$$
(120)

Figure 7 depicts computations based on Eqs. (117) and (120). They indicate that the effect of porosity expressed in the form of the Darcy number on the normalized Nusselt number is almost the same for the boundary conditions Tw = const and qw = const.

Fig. 7
figure 7

Comparisons of the numerically and analytically computed normalized Nusselt numbers \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) depending on the Darcy number and the parameter Fh

Full size image

5.2 Stationary Fluid

The solution of the dimensionless linearized Eq. (108) with the boundary conditions (13), (14) or (15), (16) has the following form

$$U = \frac{{{\text{Da}}^{ 2} }}{{{\text{Da}} + {\text{Fh}}U_{m} }}\left( {1 - \frac{{\cosh \left( {\left( {1 - 2\eta } \right)\frac{{\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } }}{{ 2 {\text{Da}}}}} \right)}}{{\cosh \left( {\frac{{\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } }}{{ 2 {\text{Da}}}}} \right)}}} \right).$$
(121)

The mean velocity Um can be expressed from Eq. (109) as

$$U_{m} = \int\limits_{0}^{1} {Ud\eta = } \frac{{{\text{Da}}^{ 2} }}{{\left( {{\text{Da}} + {\text{Fh}}U_{m} } \right)^{3/2} }}\left( {\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } - 2{\text{Da}}\tanh \left( {\frac{{\sqrt {{\text{Da}} + {\text{Fh}}U_{m} } }}{{ 2 {\text{Da}}}}} \right)} \right).$$
(122)

In frames of the model approach for high Darcy numbers, the mean velocity Um, Eq. (110), can be presented in the following form

$$U_{m} = \frac{{1 - \frac{1}{{ 1 0 {\text{Da}}}}}}{{12\left( { 1+ \frac{{1{\text{Fh}}}}{{ 1 2 0 {\text{Da}}^{ 2} }}} \right)}}.$$
(123)

A substitution of Eq. (123) into Eq. (121) yields

$$U = \frac{{ 1 2 0 {\text{Da}}^{ 2} + {\text{Fh}}}}{{10\left( { 1 2 {\text{Da}} + {\text{Fh}}} \right)}}\left( {1 - \frac{{\cosh \left( {\left( {1 - 2\eta } \right)\sqrt {\frac{{5\left( { 1 2 {\text{Da}} + {\text{Fh}}} \right)}}{{2\left( { 1 2 0 {\text{Da}}^{ 2} + {\text{Fh}}} \right)}}} } \right)}}{{\cosh \left( {\sqrt {\frac{{5\left( { 1 2 {\text{Da}} + {\text{Fh}}} \right)}}{{2\left( { 1 2 0 {\text{Da}}^{ 2} + {\text{Fh}}} \right)}}} } \right)}}} \right).$$
(124)

Results of the numerical integration of Eq. (108) and calculations by the linearized analytical solution (124) are shown in Fig. 8.

Fig. 8
figure 8

Comparisons of the numerically integrated and the analytically predicted velocity profiles in a vapor film depending on the Darcy number and the parameter Fh

Full size image

One can see from here that the agreement between the numerical solution of Eq. (108) and the approximate analytical solution (124) is very good.

An expansion of Eq. (123) into McLaren series yields a dimensional relation for the mean velocity

$$u_{m} = \frac{{g\Delta \rho \delta^{2} }}{12\mu } - \frac{{c_{F} g^{2} \rho \Delta \rho^{2} \delta^{6} }}{{1440\mu^{3} \sqrt K }} - \frac{{g\Delta \rho \delta^{4} }}{120\mu K}.$$
(125)

As above, let us again consider two separate subcases: Tw = const and qw = const.

5.2.1 T w = const

Using Eqs. (125), (17) and (18), one can obtain

$$\delta^{4} - \frac{{\delta^{6} }}{10K} - \frac{{c_{F} g\rho \Delta \rho }}{{120\mu^{2} \sqrt K }}\delta^{8} = \frac{{12k_{eff} \mu \Delta Tx}}{{gL_{\upsilon } \rho \Delta \rho }} = \delta_{0}^{4} .$$
(126)

A substitution of Eq. (54) allows recasting Eq. (126) as

$$z^{2} - \frac{1}{{10{\text{Da}}_{ 0}^{{}} }}z^{3} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}z^{4} = 1.$$
(127)

Again, we have four roots of Eq. (127).

An approximate solution of Eq. (126) at \({\text{Da}}_{ 0} \approx {\text{Da}}\) looks as

$$\delta = \frac{{\delta_{0}^{{}} }}{{\sqrt[4]{{1 - \frac{1}{{10{\text{Da}}}} - \frac{\text{Fh}}{{120{\text{Da}}_{{}}^{ 2} }}}}}} \approx \frac{{\delta_{0}^{{}} }}{{\sqrt[4]{{1 - \frac{1}{{10{\text{Da}}_{ 0}^{{}} }} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}}}}}.$$
(128)

Equation (128) leads to the solution for the Nusselt number

$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt[4]{{1 - \frac{1}{{10{\text{Da}}}} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}}} \approx \sqrt[4]{{1 - \frac{1}{{10{\text{Da}}_{ 0}^{{}} }} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}}}.$$
(129)

Figure 9 demonstrates comparisons of the calculations by Eq. (58) in combination with one of four roots of Eq. (127) and by Eq. (129). It is evident that both compared solutions agree well with each other.

Fig. 9
figure 9

Comparisons of the numerically and analytically computed normalized Nusselt numbers \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) depending on the Darcy number and the parameter Fh

Full size image

The curves depicted in Figs. 5 and 8 indicate that for a stationary fluid the effects of porosity vanish at lower values of the Darcy number than for the conditions of “no mechanical interaction”.

5.2.2 Condition q w = const

In this case we have

$$\delta^{3} - \frac{1}{10K}\delta^{5} - \frac{{c_{F} g\rho \Delta \rho }}{{120\mu^{2} \sqrt K }}\delta^{7} = \frac{{3q_{w} \mu x}}{{gL_{\upsilon } \rho \Delta \rho }} = \delta_{0}^{3} .$$
(130)

Approximate solutions for the film thickness and the Nusselt number are

$$\delta = \frac{{\delta_{0}^{{}} }}{{\sqrt[3]{{1 - \frac{1}{{10{\text{Da}}}} - \frac{\text{Fh}}{{120{\text{Da}}_{{}}^{ 2} }}}}}} \approx \frac{{\delta_{0}^{{}} }}{{\sqrt[3]{{1 - \frac{1}{{10{\text{Da}}_{ 0}^{{}} }} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}}}}},$$
(131)
$$\frac{\text{Nu}}{{{\text{Nu}}_{ 0} }} = \sqrt[3]{{1 - \frac{1}{{10{\text{Da}}}} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}}} \approx \sqrt[3]{{1 - \frac{1}{{10{\text{Da}}_{ 0}^{{}} }} - \frac{\text{Fh}}{{120{\text{Da}}_{ 0}^{ 2} }}}}.$$
(132)

Figure 10 shows curves calculated by Eqs. (129) and (132). In this case, the effect of the parameter Fh is somewhat stronger than at the condition of qw = const.

Fig. 10
figure 10

Comparisons of the numerically and analytically computed normalized Nusselt numbers \({{\text{Nu}} \mathord{\left/ {\vphantom {{\text{Nu}} {{\text{Nu}}_{ 0} }}} \right. \kern-0pt} {{\text{Nu}}_{ 0} }}\) depending on the Darcy number and the parameter Fh

Full size image

All solutions in frames of the Darcy–Forchheimer–Brinkman approach are summarized in Table 3.

Table 3 Parameters of heat transfer at film boiling of a liquid on a vertical heated wall for the Darcy–Forchheimer–Brinkman approach
Full size table

6 Conclusions

The present work focuses on a study of heat transfer during film boiling of a liquid on a vertical heated wall immersed in a porous medium subject to variation of different parameters of the porous medium and heating conditions at the wall.

Heat transfer peculiarities during film boiling in a porous medium and in the absence of a porous medium were compared against each other. Significant differences are observed for small Darcy numbers Da < 2.

The simulations of heat transfer were performed using a model of a porous medium in the Darcy–Brinkman and Darcy–Brinkman–Forchheimer approximation.

Taking into account the non-linear Forchheimer resistance causes the Nusselt number to tend to zero at large values of permeability (Da), i.e. when the properties of the porous medium tend closer to those of pure liquid.

It was shown that heat transfer intensity during film boiling in a porous medium is weaker than in a free fluid (without porosity) and decreases with the decreasing permeability of the porous medium. In the Darcy–Brinkman model depending on the interaction conditions at the vapor–liquid interface an abrupt decrease in heat transfer was observed at Da < 5.

The predicted values of the heat transfer coefficient for the boundary conditions of the constant wall temperature and constant heat flux on the wall were compared with the data of Cheng and Verma (1981). Our results for the Nusselt numbers demonstrate qualitative agreement with those of Cheng and Verma (1981), although lie below them.

The use of a porous medium model in the Darcy–Brinkman–Forchheimer approximation showed the effect of the Forchheimer parameter on heat transfer during film boiling in a porous medium. An increase in the Forchheimer parameter leads to heat transfer deterioration, which is more significant at small values of the Darcy number. Effects of different thermal boundary conditions on the heated wall on the heat transfer are insignificant.