Modeling efficient fractal features to simulate the impact of porosity and viscosity on fluid interfacial stability

Abstract

Efficient fractal theory is presented as a groundbreaking mathematical framework that precisely models the properties of porosity and viscosity in fluid flow, particularly for spinning ferrofluid columns in porous media. The novelty of this study lies in its demonstration that an inviscid fluid in fractal space replicates the behavior of a viscous fluid in traditional space, providing an entirely new perspective on fluid dynamics. This innovative approach leverages fractal theory to capture the intricate, fractal-like features within porous media and their influence on fluid flow under rotation and magnetic fields. The study’s contributions include the development of accurate characterizations for the system’s nonlinear dynamics, which were previously unattainable using traditional methods. A novel mathematical analysis establishes the stability criteria governing the behavior of viscous ferrofluids in porous media, offering key insights into their dynamics. Numerical validation further underscores the interplay between azimuthal magnetic field strength, angular rotation, and porosity, revealing that kinematic viscosity stabilizes the system when the inner fluid density is lower than the outer fluid’s. These findings not only enhance our theoretical understanding but also pave the way for practical advancements in applications involving ferrofluid stability in porous media, emphasizing the transformative potential of fractal theory in fluid dynamics.

Introduction

Fluid motion in fractal spaces is a unique intersection of fluid dynamics and fractal properties, offering insights into fluid behavior in complex, self-similar environments. Fractal spaces have non-integer dimensions and intricate patterns, leading to distinct fluid behaviors. The increased surface area and multiple roughness scales enhance the mixing and transport properties of fluids, resulting in more effective particle and energy dispersal. Anomalist diffusion is observed in fractal spaces, particularly in natural systems like water movement and airflow. Research has shown that fluids in fractal porous media have different permeability and retention characteristics, providing practical applications in environmental science and engineering.^1,2 The relationship between porosity and fractal dimensions is pivotal in understanding the behavior of fluid flow and stability within porous media. Porosity, which quantifies the void spaces in a medium, directly impacts the permeability and fluid transport properties. At the same time, fractal dimensions measure the complexity and irregularity of the porous structure.

Fluid flow in porous media is a crucial aspect of groundwater hydrology, petroleum engineering, and soil science. Characterized by porosity and permeability, these properties influence fluid behavior, dictating flow rates and distribution patterns. ³ Darcy’s law, which relates volumetric flow rate to pressure gradient, fluid viscosity, and medium permeability, is used to understand fluid movement in various applications. ⁴ However, real-world scenarios often involve more complexity due to heterogeneity in pore sizes, anisotropy in permeability, and multiple fluid phases. ⁵ Understanding fluid flow in porous rocks is essential for enhancing oil recovery and optimizing reservoir management. ⁶

The stability of surface waves plays a crucial role in fluid mechanics. These waves are influenced by the interaction between two liquid phases under the influence of electromagnetic fields, affecting their propagation and stability.^7–10 Understanding these interactions is essential for various applications, including environmental engineering, geophysics, petroleum extraction, and coastal engineering. ³ The extended Biot theory describes surface waves in porous media, considering factors like porosity, permeability, and viscosity. ¹¹ Recent computational modeling has provided deeper insights into wave behavior under various conditions, highlighting how different properties of porous media can influence wave propagation. ¹² Understanding surface wave stability in porous media is crucial for interpreting seismic data, assessing earthquake risks, and optimizing techniques in various industries, such as improving oil recovery methods and designing coastal protection structures. ¹³

The rotational azimuthal motion of the fluid column can induce various flow patterns and instabilities, including surface waves. ¹⁴ The mechanisms work together to achieve the spin-up of ferrofluids under the influence of a revolving magnetic field, demonstrating the complex interplay between magnetic, viscous, and hydrodynamic forces in such systems. ¹⁵ Ferrofluid parabolized stability equations, derived from Rosensweig and ferrofluid stability equations, describe the ellipticity of ferrofluid flow and are used to numerically test its interfacial instability in vacuum magnetic fields. ¹⁶ The study of rotating magnetic fluid columns involves analyzing the behavior of magnetic fluids, or ferrofluids, in the presence of rotating magnetic fields. These fluids are colloidal suspensions of magnetic nanoparticles that can respond to external magnetic fields due to their magnetic properties. ¹⁷ When subjected to a rotating magnetic field, these fluids can form stable columns that rotate in sync with the applied field, governed by the interplay between magnetic forces, viscous forces, and the fluid’s intrinsic properties.^9,18 The dynamics of azimuthal rotating magnetic fluid columns are influenced by factors such as the strength and frequency of the rotating magnetic field, and the viscosity of the ferrofluid. ¹⁹ The impact of a magnetic field on the fracture of ferrofluid flow is examined experimentally. ²⁰ At lower frequencies, the fluid column tends to follow the rotating magnetic field closely, while at higher frequencies, a phase lag can develop, leading to complex oscillatory behavior. ²¹

Ji-Huan He has made significant contributions to the field of fractal spaces through his extensive research and publications. His work often focuses on applying fractal calculus and variational principles to various physical and engineering problems. ²² For example, in his work, ²³ A Fractal Variational Theory for One-Dimensional Compressible Flow in Microgravity Space, he introduces a fractal variational approach to analyze compressible flow, shedding light on the distinctive behaviors that arise within fractal geometries. Another notable publication is “On the Fractal Variational Principle for the Telegraph Equation,” he extends traditional variational principles to fractal spaces, providing new insights into wave propagation and stability in these complex structures. ²² Additionally, the fractal undamped Duffing equation introduces a novel method to address nonlinear oscillators within fractal geometries, showcasing the versatility of fractal calculus in handling real-world engineering problems. ²⁴ These publications demonstrate He’s pioneering role in integrating fractal geometry with classical and modern physics to address challenging problems in science and engineering.

Fractal oscillators are complex, chaotic systems found in physical and biological systems. They generate a wide range of frequencies through period-doubling bifurcation, making them unpredictable. ²⁵ Recent advances in computational modeling have provided deeper insights into their behavior, enabling better predictions and control strategies. ²⁶ The study of nonlinear oscillations often employs a homotopy perturbation approach. ²⁷ The harmonic equivalent linearized approach (HELA) simplifies the analysis of nonlinear systems by approximating them with linear models, making it beneficial in engineering fields like structural dynamics and control systems. ²⁸

Recent publications have significantly advanced the understanding of stability in fractal spaces using HELA.^29–31 By leveraging HELA, El-Dib’s research on the damping Duffing-jerk oscillator using fractal models highlights the effectiveness of HELA in solving complex oscillatory systems. His work underscores the importance of these advanced techniques in both theoretical and practical applications, paving the way for future innovations in the field. ³¹ By employing a HELA, El-Dib has been able to delve deeply into the nonlinear characteristics of fractal spaces, bypassing the limitations imposed by small parameter expansions and linear approximations. ³⁰ This method allows for a more thorough examination of the interactions and behaviors within fractal environments, revealing phenomena such as Mathieu–Duffing oscillation that are critical for understanding complex systems. ³² Recent publications have significantly advanced the understanding of the properties of fractal space through the study of the stability of the time-delay nonlinear oscillator in fractal space. ³³

Viscous fluid flow in porous media is a complex phenomenon with significant implications in various scientific and engineering fields, including groundwater hydrology, petroleum engineering, and environmental remediation. ³⁴ The behavior of viscous fluids in porous structures is governed by key parameters such as permeability, porosity, and the viscosity of the fluid, which together influence the flow dynamics and transport properties within the medium. ³⁵ Classical models, such as Darcy’s law, provide a fundamental framework for understanding steady-state flow in homogeneous porous media. However, real-world applications often involve more complex scenarios, including heterogeneous media, non-Newtonian fluid behaviors, and multi-phase flow conditions.

