Analysis of travelling wave solutions for Eyring-Powell fluid formulated with a degenerate diffusivity and a Darcy-Forchheimer law

Nomenclature

$K_{1}$: Permeability of the medium.

$K_{2}$: Inertial permeability.

L: Fluid characteristic length.

M: Hartmann number.

MHD: Magnetohydrodynamics.

PME: Porous Medium Equation.

Re: Reynolds.

$U_0$: Fluid characteristic velocity.

$\rho$: Fluid density.

$\mu$: Fluid viscosity.

$\nu$: Fluid kinematic viscosity.

$B_{0}$: Applied magnetic field.

$\sigma$: Electrical conductivity in the medium.

1.   Introduction

Any kind of material exhibiting pores filled by a flowing fluid leads to a physical interpretation beyond the Newtonian description related with viscosity and diffusion. Currently, flow in porous medium is a source of research, not only because of purely mathematical interests, but also because of their applications in electrochemical, mechanical, geophysical, metallurgical and chemical processes. In 1856, the observations made by Henry Darcy ^[1] associated to fluid flow in porous medium led to the formulation of a physical law to describe the basic dynamic experienced. Such law was intended to describe the fluid principles for low Reynolds numbers, typical of flows in this kind of porous medium. Afterwards, Forchheimer in 1901 ^[2] and Jaeger in 1956 ^[3] provided an extension of Darcy's law to account for turbulent flows at higher Reynolds numbers. Based on these works, the Darcy's law was renamed as Darcy-Forchheimer's law by Mustak in ^[4] and ^[5].

It shall be noted that the above mentioned researchers considered the viscosity and diffusive terms as provided by the classical Newtonian description. Nonetheless and since then, several researchers have tried to extend the viscosity formulation to increase the modelling accuracy beyond the Newtonian scope. As a set of representative examples of the described non-Newtonian fluid models, the reader is referred to the works ^[6,7,8,9] together with references therein. The non-Newtonian fluids are classified into several classes. It is relevant to highlight that efforts have been made to describe diffusion phenomena based on the kinetic theory of liquids rather than on experimental or intuitive principles related with viscosity. Such efforts were conducted with success leading to a kind of non-Newtonian fluid known as Eyring-Powell. The use of the kinetic theory of liquids to describe viscosity and diffusive mechanisms makes the Eyring-Powell fluid particularly attractive to model flows given by electrically conducting particles arising in Magnetohydrodynamics (MHD). In a set of interesting analysis, Akbar et al. ^[10] obtained the numerical results for two dimensional MHD Erying-Powell fluid over stretching surface. Hina et al. ^[11] analyzed the heat transfer for MHD Erying-Powell fluid under the slip effect at the boundary. Bhatti et al. ^[12] considered two dimensional MHD Erying-Powell fluid under permeable surface and analyzed the heat transfer by converting the equations into nonlinear ordinary equations and that are solved based on numerical principles. Other interesting results in non-Newtonian fluids are given in the references [13,14,15,16,17,18], but particularly there are some recent achievements in this regards to highlight: In ^[19] a modified Eyring-Powell model is introduced to improve predictions related with the Eyring-Powell parameter on flow velocity. In ^[20], the authors discuss the adequacy of Eyring-Powell models for viscous dissipation in a flow past a convectively heated stretching sheet. Furthermore, in ^[21] and ^[22], an Eyring-Powell fluid is considered to model Coriolis effects over a non-uniform surfaces together with a discussion on the adequacy of the formulation. Finally, the reader is referred to ^[23] and ^[24] for discussions on other types of fluids and other types of analysis based on purely numerical schemes respectively.

Other techniques have been of help to describe flows under non-Newtonian approaches. For instance, some kinds of solutions known as travelling and solitary waves have provided interesting results from analytical and numerical perspective. In this regard, the reader is referred to the set of studies ^[25,26,27].

As a main research question driving the present study, it shall be mentioned the idea of searching for appropriate regularity conditions even when the Eyring-Powell viscosity term is formulated with a degenerate diffusion. Specially but not limited, the question of finding an exponential profile of solution.

