Flow Stability of Nanofluid Thin Films on Non-Uniformly Heated Porous Slopes

TL;DR

This paper studies how nanofluid films behave on heated porous slopes, showing that nanoparticle concentration can stabilize the flow.

Contribution

The study introduces a new analysis of nanofluid stability on non-uniformly heated porous slopes using nonlinear simulations and stability theory.

Findings

Higher nanoparticle volume fraction stabilizes nanofluid films by increasing the critical Reynolds number.

Porous medium permeability, density difference, and Marangoni number destabilize the flow.

Nonlinear simulations confirm that nanoparticle concentration suppresses disturbance amplitudes.

Abstract

Thin liquid film flows of nanofluids over porous surfaces are central to applications ranging from microfluidic thermal management to precision coating technologies. This study investigates the hydrodynamic and thermal stability of a nanofluid flowing down a non-uniformly heated inclined porous plane subject to the Beavers-Joseph slip boundary condition. Using the long-wave approximation, a nonlinear evolution equation governing the film thickness is derived. The stability characteristics are systematically analyzed via linear stability theory, weakly nonlinear analysis, and fast Fourier transform (FFT) numerical simulations. Quantitative results indicate that the porous medium permeability, density difference, and Marangoni number act as destabilizing factors; specifically, increasing the porous parameter β (from 0 to 0.3), the density ratio ζ0 (from 0 to 5), and the Marangoni number Mn (from 0 to 0.3) significantly reduces the critical Reynolds number and accelerates the onset of interfacial instabilities. In contrast, increasing the nanoparticle volume fraction ϕ from 0 to 0.3 exerts a dominant stabilizing effect by elevating the critical Reynolds number and shrinking the unstable wavenumber domain. Furthermore, nonlinear simulations confirm that higher nanoparticle concentrations effectively suppress the saturation amplitude of disturbances, promoting the eventual stabilization of the liquid film.

Figures

Full text

1. Introduction

Film flow is a ubiquitous phenomenon where a thin, continuous liquid layer moves along a solid surface driven by external forces. It is observed in daily scenarios, such as rainwater sliding on glass [1,2]. In industrial applications, liquid films are pivotal for gas separation, toxic metal recovery, and specifically in the coating industry [3,4]. Crucially, in thermal management systems like microelectronic cooling and thin-film heat exchangers, liquid films serve as efficient agents for heat dissipation and condensation [5,6]. One critical application is liquid film cooling in liquid rocket engines, where coolant is introduced through a porous structure to form a protective film on the combustion chamber wall. Excessive fluctuation of this film can lead to catastrophic engine damage, underscoring the importance of stability analysis.

While the fundamental stability criteria for Newtonian films were established in classical studies [7,8], modern applications often involve complex fluids and multi-physical coupling in porous media. Central to porous wall analysis is the slip boundary condition proposed by Beavers and Joseph [9], later applied by Pascal [10,11] to determine instability conditions for films on porous slopes. However, recent studies demonstrate that non-Newtonian rheology significantly alters flow dynamics in permeable layers. For instance, Mohamad et al. [12] investigated heat and mass transfers on the chemically reactive thermosolutal convective flow of Rivlin–Ericksen fluid over a porous medium, emphasizing the effect of viscous dissipation. Similarly, analytical and numerical simulations were conducted by Yadav et al. [13] to analyze the effects of temperature-reliant thermal conductivity and viscosity disparities on the onset of cellular convective motion in a viscoelastic Oldroyd-B type fluid saturated permeable layer. Furthermore, the onset of Casson fluid convection in a permeable medium layer produced by purely inner heating with a magnetic field was explored by Yadav et al. [14]. Building on these complex interactions, Liu [15] identified three unstable modes (surface, shear, and porous) in film flow, and Hill and Straughan [16] revealed distinct instabilities in fluid-porous layers.

In contrast to classical fluids, nanofluids have garnered attention for their superior thermal conductivity and utility in energy conversion systems [17,18]. While often modeled as homogeneous mixtures, extensive research has established correlations for their thermophysical properties, such as the dependence of dynamic viscosity and effective thermal conductivity on nanoparticle concentration [19,20,21]. Regarding flow dynamics, Lin et al. [22] demonstrated that nanoparticles can suppress linear instability in channel flows, while Zhang et al. [23] analyzed the coupled effects of magnetic fields and velocity slip on power-law nanofluid films. Moreover, thermal convection in a layer of micropolar nanofluid has been investigated by Chand et al. [24] to elucidate the onset of instability and the associated enhancement of heat transfer mechanisms.