Recent advancements have introduced enhanced modeling techniques, incorporating fractal and fractional calculus to better capture the heterogeneity and anomalous diffusion observed in natural porous media.^36,37 This work discusses an advanced approach to the state of viscous fluid flow in porous media, emphasizing theoretical breakthroughs and computational methods for fractal derivatives. The conclusions are presented, including the fact that an inviscid fluid in fractal space behaves similarly to a viscous fluid in traditional space and that fractal theory can result in an accurate definition of the system’s nonlinear dynamics. The scientific issues related to this study use efficient fractal theory to characterize the porosity and viscosity features of fluid flow. The use of fractal features reduces the mathematical work of solving the Navier–Stokes equation, applying appropriate boundary conditions, and calculating arbitrary constants. ³⁸

The relationship between porosity and fractal dimensions is a cornerstone of the study, as it bridges the physical properties of porous media with the mathematical framework provided by fractal theory. Highlighting this connection more prominently underscores the novelty of the work and its practical implications in modeling and analyzing complex fluid behaviors in heterogeneous porous environments. This paper introduces a pioneering approach to analyzing fluid dynamics by employing efficient fractal theory to model the nuanced properties of porosity and viscosity in fluid flow, specifically in the context of spinning ferrofluid columns in porous media. The key novelty lies in the demonstration that an inviscid fluid in fractal space behaves equivalently to a viscous fluid in traditional space, offering a transformative perspective on fluid behavior and stability analysis. By integrating fractal theory into the study of ferrofluids, the research provides a deeper understanding of how complex, fractal-like structures within porous media interact with rotational and magnetic forces. The study’s contributions include the development of a novel mathematical framework that accurately characterizes nonlinear dynamics, surpassing the limitations of traditional fluid models. Furthermore, the research identifies critical stability criteria influenced by parameters such as azimuthal magnetic field strength, angular rotation, and porosity, while highlighting the stabilizing effect of kinematic viscosity when the inner fluid density is lower than the outer fluid’s density. These findings not only enhance theoretical insights but also hold practical implications for applications requiring precise stability control in ferrofluid systems within porous environments, establishing this work as a significant advancement in the field.

The study of the azimuthal magnetic fluid column here is to apply this novel technique. In addition, we will look at how the density of the inner and outer fluids, the kinematic viscosity, the azimuthal magnetic field strength, the angular rotation parameter, and the porosity coefficient affect how stable viscous ferrofluids are in porous media. This research could be useful for understanding fluid dynamics in complex environments, with practical implications in environmental science and engineering.

A brief overview of the fractal problem

For the majority of magnetic fluid phases, the fundamental equations governing the motion of the rotating inviscid fluid column are stated in the rotating form, as shown in Ref. 39

$\frac{\partial }{\partial t}\underset{¯}{V}+\left(\underset{¯}{V}.\nabla \right)\underset{¯}{V}+2\left(\underset{¯}{\Omega }\wedge \underset{¯}{V}\right)=-\frac{1}{\rho }\nabla \pi -g{\underset{¯}{e}}_z,$

(1)

connected with the incompressibility requirement

$\nabla ·\underset{¯}{V}=0,$

(2)

where $\underset{¯}{V}$ denotes the fluid velocity relative to the rotating frame and $\frac{\partial }{\partial t}\underset{¯}{V}$ denotes the time derivative of the fluid velocity $\underset{¯}{V}=\left(u,v,w\right)$ at the fixed $\underset{¯}{r}$ position in the frame of reference, and the modified fluid pressure is defined as

$\pi =p-\frac{1}{2}\rho {\left(\underset{¯}{\Omega }\wedge \underset{¯}{r}\right)}^2-\frac{1}{2}\epsilon H^2.$

(3)

Referring to the unit vectors along the coordinate axes as $\left({\underset{¯}{e}}_r,{\underset{¯}{e}}_{\theta },{\underset{¯}{e}}_z\right)$ the angular velocity of the fluids, $\underset{¯}{\Omega }=\Omega {\underset{¯}{e}}_z\rho$ refers to fluid density, p hydrostatic pressure, $\pi$ indicates a modified pressure, $\epsilon$ represents the magnetic permeability, H is the magnetic field strength, and g is the gravitational acceleration.

The relationship between porosity and fractal dimensions is pivotal in understanding the behavior of fluid flow and stability within porous media. Porosity, which quantifies the void spaces in a medium, directly impacts the permeability and fluid transport properties. At the same time, fractal dimensions measure the complexity and irregularity of the porous structure.

When the fluid flows, in the fractal space, the governing equation of motion has the fractal derivative form as

$\frac{\partial }{\partial t^{\alpha }}\underset{¯}{V}+\left(\underset{¯}{V}.\nabla \right)\underset{¯}{V}+2\left(\underset{¯}{\Omega }\wedge \underset{¯}{V}\right)=-\frac{1}{\rho }\nabla \pi -g{\underset{¯}{e}}_z;0<\alpha <1,$

(4)

where $\alpha$ is the fractal dimension, $\frac{\partial }{\partial t^{\alpha }}\underset{¯}{V}$ represents the fractal derivative of the fluid velocity, which is He’s fractal derivative^27,40,41 and is defined as follows

$\frac{df\left(t_0\right)}{d{t}^{\alpha }}=\Gamma \left(1+\alpha \right)\underset{\begin{array}{l}t-t_0\to \Delta t,\\ \Delta t\ne 0\end{array}}{\mathit{lim}}\frac{f\left(t\right)-f\left(t_0\right)}{{\left(t-t_0\right)}^2}.$

(5)

This fractal derivative has several qualities like the following

$\underset{\alpha \to 0}{\mathit{lim}}\frac{df}{d{t}^{\alpha }}=f,\underset{\alpha \to 1}{\mathit{lim}}\frac{df}{d{t}^{\alpha }}=\frac{df}{dt},\text{and}\underset{\alpha \to 2}{\hspace{0.17em}\mathit{lim}}\frac{df}{d{t}^{\alpha }}=\frac{df}{d{t}^2}.$

(6)

Fractal calculus is local calculus and it is different from fractional calculus. The fractal derivative described in equation (4) has been successfully used in porous or hierarchical systems.^36,37 The fluid in fractal space acts like a porous medium. ⁴² It is observed that when $\alpha \to 1$ into equation (4) is translated to equation (1), the fluid flows with no porosity. If the fluid comprises porosity structures, equation (1) of fluid motion will contain the relationship of the porosity parameter η and the fluid velocity $\underset{\_}{V}$ as follows

$\frac{\partial \underset{\_}{V}}{\partial t}+\left(\underset{\_}{V}.\nabla \right)\underset{\_}{V}+2\left(\underset{\_}{\Omega }\wedge \underset{\_}{V}\right)=-\frac{1}{\rho }\nabla \pi -\eta \underset{\_}{V}-g\underset{\_}{e}.$

(7)

This article explores the relationship between the fractal dimension α and η = μ/κ, where μ is the kinematic viscosity parameter and κ is medium permeability.

El-Dib and Elgazery created the global He’s fractal derivative of the first order, which changes the original He’s fractal. A smart technique is to stick to the most recent articles,^31–36 in which the following proposition was presented:

$\frac{\partial ..}{\partial t^{\alpha }}={\delta }^{\alpha -1}\hspace{0.17em}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)\frac{\partial ..}{\partial t}+{\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)..,$

(8)

