Study on Heat Transfer of Fluid in a Porous Media by VOF Method with Fractal Reconstruction

Abstract

This paper addresses the critical gap in the existing literature regarding the combined buoyancy–Marangoni convection of power-law fluids in three-dimensional porous media with complex evaporation surfaces. Previous studies have rarely investigated the convective heat transfer mechanisms in such systems, and there is a lack of effective methods to accurately track fractal evaporation surfaces, which are ubiquitous in natural and engineering porous media (e.g., geological formations, industrial heat exchangers). This research is significant because understanding heat transfer in these complex porous media is essential for optimizing energy systems, enhancing thermal management in industrial processes, and improving the efficiency of phase-change-based technologies. For this scientific issue, a general model is designed. There is a significant temperature difference on the left and right sides of the model, which drives the internal fluid movement through the temperature difference. The upper end of the model is designed as a complex evaporation surface, and there is flowing steam above it, thus forming a coupled flow field. The VOF fractal reconstruction method is adopted to approximate the shape of the complex evaporation surface, which is a major highlight of this study. Different from previous research, this method can more accurately reflect the flow and phase change on the upper surface of the porous medium. Through numerical simulation, the influence of the evaporation coefficient on the flow and heat transfer rate can be determined. Key findings from numerical simulations reveal the following: (1) Heat transfer rates decrease with increasing fractal dimension (surface complexity) and evaporation coefficient; (2) As the thermal Rayleigh number increases, the influence of the Marangoni number on heat transfer diminishes; (3) The coupling of buoyancy and Marangoni effects in porous media with complex evaporation surfaces significantly alters flow and heat transfer patterns compared to smooth-surfaced porous media. This study provides a robust numerical framework for analyzing non-Newtonian fluid convection in complex porous media, offering insights into optimizing thermal systems involving phase changes and irregular surfaces. The findings contribute to advancing heat transfer theory and have practical implications for industries such as energy storage, chemical engineering, and environmental remediation.

1. Introduction

Observation of the Marangoni effect has long been an important perspective in the study of fluid flow. It mainly results from the local interfacial tension gradient changes caused by non-uniform heat transfer, leading to the instability of the phase interface. Qin et al. [1,2] conducted numerical studies on the Marangoni convection of fluids in special media. Zhang et al. [3] investigated the special laws of thermocapillary convection in annular pools. Liu et al. [4] carried out experimental research on the evaporation interface during the fluid phase-change process. In addition, the fractal theory, due to its rich dynamic behavior and complex flow pattern transformation process, has received extensive attention in heat and mass transfer problems. In recent years, many fruitful research results have been achieved on the heat transfer laws of fluids in fractal porous media [5,6].

Porous media are widely present in human society and nature and are applied in various engineering and technical fields. The fluid flow phenomena within them are what first attracted people’s attention. Wang et al. [7], taking the enhanced film condensation of capillary porous structures as the background, started from the spatially discrete porous medium model to study the influence of the local liquid–gas interface shape, curvature, and effective liquid thickness in the porous structure on phase-change performance, and designed the material and surface geometric structure characteristics for related applications. In the literature [8,9], many authors used the Ergun empirical correlation to derive the correlation of pressure drop in foam structures. The key to determining the correlation lies in reliably defining the structural characteristics to determine the equivalent particle diameter in the equation [10,11].

Considering the relevant influencing factors of the condensation process, it is well known that Minkowycz and Sparrow et al. [12] studied the natural convection condensation of steam on an isothermal vertical plate, analyzing the influence of non-condensable gas, interfacial resistance, superheat, temperature, and concentration gradient on the coefficients of each phase through theoretical models, and concluded that the volume concentration of non-condensable gas has a decisive influence on the heat transfer rate when it is small. Siow et al. [13] studied the effects of gas concentration, Reynolds number, inlet pressure, and inlet wall temperature difference on the laminar film condensation of mixed steam, air, and R134a in a horizontal plate channel through numerical methods, proving that the mixed gas working fluid has a significant impact on the condensate film thickness, pressure gradient, and Nusselt number. Ali et al. [14] discussed the influence of the characteristics of the condensation gas–liquid phase, such as liquid viscosity, density, thermal conductivity, and condensation latent heat, on the liquid film distribution and condensation heat transfer coefficient, aiming to improve the heat transfer performance of the natural convection condensation process in horizontal tubes.

Summarizing the previous work, theoretical (numerical or analytical) studies on the combined buoyancy–Marangoni convection of fluids in complex three-dimensional porous media with evaporation surfaces are still scarce. At the same time, there is no suitable method to track the fractal evaporation surface. Therefore, in this study, we consider a comprehensive model to study the combined buoyancy–Marangoni convection of fluids in complex evaporation surface porous media. The fractal reconstruction VOF method is adopted to track the fractal evaporation surface.

The remainder of this paper is structured as follows: Section 2 establishes the mathematical model, including the problem statement, the VOF method with fractal reconstruction for tracking complex evaporation surfaces, governing equations for both external and internal fluid regions, and boundary conditions at the liquid–gas interface. Section 3 describes the numerical methods, covering heat and mass transfer performance metrics (e.g., Nusselt number), computational procedures using the finite volume method, and validation of the model through comparisons with the existing literature. Section 4 presents and discusses the results, analyzing the combined effects of fractal dimension and evaporation coefficient, Marangoni number and thermal Rayleigh number, and Marangoni number and evaporation Biot number on flow and heat transfer. Finally, Section 5 summarizes the key conclusions, highlights the limitations of the study, and outlines potential implications for future research.

2. Mathematical Model

2.1. Problem Statement

To investigate the convection of fluid in a heterogeneous porous cavity, we have constructed a specific model (its physical structure is shown in Figure 1). The left wall of the cavity is at a high temperature, the right wall is at a relatively low temperature, and the top wall is a complex evaporation interface in contact with vapor. For the simulation of the morphology of this complex evaporation interface, the study employs the Weierstrass–Mandelbrot function and the VOF method.

In the porous medium with a complex upper surface, nanofluids are used as a medium to study the convective heat transfer process. The physical model has been shown in Figure 1. A temperature difference is applied to the two side walls of the model, with the left side set at a relatively high temperature and the right side set at a relatively low temperature. The upper surface is a fractal complex surface and remains in contact with the external gas. The method of fractal reconstruction and VOF (Volume of Fluid) is adopted to approximate this complex evaporating surface.

2.2. VOF Method with Fractal Reconstruction

The VOF method proposes the basic idea of phase interface structure, which involves the numerical simulation method of a moving phase interface. Compared with the MAC method, the VOF method saves a lot of computer storage space, requires less computer hardware, and has more obvious advantages when applied to the calculation of 3D problems. The calculation of the VOF method is relatively simple, the phase interface has a high degree of sharpness, and it can represent the structure and changes of the complex phase interface; the VOF method is far superior to other phase interface tracking methods in describing a complex phase interface and dealing with the fusion and fragmentation of a three-dimensional phase interface. For the gas–liquid two-phase interface, ${\epsilon }_l$ is the volume percentage of liquid.

Gas phase: ${\epsilon }_l=0$; liquid phase: ${\epsilon }_l=1$; gas–liquid mixed phase: $0<{\epsilon }_l<1$.

In general, the VOF model traces the gas–liquid interface by solving the volume fraction equation:

(1)${\displaystyle \frac{\partial {\epsilon }_l}{\partial t}}+\overrightarrow{u}\cdot \nabla {\epsilon }_l=0$

where $\overrightarrow{u}$ is velocity; t is time.

However, the traditional VOF method can only track relatively simple and regular interfaces. For the model in this paper, the evaporation surface is a complex parting surface. If it is directly solved by the above equation, the result is rough, inaccurate, and even distorted. The relationship between the evaporation surface and pore fractal dimension, pore size and other parameters of the porous medium cannot be described. Therefore, we must apply fractal reconstruction to meet all the needs of solving similar problems.

2.2.1. W–M Function

We observed the condition of the upper surface of the porous aluminum foam material using an electron microscope (Figure 2). Based on its characteristics, it can be analyzed that it satisfies the self-similarity of fractal geometry, that is, under different scales, there is similarity in its shape distribution [15,16].

Based on the literature research, we can use the Weierstrass–Mandelbrot function to approximately simulate this fractal surface [17]:

(2)$\begin{array}{cc}\hfill \mathrm{z}(\mathrm{x},\mathrm{y})=& L{\left({\displaystyle \frac{G}{L}}\right)}^{D-2}{\left({\displaystyle \frac{\ln \tilde{\gamma }}{M}}\right)}^{{\displaystyle \frac{1}{2}}}{\displaystyle \sum _{m=1}^M{\displaystyle \sum _{n^{*}=n_1^{*}}^{n_{\max }^{*}}{\tilde{\gamma }}^{\left(D-3\right)n^{*}}}}\hfill \\ \hfill & \times \left\{\cos {\varphi }_a-\cos \left[{\displaystyle \frac{2\pi {\tilde{\gamma }}^n{\left(x^2+y^2\right)}^{1/2}}{L}}\cos \left(\arctan \left({\displaystyle \frac{y}{x}}\right)-{\displaystyle \frac{\pi m}{M}}\right)+{\varphi }_a\right]\right\}\hfill \end{array}$