1.1.   Model formulation

Let us start by considering a one dimensional Darcy-Forchheimer flow with a diffusion term given by an Erying-Powell principle. The velocity vector is then $V = (v_1, 0, 0)$ where $v_{1}(y, t)$ is the velocity profile in the $x-$direction that is dependant on the transversal direction $y$. The geometry under study is given by an open porous domain where the flow is developed through the $x-$direction and is initially distributed by an initial condition in the $y-$direction. In addition, the diffusion is further modified to account for a nonlinear term of the form $\frac{\partial^2 v_1^m}{\partial y^2}$ (refer to ^[28] and ^[29] for a complete discussion about this kind of non-linearity and to ^[14,15,16] for different applications of particular diffusion terms). Note that the basic homogeneous problem

is well known as porous medium equation. The PME is a nonlinear parabolic equation with degenerate diffusivity. Indeed, the diffusivity is given by $D(v_1) = m \, v_1^{m-1}$ that is null whenever $v_1 = 0$. This degeneracy leads to the loss of a positive principle for null (or slightly null) solutions, and hence to the loss of regularity as understood for gaussian diffusion problems.

Based on the exposed principles, the equation under study emerges from the MHD theory and reads

where $m > 1$, $P$ is the pressure field acting along the $x-$direction and $\beta$ and $d_{1}$ are the characteristic constants of the Eyring-Powell fluid model.

To account for a single general equation, introduce the following non-dimensional quantities (see ^[30]) as follows

where $L$ is a characteristic length that can be given by the pores dimension. After replacement of the set of expressions in (1.3) into the Eq (1.2) (ignoring * for simplicity)

where $R_{e} = \frac{U_{0}L}{\nu}$ is the Reynolds number and and $M = B_{0}\sqrt{\frac{L\sigma}{\rho}}$ refers to the Hartmann number. Making the differentiation of Eq (1.4) with respect to $x$

After using the value of $\frac{d P}{dx}$ in Eq (1.4)

In addition, consider the following initial generalized condition

where $v_{0}(y)$ is the initial velocity that represents any kind of non-homogeneous distribution with the only requirement to be bounded in $L^\infty-$norm and under a finite flowing mass quantity, i.e., a bound in $L^1-$norm. In addition, note that this initial condition allows modelling any kind of irregularities or discontinuities that may exist in the porous medium involved.

2.   Regularity, existence and uniqueness of solutions

Given the degeneracy in the diffusivity introduced by the non-linear term $\frac{\partial^2 v_1^m}{\partial y^2}$ and given the generalized initial velocity distribution, solutions are analyzed based on a weak formulation. To this end, assume the existence of a test function $\Omega \in C^\infty (R)$ with compact support (or at least with a decreasing behaviour over $R$) such that for $0 < \tau < t < T$

Given a finite $r_{0},$ admit the following ball $B_{r}$ centered in $r_{0}$ and with radium $r > > r_{0}$, such that the following equation is defined

in $B_{r}\times \left[0, T\right],$ with the following boundary like and initial conditions

Firstly, the following theorem expresses the regularity properties of any solution to (1.5), i.e., solutions exist globally.

Theorem 2.1. Given $v_{0}(y)\in L^{1}(R)\cap L^{\infty}(R),$ then any non-negative solution $(v_1 \geq 0)$ to (1.2) is bounded for all$(y, t)\in B_{r}\times \left[0, T\right]$ with $r > > 1.$

Proof. Consider $\lambda \in R^{+}$, the following supporting cut-off function is defined as

such that

Multiplying Eq (2.2) by $\psi_{\lambda}$ and integrating over $B_{r}\times\left[\tau, T\right]$

Based on some results in ^[28] and ^[29], for some large $r > > r_{0} > 1$ the following holds

and

Considering the spatial variable $y$ close to $\partial B_{r},$ it can be assumed $y\sim r.$ Then

After using the expression (2.4) and after integration in the nonlinear diffusion term, the following holds