where $\delta$ is a real parameter that depends on the fractal order $\alpha$ . This represents a link between the traditional and fractal derivatives. Combining equations (8) and (4) produces

$\frac{\partial \underset{¯}{V}}{\partial t}+\frac{1}{{\delta }^{\alpha -1}\hspace{0.17em}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)}\left(\left(\underset{¯}{V}.\nabla \right)\underset{¯}{V}+2\left(\underset{¯}{\Omega }\wedge \underset{¯}{V}\right)+\frac{1}{\rho }\nabla \pi +g{\underset{¯}{e}}_z\right)=-\delta \hspace{0.17em}\mathit{cot}\left(\frac{1}{2}\pi \alpha \right)\underset{¯}{V}.$

(9)

A comparison of equations (9) and (7) indicates that

${\delta }^{\alpha -1}\hspace{0.17em}\mathit{sin}\left(\frac{1}{2}\mathit{\pi }\mathit{\alpha }\right)=1\hspace{0.17em}\text{and}\hspace{0.17em}\delta \hspace{0.17em}\mathit{cot}\left(\frac{1}{2}\pi \alpha \right)\equiv \eta .$

(10)

Combing the two equalities of (10) yields the following relation

${\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)=\eta .$

(11)

Employ (10) and (11) into the relation (8) yields

$\frac{\partial ..}{\partial t^{\alpha }}=\frac{\partial ..}{\partial t}+\eta ...$

(12)

The definition above illustrates the relationship between porosity and fractal dimensions by highlighting how fractal dimensions quantify the complexity and irregularity of porous structures, while porosity measures the void spaces that govern fluid transport. The interplay between these two properties is critical for understanding fluid flow dynamics in heterogeneous media. Fractal dimensions offer a framework to model multi-scale structures, capturing how porosity scales across different spatial resolutions and influences stability and flow behavior. This relationship bridges the physical characteristics of porous media with mathematical models, enabling precise analysis of fluid behavior under varying external forces, such as rotation and magnetic fields.

Furthermore, if the kinematic viscosity contribution is included in equation (7) the Navier–Stokes equation arises as

$\frac{\partial \underset{¯}{V}}{\partial t}+\left(\underset{¯}{V}.\nabla \right)\underset{¯}{V}+2\left(\underset{¯}{\Omega }\wedge \underset{¯}{V}\right)=-\frac{1}{\rho }\nabla \pi -\eta \underset{¯}{V}+\mu {\nabla }^2\underset{¯}{V}-g{\underset{¯}{e}}_z.$

(13)

It is convenient to enhance the fractal transformation (8) to have the form

$\frac{\partial ..}{\partial t^{\alpha }}={\delta }^{\alpha -1}\hspace{0.17em}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)\frac{\partial ..}{\partial t}+{\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)\left(1-\mathit{sin}\left(\frac{1}{2}\pi \alpha \right){\nabla }^2\right).$

(14)

Employ this formula to equation (4) yields

$\frac{\partial \underset{¯}{V}}{\partial t}+\frac{1}{{\delta }^{\alpha -1}\hspace{0.17em}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)}\left(\left(\underset{¯}{V}.\nabla \right)\underset{¯}{V}+2\left(\underset{¯}{\Omega }\wedge \underset{¯}{V}\right)+\frac{1}{\rho }\nabla \pi +g{\underset{¯}{e}}_z\right)=-\delta \hspace{0.17em}\mathit{cot}\left(\frac{1}{2}\pi \alpha \right)\underset{¯}{V}+\delta \hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right){\nabla }^2\underset{¯}{V}.$

(15)

As mentioned before, the comparison of (15) with equation (13)

${\delta }^{\alpha -1}\hspace{0.17em}\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)=1,\delta \hspace{0.17em}\mathit{cot}\left(\frac{1}{2}\pi \alpha \right)=\eta \hspace{0.17em}\text{and}\hspace{0.17em}\delta \hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)=\mu .$

(16)

The three equality of (16) can be combined to yield the following relations

$\eta \equiv {\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right),\text{and}\hspace{0.17em}\mu \equiv {\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)$

(17)

According to the above relations, the medium’s permeability is

$\kappa =\mathit{sin}\left(\frac{1}{2}\pi \alpha \right).$

(18)

With these notations the formula (14) becomes

$\frac{\partial ..}{\partial t^{\alpha }}=\frac{\partial ..}{\partial t}+\left(\eta -\mu {\nabla }^2\right)...$

(19)

It is noted that in the limit case as $\alpha \to 1$ we have $\eta =\mu =0$ so that $\frac{\partial ..}{\partial t^{\alpha }}=\frac{\partial ..}{\partial t}$ and therefore the fundamental equation (1) arises. On the other hand, as $\alpha \to 0$ we have the fact $\underset{\alpha \to 0}{\mathit{lim}}\frac{\partial \underset{\_}{V}}{\partial t^{\alpha }}=\underset{¯}{V}$ has been satisfied. This shows that the expression $\left({\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)\underset{¯}{V}\right)$ is referring to the porosity of the fluid medium. ⁴² The equation of motion for inviscid fluids in fractal space can be adapted to describe the behavior of porosity fluids in continuous space. Equation (4) can be used to simplify and expedite the procedures for the analysis of fluids in porous media. Furthermore, the quantity $\left({\delta }^{\alpha }\hspace{0.17em}\mathit{cos}\left(\frac{1}{2}\pi \alpha \right)\mathit{sin}\left(\frac{1}{2}\pi \alpha \right)\right)$ represents kinematic viscosity, as indicated in the Navier–Stokes equation (13).

The definition (19) illustrates the relationship between porosity, kinematic viscosity, and fractal dimensions, emphasizing their combined role in governing fluid dynamics in porous media. Porosity quantifies the void spaces within the medium, directly impacting the permeability and fluid transport. Kinematic viscosity characterizes the fluid’s resistance to flow and its ability to dissipate momentum, which interacts with the porous structure to influence stability and flow behavior. Fractal dimensions provide a mathematical framework to describe the complexity and scale dependence of the porous medium’s structure. Together, these parameters enable a deeper understanding of how intricate geometries, fluid properties, and multi-scale dynamics interact, particularly in systems subject to external forces like rotation and magnetic fields. This relationship forms the foundation for accurately modeling and analyzing the stability and nonlinear behavior of fluids in heterogeneous porous environments.

The physical problem of this research has been expanded to the current work ¹⁹ to provide an analytical solution for the fractal differential equation of the oscillator in a spinning medium or damping force, as in Ref. 43. After using the required boundary conditions described in the references, ¹⁹ the dispersion equation can be produced. This characteristic equation should be expressed in terms of fractal effects, followed by the application of relationship (8) to transform it into a continuous system that expresses porous effects.

Statement of the physical problem

The motion of a magnetized, incompressible fluid column with time-fractal power properties is investigated. In this case, the fluid rotates at a constant angular velocity $\underset{¯}{\Omega }=\Omega {\underset{¯}{e}}_z$ . When referring to the unit vectors along the coordinate axes $\left({\underset{¯}{e}}_r,{\underset{¯}{e}}_{\theta },{\underset{¯}{e}}_z\right)$ , the fluid’s previous angular velocity is thoroughly studied. This approach requires understanding how the rotational velocity, which has remained constant over time, interacts with the fluid’s magnetization and incompressibility. The study looks into how time-fractal power affects fluid motion while considering the unique properties that time’s fractal structure gives it. The historical background of the fluid’s angular velocity reveals information about its behavior and stability under these precise settings.

The azimuthal magnetic field could be expressed as

$\underset{¯}{H}=H_0{\underset{¯}{e}}_{\theta }.$

(20)

