Analysis of travelling wave solutions for Eyring-Powell ﬂuid formulated with a degenerate di ﬀ usivity and a Darcy-Forchheimer law

: The goal of this paper is to provide analytical assessments to a ﬂuid ﬂowing in a porous medium with a non-linear di ﬀ usion linked to a degenerate di ﬀ usivity. The viscosity term is formulated with an Eyring-Powell law, together with a non-homogeneous di ﬀ usion typical of porous medium equations (as known in the theory of partial di ﬀ erential equations). Further, the equation is supplemented with an absorptive reaction term of Darcy-Forchheimer, commonly used to model ﬂows in porous medium. The work starts by analyzing regularity, existence and uniqueness of solutions. Afterwards, the problem is transformed to study travelling wave kind of solutions. An asymptotic expansion is considered with a convergence criteria based on the geometric perturbation theory. Supported by this theory, there exists an exponential decaying rate in the travelling wave proﬁle. Such exponential behaviour is validated with a numerical assessment. This is not a trivial result given the degenerate di ﬀ usivity induced by the non-linear di ﬀ usion of porous medium type and suggests the existence of regularity that can serve as a baseline to construct numerical or energetic approaches.


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.

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 ∂ 2 v m 1 ∂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 m−1 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 β 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 = U 0 L ν is the Reynolds number and and M = B 0 Lσ ρ refers to the Hartmann number. Making the differentiation of Eq (1.4) with respect to x After using the value of dP 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 ∞ −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.

Regularity, existence and uniqueness of solutions
Given the degeneracy in the diffusivity introduced by the non-linear term ∂ 2 v m 1 ∂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 Ω ∈ C ∞ (R) with compact support (or at least with a decreasing behaviour over R) such that for 0 < τ < 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 v 1 (t) ∂Ω ∂s , 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.
Proof. Consider λ ∈ R + , the following supporting cut-off function is defined as Multiplying Eq (2.2) by ψ λ and integrating over Based on some results in [28] and [29], for some large r >> r 0 > 1 the following holds Considering the spatial variable y close to ∂B r , it can be assumed y ∼ r. Then After using the expression (2.4) and after integration in the nonlinear diffusion term, the following holds ∂Ω ∂y ∂ψ λ ∂y dyds .
Since r >> 1 and since Ω is sufficiently small and regular in the asymptotic approach, it is possible to choose Ω in such a way ∂Ω ∂y ψ λ << 1 in the proximity of the boundary ∂B r . Then, the above expression becomes (2.8) 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 Choose α in a way that the expression (2.10) is convergent, therefore For convergence, take α > m+2 m−1 and r → ∞, so that the expression (2.12) becomes where Υ is a suitable value obtained for each B r and is sufficiently small for r → ∞; hence, a finite value of Υ holds. Note that Ω and v 1 are not negative, then the integrals in the expression (2.13) are Now, it is the aim to show additional regularity results, in particular that ∂v 1 ∂y is bounded. Theorem 2.2. Assume that v 1 (y) is a non-negative solution to Eq (1.5), then ∂v 1 ∂y is bounded for all Proof. Multiplying Eq (1.5) by v 1 and applying integration by parts, the following holds 14) 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 ∂v 1 ∂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).

Travelling waves analysis
The travelling wave profiles of solutions are given as per the following definitions: v 1 (y, t) = k(ζ), where ζ = y − a 1 t, a 1 denotes travelling wave speed and the profile k : R → (0, ∞) belongs to L ∞ (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.

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 β close to I 0 in the critical point (Y 1 , 0) is defined as where β denotes a perturbation parameter close to the equilibrium (Y 1 , 0) and C 5 is a suitable constant obtained after root factorization. For the coming assessments, assume 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 β .
The related eigenvalues are both real. Consequently, I 0 is a hyperbolic manifold. Now, it is shown that the manifold I β is locally invariant under the flow (3.3), so that the manifold I 0 can be shown as an asymptotic approach to I β . For this basis, consider the functions in the proximity of the critical point (Y 1 , 0) . In this case, δ 1 is determined based on the following flows that are considered to be measurable a.e. in R where ε 0 is the infinitesimal distance between the critical point (Y 1 , 0) and the flow ψ 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 δ 1 . So the distance between the manifolds holds the normal hyperbolic condition for δ 1 ∈ (0, ∞) and β sufficiently small in the proximity of the critical point (Y 1 , 0) . Now, consider a perturbed manifold I γ associated to I 0 , but for the critical point (Y 2 , 0) and defined as where γ denotes a perturbation parameter in the proximity of the equilibrium (Y 2 , 0) and C 6 is a suitable constant obtained after root factorization. Suppose Y 4 = Y − Y 1 , then the following flow associated to I γ 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 γ is locally invariant under the flow (3.3), so that the manifold I 0 can be shown as an asymptotic approach to I γ . For this basis consider the functions in the proximity of the critical point (Y 2 , 0) . In this case, δ 1 is determined based on the following flows that are considered to be measurable a.e. in R : where ε 1 is the infinitesimal distance between the critical point (Y 2 , 0) and the flow ψ I 0 , which exists given the normal hyperbolic condition of I 0 . As the solutions are bounded, δ 1 is finite. Hence, the distance between the manifolds holds the normal hyperbolic condition for δ 1 ∈ (0, ∞) and γ sufficiently small, in the proximity of the critical point (Y 2 , 0) .

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 β and I γ .
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 1, the function H > 0 while for 0 < Y 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 1. In addition, in the assumption of a sufficiently large characteristic length, L 1, the parameter 1 (see (1.3)) for typical fluid densities. As a consequence, such is considered to construct a solution as an asymptotic expansion in terms of the mentioned parameter , understood as a perturbation and as given in (1.3). The selection of this parameter permits to account for a physical nuisance as is related with some of the constants involved in the fluid problem. Then, the following expansion holds for suitable regular functions Y 1 j , j = 0, 1, 2...
After using (3.12) into (3.11), the following holds Comparing the terms of orders 0 and 1 and The standard resolution of (3.14) leads to (3.17) Upon substitution of (3.16) and (3.17) in (3.15), it holds After solving the Eq (3.18) where . Introducing the expressions (3.16) and (3.19) into (3.12), the following holds 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 β 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 This last expression shows the conservation of the exponential profile in the proximity of the critical point kept by the asymptotic manifold I β .
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 As a consequence, a solution is given by the following expression v 1 (y, t) = 1 Again, the intention now is to show that the defined supporting manifold I γ 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 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 γ .

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 −γ √ C 6 ζ . 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 0 < A 5 1, U 0 = 1, A 1 = 1, 0 < 1, M = 1, ν = 1, K 1 = 1, K 2 = 1. 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 ζ, 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. Both solutions preserves the exponential behaviour. The picture on the right represents a zoom in the interval [3,5]. The reader can observe a certain order deficiency between the analytical and numerical assessments. Nonetheless, both solutions are close, the distance is of the order 10 −7 .

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.