Apart from fluid properties, the thermal boundary condition plays a pivotal role in film stability. Since Lin [25] pioneered the study of the Marangoni effect in non-isothermal films, this topic has been extensively explored. Samanta [26,27] investigated film flow on non-uniformly heated inclines, using long-wave perturbation to derive nonlinear surface wave equations that capture both linear and sideband instabilities. Further expanding on complex scenarios, Kirkinis and Andreev [28] examined the stabilizing role of odd viscosity, while Mukhopadhyay et al. [29] revealed that the Marangoni effect destabilizes flow over wavy bottoms, with bottom steepness exhibiting a dual effect on linear stability.

However, existing literature rarely considers the simultaneous interaction of nanofluid rheology, non-uniform heating, and porous substrates under slip conditions. This paper addresses this gap by deriving a nonlinear evolution equation for nanofluid film thickness, incorporating the Beavers-Joseph slip boundary condition and systematic asymptotic expansion. The novelty of this work lies in the integrated analysis of how nanoparticle volume fraction and density difference couple with the slip boundary to alter flow stability. By employing linear stability analysis, weakly nonlinear theory, and numerical simulation, we clarify these complex mechanisms. These findings provide theoretical guidance for optimizing thermal management in micro-electromechanical systems (MEMS) and controlling coating uniformity in precision manufacturing.

2. Mathematical Model

We have conducted a study on the flow of incompressible viscous nanofluids under the drive of gravity within a two-dimensional space. The nanofluids flow downward along an inclined slip plane with non-uniform heating. As depicted in Figure 1, a Cartesian coordinate system $[eqn]$ has been introduced. In this system, the x-axis coincides with the bottom of the plane and is oriented in the direction of the flow, while the y-axis is perpendicular to the inclined plane and points upwards towards the interior of the liquid. The inclination angle of the bottom plate is denoted as $[eqn]$ , and the acceleration resulting from gravity is indicated by g. The position of the free surface of the liquid film is expressed as $[eqn]$ , where $[eqn]$ stands for the thickness of the undisturbed film.

Within the $[eqn]$ coordinate system, the velocity field is defined as $[eqn]$ , with u and v denoting the flow velocities in the horizontal and vertical directions respectively. The dynamic characteristics of incompressible fluids can be depicted by means of the continuity equation and the momentum equation.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

among these parameters, $[eqn]$ and $[eqn]$ denote the densities of the base fluid and the nanofluid respectively. p refers to the liquid pressure, $[eqn]$ stands for the viscosity coefficient of the nanofluid, T indicates the temperature, and $[eqn]$ is the thermal diffusivity coefficient of the nanofluid, which is supposed to be a constant.

In this framework, the nanofluid is treated as a single-phase homogeneous fluid with effective thermophysical properties. While this approach is a standard first approximation in thin-film stability analysis, we acknowledge its inherent limitations. In actual nanofluid systems, effects such as nanoparticle migration, thermophoresis, and Brownian diffusion can be significant under strong thermal gradients or high shear rates. However, for the stable, well-dispersed dispersions at moderate volume fractions considered in the present study, these secondary slip mechanisms are assumed to be negligible compared to the dominant Marangoni and capillary forces. This simplification allows for a focused investigation on the fundamental interfacial dynamics while accounting for the enhanced transport properties of the nanofluid.

The thin film employs the Beavers-Joseph slip boundary condition at the location where the inclined bottom plate has $[eqn]$ . The applicable condition with respect to the tangential velocity at the interface is

$$[eqn]$$

where $[eqn]$ stands for the permeability of the porous medium, $[eqn]$ is a dimensionless parameter relying on the structure of the porous medium, and $[eqn]$ represents the Darcy’s average flow velocity in the x direction. Given that the normal velocity within the boundary layer stays constant, the condition at the interface is

$$[eqn]$$

where $[eqn]$ represents the Darcy’s average flow velocity in the y direction.

When the characteristic length of the pores in the porous medium is significantly smaller than the thickness of the liquid film, the wall condition is simplified through the application of Pascal’s scale analysis as

$$[eqn]$$

This simplification captures the kinematic effect of the porous substrate via the interface slip relationship, while the secondary internal Darcy dynamics are neglected following Pascal’s scale analysis to focus on the film’s interfacial evolution [10,11].

At the substrate with $[eqn]$ , the boundary condition regarding the temperature is

$$[eqn]$$

where $[eqn]$ denotes the ambient temperature. The parameter $[eqn]$ represents the linear rate of temperature variation, where $[eqn]$ . Here, $[eqn]$ and $[eqn]$ refer to the temperatures of the high temperature and low temperature regions of the substrate respectively.

At the free surface with $[eqn]$ , the dynamic and kinematic boundary conditions [29] are given respectively as (The detailed derivation of these equations can be found in Appendix A).

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

here, $[eqn]$ denotes the surface tension, and it is postulated that $[eqn]$ varies linearly within a narrow temperature difference range

$$[eqn]$$