Since $r > > 1$ and since $\Omega$ is sufficiently small and regular in the asymptotic approach, it is possible to choose $\Omega$ in such a way $\frac{\partial \Omega}{\partial y}\psi_{\lambda} < < 1$ in the proximity of the boundary $\partial B_r$. Then, the above expression becomes

Based on (2.7), the following integral holds

After using the expressions (2.5), (2.6), (2.8) and (2.9) into (2.4), the following holds

which implies that

Consider a test function $\Omega\in C^\infty(R)$ of the form

Choose $\alpha$ in a way that the expression (2.10) is convergent, therefore

For convergence, take $\alpha > \frac{m+2}{m-1}$ and $r\rightarrow \infty,$ so that the expression (2.12) becomes

where $\Upsilon$ is a suitable value obtained for each $B_r$ and is sufficiently small for $r\rightarrow \infty$; hence, a finite value of $\Upsilon$ holds. Note that $\Omega$ and $v_1$ are not negative, then the integrals in the expression (2.13) are finite in $\tau < s < t\leq T,$ which shows that the solution $v_1$ is bounded in $R\times\left[0, T\right].$

Now, it is the aim to show additional regularity results, in particular that $\frac{\partial v_{1}}{\partial y}$ is bounded.

Theorem 2.2. Assume that $v_{1}(y)$ is a non-negative solution to Eq (1.5), then $\frac{\partial v_{1}}{\partial y}$ is bounded for all$(y, t)\in R\times \left[0, T\right].$

Proof. Multiplying Eq (1.5) by $v_{1}$ and applying integration by parts, the following holds

which implies that

Now, applying the integration done

The Theorem 2.1 established the bound properties of non-negative solutions, as a consequence the right hand side of Eq (2.14) is bounded. Choose a finite bounding constant $A_{5}$ such that

leading to conclude on the bound of $\frac{\partial v_{1}}{\partial y}$.

The degeneracy introduce by the non-linear diffusion term may lead to the loss of uniqueness of non-negative solutions. Hence, a result about uniqueness of solutions is required to be shown.

Theorem 2.3. Let us consider two non-negative solutions to the Eq (1.5), $v_{11}(y, t)$ and$v_{12}(y, t)$. Then, if both solutions have the same initial distribution, uniqueness holds, i.e., $v_{11}(y, t) = v_{12}(y, t)$.

Proof. First of all, define a test function $\Psi(y, t)\in C^{\infty}\left(Q_{T}\right)$, where $Q_T = R \times [0, T]$. The weak formulations for both solutions $v_{11}$ and $v_{22}$ are given as

and

Assume $v_{11} \geq v_{12}$. After substraction, the following holds

For $v_{11}(y, t) > v_{12}(y, t),$ the following holds

where $C_{2} = max_{(y, t)\in Q_{T}}(v_{11}).$ Now, consider

where $C_3$ is a suitable bounding constant that can be defined supported by the regularity results shown in Theorem 2.2.

Based on the expressions (2.18) and (2.19), the Eq (2.17) becomes

Assume the following smooth shape test function defined as

where $\gamma_1$ can be chosen such that

For simplicity

For $\vert y \vert \rightarrow \infty$, the integral mass shall be null. Then, for $R > > 1$

In the asymptotic approximation $|y|\rightarrow \infty$

then

Now, define

where

The integral (2.28) is finite under the condition for $\gamma_1$ in (2.26).

Now, differentiate (2.21) with respect to $y$, so that

After integration

where $C_{4}(\gamma_1) = 4\gamma_1(\gamma_1+1).$

Now, the intention is to assess each of the integrals involved in the right hand side of Eq (2.20)

and similarly

Combining the assessments in (2.29), (2.30) and (2.31) with the Eq (2.20), the following holds

Since $\left|\chi(y)\right| < \infty$ and for $t > 0,$ with $\sup\left|v_{11}-v_{12}\right|\rightarrow 0$

which implies that $v_{11}\leq v_{12}$. This result contradicts the initial assumption $v_{11} \geq v_{12}$. As a consequence, the only compatible result is to consider $v_{11} = v_{12}$, leading to conclude on solutions uniqueness.