By substituting M = 1, ${\varphi }_a=0$ into Equation (2):

(3)$\mathrm{z}(\mathrm{x})=G^{D-2}{\left(\ln \tilde{\gamma }\right)}^{1/2}{l}^{(3-D)}\left[1-\cos \left({\displaystyle \frac{2\pi x}{l}}\right)\right]\left(0<x<l\right)$

According to fractal geometry, we express the number of protrusions with a characteristic length greater than l as [18,19]:

(4)$N_1\left(L\ge l\right)={\left({\displaystyle \frac{l_{\max }}{l}}\right)}^{D_L}$

Differentiating the above equation with respect to l, we can obtain the following distribution function:

(5)$n_l(l)=-{\displaystyle \frac{d{N}_l}{dl}}=\left(D-1\right)l_{\max }^{D-1}{l}^{-D}$

In the fractal porous media, the cumulative number of pores follows a fractal distribution. Yu [15] used the fractal scaling law to describe the pore distribution in the porous materials, which is

(6)$N_r\left(R\ge r\right)={\left({\displaystyle \frac{r_{\max }}{r}}\right)}^{D_S}$

where $D_S$: pore area fractal dimension, $r_{\max }$: maximum pore size. $N_t={\left(r_{\max }/r_{\min }\right)}^{D_S}$: the total number of pores. Then the number of pores can be expressed as follows:

(7)$-d{N}_r=D_S{r}_{\max }^{D_S}{r}^{-\left(D_S+1\right)}dr$

And,

(8)${\displaystyle \frac{-d{N}_r}{N_t}}=D_S{r}_{\min }^{D_S}{r}^{-\left(D_S+1\right)}dr=f(r)dr$

where $f(r)=D_S{r}_{\min }^{D_S}{r}^{-(D_S+1)}$ is the probability function of the pore size distribution. It also should satisfy the following:

(9)${\displaystyle {\int }_{r_{\min }}^{r_{\max }}f(r)dr}=1-{\left({\displaystyle \frac{r_{\min }}{r_{\max }}}\right)}^{D_S}$

However, based on the previous definitions and combined with probability theory, we can obtain the following:

(10)${\displaystyle {\int }_{r_{\min }}^{r_{\max }}f(r)dr}=1$

Therefore,

(11)${\left({\displaystyle \frac{r_{\min }}{r_{\max }}}\right)}^{D_S}=0$

The above formula can serve as a criterion to determine whether the pores satisfy the fractal distribution. In general, when $r_{\min }/r_{\max }<{10}^{-2}$, we can conclude that the pores of the porous medium material follow a fractal distribution.

In porous media, the effect of surface tension on evaporation surface is different, as shown in Figure 3. In a unit, the surface tension causes the evaporation interface to have an irregular shape. At the same time, the surface tension between the two units makes the evaporation interface connect them to each other. Therefore, the evaporation interface produces a self-similar surface morphology, which can be described by the W–M function.

2.2.2. Fractal Reconstruction

When analyzing the permeation of liquid moisture in porous media, we focus on a control cell located at position z with a thickness of Δz. In the interval from z to z + Δz of this control cell, the pores present different states. Specifically, the pores with a radius smaller than $r_c$ are all filled with liquid, while those with a radius greater than $r_c$ are empty.

When liquid moisture begins to permeate through the porous media, there is a clear sequence for filling the pores. Liquid water will first fill the smallest pores and then move on to fill the slightly larger ones.

It is worth noting that in porous media, the amount of liquid water content has a direct impact on the size of $r_c$. Under normal circumstances, the higher the content of liquid water in the porous media is, the larger the value of $r_c$ will be.

Here, we introduce a parameter s to represent the saturation degree of liquid moisture in porous media. Then, in the context of porous media, there is a corresponding relationship between $r_c$ and s, and this relationship can be expressed as follows [20]:

(12)$s={\displaystyle \frac{{\epsilon }_l}{\varphi }}={\displaystyle \frac{{\displaystyle {\int }_{r_{\min }}^{r_c}\pi r^{2}d{N}_r}}{{\displaystyle {\int }_{r_{\min }}^{r_{\max }}\pi r^{2}d{N}_r}}}={\displaystyle \frac{r_c^{2-D_S}-r_{\min }^{2-D_S}}{r_{\max }^{2-D_S}-r_{\min }^{2-D_S}}}$

where ${\epsilon }_l$ is the volume percentage of liquid.

From Equation (12), we have

(13)$r_c={\left[\left(r_{\max }^{2-D_S}-r_{\min }^{2-D_S}\right){\displaystyle \frac{{\epsilon }_l}{\varphi }}+r_{\min }^{2-D_S}\right]}^{{\displaystyle \frac{1}{2-D_S}}}$

And,

(14)${\displaystyle \frac{\partial r_c}{\partial z}}={\displaystyle \frac{1}{2-D_S}}{\left[\left(r_{\max }^{2-D_S}-r_{\min }^{2-D_S}\right){\displaystyle \frac{{\epsilon }_l}{\varphi }}+r_{\min }^{2-D_S}\right]}^{{\displaystyle \frac{1}{2-D_S}}-1}\left(r_{\max }^{2-D_S}-r_{\min }^{2-D_S}\right){\displaystyle \frac{1}{\varphi }}{\displaystyle \frac{\partial {\epsilon }_l}{\partial z}}$

The pressure difference in capillaries, from z to $z+\Delta z$, can be expressed as follows:

(15)$\Delta p=\left[p_0-{\displaystyle \frac{2\sigma \cos \theta }{r_c(z)+\Delta r_c(z)}}\right]-\left[p_0-{\displaystyle \frac{2\sigma \cos \theta }{r_c(z)}}\right]$

where $p_0$ is the air pressure, $\Delta p$ is the pressure difference in capillaries, $\theta$ is the contact angle. When $\Delta z\to 0$

(16)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial p}{\partial z}}& ={\displaystyle \frac{2\sigma \cos \theta }{r_c^2(z)}}{\displaystyle \frac{\partial r_c(z)}{\partial z}}\hfill \\ \hfill & ={\displaystyle \frac{2\sigma \cos \theta }{2-D_p}}{\left[\left(r_{\max }^{2-D_S}-r_{\min }^{2-D_S}\right){\displaystyle \frac{{\epsilon }_l}{\varphi }}+r_{\min }^{2-D_S}\right]}^{{\displaystyle \frac{1}{2-D_S}}-1}\left(r_{\max }^{2-D_S}-r_{\min }^{2-D_S}\right){\displaystyle \frac{1}{\varphi }}{\displaystyle \frac{\partial {\epsilon }_l}{\partial z}}\hfill \end{array}$

From the above equations, we can get the approximate shape of a gas–liquid interface. Next, with the help of the W–M function, we perform fractal reconstruction of the gas–liquid interface. This process occurs within a time step t, so a virtual time step $t^{*}$ is introduced here.

Considering the evaporation, the function of mass flow rate for liquid film can be expressed as

(17)${\displaystyle \frac{\partial \rho {\epsilon }_l}{\partial t^{*}}}=-{\displaystyle \frac{\partial \left(\dot{m}{\epsilon }_l\right)}{\partial z}}-h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right)$

where $S_k$: the specific surface area of the evaporation surface, $h_{l\leftrightarrow g}$: evaporation coefficient, $C^{*}\left(T\right)$: the saturated vapor concentration, $C_{ev}$: the ambient vapor concentration.

The Hagen–Poiseuille’s law can be used to derive the mean velocity [20]:

(18)$\overline{u}={\displaystyle \frac{n}{3n+1}}{\left(-{\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\left({\displaystyle \frac{\partial p}{\partial z}}\right)}^{{\displaystyle \frac{1}{n}}}\cdot r^{{\displaystyle \frac{1}{n}}+1}$

where n: power law index, $\overline{\mu }$: the fluid dynamic viscosity.

Based on Equation (18), the total mass flow rate of liquid can be expressed as follows:

(19)$\begin{array}{cc}\hfill \dot{M}& ={\displaystyle {\int }_{r_{\min }}^{r_c}\pi r^2\cdot {\rho }_{nf}}\overline{u}\cdot f(r)dr\hfill \\ \hfill & =\rho \cdot \pi {\displaystyle \frac{n}{3n+1}}{\left(-{\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{D}_S{r}_{\max }^{D_S}{\left({\displaystyle \frac{\partial p}{\partial z}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{r_c^{\frac{1}{n}-D_S+3}-r_{\min }^{\frac{1}{n}-D_S+3}}{\frac{1}{n}-D_S+3}}\hfill \end{array}$

The cross-sectional area of the channel $A_c$ can be expressed as

(20)$A_c={\displaystyle {\int }_{r_{\min }}^{r_c}\pi r^2}\cdot f(r)dr=\pi D_S{r}_{\max }^{D_S}{\displaystyle \frac{r_c^{2-D_S}-r_{\min }^{2-D_S}}{2-D_S}}$

The mass flux of the liquid in the vertical direction can be expressed as

(21)$\dot{m}={\displaystyle \frac{\dot{M}}{A_c}}=\rho \cdot {\displaystyle \frac{n}{3n+1}}{\left(-{\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{2-D_S}{\frac{1}{n}-D_S+3}}{\left({\displaystyle \frac{\partial p}{\partial z}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{r_c^{\frac{1}{n}-D_S+3}-r_{\min }^{\frac{1}{n}-D_S+3}}{r_c^{2-D_S}-r_{\min }^{2-D_S}}}$

Combining Equations (17) and (21), we can get

(22)${\displaystyle \frac{\partial {\epsilon }_l}{\partial t^{*}}}={\displaystyle \frac{n}{3n+1}}{\left(-{\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{2-D_S}{\frac{1}{n}-D_S+3}}{\displaystyle \frac{\partial }{\partial z}}\left[{\left({\displaystyle \frac{\partial p}{\partial z}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{r_c^{\frac{1}{n}-D_S+3}-r_{\min }^{\frac{1}{n}-D_S+3}}{r_c^{2-D_S}-r_{\min }^{2-D_S}}}{\epsilon }_l\right]-h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right){\displaystyle \frac{1}{\rho }}$

Due to the $r_{\min }\ll r_{\max }$ Equation (23) can be written as

(23)${\displaystyle \frac{\partial {\epsilon }_l}{\partial t^{*}}}={\displaystyle \frac{n}{3n+1}}{\left(-{\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{2-D_S}{\frac{1}{n}-D_S+3}}{\displaystyle \frac{\partial }{\partial z}}\left[{\left({\displaystyle \frac{\partial p}{\partial z}}\right)}^{{\displaystyle \frac{1}{n}}}{r}_c^{{\displaystyle \frac{1}{n}}+1}{\epsilon }_l\right]-h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right){\displaystyle \frac{1}{\rho }}$

Substituting Equation (16) into Equation (23), we have

(24)${\displaystyle \frac{\partial {\epsilon }_l}{\partial t^{*}}}=D_{lz}\left({\epsilon }_l\right){\displaystyle \frac{{\partial }^2{\epsilon }_l}{\partial z^2}}-h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right){\displaystyle \frac{1}{\rho }}$

where

(25)$D_{lz}\left({\epsilon }_l\right)={\displaystyle \frac{1}{3n+1}}{\left({\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{2-D_S}{\frac{1}{n}-D_S+3}}{\displaystyle \frac{2\sigma \cos \theta r_{\max }^{\frac{1}{n}}}{2-D_S}}{\left(-{\displaystyle \frac{\partial p}{\partial z}}\right)}^{{\displaystyle \frac{1}{n}}-1}{\left({\displaystyle \frac{{\epsilon }_l}{\varphi }}\right)}^{{\displaystyle \frac{\frac{1}{n}}{2-D_S}}}$

Similarly, we have

(26)$D_{lx}\left({\epsilon }_l\right)={\displaystyle \frac{1}{3n+1}}{\left({\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{2-D_S}{\frac{1}{n}-D_S+3}}{\displaystyle \frac{2\sigma \cos \theta r_{\max }^{\frac{1}{n}}}{2-D_S}}{\left(-{\displaystyle \frac{\partial p}{\partial x}}\right)}^{{\displaystyle \frac{1}{n}}-1}{\left({\displaystyle \frac{{\epsilon }_l}{\varphi }}\right)}^{{\displaystyle \frac{\frac{1}{n}}{2-D_S}}}$

(27)$D_{ly}\left({\epsilon }_l\right)={\displaystyle \frac{1}{3n+1}}{\left({\displaystyle \frac{1}{2\overline{\mu }}}\right)}^{{\displaystyle \frac{1}{n}}}{\displaystyle \frac{2-D_S}{\frac{1}{n}-D_S+3}}{\displaystyle \frac{2\sigma \cos \theta r_{\max }^{\frac{1}{n}}}{2-D_S}}{\left(-{\displaystyle \frac{\partial p}{\partial y}}\right)}^{{\displaystyle \frac{1}{n}}-1}{\left({\displaystyle \frac{{\epsilon }_l}{\varphi }}\right)}^{{\displaystyle \frac{\frac{1}{n}}{2-D_S}}}$

In order to find $S_k$, integrating Equation (4), the length of profile $L_p$ is as follows:

(28)$\begin{array}{cc}\hfill L_p& ={\displaystyle {\int }_0^l{G}^{D-2}{\left(\ln \tilde{\gamma }\right)}^{1/2}{l}^{(3-D)}\left[1-\cos \left(\frac{2\pi x}{l}\right)\right]} ds\hfill \\ \hfill & ={\displaystyle {\int }_0^l{G}^{D-2}{\left(\ln \tilde{\gamma }\right)}^{1/2}{l}^{(3-D)}\left[1-\cos \left(\frac{2\pi x}{l}\right)\right]}\sqrt{1+{\sin }^2\left({\displaystyle \frac{2\pi x}{l}}\right)} dx\hfill \\ \hfill & =G^{D-2}{\left(\ln \tilde{\gamma }\right)}^{1/2}{l}^{(3-D)}\cdot {\displaystyle \frac{l}{\pi }}\left[\sqrt{2}+\ln \left(1+\sqrt{2}\right)\right]\hfill \end{array}$

then $S_k$ can be calculated by the following:

(29)$\begin{array}{cc}\hfill S_k& ={{\displaystyle {\int }_0^{l_{\max }}{n}_l(l)\pi \left(\frac{L_p}{2}\right)}}^{2}dl={\displaystyle {\int }_0^{l_{\max }}\left(D-1\right)l_{\max }^{D-1}{l}^{-D}\frac{{\left[\sqrt{2}+\ln \left(1+\sqrt{2}\right)\right]}^2}{4\pi }{G}^{2D-4}\left(\ln \tilde{\gamma }\right)l^{8-2D}}dl\hfill \\ \hfill & ={\displaystyle \frac{{\left[\sqrt{2}+\ln \left(1+\sqrt{2}\right)\right]}^2{G}^{2D-4}\left(\ln \tilde{\gamma }\right)\left(D-1\right)l_{\max }^{8-2D}}{4\pi \cdot \left(9-3D\right)}}\hfill \end{array}$

In this model, the evaporation interface can be regarded as a thin water film whose shape is determined by the fractal structure of the upper surface. Then, we have $l=2r,l_{\max }=2r_c$. Equation (29) can be written as

(30)$\begin{array}{cc}\hfill & S_k={\displaystyle \frac{{\left[\sqrt{2}+\ln \left(1+\sqrt{2}\right)\right]}^2{G}^{2D-4}\left(\ln \tilde{\gamma }\right)\left(D-1\right){\left(2r_c\right)}^{8-2D}}{4\pi \cdot \left(9-3D\right)}}\hfill \\ \hfill & ={\displaystyle \frac{{\left[\sqrt{2}+\ln \left(1+\sqrt{2}\right)\right]}^2{G}^{2D-4}\left(\ln \tilde{\gamma }\right)\left(D-1\right){\left(2r_{\max }\right)}^{8-2D}{\left(\frac{{\epsilon }_l}{\varphi }\right)}^{\frac{8-2D}{2-D_s}}}{4\pi \cdot \left(9-3D\right)}}\hfill \end{array}$

The Equation (5) can be written as

(31)$N_1\left(R\ge r\right)={\left({\displaystyle \frac{2r_c}{2r}}\right)}^{D-1}={\left({\displaystyle \frac{r_c}{r}}\right)}^{D-1}$

In fractal porous materials, the relationship between D and Ds is as follows:

(32)${\left({\displaystyle \frac{r_{\max }}{r_{\min }}}\right)}^{D_S}-{\left({\displaystyle \frac{r_{\max }}{r_c}}\right)}^{D_S}={\left({\displaystyle \frac{r_c}{r_{\min }}}\right)}^{D-1}-1$

where

(33)${\left({\displaystyle \frac{r_{\max }}{r_{\min }}}\right)}^{D_S}-{\left({\displaystyle \frac{{\epsilon }_l}{\varphi }}\right)}^{-{\displaystyle \frac{D_S}{2-D_S}}}={\left({\displaystyle \frac{r_{\max }}{r_{\min }}}{\left({\displaystyle \frac{{\epsilon }_l}{\varphi }}\right)}^{{\displaystyle \frac{1}{2-D_S}}}\right)}^{D-1}-1$

After fractal reconstruction, ${\epsilon }_l$ is the volume percentage of liquid which should satisfy

(34)${\displaystyle \frac{\partial \left(\rho {\epsilon }_l\right)}{\partial t^{*}}}={\displaystyle \frac{\partial }{\partial x}}\left[D_{lx}\left({\epsilon }_l\right){\displaystyle \frac{\partial \left(\rho {\epsilon }_l\right)}{\partial x}}\right]+{\displaystyle \frac{\partial }{\partial y}}\left[D_{ly}\left({\epsilon }_l\right){\displaystyle \frac{\partial \left(\rho {\epsilon }_l\right)}{\partial y}}\right]+{\displaystyle \frac{\partial }{\partial z}}\left[D_{lz}\left({\epsilon }_l\right){\displaystyle \frac{\partial \left(\rho {\epsilon }_l\right)}{\partial z}}\right]-h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right)$

When $t^{*}$ approaches infinity, the fractal gas–liquid interface is generated.

Traditional surface modeling methods fail to capture the self-similar, multi-scale roughness of natural porous media evaporation surfaces, leading to oversimplified phase-change and flow behaviors. In contrast, the Weierstrass–Mandelbrot (W–M) function, as a fractal tool, inherently describes the self-similarity of complex surfaces (verified by SEM images in Figure 2), enabling accurate quantification of relationships between surface complexity (fractal dimension Ds), pore structure, and specific surface area (via Equation (28)). This innovation allows (1) precise tracking of gas–liquid interfaces with irregular morphologies, which traditional VOF methods cannot resolve; (2) direct linkage of surface roughness to evaporation efficiency; and (3) more realistic simulation of coupled buoyancy–Marangoni convection by capturing local temperature gradients induced by surface asperities.

For the VOF model, the physical parameters of the fluid in the governing equations are jointly determined by the phases in each control body. For the gas–liquid interface, the density, viscosity, specific heat capacity, and thermal conductivity in each unit can be given by the following equations:

(35)$\begin{array}{c}{\rho }_m={\rho }_g(1-{\epsilon }_l)+{\rho }_l{\epsilon }_l\hspace{1em}\hspace{1em}{C}_{pm}=C_{pg}(1-{\epsilon }_l)+C_{pl}{\epsilon }_l\\ K_m=K_g(1-{\epsilon }_l)+K_l{\epsilon }_l\hspace{1em}\hspace{1em}{\mu }_m={\mu }_g(1-{\epsilon }_l)+{\mu }_l{\epsilon }_l\end{array}$

2.3. Governing Equations and Boundary Conditions

2.3.1. Governing Equations in the External Fluid Region

(36)$\begin{array}{c}\nabla \cdot \overrightarrow{u}=0\\ {\rho }_g({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \nabla \overrightarrow{u})=-(\nabla p+{\sigma }^{*}\widehat{\kappa }\nabla {\epsilon }_l)+\nabla {\mu }_g(\nabla \overrightarrow{u}+{(\nabla \overrightarrow{u})}^T)-{\rho }_{g}g\overline{Q}\\ {\displaystyle \frac{\partial T}{\partial t}}+{\rho }_g{c}_{pg}(\overrightarrow{v}\cdot \nabla )T=\nabla \cdot [K_g(x,y,z)\nabla T]\end{array}$

where, $\overrightarrow{u}$ is the velocity vector, p is the pressure, t is the time, $\overline{Q}$ is the unit vector in the z direction, g is the acceleration of gravity, ${\rho }_g$ the density of vapor, T is the temperature, and $K_g$ is the effective thermal conductivity of water vapor.

The specific heat of pure water vapor in the above formula is

(37)$c_{pg}=1877.8-0.4417T+1.568\times {10}^{-3}{T}^2-7.286\times {10}^{-7}{T}^3$

2.3.2. Governing Equations in the Internal Fluid Region

(38)$\begin{array}{c}\nabla \cdot \overrightarrow{u}=0\\ {\rho }_l({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \nabla \overrightarrow{u})=-(\nabla p+{\sigma }^{*}\widehat{\kappa }\nabla {\epsilon }_l)+\nabla {\mu }_l(\nabla \overrightarrow{u}+{(\nabla \overrightarrow{u})}^T)-{\displaystyle \frac{{\mu }_l}{K(x,y,z)}}\overrightarrow{u}-\rho g\overline{Q}\\ [(1-\varphi ){(\rho C_p)}_s+\varphi {(\rho C_p)}_l]{\displaystyle \frac{\partial T}{\partial t}}+{(\rho C_p)}_l\overrightarrow{u}\cdot \nabla T=\nabla \cdot [K_l(x,y,z)\nabla T]\end{array}$

where $\overrightarrow{u}$ is the velocity vector, p is the pressure, t is the time, $\overline{Q}$ is the unit vector in the z direction, g is the gravitational acceleration, and $C_{ps}$, $C_{pl}$ are the volumetric specific heat capacities of the porous medium and fluid, respectively, $\varphi$ is the porosity of the porous medium, and T is the temperature, the surface tension coefficient considers the Marangoni effect:

$\overrightarrow{n}=\nabla {\epsilon }_l/\left|{\epsilon }_l\right|,\hspace{1em}\hspace{1em}\widehat{\kappa }=\nabla \cdot \overrightarrow{n}\hspace{1em}\hspace{1em}{\sigma }^{*}={\sigma }_0-{\beta }_T(T-T_0)$

where ${\mathit{\beta }}_T=-{\left.\mathit{\partial }\mathit{\sigma }/\mathit{\partial }T\right|}_C$, ${\sigma }_0$ is the initial surface tension.

Since the permeability is anisotropic, the actual permeability in each direction can be calculated by the following formula [21] (taking the x direction as an example):

(39)$K_x={\displaystyle \frac{1}{32}}{\left({\displaystyle \frac{N_{\max }}{N_{total}}}\right)}^{{\displaystyle \frac{1+D_t}{2}}}{\displaystyle \frac{D_s}{3+D_t-D_s}}{N}_{\max }{A}_0$

where $D_t$ is the tortuosity fractal dimension of the yoz section; image enhancement and the box-counting method can be used to calculate it [21,22,23].

2.3.3. Liquid–Gas Interface Matching Conditions

Considering that the volume flow through the interface is continuous, most previous studies have adopted the following interface velocity and pressure matching conditions:

(40)${\overrightarrow{u}}_l={\overrightarrow{u}}_g{p}_l=p_g$

And consider the evaporation endothermic effect, the energy exchange between the two phases:

(41)$K_m{\displaystyle \frac{\partial T}{\partial z}}=-h_t\left(T-T_{ev}\right)-\lambda h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right)$

$j=h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right)$ is the evaporative mass flow

2.3.4. The Governing Equations of the Flow Field in the Whole Area

(42)$\begin{array}{c}\nabla \cdot \overrightarrow{u}=0\\ {\rho }_m({\displaystyle \frac{\partial \overrightarrow{u}}{\partial t}}+\overrightarrow{u}\cdot \nabla \overrightarrow{u})=-(\nabla p+{\sigma }^{*}\widehat{\kappa }\nabla {\epsilon }_l)-{\alpha }_{H^{*}}{\displaystyle \frac{{\mu }_f}{K_f(x,y,z)}}\overrightarrow{u}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\nabla \cdot {\mu }_m(\nabla \overrightarrow{u}+{(\nabla \overrightarrow{u})}^T)-{\rho }_{m}g\overline{Q}\\ \left([(1-\varphi ){(\rho C_p)}_p+\varphi {(\rho C_p)}_l](1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}\right){\displaystyle \frac{\partial T}{\partial t}}+{\rho }_m{C}_{pm}\overrightarrow{u}\cdot \nabla T\\ =\nabla \cdot [K_m(x,y,z)\nabla T]\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}$

where, $\overrightarrow{u}$ is the velocity vector, p is the pressure, t is the time, $\overline{Q}$ is the unit vector in the z direction, g is the acceleration of gravity, $0<{\epsilon }_l\le 1$: ${\alpha }_{H^{*}}=1$, ${\epsilon }_l=0$: ${\alpha }_{H^{*}}=0$.

The boundary conditions used in this study are as follows:

Left wall: x = 0, T = T_h (high temperature), u= 0 (no-slip condition);

Right wall: Left wall: x = 1, T = T_l (low temperature), u= 0;

Bottom and side walls: z = 0, y = 0, and y = 1 are adiabatic $(\partial T/\partial n=0)$ with u= 0;

Upper evaporation surface: z = 0.3 is coupled with external vapor, satisfying the evaporative mass flux equation $j=h_{l\leftrightarrow g}{S}_k\left(C^{*}\left(T\right)-C_{ev}\right)$ where $C^{*}\left(T\right)$ is the saturated vapor concentration (determined solely by ambient temperature) and $C_{ev}$ is the vapor concentration in the void space.

Physical Interpretation:

The above boundary conditions drive buoyancy convection via temperature differences (left/right walls) and Marangoni convection via surface tension gradients (gas–liquid interface). Evaporation–condensation processes enable coupling between gas and liquid phases. Flow in the porous medium is jointly regulated by pore structure (permeability) and fluid properties (power-law index). The fractal-reconstructed complex evaporation surface enhances gas–liquid exchange by increasing specific surface area, making the model more representative of real-world heat transfer scenarios in natural porous media (e.g., soil, foam materials).

The following dimensionless variables are used:

(43)$\begin{array}{c} (\mathrm{U},\mathrm{V},\mathrm{W})={\displaystyle \frac{(\mathrm{u},\mathrm{v},\mathrm{w})\mathrm{H}}{\mathit{\nu }}},\tau ={\displaystyle \frac{vt}{H^2}},\Theta ={\displaystyle \frac{T-T_0}{T_h-T_l}},\alpha ={\displaystyle \frac{k_0}{{(\rho C_p)}_{nf}}},\mathrm{Pr}={\displaystyle \frac{\mathit{\nu }}{\alpha }},R{a}_T={\displaystyle \frac{g{({\beta }_T)}_f{H}^3(T_h-T_l)}{\mathit{\nu }\alpha }}\\ P={\displaystyle \frac{p{H}^2}{{\rho }_f{v}^2}},\mathit{\nu }={({\displaystyle \frac{u_f}{{\rho }_f}})}^{{\displaystyle \frac{1}{2-n}}}{H}^{{\displaystyle \frac{2(1-n)}{2-n}}},\overline{\mu }={\displaystyle \frac{{\mu }_{eff}}{{\mu }_f{(\frac{v}{H^2})}^{n-1}}} ,Da={\displaystyle \frac{K_0}{H^2}},\end{array}$

where Da is the Darcy number, Pr is the Prandtl number, $R{a}_T$ is the thermal Rayleigh number.

The dimensionless governing equations of the fluid region in porous media are specially introduced here:

(44)${\displaystyle \frac{\partial U}{\partial X}}+{\displaystyle \frac{\partial V}{\partial Y}}+{\displaystyle \frac{\partial W}{\partial Z}}=0$

(45)$\begin{array}{c}\left({\displaystyle \frac{\partial U}{\partial \tau }}+U{\displaystyle \frac{\partial U}{\partial X}}+V{\displaystyle \frac{\partial U}{\partial Y}}+W{\displaystyle \frac{\partial U}{\partial Z}}\right)((1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}E)=-\left({\displaystyle \frac{\partial P}{\partial X}}+{\displaystyle \frac{\partial {\epsilon }_l}{\partial X}}\right)\hfill \\ +((1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}F)\left[\left({\displaystyle \frac{\partial {\tau }_{XX}}{\partial X}}+{\displaystyle \frac{\partial {\tau }_{XY}}{\partial Y}}+{\displaystyle \frac{\partial {\tau }_{XZ}}{\partial Z}}\right)\right.\left.-\overline{\mu }{\displaystyle \frac{U}{Da}}\right]\hfill \end{array}$

(46)$\begin{array}{c}\left({\displaystyle \frac{\partial V}{\partial \tau }}+U{\displaystyle \frac{\partial V}{\partial X}}+V{\displaystyle \frac{\partial V}{\partial Y}}+W{\displaystyle \frac{\partial V}{\partial Z}}\right)((1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}E)=-\left({\displaystyle \frac{\partial P}{\partial Y}}+{\displaystyle \frac{\partial {\epsilon }_l}{\partial Y}}\right)\hfill \\ +((1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}F)\left[\left({\displaystyle \frac{\partial {\tau }_{YX}}{\partial X}}+{\displaystyle \frac{\partial {\tau }_{YY}}{\partial Y}}+{\displaystyle \frac{\partial {\tau }_{YZ}}{\partial Z}}\right)\right.\left.-\overline{\mu }{\displaystyle \frac{V}{Da}}\right]\hfill \end{array}$

(47)$\begin{array}{c}\left({\displaystyle \frac{\partial W}{\partial \tau }}+U{\displaystyle \frac{\partial W}{\partial X}}+V{\displaystyle \frac{\partial W}{\partial Y}}+W{\displaystyle \frac{\partial W}{\partial Z}}\right)((1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}E)=-\left({\displaystyle \frac{\partial P}{\partial Z}}+{\displaystyle \frac{\partial {\epsilon }_l}{\partial Z}}\right)\hfill \\ +((1-{\alpha }_{H^{*}})+{\alpha }_{H^{*}}F)\left[\left({\displaystyle \frac{\partial {\tau }_{ZX}}{\partial X}}+{\displaystyle \frac{\partial {\tau }_{ZY}}{\partial Y}}+{\displaystyle \frac{\partial {\tau }_{ZZ}}{\partial Z}}\right)\right.\left.-\overline{\mu }{\displaystyle \frac{W}{Da}}\right]+{\displaystyle \frac{R{a}_T}{\mathrm{Pr}}}\Theta \hfill \end{array}$

(48)$\gamma {\displaystyle \frac{\partial \Theta }{\partial \tau }}+\left(U{\displaystyle \frac{\partial \Theta }{\partial X}}+V{\displaystyle \frac{\partial \Theta }{\partial Y}}+W{\displaystyle \frac{\partial \Theta }{\partial Z}}\right)={\displaystyle \frac{1}{\mathrm{Pr}}}\left[{\displaystyle \frac{\partial }{\partial X}}\left(K{\displaystyle \frac{\partial \Theta }{\partial X}}\right)+{\displaystyle \frac{\partial }{\partial Y}}\left(K{\displaystyle \frac{\partial \Theta }{\partial Y}}\right)+{\displaystyle \frac{\partial }{\partial Z}}\left(K{\displaystyle \frac{\partial \Theta }{\partial Z}}\right)\right]$

(49)$E=1-\mathit{\delta }+\mathit{\delta }{\displaystyle \frac{{\mathit{\rho }}_s}{{\mathit{\rho }}_f}},F={\displaystyle \frac{1}{{(1-\mathit{\delta })}^{2.5}}}$

where $\mathit{\gamma }={\displaystyle \frac{(1-\mathit{\varphi }){(\mathit{\rho }{C}_p)}_p\mathit{+}\mathit{\varphi }{(\mathit{\rho }{C}_p)}_f}{{(\mathit{\rho }{C}_p)}_f}}$ is heat capacity ratio. For the power-law nanofluids, the stress tensors are expressed as

(50)${\mathit{\tau }}_{XX}=2\overline{\mathit{\mu }}{\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }X}},{\mathit{\tau }}_{YY}=2\overline{\mathit{\mu }}{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Y}},{\mathit{\tau }}_{ZZ}=2\overline{\mathit{\mu }}{\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Z}}$

(51)${\mathit{\tau }}_{XY}={\mathit{\tau }}_{YX}=\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Y}}+{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }X}}\right),{\mathit{\tau }}_{XZ}={\mathit{\tau }}_{ZX}=\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Z}}+{\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }X}}\right),{\mathit{\tau }}_{ZY}={\mathit{\tau }}_{YZ}=\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Y}}+{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Z}}\right)$