where $[eqn]$ represents the surface tension at the reference temperature $[eqn]$ , and $[eqn]$ is defined as $[eqn]$ . $[eqn]$ is the coefficient of thermal surface tension.

The thermal balance at the deformable free surface $[eqn]$ is governed by Newton’s law of cooling, which characterizes the convective heat exchange between the liquid film and the ambient environment [29]. Specifically, Equation (13) defines a Robin-type boundary condition:

$$[eqn]$$

where $[eqn]$ denotes the thermal conductivity coefficient and $[eqn]$ represents the heat transfer coefficient between the fluid and the air.

The following dimensionless parameters are introduced to nondimensionalize the equations (indicated by the superscript asterisk)

$$[eqn]$$

where $[eqn]$ represents the Nusselt velocity. $[eqn]$ denotes the undisturbed film thickness as well as the transverse length. $[eqn]$ is the characteristic longitudinal length scale, being significantly longer than the film thickness. $[eqn]$ is defined as the aspect ratio, and $[eqn]$ is the average volume fraction of nanoparticles over a certain time interval.

By substituting the dimensionless parameters (14) into the governing equations and boundary conditions and removing the asterisks, we obtain

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where the dimensionless parameters are defined as the Reynolds number $[eqn]$ , the Prandtl number $[eqn]$ , the relative density $[eqn]$ , and the relative viscosity $[eqn]$ :

$$[eqn]$$

The Beavers-Joseph slip boundary condition of the thin film on the interface is expressed as

$$[eqn]$$

where $[eqn]$ is used to denote the dimensionless parameter for the permeability of the porous medium.

The thermal boundary condition at the substrate ( $[eqn]$ ) is simplified as:

$$[eqn]$$

At the free surface where $[eqn]$

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ , and the dimensionless parameters are the Marangoni number $[eqn]$ , the capillary number $[eqn]$ and the Biot number $[eqn]$ , whose forms are respectively

$$[eqn]$$

Expand the physical quantities u, v, p and T into the form of power series of the small parameter $[eqn]$

$$[eqn]$$

Substitute the asymptotic expansion (25) into the dimensionless equations, ignoring the $[eqn]$ and higher-order terms, and the resulting zeroth order governing equation is given by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The boundary conditions at the zeroth order are as follows

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

By solving the system of equations, we obtain

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

This result arises from the assumption that the Biot number is of order $[eqn]$ (i.e., $[eqn]$ ), representing a regime where heat transfer to the ambient gas is negligible compared to the streamwise thermal transport. This scaling argument is consistent with the established framework for non-isothermal film flows on linearly heated substrates, as detailed by Mukhopadhyay and Haldar [30].

Similar to the derivation of the zero-order equation, we substitute the asymptotic expansion into the dimensionless equations again, and thus the first-order governing equation is given by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The boundary conditions at the first order are as follows

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

By solving the first-order equation, the expression for $[eqn]$ is obtained as

$$[eqn]$$

Define the local flow rate $[eqn]$ as

$$[eqn]$$

where $[eqn]$ . Substitute u into Equation (52) and solve the integral to obtain

$$[eqn]$$

Employing an alternative form of the kinematic boundary condition

$$[eqn]$$

We derive the nonlinear evolution equation for the film thickness

$$[eqn]$$

To validate the derivation of the nonlinear evolution equation, we consider the limiting case of a pure Newtonian fluid flowing down a solid, isothermal inclined plane. By setting the nanoparticle volume fraction $[eqn]$ (yielding $[eqn]$ ), the porous permeability parameter $[eqn]$ , and the Marangoni number $[eqn]$ , Equation (55) strictly reduces to the classical Benney equation describing the evolution of a thin liquid film [31,32]. This consistency confirms that the present model correctly captures the fundamental nonlinear hydrodynamics in the Newtonian limit.

3. Linear Stability Analysis

For the investigation of the stability of film flow, the film thickness under the perturbed state can be expressed as $[eqn]$ , where $[eqn]$ denotes the dimensionless unperturbed film thickness. Upon substituting it into Equation (54) and neglecting the higher-order terms of $[eqn]$ and retaining up to $[eqn]$ , we obtain

$$[eqn]$$

where

$$[eqn]$$

among them, A, B, and C as well as their corresponding derivatives denote the values when $[eqn]$ .

The physical property parameters of nanofluids are provided by the subsequent formula, where $[eqn]$ denotes the nanoparticle volume fraction.

$$[eqn]$$

in this paper, the formulas for the relative viscosity and relative density of nanofluids are estimated values derived from theoretical models, which are used to characterize the physical properties of nanofluids under the conditions of this study. Specifically, $[eqn]$ adopts the formula extended by Brinkman [33], which is based on the well-known Einstein equation [34]. Meanwhile, $[eqn]$ has been proven to be appropriate through the experiments of Pak and Cho [35] and is widely accepted. The coefficient $[eqn]$ represents the normal density difference between nanoparticles and the base fluid. Since most types of nanoparticles are heavier than common base fluids, in this study, the value of $[eqn]$ is considered to be in the range of 0 to 5.