3.   Travelling waves analysis

The travelling wave profiles of solutions are given as per the following definitions: $v_{1}(y, t) = k(\zeta),$ where $\zeta = y-a_{1}t,$ $a_{1}$ denotes travelling wave speed and the profile $k: R\rightarrow (0, \infty)$ belongs to $L^{\infty}(R)$. The Eq (1.5) is then transformed in terms of the travelling wave profile as

Let us introduce the following new variables

such that the Eq (3.1) can be analyzed in a phase space driven by the following set of equations

Making $Y' = 0$ and $Z' = 0$ to determine the critical points, the following is obtained

with solutions

and

Consequently, the system critical points are given by $(Y_{1}, 0)$ and $(Y_{2}, 0)$. The next intention is to characterize such critical points together with the bundles of orbits in their proximity.

3.1.   Geometric perturbation theory

The Geometric perturbation theory is used to show the asymptotic behavior of specifically defined manifolds in the proximity of the system critical points. Firstly, assume the following manifold $I_0$ defined as

with the system critical points $(Y_{1}, 0)$ and $(Y_{2}, 0).$

A perturbed manifold $I_{\beta}$ close to $I_{0}$ in the critical point $(Y_{1}, 0)$ is defined as

where $\beta$ denotes a perturbation parameter close to the equilibrium $\left(Y_{1}, 0\right)$ and $C_{5}$ is a suitable constant obtained after root factorization. For the coming assessments, assume $Y_{3} = Y-Y_{2}$.

The intention is to apply the Fenichel invariant manifold theorem ^[31] as formulated in ^[32] and ^[33]. To this end, it shall be shown that $I_{0}$ is a normally hyperbolic manifold, i.e., the eigenvalues associated to $I_{0}$ in the linear frame proximal to the critical points, and transversal to the tangent space, have non-zero real part. This is shown based on the following equivalent linear flow associated to $I_{\beta}.$

The related eigenvalues are both real. Consequently, $I_{0}$ is a hyperbolic manifold. Now, it is shown that the manifold $I_{\beta}$ is locally invariant under the flow (3.3), so that the manifold $I_{0}$ can be shown as an asymptotic approach to $I_{\beta}.$ For this basis, consider the functions

which are $C^{i}\left(R\times\left[0, \delta_{1}\right]\right), \, i \geq 1$, in the proximity of the critical point $\left(Y_{1}, 0\right).$ In this case, $\delta_{1}$ is determined based on the following flows that are considered to be measurable a.e. in $R$

where $\varepsilon_0$ is the infinitesimal distance between the critical point $(Y_1, 0)$ and the flow $\psi^{I_{0}}$, which exists given the normal hyperbolic condition of $I_0$. As the solutions have been shown to be bounded, it is possible to conclude on a finite $\delta_{1}$. So the distance between the manifolds holds the normal hyperbolic condition for $\delta_{1}\in\left(0, \infty\right)$ and $\beta$ sufficiently small in the proximity of the critical point $\left(Y_{1}, 0\right).$

Now, consider a perturbed manifold $I_{\gamma}$ associated to $I_{0}$, but for the critical point $(Y_{2}, 0)$ and defined as

where $\gamma$ denotes a perturbation parameter in the proximity of the equilibrium $\left(Y_{2}, 0\right)$ and $C_{6}$ is a suitable constant obtained after root factorization. Suppose $Y_{4} = Y-Y_{1}$, then the following flow associated to $I_{\gamma}$ is given by

The associated eigenvalues are both real, leading to state that $I_{0}$ is, indeed, a hyperbolic manifold. Now, it is shown that the manifold $I_{\gamma}$ is locally invariant under the flow (3.3), so that the manifold $I_{0}$ can be shown as an asymptotic approach to $I_{\gamma}.$ For this basis consider the functions

which are $C^{i}\left(R\times\left[0, \delta_{1}\right]\right), \, i \geq 1$, in the proximity of the critical point $\left(Y_{2}, 0\right).$ In this case, $\delta_{1}$ is determined based on the following flows that are considered to be measurable a.e. in $R:$