(52)$\overline{\mathit{\mu }}=\left[2{\left({\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }X}}\right)}^2+\right.2{\left({\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Y}}\right)}^2+2{\left({\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Z}}\right)}^2{\left.+{\left({\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Y}}+{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }X}}\right)}^2+{\left({\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Z}}+{\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Y}}\right)}^2+{\left({\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }X}}+{\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Z}}\right)}^2\right]}^{(n-1)/2}$

3. Numerical Methods

3.1. Heat and Mass Transfer Performance

The Nusselt number (Nu) on the left wall of the model is used to reflect the heat transfer efficiency [24]:

(53)$\mathrm{N}\mathrm{u}=-{\displaystyle \frac{k_{eff}(x,y,z)}{k_c(x,y,z)}}{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial X}}{|}_{X=0}$

where $k_{eff}(x,y,z)=\mathit{\varphi }{k}_f+(1-\mathit{\varphi })k_p(x,y,z),$ $k_c(x,y,z)=\mathit{\varphi }{k}_f+(1-\mathit{\varphi })k_p(x,y,z)$

The definitions of average Nusselt number $(\overline{Nu})$ is

(54)$\overline{Nu}={\displaystyle {\mathbf{\int }}_{Z=0}^{H^{*}}{\displaystyle {\mathbf{\int }}_{Y=0}^1\left|Nu\right|dZdY,}}$