The homogeneous azimuthal magnetic field intensity $H_0$ is represented by the unit vector ${\underset{\_}{e}}_{\theta }$ . Accordingly, Maxwell’s equations can be reduced to

$\nabla ·\epsilon \underset{¯}{H}=0,\text{and}\nabla \times \underset{¯}{H}=0.$

(21)

In developing Maxwell’s equations for the situation, we assume that the magneto-quasistatic approximation is valid.^10,44 With a quasi-static model, it is recognized that important temporal rates of change are low enough to ignore contributions from a specific dynamical process.

After justifying the quasi-static approximation, a scalar function $\psi$ representing the magnetic potential may develop, as shown below:

$\underset{¯}{H}=H_0{\underset{¯}{e}}_{\theta }-\nabla \psi .$

(22)

The scalar function $\psi$ fulfills the following Laplace’s equation:

${\nabla }^2\psi =0.$

(23)

In the equilibrium configuration, the interface surface is assumed to have a circular cross-section with a radius R. A schematic sketch of the physical model is provided in Figure 1, illustrating this configuration. This figure helps visualize the system, showing the assumed circular interface which is critical for understanding the equilibrium properties and behavior of the magnetized, incompressible fluid column under study. Around the axis of symmetry, the inner column rotates with a constant angular velocity ${\Omega }_1$ . This fluid column has both density ${\rho }_1$ and magnetic permeability ${\epsilon }_1$ . Surrounding this column is an infinite rotating fluid with density ${\rho }_2$ , magnetic permeability ${\epsilon }_2$ , and a constant angular speed ${\Omega }_2$ . Both fluids are incompressible, and have surface tension at the interface.

Figure 1.

A schematic sketch of the physical model.

Next, assuming a minor variation from the fundamental configuration, the surface deformation is represented by the methodology reported in Ref. 19. According to this analysis, the interfacial wave can be established with a limited amplitude. According to the theory of tiny disturbances, the separation surface can be represented as:

$r=R+\xi \left(\theta ,t\right);\hspace{0.17em}\left(0\le \theta \le 2\pi \right).$

(24)

This representation accounts for the slight perturbations in the system, allowing for an accurate depiction of the interfacial wave’s behavior with a limited amplitude. The initial conditions can be described as:

$\xi \left(\theta ,0\right)=a{e}^{im\theta }+c.c.\&\dot{\xi }\left(\theta ,0\right)=b{e}^{im\theta }+c.c.$

(25)

The integer m is used to denote the azimuthal wavenumber. It is believed to be genuine and positive. The c.c. refers to the complex conjugate of the previous term.

Derivation of the nonlinear characteristic equation using fractal derivatives

When analyzing the equation of motion (4) under the fractal suitable boundary conditions, the following nonlinear characteristic fractional differential equation is produced, as stated in Ref. 19

$\left(1-\frac{m}{R}\frac{\left({\rho }_1-{\rho }_2\right)}{\left({\rho }_1+{\rho }_2\right)}\xi \right)\frac{\partial \xi }{\partial t^{2\alpha }}-\frac{2i}{\left({\rho }_1+{\rho }_2\right)}\left(\left(m-1\right){\rho }_1{\Omega }_1-{\rho }_2{\Omega }_2\right)\frac{\partial \xi }{\partial t^{\alpha }}$

$+\frac{m}{R\left({\rho }_1+{\rho }_2\right)}\left[\frac{m}{R}{H}^2\frac{{\left({\epsilon }_1-{\epsilon }_2\right)}^2}{{\epsilon }_1+{\epsilon }_2}+\frac{mS}{R^2}-R\left(\left(m-2\right){\rho }_1{\Omega }_1^2+\left(m+2\right){\rho }_2{\Omega }_2^2\right)\right]\xi$

$-\frac{m}{R\left({\rho }_1+{\rho }_2\right)}\left[\frac{m^2}{R^2}{H}^2\frac{{\left({\epsilon }_1-{\epsilon }_2\right)}^3}{{\left({\epsilon }_1+{\epsilon }_2\right)}^2}-\frac{3mS}{2R^3}-\left(\left(m-2\right){\rho }_1{\Omega }_1^2-\left(m+2\right){\rho }_2{\Omega }_2^2\right)\right]{\xi }^2-\frac{2im}{R\left({\rho }_1+{\rho }_2\right)}\left({\rho }_1{\Omega }_1-\left(m+1\right){\rho }_2{\Omega }_2\right)\xi \frac{\partial \xi }{\partial t^{\alpha }}-\frac{m^{4}S}{2R^5\left({\rho }_1+{\rho }_2\right)}{\xi }^3=0,$

(26)

with its complex conjugate formula, where S is the amount of surface tension. This equation is a partial differential oscillator constructed by combining the time-fractal derivative’s effects into the equation of motion and appropriate boundary conditions, encapsulating the system’s nonlinear dynamics and complicated behavior.

To solve equation (26) in simple techniques, apply HELA for nonlinear oscillation. ⁴⁵ This method is focused on obtaining the equivalent linearized form of equation (26). To accomplish this, the typical fractional differential equation (26) can be reduced to its simplest form. This is the topic of the next section.

Transformation of the fractional derivative into the traditional derivative

The study of the nonlinear characteristic equation (26) can be streamlined by converting the fractal derivative to a derivative with an integer order form, making analytical solutions easier to implement. Utilizing the relation (12):

$\frac{\partial \xi }{\partial t^{\alpha }}={\xi }_t+\eta \xi ;0<\alpha <1.$

(27)

Therefore, it has

$\frac{\partial \xi }{\partial t^{2\alpha }}={\xi }_{tt}+2\eta {\xi }_t+{\eta }^2\xi .$

(28)

This method expresses the fractal nonlinear oscillator in a more manageable integer-order form, making it easier to analyze and solve. The aforementioned strategy is applied to the fractal nonlinear characteristic problem (26). Introducing a dimensionless formula for private measurements of length, mass, and time such that in private:

$R=R^{*}L,{\rho }_1=\rho {\rho }_2,{\Omega }_1=\Omega {\Omega }_2,{\epsilon }_1=\epsilon {\epsilon }_2,H^2=H^{*2}{\rho }_2{\Omega }_2^2{L}^2/{\epsilon }_2,\eta ={\eta }^{*}{\Omega }_2,\mu ={\mu }^{*}{L}^2{\Omega }_2,\text{and}\hspace{0.17em}L={\left(S/{\rho }_2{\Omega }_2^2\right)}^{1/3}.$

As a result, the fractal nonlinear oscillator (26) should be converted to

${\xi }_{tt}+k{e}^{i{\sigma }_k}{\xi }_t+w{e}^{i{\sigma }_w}\xi +q{e}^{i{\sigma }_q}\xi {\xi }_t+\lambda e^{i{\sigma }_{\lambda }}{\xi }^2+C_0\xi {\xi }_{tt}-C_1{\xi }^3=0,$

(29)

where the real constants, $C_0$ and, $C_1$ represent the coefficients, and the polar form represents the remaining complex coefficients. At this stage, $k,w,q$ and λ is the amplitudes of the polar representation while ${\sigma }^{\prime }s$ referring to the arguments of these representations. These constants are

$k=\sqrt{k_r^2+k_i^2},{\sigma }_k={\mathit{tan}}^{-1}\left(\frac{k_i}{k_r}\right),w=\sqrt{w_r^2+w_i^2},{\sigma }_w={\mathit{tan}}^{-1}\left(\frac{w_i}{w_r}\right),q=\sqrt{q_r^2+q_i^2},{\sigma }_q={\mathit{tan}}^{-1}\left(\frac{q_i}{q_r}\right),\lambda =\sqrt{{\lambda }_r^2+{\lambda }_i^2},{\sigma }_{\lambda }={\mathit{tan}}^{-1}\left(\frac{{\lambda }_i}{{\lambda }_r}\right),$