where $\varepsilon_1$ is the infinitesimal distance between the critical point $(Y_2, 0)$ and the flow $\psi^{I_{0}}$, which exists given the normal hyperbolic condition of $I_0$. As the solutions are bounded, $\delta_{1}$ is finite. Hence, the distance between the manifolds holds the normal hyperbolic condition for $\delta_{1}\in\left(0, \infty\right)$ and $\gamma$ sufficiently small, in the proximity of the critical point $\left(Y_{2}, 0\right).$

3.2.   Construction of travelling wave profiles

Since the manifold $I_{0}$ has been shown to satisfy the normal hyperbolic condition under the flow (3.3), it is possible to obtain asymptotic profiles operating in the linearized flows $I_\beta$ and $I_\gamma$.

For this purpose, consider the Eq (3.3) such that a family of trajectories in the phase plane $(Y, Z)$ can be obtained by

For $Y \gg 1$, the function $H > 0$ while for $0 < Y \ll 1,$ the function $H < 0$. In addition, as the function $H$ is singular for $Z = 0$, a trajectory in the proximity of the critical condition $(Y_{1}, 0)$ shall require

After using the Eq (3.3) into Eq (3.10), the following holds

Typically, the flow velocity in a porous medium permits the hypothesis $U_0 \ll 1$. In addition, in the assumption of a sufficiently large characteristic length, $L \gg 1$, the parameter $\epsilon \ll 1$ (see (1.3)) for typical fluid densities. As a consequence, such $\epsilon$ is considered to construct a solution as an asymptotic expansion in terms of the mentioned parameter $\epsilon$, understood as a perturbation and as given in (1.3). The selection of this parameter permits to account for a physical nuisance as $\epsilon$ is related with some of the constants involved in the fluid problem. Then, the following expansion holds for suitable regular functions $Y_{1j}, \, j = 0, 1, 2...$

After using (3.12) into (3.11), the following holds

Comparing the terms of orders $\epsilon^{0}$ and $\epsilon^{1}$

and

The standard resolution of (3.14) leads to

Differentiate (3.16) to get

Upon substitution of (3.16) and (3.17) in (3.15), it holds

where $a = \left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right), \, \, \, \, \, \,$ $b = \frac{2A_{3}\sqrt{\left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)^{2}+\frac{A_{5}}{A_{3}}}}{a_{1}}, \, \, \, \, \, \, \,$ $c = \left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)+ \sqrt{\left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)^{2}+\frac{A_{5}}{A_{3}}}, \, \, \, \, \, \,$ $d = \left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)- \sqrt{\left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)^{2}+\frac{A_{5}}{A_{3}}} \, \,$ and $f = -\frac{2M^{2}}{R_{e}a_{1}}+\frac{2A_{2}}{a_{1}}.$

After solving the Eq (3.18)

where $C_{1} = \frac{2a}{d-c}- \frac{2A_{3}A_{5}(c+d)}{a_{1}bcd}$, $C_{2} = \frac{A_{5}}{c}+\frac{3(2A_{3}a+a_{1}b)}{2a_{1}}+f$ and $C_{3} = \frac{2a}{d-c}-\frac{12aA_{3}}{a_{1}^{2}b}$.

Introducing the expressions (3.16) and (3.19) into (3.12), the following holds

which implies that

The Eq (3.21) shows the existence of an exponential decaying under the travelling wave profile asymptotic expansion. This results is not trivial in this case given the degeneracy in the diffusivity under the non-linear diffusion of porous medium kind.

Now, the intention is to show that the defined supporting manifold $I_{\beta}$ preserves the exponential behavior shown in the proximity of the associated critical point. For this purpose, consider the expression (3.7), so that

After solving (3.22)

Introducing the expression (3.7) into (3.23), the following holds

such that a solution is

which implies that

This last expression shows the conservation of the exponential profile in the proximity of the critical point kept by the asymptotic manifold $I_{\beta}.$

The existence and exact assessments of an orbit in the proximity of the critical point $(Y_{2}, 0)$ follows the same procedure as discussed above for the point $(Y_{1}, 0)$. Indeed, a trajectory in the phase plane is given as