3.2. Computing Procedure

Finite volume method is adopted to solve numerically the governing equations with the 3D staggered grid arrangement. The spatial discretizations of the governing are

(55)$\mathbf{\left(}{U}_{i,j,k}^{M+1}-U_{i-1,j,k}^{M+1}\mathbf{\right)}\mathbf{\Delta }Y\mathbf{\Delta }Z+\mathbf{\left(}{V}_{i,j,k}^{M+1}-V_{i,j-1,k}^{M+1}\mathbf{\right)}\mathbf{\Delta }X\mathbf{\Delta }Z+\mathbf{\left(}{W}_{i,j,k}^{M+1}-W_{i,j,k-1}^{M+1}\mathbf{\right)}\mathbf{\Delta }X\mathbf{\Delta }Y=0$

(56)$\begin{array}{c}{\displaystyle \frac{\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z}{\mathbf{\Delta }\mathit{\tau }}}(U_{i,j,k}^{M+1}-U_{i,j,k}^M)+{\displaystyle \frac{\overline{\mathit{\mu }}\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z}{Da}}{U}_{i,j,k}^M+(F_{i,j,k}^{(1)}-F_{i-1,j,k}^{(1)})\mathbf{\Delta }Y\mathbf{\Delta }Z\hfill \\ +(G_{i,j,k}^{(1)}-G_{i,j-1,k}^{(1)})\mathbf{\Delta }X\mathbf{\Delta }Z+(H_{i,j,k}^{(1)}-H_{i,j,k-1}^{(1)})\mathbf{\Delta }X\mathbf{\Delta }Y+(P_{i+1,j,k}^{M+1}-P_{i,j,k}^{M+1})\mathbf{\Delta }Y\mathbf{\Delta }Z\hfill \\ +{\mathit{\alpha }}_{H^{*}}({{\widehat{H}}^{*}}_{i+1,j,k}^{M+1}-{{\widehat{H}}^{*}}_{i,j,k}^{M+1})\mathbf{\Delta }Y\mathbf{\Delta }Z=0\hfill \end{array}$