(30)

where

$k_r=2\eta ,\text{and}{k}_i=-\frac{2\left(\left(m-1\right)\rho \Omega -1\right)}{\left(1+\rho \right)},$

$w_r={\eta }^2+\frac{m}{R\left(\rho +1\right)}\left[\frac{m}{R}{H}^2\frac{{\left(\epsilon -1\right)}^2}{\epsilon +1}+\frac{m}{R^2}-R\left(\left(m-2\right)\rho {\Omega }^2+\left(m+2\right)\right)\right],\&w_i=-\frac{2\eta }{\left(\rho +1\right)}\left(\left(m-1\right)\rho \Omega -1\right),$

$q_r=\frac{2m\left(1-\rho \right)\eta }{R\left(\rho +1\right)},\&q_i=\frac{2m\left(1+m-\rho \Omega \right)}{R\left(\rho +1\right)},$

${\lambda }_r=-\frac{m}{R\left(\rho +1\right)}\left(\frac{m^2}{R^2}{H}^2\frac{{\left(\epsilon -1\right)}^3}{{\left(\epsilon +1\right)}^2}-\frac{3m}{2R^3}-\left(\left(m-2\right)\rho {\Omega }^2-\left(m+2\right)\right)\right)-\frac{m}{R}\frac{\left(\rho -1\right)}{\left(\rho +1\right)}{\eta }^2,{\lambda }_i=\frac{2m}{R\left(\rho +1\right)}\left(1+m-\rho \Omega \right)\eta ,$

$C_0=-\frac{m\left(\rho -1\right)}{R\left(\rho +1\right)},\&C_1=\frac{m^4}{2R^5\left(\rho +1\right)}.$

It is noted that the “*” has been omitted for simplicity.

Adding the complex conjugate form of (29) to the real parameter ξ, which represents wave elevation, produces

${\xi }_{tt}+K{\xi }_t+W\xi +Q\xi {\xi }_t+\Lambda {\xi }^2+C_0\xi {\xi }_{tt}-C_1{\xi }^3=0,$

(31)

where

$K=k\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\sigma }_k,W=w\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\sigma }_w,Q=q\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\sigma }_q,\Lambda =\lambda \hspace{0.17em}\mathit{cos}\hspace{0.17em}{\sigma }_{\lambda }.$

(32)

It is better to simplify the partial differential equation (31) into a basic nonlinear equation in one variable. This can be performed by applying the traveling wave transformation. Introducing a new variable $\tau \left(\theta ,t\right)=m\theta +t$ involves the following transformations:

$\frac{\partial \xi \left(\theta ,t\right)}{\partial t}=\frac{\partial \xi }{\partial \tau }\frac{\partial \tau }{\partial t}=\frac{\partial \xi \left(\tau \right)}{\partial \tau },$

(33)

$\frac{\partial \xi \left(\theta ,t\right)}{\partial \theta }=\frac{\partial \xi }{\partial \tau }\frac{\partial \tau }{\partial \theta }=m\frac{\partial \xi \left(\tau \right)}{\partial \tau }.$

(34)

According to (33) and (34), the nonlinear characteristic equation (31) should be translated into the following ordinary differential equation

${\xi }^{″}+K{\xi }^{\prime }+W\xi +Q\xi {\xi }^{\prime }+\Lambda {\xi }^2+C_0\xi {\xi }^{″}-C_1{\xi }^3=0.$

(35)

The overhead dash shows the derivative of the variable $\tau$ . It should be emphasized that equation (35) describes the Helmholtz–Duffing oscillator with damping effects. It should be emphasized that $\tau \to 0$ it accomplished both $\theta \to 0$ and $t\to 0$ . Therefore, we accept the following initial conditions:

$\xi \left(\tau \to 0\right)=a\hspace{0.17em}\text{and}\hspace{0.17em}{\xi }^{\prime }\left(\tau \to 0\right)=b.$

(36)

The constants a and b consist of the non-zero initial conditions defined in (26).

It is observed that equation (35) can be obtained in the form of the Van der Pol oscillator as

${\xi }^{″}+G\left({\xi }^2\right){\xi }^{\prime }+F\left(\xi \right)=0,$

(37)

where the even function $G\left({\xi }^2\right)$ denotes a nonlinear damping component of the Van der Pol type and $F\left(\xi \right)$ is the restoring force, which is made up of all odd terms and an enhanced quadratic term. Following El-Dib’s publication,^32,45 the coefficients of equation (37) can be found as

$G\left({\xi }^2\right)=K+Q\underset{0}{\overset{\xi }{{\displaystyle \int }}}xdx=K+\frac{1}{2}Q{\xi }^2,$

(38)

$F\left(\xi \right)=W\xi -C_1{\xi }^3+\underset{0}{\overset{\xi }{{\displaystyle \int }}}\left(\Lambda x^2+C_{0}x{x}^{″}\right)dx=W\xi +\left(\frac{1}{3}\Lambda -C_1\right){\xi }^3+\frac{1}{2}{C}_0{\xi }^2{\xi }^{″}.$

(39)

Postulate that the nonlinear oscillation at hand has a complete frequency denoted by $\omega$ and amplitude represented by A.

Representation using the harmonic equivalent linearization approach

The application of HELA simplifies the original nonlinear equation by transforming it into a form that is easier to analyze and solve, taking advantage of the properties of linear damping oscillators to approximate the behavior of the system. Given the HELA, equation (35) can be rearranged into the damped linear harmonic equation. As mentioned in Ref. 45, the equivalent linearized form of equation (37) can be sought in the following linear damping oscillator

${\xi }^{″}+\chi \left(A\right){\xi }^{\prime }+{\varpi }^2\left(A\right)\xi =0.$

(40)

This transformation leverages the simplicity and well-understood nature of linear damping oscillators to provide insights into the dynamics of the original nonlinear system.

The damping coefficient is evaluated ⁴⁵ as

$\chi \left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}{{\xi }^{\prime }}_0^{2}G\left({\xi }_0^2\right)d\tau }{\underset{0}{\overset{T}{{\displaystyle \int }}}{{\xi }^{\prime }}_0^{2}d\tau };T=\frac{2\pi }{\omega }$

(41)

The auxiliary part of equation (40) can be estimated ⁴⁵ as

${\varpi }^2\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}{\xi }_{0}F\left({\xi }_0\right)d\tau }{\underset{0}{\overset{T}{{\displaystyle \int }}}{\xi }_0^{2}d\tau },$

(42)

where ${\xi }_0$ denoted the initial periodic solution, which may be introduced in the form

${\xi }_0\left(\tau \right)=a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\omega \tau +\frac{b}{\omega }\mathit{sin}\hspace{0.17em}\omega \tau ,$

(43)

where $\omega$ signifies the frequency of the oscillating system that needs to be determined. The trial solution (43) has the following simple representation

${\xi }_0\left(\tau \right)=A\hspace{0.17em}\mathit{cos}\left(\omega \tau -\varphi \right),$

(44)

where the amplitude and the phase $\varphi$ of this solution are given by

$A^2=a^2+\frac{b^2}{{\omega }^2}.$

(45)

$\varphi ={\mathit{tan}}^{-1}\left(\frac{b}{a\omega }\right).$

(46)

Combining (44) with (41) and (42) yields

$\chi \left(A\right)=K+\frac{1}{8}Q{A}^2,$

(47)