By assuming that the perturbation presents in the form of a sine wave, the linear response of the thin film was studied, that is

$$[eqn]$$

where $[eqn]$ is the amplitude of the perturbation, and $[eqn]$ represents the complex conjugate. The real number k denotes the wave number, and $[eqn]$ is the complex frequency. Substituting Equation (59) into Equation (56) and considering its linear part, we obtain the dispersion relation in the sinusoidal mode as

$$[eqn]$$

In Equation (60), the real part and the imaginary part of $[eqn]$ are respectively expressed as

$$[eqn]$$

where $[eqn]$ is the oscillation frequency, and $[eqn]$ is the linear growth rate of the amplitude. If the linear growth rate of the amplitude $[eqn]$ , the flow is linearly unstable; on the contrary, the flow is stable. If $[eqn]$ , the flow is said to be neutrally stable, and the critical Reynolds number at this time can be obtained as

$$[eqn]$$

When we assume that both the relative density and the relative viscosity of the nanofluid are 1, and the porous medium parameter and the Marangoni number are 0, that is, $[eqn]$ and $[eqn]$ , $[eqn]$ , the critical Reynolds number obtained at this time is consistent with the results of Benjamin and Yih.

4. Weakly Non-Linear Analysis

In order to predict the characteristics of thin film flow more precisely, we adopt a multiscale method and expand the interface perturbation $[eqn]$ in the following form:

$$[eqn]$$

where

$$[eqn]$$

here t and x represent the rapidly varying scales, while $[eqn]$ , $[eqn]$ represent the slowly varying scales. Assuming that these variables are independent of each other, the time and space derivatives become

$$[eqn]$$

Substituting Equations (63)–(65) into Equation (56), we can obtain

$$[eqn]$$

where the expressions of the operators $[eqn]$ , $[eqn]$ , $[eqn]$ and the nonlinear terms $[eqn]$ , $[eqn]$ in Equation (66) are as follows

$$[eqn]$$

For the first-order problem $[eqn]$ of $[eqn]$ , we have:

$$[eqn]$$

assume that the solution of this equation has the following form

$$[eqn]$$

where $[eqn]$ represents the complex amplitude. It should be noted that in the linear stability analysis, we have already obtained the solution (67), except that $[eqn]$ is replaced by $[eqn]$ , because near the neutral stability curve, the order of $[eqn]$ is $[eqn]$ . Therefore, the function $[eqn]$ changes slowly and can be absorbed by $[eqn]$ .

For the second-order equation $[eqn]$

$$[eqn]$$

by substituting Equation (68) into Equation (69), we can obtain

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$

To ensure the validity and accuracy of the calculation results, we eliminate the secular term generated by the first term on the right side of Equation (70), thereby obtaining the expression of the solution $[eqn]$ :

$$[eqn]$$

by introducing the coordinate transformation $[eqn]$ , where $[eqn]$ , and utilizing the solvability condition of the third-order equation, the amplitude satisfies the following equation

$$[eqn]$$

where the coefficients $[eqn]$ , $[eqn]$ , and $[eqn]$ are expressed as

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Ignoring the diffusion effect of the second term in Equation (72), we obtain

$$[eqn]$$

Assume that the solution of Equation (73) can be written in the following form

$$[eqn]$$

Substituting Equation (74) into Equation (73), we get

$$[eqn]$$ $$[eqn]$$

The second term on the right side of Equation (75) is a nonlinear term, and this term can adjust the exponential change of the linear perturbation according to the signs of $[eqn]$ and $[eqn]$ . When $[eqn]$ , Equation (75) simplifies to the equation given in the linear theory. When the right side of Equation (75) is zero, the equilibrium amplitude can be obtained as

$$[eqn]$$

When $[eqn]$ , the amplitude reaches saturation; while when $[eqn]$ , the amplitude does not saturate, that is, $[eqn]$ will lead to the instability of the system. According to the signs of $[eqn]$ and $[eqn]$ , the following four nonlinear regions are defined: the supercritical stable region $[eqn]$ ( $[eqn]$ ), the subcritical unstable region $[eqn]$ ( $[eqn]$ ), the unconditionally stable region $[eqn]$ ( $[eqn]$ ), and the explosive state region $[eqn]$ ( $[eqn]$ ).

5. Numerical Simulations

Nonlinear evolution equations are vital mathematical tools for describing complex film dynamics. To solve the governing Equation (55) in a periodic domain, we employ a pseudospectral method based on the Fast Fourier Transform (FFT) algorithm. This approach transforms the spatial variables into a series of wavenumbers, converting the partial differential equation (PDE) into a coupled system of ordinary differential equations (ODEs).