(57)$\begin{array}{c}{\displaystyle \frac{\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z}{\mathbf{\Delta }\mathit{\tau }}}(V_{i,j,k}^{M+1}-V_{i,j,k}^M)+{\displaystyle \frac{\overline{\mathit{\mu }}\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z}{Da}}{V}_{i,j,k}^M+(F_{i,j,k}^{(2)}-F_{i,j-1,k}^{(2)})\mathbf{\Delta }X\mathbf{\Delta }Z\hfill \\ +(G_{i,j,k}^{(2)}-G_{i,j,k-1}^{(2)})\mathbf{\Delta }X\mathbf{\Delta }Y+(H_{i,j,k}^{(2)}-H_{i-1,j,k}^{(2)})\mathbf{\Delta }Y\mathbf{\Delta }Z+(P_{i,j+1,k}^{M+1}-P_{i,j,k}^{M+1})\mathbf{\Delta }X\mathbf{\Delta }Z\hfill \\ +{\mathit{\alpha }}_{H^{*}}({{\widehat{H}}^{*}}_{i,j+1,k}^{M+1}-{{\widehat{H}}^{*}}_{i,j,k}^{M+1})\mathbf{\Delta }X\mathbf{\Delta }Z=0\hfill \end{array}$

(58)$\begin{array}{c}{\displaystyle \frac{\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z}{\mathbf{\Delta }\mathit{\tau }}}(W_{i,j,k}^{M+1}-W_{i,j,k}^M)+{\displaystyle \frac{\overline{\mathit{\mu }}\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z}{Da}}{W}_{i,j,k}^M+(F_{i,j,k}^{(3)}-F_{i,j,k-1}^{(3)})\mathbf{\Delta }X\mathbf{\Delta }Y\hfill \\ +(G_{i-1,j,k}^{(3)}-G_{i,j,k}^{(3)})\mathbf{\Delta }Y\mathbf{\Delta }Z+(H_{i,j,k}^{(3)}-H_{i,j-1,k}^{(3)})\mathbf{\Delta }X\mathbf{\Delta }Z+(P_{i,j,k+1}^{M+1}-P_{i,j,k}^{M+1})\mathbf{\Delta }X\mathbf{\Delta }Y\hfill \\ +{\mathit{\alpha }}_{H^{*}}({{\widehat{H}}^{*}}_{i,j,k+1}^{M+1}-{{\widehat{H}}^{*}}_{i,j,k}^{M+1})\mathbf{\Delta }X\mathbf{\Delta }Y+{\displaystyle \frac{R{a}_T^{*}}{\mathrm{Pr}}}\left({\Theta }_{i,j,k}^{M+1}\right)\mathbf{\Delta }X\mathbf{\Delta }Y\mathbf{\Delta }Z=0\hfill \end{array}$

where

$F^{(1)}=U^2-2\overline{\mathit{\mu }}{\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }X}},G^{(1)}=UV-\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Y}}+{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }X}}\right),H^{(1)}=UW-\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Z}}+{\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }X}}\right)$

$F^{(2)}=V^2-2\overline{\mathit{\mu }}{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Y}},G^{(2)}=WV-\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Z}}+{\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Y}}\right),H^{(2)}=UV-\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }X}}+{\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Y}}\right)$

$F^{(3)}=W^2-2\overline{\mathit{\mu }}{\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Z}},G^{(3)}=WU-\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }X}}+{\displaystyle \frac{\mathbf{\partial }U}{\mathbf{\partial }Z}}\right),H^{(3)}=WV-\overline{\mathit{\mu }}\left({\displaystyle \frac{\mathbf{\partial }W}{\mathbf{\partial }Y}}+{\displaystyle \frac{\mathbf{\partial }V}{\mathbf{\partial }Z}}\right)$

The time terms of the governing equations are discretized by the third-order Runge-Kutta method, i.e.,

(59)${\displaystyle \frac{dU}{d\tau }}=H(U)$

Here

$\begin{array}{c}{U}^{(1)}=U^n+\Delta \tau H(U^n)\hfill \\ U^{(2)}={\displaystyle \frac{3}{4}}{U}^n+{\displaystyle \frac{1}{4}}{U}^{(1)}+{\displaystyle \frac{1}{4}}\Delta \tau H(U^{(1)})\hfill \\ U^{n+1}={\displaystyle \frac{1}{3}}{U}^n+{\displaystyle \frac{2}{3}}{U}^{(1)}+{\displaystyle \frac{2}{3}}\Delta \tau H(U^{(2)})\hfill \end{array}$

We define the situation that meets the following criteria as the achievement of convergence:

(60)${\displaystyle \frac{{\displaystyle \mathbf{\sum }\left|{\phi }_{i,j,k}^m-{\phi }_{i,j,k}^{m-1}\right|}}{{\displaystyle \mathbf{\sum }\left|{\phi }_{i,j,k}^m\right|}}}\le {10}^{-6}$

where m: current iteration, φ: (U, V, W, Θ), (i, j, k): grid points. All the codes are written by the authors.

3.3. Grid Independency and Validation

Grid independence was verified by comparing the average Nusselt number (Nu) on the left wall across three grid densities (Table 1):

The relative error between the medium and fine grids is <2%, indicating that the 80 × 80 × 80 grid sufficiently resolves the flow and temperature fields. Thus, the medium grid was adopted for all simulations to balance accuracy and computational efficiency.

To verify the grid independence, previous research results are selected for comparison. Figure 4 and Figure 5 present the comparison between the isotherms in the X–Z plane obtained by the current computation and the results obtained by Liu et al. [25]. By comparing the isotherms obtained from the test data and the literature data, it can be found that the distribution trends of the temperature gradients are basically consistent.

To quantify the agreement between the current results and those of Liu et al. [25], quantitative error analysis was performed on the isotherms in the X–Z plane. The average absolute error (MAE) between the two datasets was calculated as 2.7%, and the root-mean-square error (RMSE) was 3.1%. These values are well within the typical acceptable error range (<5%) for numerical simulations of convective heat transfer in porous media, confirming the reliability of the present model. The L2 norm of isotherm differences between present results and Liu et al. [25] is 2.9%. In summary, the error stays within acceptable limits, and the actual conditions are consistent with what was expected.

4. Results and Discussion

This section presents the outcomes of numerical simulations concerning the temperature field and Nusselt number. Additionally, an analysis has been conducted on how factors such as the Marangoni number, thermal Rayleigh number, evaporation coefficient, fractal dimension of the upper surface, and temperature field in the porous medium influence power-law fluids. The values of Pr, Da are fixed at 10, ${10}^{-3}$, δ = 0.02, $n=0.86$ respectively. The other parameters are in the range $0\le M{a}_T\le 1000$, ${10}^4\le R{a}_T\le {10}^5$, $1.0\le D_s\le 1.6$, $0.1\le h_{l\leftrightarrow g}\le 0.3$ [24,25].

4.1. Combined Effects of Fractal Dimension and Evaporation Coefficient

Firstly, we investigated the influence of fractal dimension and evaporation coefficient on convective heat transfer of fluid in porous media. The thermal Rayleigh number $(R{a}_T)$ is fixed at ${10}^4$. The Marangoni number $(M{a}_T)$ is set to 100.

Figure 6 shows the heat transfer situation in the X–Z plane of the model’s numerical simulation results. From top to bottom, the evaporation coefficient increases continuously $(h_{l\leftrightarrow g}=0.1,0.2,0.3)$; from left to right, the fractal dimension increases $(\mathrm{D}\mathrm{s}=1.0,1.2,1.6)$, which also means that the unevenness of the upper surface increases. According to the image, it is obvious that for the same fractal dimension, as the evaporation coefficient gradually increases, the blue gradient layers in the porous medium increase, which also means that the heat transfer efficiency decreases. This is because the larger the evaporation coefficient, the more heat from the internal fluid is transferred to the outside, affecting the convective heat transfer of the hot fluid. At the same time, we also found that this trend is more pronounced in the images with higher fractal dimensions. This also means that the fractal dimension has a significant impact on the heat transfer performance.

From another perspective, for the same evaporation coefficient, as the fractal dimension gradually increases (the complexity of the evaporation surface increases), the proportion of the blue gradient layer representing low temperature increases. This is because a higher fractal dimension means a larger gas–liquid interface specific surface area, and the more complex the surface, the more conducive it is to the evaporation and heat absorption process.