Solving the Eq (3.27) leads to obtain

where $a = \left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right),$ $b = \frac{2A_{3}\sqrt{\left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)^{2}+\frac{A_{5}}{A_{3}}}}{a_{1}},$

$c = \left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)+ \sqrt{\left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)^{2}+\frac{A_{5}}{A_{3}}},$

$d = \left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)- \sqrt{\left(\frac{M^{2}}{2A_{3}R_{e}}+\frac{A_{2}}{2A_{3}}\right)^{2}+\frac{A_{5}}{A_{3}}}$,

$f = -\frac{2M^{2}}{R_{e}a_{1}}+\frac{2A_{2}}{a_{1}},$ $C_{1} = \frac{2a}{d-c}- \frac{2A_{3}A_{5}(c+d)}{a_{1}bcd}$, $C_{2} = \frac{A_{5}}{c}+\frac{3(2A_{3}a+a_{1}b)}{2a_{1}}+f$ and $C_{3} = \frac{2a}{d-c}-\frac{12aA_{3}}{a_{1}^{2}b}$.

As a consequence, a solution is given by the following expression

Again, the intention now is to show that the defined supporting manifold $I_{\gamma}$ preserves the exponential behavior in the proximity of the critical point $(Y_2, 0)$. For this purpose, consider the expression (3.8) so that

After solving (3.30), the following holds

Introducing the expression (3.8) into (3.31)

such that a solution is given by

which implies that

This last expression permits to show the conservation of the exponential profile in the proximity of the critical point under study by the asymptotic manifold $I_{\gamma}.$

4.   Numerical validation

According to expressions (3.25) and (3.33), the analytical assessments have provided an exponential tail in the proximity of the system critical points. To validate the assessments, it is the intention to compare such analytical exponential tail with a numerical simulation of the Eq (2.1). Note that the numerical exploration is not intended to introduce a wide parametrical analysis with each of the involved fluid constants, on the contrary the objective is to provide evidences on the accuracy of the analytical process followed. In addition, for the validation exercise the expression (3.33) has been considered, under the form $Y-Y_{1} = e^{-\gamma\sqrt{C_{6}}\zeta}.$ The numerical assessment is executed with the bvp4c in Matlab (refer to ^[34] for further details on this software tool).

The $R_e$ number has been considered as a free parameter to establish the accuracy of solutions. In addition, the following particular values in the other fluid constant have been considered for the sake of simplicity, but with no impact in the conclusions

Note that $L = 1000$ as the domain of integration has been considered as $(-500,500)$. This domain has been selected after some numerical trials and provides an efficient compromise between low impact of the collocation methods at the borders and a moderate computational time. In addition, the global absolute error allowed is of $10^{-4}$.

Note that the $R_e$ in porous media is typically low, hence the considered values are $R_e = 1; 100$.

As it can be noticed from Figures 1 and 2, the analytical asymptotic solution (3.33) based on the geometric perturbation theory and the numerical solution to (2.1) fit an exponential behaviour with an order deficiency in the decreasing rate. Nonetheless, solutions are very close for increasing values of $\zeta$, this is in the asymptotic approach where the analytical solution has been obtained. This is not a trivial conclusion for the degenerate diffusion given in (2.1) and represents a finding to remark.

5.   Conclusions

The set of assessments and discussions presented in this paper have provided evidences on the regularity of solutions to a kind of fluid flow with degenerate diffusivity and formulated with Darcy-Forchheimer porosity principles and Eyring-Powell viscosity terms. Afterwards, asymptotic travelling waves solutions have been explored based on a perturbation technique. Particularly, the Geometric Perturbation Theory has shown that an exponential decreasing rate holds even when the equation discussed presents a degenerate diffusivity. Furthermore, two linearized manifolds in the proximity of the critical points have been shown to hold and an exponential behaviour has been obtained accordingly. Eventually, a numerical exercise has been introduced to account for the validation of the analytical assessments done.

Conflict of interest

The authors declare that there are not any conflict of interests.