Specifically, the film thickness $[eqn]$ is mapped to the spectral domain as $[eqn]$ , where F and $[eqn]$ denote the forward and inverse FFT operators, respectively. This technique allows spatial derivatives to be computed algebraically as $[eqn]$ . To maintain computational efficiency and avoid expensive convolution sums, nonlinear terms are evaluated in the physical domain through point-wise multiplication after applying $[eqn]$ , before being transformed back to the spectral domain. Applying this procedure to Equation (55) yields the following evolution system:

$$[eqn]$$

For the numerical setup, the computational domain is set to $[eqn]$ , ensuring that multiple wavelengths of the most unstable mode are captured. The domain is discretized with $[eqn]$ grid points, a resolution verified through grid-independence tests to ensure that truncation errors remain negligible. The initial condition is prescribed as a small harmonic perturbation $[eqn]$ , where $[eqn]$ is the most unstable wavenumber determined by linear stability analysis.

Temporal integration of Equation (78) is performed using the MATLAB ode45 solver, an adaptive-step scheme based on the explicit fourth- and fifth-order Runge-Kutta method. Unlike fixed-step methods, ode45 automatically adjusts the internal integration step to satisfy a relative error tolerance of $[eqn]$ , ensuring numerical stability during rapid nonlinear deformation. The time increment for data output is fixed at 0.1, which is sufficient to resolve the intricate interfacial features of the film.

6. Results and Discussion

6.1. Linear Stability Analysis

Figure 2 illustrates the impact of the porous medium’s permeability parameter, $[eqn]$ , on the stability of the nanofluid flow, presenting both the dispersion relation (Figure 2a) and the neutral stability curves (Figure 2b). In Figure 2a, the growth rate of linear perturbations, $[eqn]$ , is plotted against the wave number k. The flow is considered stable when $[eqn]$ , unstable when $[eqn]$ , and neutrally stable at $[eqn]$ . As $[eqn]$ increases, the dispersion curves shift upward, signifying an enhancement of flow instability. Physically, a higher $[eqn]$ corresponds to increased permeability of the porous substrate, which reduces flow resistance and weakens the stabilizing interaction between the fluid and the wall, thereby making the flow more susceptible to disturbances. This destabilizing effect is further corroborated in Figure 2b, where neutral stability curves are plotted in the ( $[eqn]$ ) plane. Here, the region below each curve represents unstable flow conditions. The intersection of each curve with the vertical axis ( $[eqn]$ ) defines the critical Reynolds number ( $[eqn]$ ). The results indicate that an increase in $[eqn]$ leads to a decrease in $[eqn]$ , confirming that higher permeability destabilizes the film flow. This finding aligns with the conclusions reported by Sadiq and Usha [36], reinforcing the validity of the current model. In practical applications such as liquid film cooling, excessive permeability could therefore intensify fluctuations, potentially compromising system performance.

Figure 3 delineates the influence of the nanoparticle volume fraction, $[eqn]$ , on the hydrodynamic stability of the film flow. Figure 3a reveals that the peak growth rate $[eqn]$ decreases monotonically as $[eqn]$ increases, indicating a strong stabilizing effect. This trend is corroborated by the neutral stability curves in Figure 3b, where the critical Reynolds number ( $[eqn]$ ) shifts markedly towards higher values with increasing $[eqn]$ . Physically, this stabilization stems from the modified transport properties of the nanofluid. The addition of nanoparticles elevates the effective viscosity, intensifying the dissipation of kinetic energy within disturbance waves [33,37]. Consequently, the nanofluid system requires significantly larger inertial forces to trigger instability compared to the base fluid.

Figure 4 elucidates the impact of the density ratio parameter, $[eqn]$ , on the hydrodynamic stability of the nanofluid film. Figure 4a presents the dispersion relation, revealing that the peak growth rate $[eqn]$ increases monotonically as $[eqn]$ rises from 0 to 5. This trend indicates that a larger density difference between the nanoparticles and the base fluid exerts a destabilizing effect on the flow. Physically, this phenomenon stems from the alteration of the gravitational body force acting on the fluid. The parameter $[eqn]$ characterizes the relative density contribution of the nanoparticles; a higher $[eqn]$ implies a denser nanofluid mixture. Since the film flow is gravity-driven, an increased effective density amplifies the streamwise component of gravity [7]. This enhanced driving force augments the fluid’s momentum and the associated inertial forces. According to classical stability theory, inertia promotes the amplification of surface waves, while viscosity acts as a stabilizing damping force [8]. Consequently, when the density-induced inertial effects overwhelm viscous damping, the flow becomes more susceptible to disturbances. This conclusion is quantitatively corroborated by the neutral stability curves in Figure 4b, where the critical Reynolds number ( $[eqn]$ ) decreases significantly as $[eqn]$ increases. This reduction in $[eqn]$ confirms that the presence of heavier nanoparticles lowers the threshold for the onset of instability, thereby expanding the unstable domain.