Figure 7 shows the three-dimensional distribution diagram of the model. The morphology of 3D isothermal surfaces changes significantly with $h_{l\leftrightarrow g}$ and $Ds$. Under high $Ds$ (e.g., 1.6) and high $h_{l\leftrightarrow g}$ (e.g., 0.3), low-temperature isothermal surfaces are present deeper into the model (maximum indentation depth reaches 30% of the model height), reflecting enhanced local evaporation on complex surfaces. For low $Ds$ (e.g., 1.0) and low $h_{l\leftrightarrow g}$ (e.g., 0.1), isothermal surfaces are flatter, with high-temperature regions occupying 60% of the left model volume, indicating weaker resistance to heat transfer. On the one hand, it confirms the above view; on the other hand, we can see more comprehensively and clearly that the fractal dimension and evaporation coefficient have a significant impact on the convective heat transfer of fluid in porous media.

Figure 8 specifically presents the data visualization. Through the chart, we can find that as the evaporation coefficient increases, the average Nusselt number gradually decreases, which means that the heat transfer efficiency gradually decreases; as the fractal dimension gradually decreases, the heat transfer efficiency becomes higher and higher. For the largest evaporation coefficient and the largest fractal dimension in the test, the average Nusselt number reaches the lowest point. In conclusion, the trend of this chart is completely consistent with the results of the above image analysis. This indicates that the combination of evaporation coefficient and fractal dimension has a huge impact on the convective heat transfer of fluid in porous media.

4.2. Combined Effects of the Marangoni Number and Thermal Rayleigh Number

In this section, we perform a fractal reconstruction of the complex evaporative surface of the model by the VOF method. By observing the change of the interface at different times and different evaporation coefficients, we explore the change of the convective diffusion heat transfer of the fluid inside the porous medium. The effect of the evaporation coefficient on the fractal reconstruction of the evaporation interface is studied.

In order to better observe the changes in the interface as time progresses, we select the (X, Z) plane observation surface of the model for analysis. It can be seen that, in the case of either $h_{l\leftrightarrow g}=0.1$ or $h_{l\leftrightarrow g}=0.2$, the evaporation interface tends to be flat as time goes on. In addition, through the comparison of Figure 9 and Figure 10, it can be found that when $h_{l\leftrightarrow g}=0.1$, the low-temperature distribution on the evaporation interface is narrower and the heat transfer efficiency is higher; when $h_{l\leftrightarrow g}=0.2$, the low-temperature region expands, accounting for a higher proportion of the evaporation surface, the heat transfer efficiency is lower. This is because when the evaporation coefficient is high, the evaporation effect is more obvious, and more heat absorption leads to a greater drop in the temperature of the evaporation interface, which affects the heat and mass transfer efficiency inside the porous medium.

This situation can also be clearly seen on the map of the three-dimensional isothermal surface distribution. As time progresses, in the process of convection, heat is first rapidly introduced into the system from the thermal boundary layer of the high-temperature wall, and then horizontally transferred to the low-temperature wall along the interior of the system, where the heat is rapidly released at the boundary layer. The temperature gradient, concentration gradient and velocity near the hot and cold walls are much larger than those in other areas of the cavity.

4.3. Combined Effects of the Marangoni Number and Evaporation Biot Number

Thirdly, we vary the values of $\mathit{M}{\mathit{a}}_{\mathit{T}}$ from 0 to 1000 and $h_{l\leftrightarrow g}$ from 0.1 to 0.3 to investigate the combined effects of Marangoni number and evaporation coefficient on the gas–liquid interface.

The tangential change of the interface temperature along the interface is presented in Figure 11, and it can be seen that thermal capillary convection increases the gas–liquid interface temperature. Under the same evaporation coefficient, the interface temperature augments with the increase in convection intensity. As mentioned above, the increase in the evaporation coefficient is able to not only lower the temperature of the gas–liquid interface but also reduce the temperature gradient of the interface near the cold wall, thereby reducing the influence of thermal convection on the interface distribution.

The variation of local mass flow j at the interface with Ma and $h_{l\leftrightarrow g}$ is presented in Figure 12. When the thermal capillary convection is extremely weak, most of the interface mass flow augments with the increase of $h_{l\leftrightarrow g}$, such as the curve of Ma = 100. As the thermal capillary convection reaches a certain intensity (Ma = 1000), the variation law of local mass flow j is divided into three regions: Near the cold wall, j increases as $h_{l\leftrightarrow g}$ decreases. This is because under the small $h_{l\leftrightarrow g}$ the thermal capillary convection intensity is the largest, which can not only transport the heat of the interface hot section to the cold end, but also increase the evaporation flow by increasing the interface temperature.

Near the hot wall, the evaporation capacity of the fluid (non-equilibrium effect) plays a major role in the contribution of evaporation flow. Due to the hot wall having a strong heating capacity, the interface temperature of different evaporation modes is extremely close and the liquid evaporation flow with a big $h_{l\leftrightarrow g}$ number is large. In the middle of the interface, the local evaporation flow augments first and then decreases as $h_{l\leftrightarrow g}$ increases under the interaction of thermal capillary convection and the non-equilibrium effect of evaporation. This is because the thermal capillary convection under the small $h_{l\leftrightarrow g}$ plays a leading role in the contribution of the evaporative flow, which raises the interface temperature. Meanwhile, the non-equilibrium effect of evaporation under the large $h_{l\leftrightarrow g}$ number plays a dominant role in the evaporative flow, which reduces the interface temperature. Therefore, there is a moderate evaporating coefficient, which can fully reflect the contribution of them both to the evaporative flow, thereby generating the maximum value of the local evaporative flow.

In addition, the area between each curve and the horizontal axis represents the total evaporation flow rate J (integrated evaporation flow rate) on the phase-change interface. It can be seen from the above analysis that in the weakly heated capillary convective liquid layer, the local mass flow rate at the interface is mostly augmented with the increase in $h_{l\leftrightarrow g}$; therefore, the total evaporation J also increases monotonically with $h_{l\leftrightarrow g}$. In the liquid layer with a certain intensity of thermal capillary convection, the change law of local mass flow j is divided into three different areas, and the local mass flow J in two areas no longer monotonically increases with $h_{l\leftrightarrow g}$. There may exist a large $h_{l\leftrightarrow g}$ liquid layer where the total evaporation flow rate J is smaller than that of the small $h_{l\leftrightarrow g}$ liquid layer. The contribution of thermal capillary convection in the evaporation layer to the total evaporation flow is shown in Figure 12.

5. Conclusions

A three-dimensional numerical investigation is conducted on the combined buoyancy–Marangoni convection of power-law fluids within porous media featuring a complex evaporating surface. The influences of the fractal dimension, evaporation coefficient, Marangoni number, and thermal Rayleigh number are systematically explored. It is found that, as anticipated, the heat transfer rate tends to decline as the fractal dimension and evaporation coefficient increase. Moreover, as the thermal Rayleigh number grows, the impact of the Marangoni number on heat transfer becomes less significant. In brief, when it comes to the effect on heat transfer, the fluid flowing in porous media with a complex evaporating surface is more notable compared to that flowing in porous media without such a surface. This implies that the findings are typical of general non-Newtonian convection problems. However, the model in this study has certain limitations: it assumes that the porous medium is homogeneous, while actual media are usually heterogeneous; numerical diffusion in the VOF method may slightly smooth sharp interfaces; and the model is only applicable to power-law fluids. In short, this research provides a validated numerical tool for analyzing heat transfer in complex porous media, aiding the design of efficient thermal systems such as phase-change heat exchangers and porous media-based energy storage devices. The fractal reconstruction method offers a new approach to modeling irregular surfaces, paving the way for more accurate simulations of natural and industrial processes involving phase changes. Future research can be expanded in several aspects: extending the model to heterogeneous porous media with variable porosity; introducing turbulence models to cover high-Rayleigh-number regions; validating with experimental data on fractal surface evaporation; and expanding the research objects to viscoelastic fluids and reactive flows, so as to further improve the research system and enhance its applicability.

Author Contributions

Conceptualization, S.L. and Q.Z.; methodology, S.L.; software, S.L.; validation, S.L. and W.X.; formal analysis, S.L. and Q.Z.; investigation, S.L.; resources S.L.; data curation, S.L.; writing—original draft preparation, S.L.; writing—review and editing, S.L.; visualization, S.L.; supervision, S.L.; project administration, W.X.; funding acquisition, S.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article as no datasets were generated or used during the research.

Acknowledgments

We would like to thank the School of Aeronautics and Astronautics, Sun Yat-sen University for the various assistance provided in this research.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

$C$	concentration	$C_p$	heat capacity of porous medium
$D_{eff}$	effective mass diffusivity	$\phi$	porosity of the porous media
$g$	gravitational acceleration	$\overline{\mu }$	dimensionless dynamic viscosity
${\beta }_C$	solutal coefficients of surface tension	${\beta }_T$	thermal coefficients of surface tension
$Nu$	local Nusselt number	$\mathrm{\Theta }$	dimensionless temperature
$P$	dimensionless pressure	$Pr$	Prandtl number
$\overline{Q}$	vertical unit vector	$Ra$	porous Rayleigh number
$T$	temperature	$\rho$	density
$p$	pressure	$\alpha$	thermal diffusion coefficient
$\overrightarrow{u}$	velocity vector	$\gamma$	heat capacity ratio
$u$	velocity in x direction	$Da$	Darcy number
$U$	dimensionless velocity in X direction	$x$	x-coordinate
$v$	velocity in y direction	$t$	time
$w$	velocity in z direction	$K$	permeability
$k$	thermal diffusivity	$M$	Marangoni number

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1 Schematic of the physical configuration.