${\varpi }^2\left(A\right)=W+\frac{1}{8}{A}^2\left(2\Lambda -6C_1-3{\omega }^2{C}_0\right).$

(48)

Equation (40) describes a linear second-order equation with damping effects. The solution can be found in the form

$\xi \left(\tau \right)=\left(a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\omega \tau +\frac{1}{2\omega }\left(2b+a\chi \right)\mathit{sin}\hspace{0.17em}\omega \tau \right)e^{-\frac{1}{2}\chi \tau },$

(49)

where the frequency $\omega$ is evaluated to be

${\omega }^2={\varpi }^2\left(A\right)-\frac{1}{4}{\chi }^2\left(A\right).$

(50)

By combining (47) and (48) into (50), the frequency formula in terms of A is as follows

${\omega }^2=\frac{1}{32\left(8+3A^2{C}_0\right)}\left(256W-{\left(8K+A^{2}Q\right)}^2-64A^2\left(3C_1-\Lambda \right)\right).$

(51)

The aforementioned frequency equation, which contains the parameter A, does not allow for the analysis of the stability behavior. Thus, using (45) and rearranging (51) as a polynomial in ${\omega }^2$ yields

${\omega }^6+{\Gamma }_4{\omega }^4+{\Gamma }_2{\omega }^2+{\Gamma }_0=0,$

(52)

where

${\Gamma }_4=\frac{1}{32\left(8+3a^2{C}_0\right)}\left({\left(8K+a^{2}Q\right)}^2-256W+96b^2{C}_0+64a^2\left(3C_1-\Lambda \right)\right),$

(53)

${\Gamma }_2=\frac{b^2}{16\left(8+3a^2{C}_0\right)}\left(8KQ+a^2{Q}^2-32\Lambda +96C_1\right),$

(54)

${\Gamma }_0=\frac{b^4{Q}^2}{32\left(8+3a^2{C}_0\right)}.$

(55)

This is an explicit frequency equation that can be used to determine the stability criterion.

Stability analysis

The stability criteria require that the characteristic equation (52) have positive roots ${\omega }_i^2>0$ (necessary condition) with real quantities as a sufficient condition.

The necessary condition can be met when

${\Gamma }_4<0,{\Gamma }_2>0\&{\Gamma }_0<0.$

(56)

The adequate condition is met when all possible discernments of the frequency equation (52) are positive. This constraint requires

${\Gamma }_4^2-3{\Gamma }_2>0\hspace{0.17em}\text{and}\hspace{0.17em}{\Gamma }_2^2-3{\Gamma }_0{\Gamma }_4>0.$

(57)

The condition ${\Gamma }_0<0$ can be satisfied when $8+3a^2{C}_0<0$ . Insert the value of $C_0$ yields that this condition can be satisfied when the fluid’s density ratio restricted to

$\rho <\frac{3a^{2}m-8R}{3a^{2}m+8R}.$

(58)

This condition states that the stability behavior demands that the inner fluid’s density be less than that of the outer fluid. This is consistent with the findings reported in Ref. 18. If the last requirement of (56) is met, the stability behavior owing to the middle condition ${\Gamma }_2>0$ and the first condition ${\Gamma }_4<0$ requires that

$8KQ+a^2{Q}^2-32\Lambda +96C_1<0,\text{and}$

(59)

${\left(8K+a^{2}Q\right)}^2-256W+96b^2{C}_0+64a^2\left(3C_1-\Lambda \right)>0.$

(60)

The first and second conditions of (57) revealed that

${\left({\left(8K+a^{2}Q\right)}^2-256W+96b^2{C}_0+64a^2\left(3C_1-\Lambda \right)\right)}^2-192b^2\left(8+3a^2{C}_0\right)\left(8KQ+a^2{Q}^2-32\Lambda +96C_1\right)>0,$

(61)

${\left(8KQ+a^2{Q}^2-32\Lambda +96C_1\right)}^2-\frac{3Q^2}{4}\left({\left(8K+a^{2}Q\right)}^2-256W+96b^2{C}_0+64a^2\left(3C_1-\Lambda \right)\right)>0.$

(62)

Fractal derivative modifications with viscous contributions

We should use enhanced modeling tools that use fractal calculus to simply show the heterogeneity and viscous contribution in natural porous media, as indicated in equations (14) and (19). Currently, the study of the nonlinear characteristic equation (26) can be reformatted to incorporate the viscous contribution by changing the fractal oscillator into a derivative with an integer-order form, using the following modified relations

$\frac{\partial \xi }{\partial t^{\alpha }}=\left(\frac{\partial }{\partial t}+\left(\eta -\mu {\nabla }^2\right)\right)\xi \left(\theta ,t\right)={\xi }_t+\eta \xi +\mu \left(m^2-1\right){\xi }^{-1},$

(63)

$\frac{\partial \xi }{\partial t^{2\alpha }}={\left(\frac{\partial }{\partial t}+\left(\eta -\mu {\nabla }^2\right)\right)}^2\xi \left(\theta ,t\right)={\xi }_{tt}+2\eta {\xi }_t+{\eta }^2\xi +\mu \left({\xi }_t\xi +2\eta \left(m^2-1\right){\xi }^2+\mu {\left(m^2-1\right)}^2\right){\xi }^{-3}.$

(64)

Putting (63) and (64) into equation (26) with the dimensionless formulas described before gets the following singular nonlinear oscillator and its complex conjugate

${\xi }_{tt}+k{e}^{i{\sigma }_k}{\xi }_t+w{e}^{i{\sigma }_w}\xi +\lambda e^{i{\sigma }_{\lambda }}{\xi }^2+q{e}^{i{\sigma }_q}\xi {\xi }_t+C_0\xi {\xi }_{tt}-C_1{\xi }^3+C_0{\mu }^2{\left(m^2-1\right)}^2{\xi }^{-2}+\left(m^2-1\right)k{e}^{i{\sigma }_k}{\xi }^{-1}+\mu {\xi }_t{\xi }^{-2}+C_0\mu {\xi }_t{\xi }^{-1}+\mu \left(m^2-1\right)q{e}^{i{\sigma }_q}+{\mu }^2{\left(m^2-1\right)}^2{\xi }^{-3}=0.$

(65)

Combining the previous equation with its complex conjugate formula provides

${\xi }_{tt}+K{\xi }_t+W\xi +\Lambda {\xi }^2+Q\xi {\xi }_t+C_0\xi {\xi }_{tt}-C_1{\xi }^3+C_0{\mu }^2{\left(m^2-1\right)}^2{\xi }^{-2}+\left(m^2-1\right)K{\xi }^{-1}+\mu {\xi }_t{\xi }^{-2}+C_0\mu {\xi }_t{\xi }^{-1}+\mu \left(m^2-1\right)Q+{\mu }^2{\left(m^2-1\right)}^2{\xi }^{-3}=0.$

(66)

Converting quadratic nonlinearity to odd nonlinearity, as described in Ref. 46, provides

${\xi }_{tt}+\left(\mu {\xi }^{-2}+K+\frac{1}{2}Q{\xi }^2\right){\xi }_t+\left(W+\mu \left(m^2-1\right)Q\right)\xi +\left(\frac{1}{3}\Lambda -C_1\right){\xi }^3+\frac{1}{2}{C}_0{\xi }^2{\xi }_{tt}+\left(m^2-1\right)\left(K-C_0{\mu }^2\left(m^2-1\right)\right){\xi }^{-1}+{\mu }^2{\left(m^2-1\right)}^2{\xi }^{-3}+\frac{1}{2}{C}_0\mu {{\xi }_t}^2{\xi }^{-1}=0.$

(67)