Figure 5 illustrates the influence of the Marangoni number, $[eqn]$ , on the hydrodynamic stability of the film flow. Figure 5a presents the dispersion relation between the growth rate $[eqn]$ and the wavenumber k. It is observed that as $[eqn]$ increases, the dispersion curve shifts noticeably upward, resulting in a higher peak growth rate. This trend indicates that the Marangoni effect acts as a destabilizing factor in the flow system. This conclusion is further corroborated by the neutral stability curves in Figure 5b, where an increase in $[eqn]$ leads to a significant reduction in the stable region and a decrease in the critical Reynolds number ( $[eqn]$ ). Physically, this destabilization is driven by the thermocapillary mechanism [38]. In a heated falling film, a local thinning of the liquid layer (wave trough) results in a higher surface temperature due to its proximity to the heated wall, while thicker regions (wave crests) remain relatively cooler. Since surface tension decreases with increasing temperature, a surface tension gradient is established along the interface. This gradient generates a thermocapillary stress that drags the fluid from the hotter, low-surface-tension troughs toward the cooler, high-surface-tension crests [39]. This mass transport exacerbates the surface deformation, thereby amplifying the disturbances and promoting flow instability [36].

6.2. Weakly Non-Linear Stability Analysis

Figure 6 delineates the stability boundaries in the ( $[eqn]$ ) plane for varying porous medium permeabilities $[eqn]$ . The domain is partitioned into four distinct regimes by the neutral stability curve ( $[eqn]$ ) and the threshold curve where the Landau coefficient vanishes ( $[eqn]$ ). These regimes are identified as follows: the supercritical stable region $[eqn]$ ( $[eqn]$ ), where finite-amplitude waves saturate to a stable equilibrium; the subcritical unstable region $[eqn]$ ( $[eqn]$ ), where the flow is linearly stable but unstable to finite-amplitude disturbances; the unconditionally stable region $[eqn]$ ( $[eqn]$ ); and the explosive instability region $[eqn]$ ( $[eqn]$ ). Comparing Figure 6a,b, it is evident that an increase in permeability $[eqn]$ leads to a reduction in the critical Reynolds number ( $[eqn]$ ). Crucially, increasing $[eqn]$ significantly expands the explosive region $[eqn]$ while shrinking the supercritical stable region $[eqn]$ . Physically, this implies that higher permeability not only lowers the threshold for linear instability but also promotes a more dangerous type of instability where disturbance amplitudes may grow unbounded, potentially accelerating the transition to turbulence. This confirms that the permeability of the inclined plane has a profound destabilizing effect on the flow field in both linear and nonlinear regimes.

Figure 7 presents the weakly nonlinear stability diagrams in the ( $[eqn]$ ) plane for different nanoparticle volume fractions $[eqn]$ . A comparative analysis of Figure 7a ( $[eqn]$ ) and Figure 7b ( $[eqn]$ ) reveals a marked stabilizing trend: the critical Reynolds number ( $[eqn]$ ) increases from approximately 0.24 to 0.32. This shift significantly expands the unconditionally stable region ( $[eqn]$ ). Crucially, from a nonlinear perspective, the explosive instability region ( $[eqn]$ ), where disturbances are predicted to grow unboundedly potentially leading to film rupture, is notably compressed. Physically, this suppression of nonlinear instability is governed by the enhanced energy dissipation mechanism inherent to nanofluids. In the weakly nonlinear regime, the evolution of the disturbance amplitude is dictated by the balance between energy supply and energy dissipation. The addition of nanoparticles elevates the effective viscosity of the fluid [32]. This enhancement intensifies the viscous damping of finite-amplitude waves, effectively counteracting the nonlinear energy transfer that typically drives the system towards explosive growth. Consequently, the threshold for subcritical instability is raised, and the nanofluid film is more likely to saturate into stable finite-amplitude traveling waves rather than undergoing catastrophic rupture, exhibiting superior global stability compared to the pure base fluid [40].