Figure 2 Electron micrographs of the upper surface of porous media at different scales.

Figure 3 The process of fractal reconstruction.

Figure 4 The isotherms (t = 10 s) in X–Z plane obtained by present computation and Liu et al. [25].

Figure 5 The isotherms (t = 30 s) in X–Z plane obtained by present computation and Liu et al. [25].

Figure 6 Effects of fractal dimension and evaporation coefficient on isotherms on the (X, Z) plane.

Figure 7 Effects of fractal dimension and evaporation coefficient on iso-surface of isotherms.

Figure 8 Variations of the average Nusselt number on the left wall with evaporation coefficient.

Figure 9 Temperature distribution of fluid on the mid of (X–Z) plane at different times.

Figure 10 Temperature distribution of fluid inside porous media at different times.

Figure 11 Temperature profiles at vapor–liquid interface at different $h_{l\leftrightarrow g}$ and $M{a}_T$.

Figure 12 Local flux of evaporation at vapor–liquid interface at different $h_{l\leftrightarrow g}$ and $M{a}_T$.

Table 1

Data for grid independence verification.

Grid Density (X × Y × Z)	Average Nusselt Number (Nu)	Relative Error Compared to Finest Grid (%)
40 × 40 × 40 (coarse)	3.12	5.8
80 × 80× 80 (medium)	3.28	1.2
120 × 120 × 120 (fine)	3.32	-

References

1. Qin, T.R.; Grigoriev, R.O. A numerical study of buoyancy-Marangoni convection of volatile binary fluids in confined geometries. Int. J. Heat Mass Transf.; 2018; 127, pp. 308-320. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.06.141]

2. Qin, T.R.; Tukovic, Z.; Grigoriev, R.O. Buoyancy-thermocapillary convection of volatile fluids under atmospheric conditions. Int. J. Heat Mass Transf.; 2014; 75, pp. 284-301. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.03.027]

3. Zhang, L.; Li, Y.R.; Wu, C.M. Effect of Surface Evaporation on Steady Thermocapillary Convection in an Annular Pool, Microgravity. Sci. Technol.; 2016; 28, pp. 499-509.

4. Zhu, Z.Q.; Liu, Q.S. Experimental investigation of thermocapillary convection in a liquid layer with evaporating interface. Chin. Phys. Lett.; 2008; 5, pp. 4046-4049. [DOI: https://dx.doi.org/10.1088/0256-307x/25/11/057]

5. Zhang, Q.; Wu, X.Y.; Meng, Q.B.; Cai, J.C. Fractal models for gas-water transport in shale porous media considering wetting characteristics. Fractals; 2020; 28, 2050138. [DOI: https://dx.doi.org/10.1142/S0218348X20501388]

6. Song, W.H.; Yao, J.; Li, Y.; Sun, H.; Wang, D.Y.; Yan, X. Gas-water relative permeabilities fractal model in dual-wettability multiscale shale porous media during injected water spontaneous imbibition and flow back process. Fractals; 2020; 28, 2050103. [DOI: https://dx.doi.org/10.1142/S0218348X20501030]

7. Wang, R.; Dion, S.A. Effect of meniscus curvature on phase-change performance during capillary-enhanced filmwise condensation in porous media. Front. Therm. Eng.; 2023; 3, 1131363. [DOI: https://dx.doi.org/10.3389/fther.2023.1131363]

8. Lacroix, M.; Nguyen, P.; Schweich, D.; Huu, C.P.; Savin-Poncet, S.; Edouard, D. Pressure drop measurements and modeling on SiC foams. Chem. Eng. Sci.; 2007; 62, pp. 3259-3267. [DOI: https://dx.doi.org/10.1016/j.ces.2007.03.027]

9. Dietrich, B.; Schabel, W.; Kind, M.; Martin, H. Pressure drop measurements of ceramic sponges-Determining the hydraulic diameter. Chem. Eng. Sci.; 2009; 64, pp. 3633-3640. [DOI: https://dx.doi.org/10.1016/j.ces.2009.05.005]

10. Jafarizadeh, A.; Ahmadzadeh, M.; Mahmoudzadeh, S.; Panjepour, M. A New Approach for Predicting the Pressure Drop in Various Types of Metal Foams Using a Combination of CFD and Machine Learning Regression Models. Transp. Porous Media; 2023; 147, pp. 59-91. [DOI: https://dx.doi.org/10.1007/s11242-022-01895-0]

11. Inayat, A.; Schwerdtfeger, J.; Freund, H.; Körner, C.; Singer, R.F.; Schwieger, W. Periodic open-cell foams: Pressure drop measurements and modeling of an ideal tetrakaidecahedra packing. Chem. Eng. Sci.; 2011; 66, pp. 2758-2763. [DOI: https://dx.doi.org/10.1016/j.ces.2011.03.031]

12. Minkowycz, W.J.; Sparrow, E.M. Condensation Heat Transfer in the Presence of Noncondensables, Interfacial Resistance, Superheating, Variable Properties, and Diffusion. Int. J. Heat Mass Transf.; 1966; 9, pp. 1125-1144. [DOI: https://dx.doi.org/10.1016/0017-9310(66)90035-4]

13. Siow, E.C.; Ormiston, S.J.; Soliman, H.M. Fully Coupled Solution of a Two-Phase Model for Laminar Film Condensation of Vapor–Gas Mixtures in Horizontal Channels. Int. J. Heat Mass Transf.; 2002; 45, pp. 3689-3702. [DOI: https://dx.doi.org/10.1016/S0017-9310(02)00094-7]

14. Ali, H.M.; Qasim, M.Z.; Ali, M. Free Convection Condensation Heat Transfer of Steam on Horizontal Square Wire Wrapped Tubes. Int. J. Heat Mass Transf.; 2016; 98, pp. 350-358. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.03.053]

15. Yu, B.M.; Cheng, P. A fractal permeability model for bi-dispersed porous media. Int. J. Heat Mass Transf.; 2002; 45, pp. 2983-2993. [DOI: https://dx.doi.org/10.1016/S0017-9310(02)00014-5]

16. Mandelbrot, B.B. The Fractal Geometry of Nature; Freeman: San Francisco, CA, USA, 1982.

17. Zhang, Q.; Chen, X.Q.; Huang, Y.Y.; Chen, Y. Fractal modeling of fluidic leakage through metal sealing surfaces. AIP Adv.; 2018; 8, 045310. [DOI: https://dx.doi.org/10.1063/1.5023708]

18. Zou, M.Q.; Yu, B.M.; Fu, Y.J.; Xu, P. A Monte Carlo method for simulating fractal surfaces. Physica A; 2007; 386, pp. 176-186. [DOI: https://dx.doi.org/10.1016/j.physa.2007.07.058]

19. Zhu, Q.Y.; Yang, W.B.; Wang, G.H.; Yu, H.Z. A study on capillary actions of power-law fluids in porous fibrous materials via W-M function. Int. J. Heat Mass Transfer.; 2019; 129, pp. 255-264. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.09.080]

20. Abu-Nada, E. Application of nanofluids for heat transfer enhancement of separated flows encountered in a backward facing step. Int. J. Heat Fluid Flow; 2008; 29, pp. 242-249. [DOI: https://dx.doi.org/10.1016/j.ijheatfluidflow.2007.07.001]

21. Zhu, Q.Y.; Yang, W.B.; Yu, H.Z. Study on the permeability of red sandstone via image enhancement. Fractals; 2017; 25, 1750055. [DOI: https://dx.doi.org/10.1142/S0218348X17500554]

22. Zhao, W.D.; Xu, Z.J.; Zhao, J.; Zhao, F.; Han, X.Z. Variational infrared image enhancement based on adaptive dual-threshold gradient field equalization. Infrared Phys. Technol.; 2014; 66, pp. 152-159. [DOI: https://dx.doi.org/10.1016/j.infrared.2014.05.018]

23. Li, J.; Du, Q.; Sun, C.X. An improved boxcounting method for image fractal dimension estimation. Pattern Recognit.; 2009; 42, pp. 2460-2469. [DOI: https://dx.doi.org/10.1016/j.patcog.2009.03.001]

24. Zhuang, Y.J.; Zhu, Q.Y. Numerical study on combined buoyancy–Marangoni convection heat and mass transfer of power-law nanofluids in a cubic cavity filled with a heterogeneous porous medium. Int. J. Heat Fluid Flow.; 2018; 71, pp. 39-54. [DOI: https://dx.doi.org/10.1016/j.ijheatfluidflow.2018.03.006]

25. Liu, S.; Zhu, Q.Y. A Study on Heat and Mass Transfer of Power-Law Nanofluids. In A Fractal Porous Medium with Complex Evaporating Surface. Fractals; 2021; 7, 2150211. [DOI: https://dx.doi.org/10.1142/S0218348X2150211X]