This is referred to as a singular Van-der-Pol oscillator. This equation, using the transformation given above in (33) and (34) can be expressed symbolically as

${\xi }^{″}+\frac{1}{{\xi }^2}{F}_1\left({\xi }^2\right){\xi }^{\prime }+\frac{1}{{\xi }^3}{F}_0\left(\xi \right)=0,$

(68)

where the nonlinear functions $F_i\left(\xi \right)$ are listed below

$F_1\left({\xi }^2\right)=\mu +K{\xi }^2+\frac{1}{2}Q{\xi }^4,$

(69)

$F_0\left(\xi \right)=\left(\frac{1}{3}\Lambda -C_1\right){\xi }^6+\frac{1}{2}{C}_0{\xi }^5{\xi }^{″}+\left(W+\mu \left(m^2-1\right)Q\right){\xi }^4+\frac{1}{2}{C}_0\mu {\xi }^{\prime 2}{\xi }^2+\mu \left(m^2-1\right)\left(K-C_0\mu \left(m^2-1\right)\right){\xi }^2+{\mu }^2{\left(m^2-1\right)}^2.$

(70)

Using the approach described in Ref. 45, the nonlinear equation (68) should be reduced to the following equivalent damping linear equation

${\xi }^{″}+{\chi }_1\left(A\right){\xi }^{\prime }+{\varpi }_1^2\left(A\right)\xi =0,$

(71)

where the real coefficient ${\chi }_1\left(A\right)$ and ${\varpi }_1^2\left(A\right)$ is given by

${\chi }_1\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}{\xi }^{\prime 2}{F}_1\left({\xi }^2\right)d\tau }{\underset{0}{\overset{T}{{\displaystyle \int }}}{\xi }^2{\xi }^{\prime 2}d\tau },\text{and}{\varpi }_1^2\left(A\right)=\frac{\underset{0}{\overset{T}{{\displaystyle \int }}}{\xi }^2{F}_0\left(\xi \right)d\tau }{\underset{0}{\overset{T}{{\displaystyle \int }}}{\xi }^{6}d\tau }.$

(72)

Apply the trial solution (44–72) and compute the integral returns

${\chi }_1\left(A\right)=\frac{1}{4A^2}\left(A^{4}Q+4K{A}^2+16\mu \right),$

(73)

${\varpi }_1^2\left(A\right)=\frac{1}{5A^4}\left(\begin{array}{l}\frac{35}{8}\left(\frac{1}{3}\Lambda -C_1-\frac{1}{2}{\omega }^2{C}_0\right)A^6+\left(5W+5\left(m^2-1\right)Q\mu +\frac{1}{2}\mu {\omega }^2{C}_0\right)A^4\\ +6\mu \left(m^2-1\right)\left(K-\left(m^2-1\right)\mu C_0\right)A^2+8{\left(m^2-1\right)}^2{\mu }^2\end{array}\right).$

(74)

Equation (71) has the following solution

$\xi \left(\tau \right)=\left(a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\omega \tau +\frac{1}{2\omega }\left(2b+a{\chi }_1\left(A\right)\right)\mathit{sin}\hspace{0.17em}\omega \tau \right)e^{-\frac{1}{2}{\chi }_1\left(A\right)\tau },$

(75)

where the frequency $\omega ,$ in the current case, has a formula in the form

${\omega }^2\left(A\right)={\varpi }_1^2\left(A\right)-\frac{1}{4}{\chi }_1^2\left(A\right).$

(76)

Inputting (73) and (74) into (76) produces

${\omega }^2\left(A\right)=16\left(\frac{35}{24}\left(\Lambda -3C_1\right)A^6+5\left(W+\left(m^2-1\right)Q\mu \right)A^4+6\mu \left(m^2-1\right)\left(K-\left(m^2-1\right)\mu C_0\right)A^2+8{\left(m^2-1\right)}^2{\mu }^2-\frac{5}{64}{\left(A^{4}Q+4K{A}^2+16\mu \right)}^2\right)/A^4\left(80+\left(35A^2-8\mu \right)C_0\right).$

(77)

It should be observed that in the limit case, $\mu \to 0$ , the frequency equation (77) reduces to equation (51). As previously stated, the stability condition can be described by using relation (45) to produce a polynomial of power 5 in ω². So, we have

${\iota }_5{\omega }^{10}+{\iota }_4{\omega }^8+{\iota }_3{\omega }^6+{\iota }_2{\omega }^4+{\iota }_1{\omega }^2+{\iota }_0=0.$

(78)

As previously stated, the stability criteria involve meeting the following conditions

${\text{Necessary}\hspace{0.17em}\text{conditions}:\iota }_i<0;i=1,3,5\hspace{0.17em}\text{and}\hspace{0.17em}{\iota }_i>0;i=0,2,4,$

(79)

${\text{Sufficient}\hspace{0.17em}\text{conditions}:2\iota }_1^2-5{\iota }_0{\iota }_2>0,{\iota }_2^2-2{\iota }_1{\iota }_3>0,{\iota }_3^2-2{\iota }_2{\iota }_4>0,2{\iota }_4^2-5{\iota }_3{\iota }_5>0.$

(80)

where

${\iota }_5=a^4\left(80+\left(35a^2-8\mu \right)C_0\right),$

(81)

${\iota }_4=\frac{5a^8{Q}^2}{4}+\frac{10}{3}{a}^6\left(3KQ-7\Lambda \right)-64\left(-3-4m^2+2m^4\right){\mu }^2$

$+32a^2\left(5b^2+K\left(8-3m^2\right)\mu \right)+20a^4\left(K^2-4W+2\left(3-2m^2\right)Q\mu \right)+a^2\left(105a^2{b}^2+16\mu \left(-b^2+6{\left(-1+m^2\right)}^2\mu \right)\right)C_0+70a^6{C}_1,$

(82)

${\iota }_3=b^2\{80b^2+5a^6{Q}^2+10a^4\left(3KQ-7\Lambda \right)+32K\left(8-3m^2\right)\mu +40a^2\left(K^2-4W+6Q\mu -4m^{2}Q\mu \right)+\left(105a^2{b}^2-8b^2\mu +96{\left(-1+m^2\right)}^2{\mu }^2\right)C_0+210a^4{C}_1\},$

(83)

${\iota }_2=\frac{5}{2}{b}^4\left(8K^2+12a^{2}KQ+3a^4{Q}^2-32W-28a^2\Lambda +16\left(3-2m^2\right)Q\mu +14b^2{C}_0+84a^2{C}_1\right),$

(84)

${\iota }_1=\frac{5}{3}{b}^6\left(6KQ+3a^2{Q}^2-14\Lambda +42C_1\right),$

(85)

${\iota }_0=\frac{5b^8{Q}^2}{4}.$

(86)

Numerical illustration

This section presents graphical representations of the stability behavior influenced by the porosity factor η and discusses the results concerning the kinematic viscosity coefficient μ. The study examines the stability conditions (59–62) in the absence of the kinematic viscosity factor μ. The stability diagrams were generated using the Mathematica 13.2 software package. Figures 2–8 illustrate the resulting curves under these conditions. The numerical findings revealed that the most effective stability condition is a condition (59). Figures 9 and 10 evaluate the role of kinematic viscosity in porous media, focusing on the stability criteria (79) and (80). Figures 8 and 10 display graphs that depict the effects of both the porosity coefficient and the kinematic viscosity on the stability results. The numerical findings revealed that the most effective stability conditions are conditions ${\iota }_i<0;i=1,5$ . Additionally, the computations indicated that the stability condition ${\iota }_1<0$ produces the same behavior as the condition (59). However, the condition ${\iota }_5<0$ specifically addresses the impact of kinematic viscosity on the stability. It was discovered that meeting the stability criterion requires specific consideration of the kinematic viscosity effects