Figure 8 presents the weakly nonlinear stability diagrams in the ( $[eqn]$ ) plane for varying density ratio parameters $[eqn]$ . A comparison between Figure 8a ( $[eqn]$ ) and Figure 8b ( $[eqn]$ ) reveals that the critical Reynolds number ( $[eqn]$ ) decreases noticeably from 0.30 to 0.26. Concurrently, the topology of the stability map undergoes a significant transformation: the unconditionally stable region ( $[eqn]$ ) shrinks, whereas the explosive instability region ( $[eqn]$ ) expands substantially. Physically, this global destabilization is driven by the enhanced inertial effects associated with heavier nanoparticles. An increase in $[eqn]$ signifies a larger density difference between the particles and the base fluid, which amplifies the effective gravitational component driving the flow [7]. This increased driving force results in higher local momentum and stronger inertial forces. In the nonlinear regime, these inertial forces facilitate the energy transfer from the mean flow to the disturbance waves. When this energy input overwhelms the viscous dissipation, finite-amplitude disturbances grow unboundedly. The expansion of the explosive region ( $[eqn]$ ) specifically indicates that the nanofluid film becomes increasingly susceptible to catastrophic rupture even at lower Reynolds numbers, as the stabilizing viscous forces are outmatched by the density-induced inertia [32].

Figure 9 delineates the weakly nonlinear stability diagrams in the ( $[eqn]$ ) plane for varying Marangoni numbers $[eqn]$ . A comparative assessment reveals that as $[eqn]$ increases, the critical Reynolds number ( $[eqn]$ ) decreases significantly. Concurrently, the topology of the stability map transforms drastically: the unconditionally stable region ( $[eqn]$ ) shrinks, while the explosive instability region ( $[eqn]$ ) expands substantially. Physically, this profound destabilization is driven by the intensification of thermocapillary stresses at the free surface. The Marangoni effect induces a flow from the hotter, thinner regions (wave troughs) to the cooler, thicker regions (wave crests) due to surface tension gradients [38]. In the weakly nonlinear regime, this mass transport reinforces surface deformations. When $[eqn]$ is high, the energy input from thermocapillary forces overwhelms the viscous dissipation that typically saturates the wave amplitude. The significant expansion of the explosive region ( $[eqn]$ ) indicates that strong Marangoni effects create a self-amplifying feedback loop, rendering the nanofluid film highly susceptible to unbounded disturbance growth and potential rupture, even at relatively low Reynolds numbers [32].

6.3. Numerical Simulations

Figure 10 illustrates the spatiotemporal evolution of the interfacial waves obtained from the direct numerical simulation for $[eqn]$ within the time span $[eqn]$ . The color gradient represents the progression of time, revealing that the amplitude of the surface elevation $[eqn]$ grows significantly as time evolves. This confirms that at $[eqn]$ , the flow is in the unstable regime, where small perturbations evolve into finite-amplitude traveling waves. Comparing the wave profiles in Figure 10a ( $[eqn]$ ) and Figure 10b ( $[eqn]$ ), a distinct difference in the nonlinear growth is observed. For the solid substrate ( $[eqn]$ ), the wave amplitude reaches a maximum of approximately 1.05. In contrast, for the porous substrate with $[eqn]$ , the wave crests are noticeably amplified, reaching a dimensionless height of roughly 1.10. This numerical evidence aligns perfectly with the predictions from the linear and weakly nonlinear stability analyses presented earlier. It demonstrates that an increase in permeability $[eqn]$ reduces the wall friction and enhances the momentum transfer near the interface, thereby accelerating the growth of disturbances and intensifying the flow instability in the nonlinear regime. This finding is consistent with the results of Sadiq et al. [41], who reported that the presence of a porous substrate and higher permeability significantly increases the amplitude of disturbance waves in heated liquid films.

Figure 11 depicts the spatiotemporal evolution of the interfacial waves for different nanoparticle volume fractions $[eqn]$ at $[eqn]$ . Similar to the observations in Figure 10, the temporal color gradient confirms that the flow is in the unstable nonlinear regime, characterized by the growth of surface disturbances into finite-amplitude waves. A comparative analysis of Figure 11a ( $[eqn]$ ) and Figure 11b ( $[eqn]$ ) reveals a notable impact of the nanoparticle concentration on the wave dynamics. It is observed that the wave amplitude in the nanofluid ( $[eqn]$ ) is significantly larger than that in the pure fluid ( $[eqn]$ ). Specifically, the wave crests in Figure 11b reach higher dimensionless values, indicating a more intense nonlinear interaction. Physically, this phenomenon can be attributed to the competition between viscous damping and inertial forces. Although increasing $[eqn]$ enhances the effective viscosity (a stabilizing factor), the presence of heavy nanoparticles (indicated by a high density ratio $[eqn]$ ) concurrently increases the effective density of the mixture. At the relatively high Reynolds number of $[eqn]$ , the density-induced inertial destabilization dominates over the viscous stabilization. Consequently, a higher concentration of heavy nanoparticles fuels the momentum of the disturbance waves, leading to larger amplitudes and a more unstable interface [37,41].

Figure 12 illustrates the spatiotemporal evolution of the interfacial waves for different density ratio parameters $[eqn]$ at $[eqn]$ . Consistent with the previous cases, the temporal color gradient indicates that the flow is in the unstable nonlinear regime, where disturbances evolve into sustained finite-amplitude traveling waves. A comparison between Figure 12a ( $[eqn]$ ) and Figure 12b ( $[eqn]$ ) reveals a positive correlation between the transient wave amplitude and the density ratio. Specifically, the wave crests for $[eqn]$ exhibit significantly larger dimensionless amplitudes compared to the case with $[eqn]$ . Physically, this amplification is governed by the density-induced inertial mechanism. The parameter $[eqn]$ characterizes the relative density of the nanoparticles; a higher $[eqn]$ corresponds to a heavier nanofluid mixture. Since the flow is gravity-driven, an increase in effective density amplifies the streamwise component of the gravitational body force [7]. This enhanced driving force imparts greater momentum to the fluid, thereby intensifying the inertial forces within the wave crests. In the nonlinear regime, this elevated inertia effectively overcomes the viscous damping, fueling the energy transfer from the mean flow to the disturbance and resulting in larger wave amplitudes [32]. This confirms that the density difference exerts a profound destabilizing influence on the film flow evolution.

Figure 13 depicts the spatiotemporal evolution of the interfacial waves for different Marangoni numbers $[eqn]$ at a Reynolds number of $[eqn]$ . The simulation results clearly indicate that the flow is in the highly nonlinear unstable regime, where infinitesimal disturbances have evolved into large-amplitude traveling waves. A comparison between Figure 13a ( $[eqn]$ ) and Figure 13b ( $[eqn]$ ) reveals the profound destabilizing impact of the Marangoni effect. In the absence of thermocapillary forces ( $[eqn]$ ), the wave maintains a moderate amplitude. However, when $[eqn]$ is increased to 0.8, the wave profiles undergo a dramatic amplification, characterized by sharper crests and significantly larger peak-to-trough heights. Physically, this intensification is driven by the thermocapillary stress induced by surface temperature gradients. In a heated film, the wave troughs are closer to the heated wall and thus hotter than the crests. Consequently, surface tension is lower in the troughs and higher at the crests. This gradient generates a tangential stress that pulls liquid from the troughs toward the crests, reinforcing the surface deformation [38]. At high Reynolds numbers ( $[eqn]$ ), this mechanism acts in concert with inertial forces to overcome viscous damping, leading to the severe amplification of disturbance waves observed in the simulation [39,41].

7. Conclusions

This study provides a comprehensive investigation into the stability of nanofluid thin film flow over a non-uniformly heated, inclined porous substrate subject to a slip boundary condition. A nonlinear evolution equation governing the film thickness was derived utilizing the long-wave approximation and the Beavers-Joseph slip model. By systematically employing linear stability analysis, weakly nonlinear analysis, and direct numerical simulations via the Fast Fourier Transform (FFT) method, we dissected the specific influences of porous medium permeability, nanoparticle volume fraction, density ratio, and the Marangoni number on flow dynamics.

The results from the linear stability analysis reveal distinct competitive mechanisms governing film stability. The porous medium permeability ( $[eqn]$ ), nanoparticle-to-base-fluid density ratio ( $[eqn]$ ), and Marangoni number ( $[eqn]$ ) were identified as destabilizing factors. Specifically, increasing the porous parameter $[eqn]$ from 0 to 0.3 and the Marangoni number $[eqn]$ from 0 to 0.3 significantly amplifies the maximum temporal growth rate of perturbations and expands the unstable wavenumber domain, effectively reducing the critical Reynolds number ( $[eqn]$ ). Similarly, a higher density ratio ( $[eqn]$ increasing from 0 to 5) enhances the instability growth rate. In stark contrast, the nanoparticle volume fraction ( $[eqn]$ ) acts as a dominant stabilizing agent; increasing $[eqn]$ from 0 to 0.3 markedly suppresses the growth rate and shifts the neutral stability curve towards higher Reynolds numbers due to increased effective viscosity. Furthermore, both the weakly nonlinear analysis and the FFT-based numerical simulations show excellent agreement with these linear predictions.

These findings offer critical theoretical guidance for engineering applications dependent on stable thin films, such as micro-scale liquid cooling and precision coating techniques (e.g., CVD or PVD processes). To enhance film uniformity, our model suggests optimizing the nanoparticle concentration and minimizing the density mismatch between particles and the base fluid. Conversely, destabilizing effects can be mitigated by controlling the substrate permeability and managing thermal gradients to limit thermocapillary forces. In summary, this work demonstrates that the stability of nanofluid films on complex substrates can be precisely tailored by balancing these competing physical mechanisms. Future research could extend this framework to consider particle aggregation kinetics, non-Newtonian rheological behaviors, and three-dimensional instability patterns.