$\mu >\frac{35a^{2}m\left(1-\rho \right)+80R\left(\rho +1\right)}{8m\left(1-\rho \right)}\hspace{0.17em}\text{provided}\hspace{0.17em}\text{that}\hspace{0.17em}\rho <1$

(87)

In all graphs illustrating critical stability conditions, the fluctuation of the constant radius R is displayed against H². In these planes, three stable regions and one unstable region were identified. Both conditions (59) and ${\iota }_1<0$ create stable regions 1 and 2. While condition ${\iota }_5<0$ yields stable Region-3.

Figure 2.

Demonstrates stability plane at different amounts of η <1.

Figure 3.

Demonstrates stability plane at different amounts of η >1.

Figure 4.

The stability regions for the variation of the rotation parameter $\mathrm{\Omega }$ when m = 1.

Figure 5.

The stability regions for the variation of the rotation parameter $\mathrm{\Omega }$ when m = 3.

Figure 6.

The stability regions for the variation of the density ratio $\rho$ when m = 1.

Figure 7.

The stability regions for the variation of the density ratio $\rho$ when m = 3.

Figure 8.

The stability areas for the identical system are shown in Figure 2, but with the azimuthal wavenumber m variation considered.

Figure 9.

The stability regions for the variation of the kinematic viscosity $\mu$ when m = 1.

Figure 10.

The stability regions for the variation of the kinematic viscosity $\mu$ when m = 3.

The numerical computation is performed using the parameters shown below

$m=2,\rho =0.1,\mathrm{\Omega }=2,a=1,b=1.$

The graphic illustration shows that instability occurs at small values of the radius R and magnetic field strength H². Stability, however, is achieved at relatively high concentrations of H². On the other hand, relatively large values of R have a stabilizing effect even in the absence of H². When H² levels are sufficiently high, a stable region (region 1) is formed. As the concentration of H² increases, the extent of this stable zone also expands. Conversely, the stable zone in region 2, which is sustained by unusually large values of R, decreases as the H² concentration rises. This suggests that increasing magnetic field intensity (H²) has a dual role: it exerts a stabilizing influence in region 1 and a destabilizing effect in region 2.

Figures 2 and 3 show the stability condition (59) when adjusting the porosity parameter η. Figure 2 illustrates the impact of η <1, while Figure 3 depicts the effect of η >1. Figure 1 shows that increasing small η decreases stable region -2 while not affecting region 1. In Figure 3, increasing η causes a decrease in stable region-1 and a larger decrease in region-2, which disappears for η ≥2. The rise in porosity parameter η has a destabilizing effect. This conclusion has been observed previously in Ref. 19.

Increasing the rotation parameter Ω significantly affects the stability profile and alters the azimuthal wavenumber m. The computations revealed that when m = 2, there is no noticeable impact on the stability profile. However, for m = 1 and m ≥ 3, the stability behavior exhibits two opposing effects. Figures 4 and 5 illustrate this tendency for m = 1 and m = 3, respectively. In Figure 4, with m = 1, increasing the rotation parameter Ω decreases the width of stable zone-1 while expanding the area of stable region-2. This pattern indicates that increasing Ω has a dual effect on the stability profile: it exerts a destabilizing influence in stable area 1 and a stabilizing influence in stable region 2. Figure 5 depicts the scenario when m reaches a value of 3. The graph shows that as Ω increases, stable region-1 broadens while stable region-2 narrows. This observation suggests that the dual role of Ω can be observed in both stability behavior and contrast profiles, highlighting its complex influence on the system’s stability.

Figures 6 and 7 illustrate the impact of increasing the density ratio ρ on stability behavior. When comparing these graphs to those in Figures 4 and 5, it becomes evident that increasing the density ratio ρ produces effects similar to those observed when Ω is increased. This comparison highlights the analogous influence of both parameters on the stability profile, suggesting that changes in density ratio ρ and rotation parameter Ω can lead to comparable modifications in the system’s stability behavior.

It is customary to investigate the influence of the azimuthal wavenumber m on the stability profile, and thus Figure 8 has been established for this purpose. This graph indicates that as mmm increases, the stable region-1 expands in breadth while the stable region-2 contracts in the area. This suggests that increasing m serves two purposes: it enhances the stability in region 1 by broadening it, while simultaneously reducing the stability in region 2 by causing it to shrink.

When stability conditions ${\iota }_i<0;i=1,5$ are activated, the influence of kinematic viscosity μ is examined. As previously mentioned, the critical situation resulted in the formation of a third stable region. Figures 9 and 10 illustrate the variation of the parameter μ for m = 1 and m = 3, respectively. The same configuration as in Figures 2–8 is employed, with the exception that a = b = 0.1. It is observed that large values of amplitude (a) tend to diminish the influence of viscosity. These graphs demonstrate that increasing μ leads to an expansion of the stable region-3. Additionally, higher values of mmm further enhance the stabilizing effect of μ.

Conclusion

This paper introduces an innovative approach to fluid dynamics by employing fractal theory to model the intricate properties of porosity and viscosity in fluid flow, specifically focusing on spinning ferrofluid columns in porous media. The study’s novelty lies in demonstrating that an inviscid fluid in fractal space behaves analogously to a viscous fluid in traditional space, offering a transformative framework for stability analysis and nonlinear dynamics. By integrating fractal derivatives, the research bridges the gap between the simplified mathematical treatment of inviscid fluids and the complexities inherent in viscous fluid dynamics, which often require solving the challenging Navier–Stokes equations with additional boundary conditions. The contributions of this study are multifaceted. First, it provides a novel mathematical framework that simplifies the analysis of viscous fluids without compromising the depth or accuracy of the results. This is achieved by transforming fractal derivatives into traditional derivatives, capturing the effects of viscosity in a computationally efficient manner. Second, the research delivers critical insights into stability criteria influenced by parameters such as azimuthal magnetic field strength, angular rotation, porosity, and kinematic viscosity. Notably, the findings highlight the stabilizing role of kinematic viscosity when the inner fluid density is lower than that of the outer fluid, a crucial factor for practical applications. The implications of these findings extend across various fields. The results deepen our understanding of the behavior of ferrofluids in porous media under rotational and magnetic influences, with applications ranging from advanced cooling systems and enhanced oil recovery to biomedical technologies like magnetic targeting. The numerical validations provided in the study further reinforce the reliability and applicability of the proposed approach. The study’s methodology offers a significant advancement in the field by reducing the computational effort and time traditionally required for analyzing viscous fluid systems. This efficiency enables more robust modeling of ferrofluid systems within porous environments, paving the way for practical implementations in engineering and applied physics.

Reduces the complexity of analyzing viscous fluids by using fractal derivatives, which emulate the behavior of viscous fluids in traditional space. Avoids the need for solving the full Navier–Stokes equations, thereby streamlining the computational process. This research not only enriches the theoretical landscape of fluid dynamics but also provides actionable insights for practical applications, marking a significant step forward in the study of ferrofluids and porous media.

Future Directions

While the findings of this paper are comprehensive and align with theoretical expectations, future research could build upon this foundation by exploring:

• Multi-scale Interactions: Extending fractal theory to address the dynamics of multi-scale porous media, considering both micro- and macro-scale interactions.

• Additional Parameters: Investigating the effects of temperature gradients, variable magnetic field strengths, or other external forces on fluid stability and nonlinear dynamics.

• Real-world Applications: Applying the results to practical systems such as magnetic cooling technologies, ferrofluid-based sensors, and advanced material processing techniques.