Regularity and Travelling Wave Proﬁles for a Porous Eyring–Powell Fluid with Darcy–Forchheimer Law

: The goal of this study is to provide analytical and numerical assessments to a ﬂuid ﬂow based on an Eyring–Powell viscosity term and a Darcy–Forchheimer law in a porous media. The analysis provides results about regularity, existence and uniqueness of solutions. Travelling wave solutions are explored, supported by the Geometric Perturbation Theory to build proﬁles in the proximity of the equation critical points. Finally, a numerical routine is provided as a baseline for the validity of the analytical approach presented for low Reynolds numbers typical in a porous medium. writing—review J.L.D.P., S.u.R., A.N.R. and J.R.G.; visualization, J.L.D.P., S.u.R., A.N.R. and J.R.G.; supervision, J.L.D.P.; project administration, J.L.D.P., S.u.R., A.N.R. and J.R.G.; funding acquisition, J.L.D.P. and A.N.R.


Introduction
Transport processes in porous media are typically given in several applications: mechanics, electrochemistry, metallurgic and geophysics, to cite some. In 1856, Darcy [1] provided a mathematical relation to study the fluid dynamics involved in porous media, particularly, water flows in sand beds.
The flow of a fluid in a porous medium typically occurs under low Reynolds numbers due to the interstitial structure of these kind of media. To increase the Reynolds numbers, it is required to enlarge the values of pressure gradients in the medium, i.e., to account for a steep pressure distribution. This fact improves turbulent mechanisms in the fluid, making the Darcy's law vague. Consequently, and to increase modelling accuracy in the mentioned circumstances, first Forchheimer in 1901 [2], and later Jaeger in 1956 [3], proposed empirical relations between the velocity and the hydraulic gradient in a turbulent flow. Currently, other modelling strategies are available. For example, in [4], analytical and numerical solutions were provided to model in shallow water. To this end, the (G /G)−expansion method was considered to obtain solutions to the classical Gardner equation, modified Korteweg-de Vries equation, and generalized Korteweg-de Vries equation. Afterward, the numerical approach allowed us to compare with the analytical findings.
Based on the studies in [2,3], the Darcy's law was renamed as Darcy-Forchheimer's principle by Mustak in [5] and Ward in [6]. Note that the mentioned formulations were considered for Newtonian fluids. Nonetheless, the Darcy-Forchheimer's law has been employed for non-Newtonian fluids as well (the reader is referred to [7][8][9][10] and references there for a wider scope).
There exist plenty of non-Newtonian fluids in the literature, but particularly, we focus our analysis on the Eyring-Powell fluid flow. This model allows us to describe the shear stress and the associated viscosity based on the kinetic theory of liquids. This underlying theory is robust enough to provide a kinematic condition describing the viscosity terms. In addition, and as it will be described afterwards in the governing equation, the fluid is assumed to flow in a porous medium given by a Darcy-Forchheimer law in independent term. A similar hybrid Eyring-Powell and Darcy-Forchheimer flow has been explored in [11] to model a fluid in a stratified stretching sheet embedded in a porous medium. In this mentioned study, numerical solutions, supported by a Chebyshev spectral method, are obtained to the set of non-linear equations describing the flow under study.
The Eyring-Powell fluid model has been considered to model flow profiles in magnetohydrodynamics (MHD). Akbar [12] studied solutions for a two dimensional MHD fluid. In addition, Hina [13] analysed a MHD Eyring-Powell fluid with heat transfer. Later, Bhatti et al. [14] considered a permeable stretching surface for an MHD flow and made an analysis involving heat transfer. Finally, the reader is referred to the set of papers [15][16][17][18][19][20] for further analysis related with Eyring-Powell fluids. Note that other kinds of non-homogeneous diffusion (not supported by the Eyring-Powell conceptions) have been considered in [21] for a Darcy-Forchheimer flow.
As it will be described, the present study provides travelling wave solutions. In this regard, and in the search of analytical profiles to non-Newtonian fluids, the studies [22][23][24] provide important achievements. In our case, the convergence and stability results in the travelling wave profiles are provided by a Geometric Perturbation approach and are validated by a numerical exercise using the module bvp4c in Matlab. Note that the numerical approach and the use of the function bvp4c have been previously employed in fluid modelling. In [25], the authors provided a modelling analysis about the impact of thermophoretic particle deposition and magnetic dipole in the flow of a Maxwell liquid over a stretching sheet. To this end, they used a numerical scheme based on a Runge-Kutta-Fehlberg 45 (RKF 45) process with shooting technique in ODEs. The authors in [26] employed the bvp4c function to solve a bioconvection flow of Sisko nanofluid confined by a stretched surface. Some other studies in MHD flows made use of advanced numerical analysis: the authors in [27] analysed a mesh-free weak-strong (MWS) method to provide solutions for a MHD flow in a pipe with different geometries and with arbitrary conducting walls. In [28], a meshless local Petrov-Galerkin method was used to search for solutions in an unsteady MHD flow for different values in the wall conductivity and for different Hartmann numbers. Further, in [29], an MHD flow was solved by a combination of finite volume method and spectral techniques. Another novel numerical analysis was provided in [30] where Crank--Nicolson schemes supported by energy methods were employed to study the convergence and stability of the solutions.
The article structure is as follows: Firstly, a discussion about the fluid model proposed is introduced. Afterward, regularity of solutions is explored. Profiles of solutions are obtained based on the theory of travelling waves and the Geometric Perturbation Theory. Finally, a numerical process aims at validating the analytical approach proposed.

Fluid Model Principles
To start with, assume a one-dimensional, incompressible, unsteady and electrically conducting Eyring-Powell fluid. The velocity field is given by is the first velocity component depending of a transverse direction y. The geometry under study can be described as a big porous media with volume V and porous structures of typical size l, such that V 1 3 l. Given any interstitial distribution in the medium, a certain initial velocity is given as it will be introduced. The intention is to characterize the evolution of the velocity profile flowing through the x−axis and while varying with the transversal direction y. The hypothesis of a one dimensional flow is typical in flowing principles of pipes of stretching surfaces (see [31] for a MHD fluid flowing along rectangular and circular pipes).
The continuity equation and constitutive equation for Eyring-Powell fluid are given by: and: where ρ is the fluid density, − → J is the current density, − → B is the magnetic field and − → A is given by where p is pressure the pressure distribution, τ ij is the shear stress tensor, − → I is the identity tensor, µ 1 the magnetic permeability, − → E is the electric field and σ is the electrical conductivity. The magnetic field can be expressed as: are the imposed and induced magnetic fields, respectively. In addition, the shear stress tensor τ ij for an Eyring-Powell fluid model is given by (refer to [16]): where µ dynamic viscosity, β and d 1 are the characteristics of Powell-Eyring model and v i the velocity component for i = 1, 2, 3. Considering: The governing equation for i = 1 in the absent of induced magnetic field can be written as (refer to [15][16][17]): where ν = µ ρ is the kinematic viscosity, K is the permeability (hydraulic conductivity) of the porous medium, K 1 is the inertial permeability, C 1 is a fluid characteristic constant and B 0 refers to the magnitude of the imposed magnetic field as describe above. Let us introduce non-dimensional quantities (refer to [32]) as follows: Substitute (8) into into Equation (7) (ignoring * for simplicity): where R e = U 0 L ν is the Reynolds number and M = B 0 L σ ρν is the Hartmann number. Now, differentiate the Equation (9) with regards to x, so that: After using the value of − dP dx in Equation (9): As pointed previously, it is required to introduce a general initial velocity distribution that is related with the interstitial structure of the medium. Mathematically, this is considered as follows: where v 1 → 0 + at the pseudo-boundary in |y| → ∞ and v 10 (y) refers to the mentioned initial velocity distribution.

Existence and Uniqueness of Solutions
Let us consider a test function Assume a finite r 0 ∈ R + and define a ball B r accordingly, such that r r 0 . The with the following boundary and initial like conditions: The following Theorem shows the existence of solutions to the defined problem (13).
Proof. Consider ζ ∈ R + . The following cut-off function is defined: Multiplying the Equation (13) by ψ ζ and integrating in B r × [τ, T], we obtain: As the intention is to prove uniqueness and regularity under a global time evolution for m > 1, assume some large r r 0 > 1 and with t 1 ( [33,34]): Considering the spatial variable y close to ∂B r , it can be assumed y ∼ r. Then, for m = 2, it holds: Then, the following integral can be bounded as: As r 1 and taking φ 1 sufficiently small such that ∂φ 1 ∂y ψ ζ 1 over ∂B r : and: Integrating, we get: Using (17) and (18) in (14), we obtain: Next, consider a test function φ 1 of the form: we can choose a in such a way that (18) is convergent, therefore: For a > 4 and r → ∞: As the previous integrals are positive, we conclude on the null values of the involved integrals for r → ∞. As a consequence, the previous integrals are regular and finite for any other r < ∞. Considering (21) and the expression (18) for any r < ∞: where Υ is a suitable constant that exists in accordance with the integrals shown in (21). As both integrals are finite in τ < s < t < T, it is possible to conclude on the Theorem postulation about the boundness of solutions in B r ×(0, T]. The next intention is to continue exploring the regularity of the proposed Equation (10) by showing the boundness of ∂v 1 ∂y .
Proof. Multiplying Equation (10) by v 1 and using integration by parts: which implies that: After integration in both sides, the following holds: Based on Theorem 2.1 results, the right hand side of Equation (23) is bounded, therefore we can choose A 2 such that: which implies that ∂v 1 ∂y is bounded.
The coming intention is to show the uniqueness of solutions. To this end, the following Gronwall's inequality is required: Proof. Suppose that v 1 is the maximal solution to Equation (10) in R × (0, T) such that: with α > 0 arbitrary small. In addition, a minimal solution to Equation (10) is given for the initial condition: v 1 (y, 0) = v 10 (y).
The maximal and minimal solutions satisfy the following equations: Then, for every test function φ 1 ∈ C ∞ (R), and after subtraction, the following expressions hold:

Based on Theorems 1 and 2 results, we can choose A 3 and A 4 such that
Consider the test function: where A 4 and b are constants. Differentiate φ 1 with regards to s and y, to obtain: Then: Using (31) into (28), the following holds: After making the differentiation with regards to t: Define: Putting (32) into (31): with with α positive and sufficiently small. After applying Gronwall's inequality (see Proposition 1) on (33), the following holds for t ∈ [0, T]. As α is sufficiently small then h(t) = 0, i.e., v 1 = v 1 , which shows the uniqueness of solutions as stated.

Travelling Waves Existence and Regularity
The travelling wave profiles are described as v 1 (y, t) = g(ζ) where ζ = y − ct ∈ R, c refer to the travelling wave speed and g : R → (0, ∞) is the travelling wave profile that is requested to satisfy g ∈ L ∞ (R).
Equation (10) is then transformed in accordance with the travelling wave change as: (34) with g (ζ) < 0, by hypothesis. Now, consider the new variables: so that the following system holds: To analyse the propose system nearness the critical points, establish X = 0 and Y = 0, yielding: G , Consequently, (X 1 , 0) and (X 2 , 0) are the system critical points. The intention, now, is to make use of the Geometric Perturbation Theory to characterize the obtained critical points and to determine any solution profile in the vicinity of such critical conditions.

Geometric Perturbation Theory
In this section, we use the Singular Geometric Perturbation Theory to show the asymptotic behaviour of solutions associated to hyperbolic manifolds to be defined.
For this purpose, assume the following manifold: under the flow (36) and with the same critical points (X 1 , 0) and (X 2 , 0). Now, the following perturbed manifold, N α , close to N 0 in the critical point (X 1 , 0) is defined as: where α denotes a perturbation parameter close to equilibrium point (X 1 , 0) and F 1 is a suitable constant, which is found after root factorization. Firstly, consider X 2 = X − X 2 . Our intention is to apply the Fenichel invariant manifold theorem in [35] as formulated in [36,37].
For this, we shall show that N 0 is a normal hyperbolic manifold, i.e., the eigenvalues of N 0 in the linearized frame close to the critical point, and transversal to the tangent space, have non-zero real part. This is shown based on the following equivalent flow associated to N 0 : The eigenvalues are both real ± √ F 1 α , which show that N 0 is hyperbolic. Now, we want to show that the manifold N α is locally invariant under the flow (36). As a consequence of this, the manifold N 0 can be regarded as an asymptotic manifold to N α . For this, we consider the functions: , in the proximity of the critical point (X 1 , 0). In this case, δ is determined based on the following flows that are considered to be measurable, a.e., in R : Since the solutions are bounded, we can call δ = F 1 X 2 that represents a finite value. Then, the distance between the manifolds holds the normal hyperbolic condition for δ ∈ (0, ∞) and for α sufficiently small close the critical point (X 1 , 0). Now, we consider the following perturbed manifold N γ , close to N 0 in the critical point (X 2 , 0), defined as: where γ denotes a perturbation parameter close to the equilibrium (X 2 , 0) and B is a suitable constant found after root factorization. Assume X 2 = X − X 2 , then the Fenichel invariant manifold theorem can apply in the same manner as for the critical point (X 1 , 0). Note the following equivalent flow associated to N 0 : The eigenvalues are both real ± √ Bγ . This shows that N 0 is a hyperbolic manifold. Now, the manifold N γ is locally invariant/under the flow (36), so that the manifold N 0 can be expressed as an asymptotic manifold to N γ . Consider the functions: , in the proximity of the critical point (X 2 , 0). In this case, β is determined based on the following flows that are considered to be measurable, a.e., in R : Since the solutions are bounded, β = B X 1 is finite. Then, the distance between the manifolds holds the normal hyperbolic condition for β ∈ (0, ∞) and γ sufficiently small close the critical point (X 2 , 0).

Travelling Waves Profiles
Based on the show normal hyperbolic condition in the manifold N 0 under the flow (36), convergent travelling wave profiles can be obtained. For this purpose, let us consider firstly the Equation (36) such that the following family of trajectories in the phase plane (X, Y) holds The intention is to determine a trajectory in the phase plane in the proximity of the equilibrium (X 1 , 0). This is shown based on a comparison with subsolutions for a sufficiently small travelling wave speed and supersolutions for a sufficiently large speed together with a topological argument and the continuity of H. Assume c → 0, then it is possible to find a suitable value of A 1 such that dY 1 dX 1 > 0 while when c 0, it is possible to conclude on a condition of the form dY 1 dX 1 < 0 for suitable values in the involved constants. Given the continuity of H, it is possible, hence to conclude on the existence of a critical trajectory close the critical point (X 1 , 0) of the form The last expression shows the existence of an exponential profile in the travelling wave frame. Now, we want to show that the defined supporting manifolds preserve the exponential behaviour close the critical points. For this purpose, consider the Equation (38), so that: After solving Equation (45), we get: From the (38), the Equation (45) becomes: After solving (46): Then, the Equation (45) becomes: , which shows that the existence of an exponential profile along the travelling wave frame holds. Now, we want to show that the defined supporting manifold N γ preserves the exponential behaviour close to the critical points. For this purpose: After solving Equation (49), we get: From the expression (40), the Equation (50) becomes: After solving Equation (51), we have: From (35), the Equation (45) becomes: This last expression permits to show the conservation of the exponential profile in the critical points defined by the asymptotic manifold N γ .
It shall be noted that the travelling wave profiles (44) and (48) defined in the proximity of the critical points (X 1 , 0) and (X 2 , 0), respectively, follow the same behaviour.

Numerical Validation
The idea in this section is to show that the analytical assessments done in previous sections are accurate to support asymptotic solutions to the Equation (10). To this end, a numerical exercise is provided to solve Equation (10) and then compared with the analytical solutions in expressions (44) and (48).
The numerical routine introduced in this section is provided for moderate values in the Reynolds number (to be representative of a porous medium) and in the assumption of: The idea is to provide a validating exercise of the analytical assessments for certain values in the fluid parameters. The numerical exercise is not intended to provide solutions and to discuss them for different combinations in the parameters.
The following properties are of relevance: • The numerical approach is done with the Matlab function bvp4c. This function is based on a Runge-Kutta implicit approach with interpolant extensions [38]. The bvp4c has a collocation method at the pseudo-boundary conditions given by g(−∞) = 1 and g(∞) = 0. • The influence of the pseudo-boundary conditions and the collocation method shall be minimized. To this end, the integration domain is sufficiently large (−1000, 1000). • To make the problem tractable in terms of computational efforts, the integration domain has been split into 100,000 nodes with an absolute accumulated error of 10 −4 . With this level of discretization, the problem has been simulated with standard computers and with reasonable computational lead times.
The results are provided in Figures 1 and 2. The solutions are close for moderate values in the Reynolds number independently of the value of a. Note that the asymptotic convergence of travelling wave profiles to the numerical solution is given for ζ as sufficiently big. As described in each figure footprint, the criteria for asymptotic convergence between the analytical and numerical solutions is considered as the global distance ≤ 10 −3 . For this value and making numerical explorations, a dedicated value of ζ min is found, such that for ζ ≥ ζ min , the global error is preserved. is the asymptotic solution as per expressions (44) and (48) and f (ζ) is the numerical resolution of Equation (36). For Re = 10 (left), the global error between the analytical and numerical solution is ≤10 −3 for ζ ≥ 1.73, while for Re = 100 (right), the global error is ≤10 −3 for ζ ≥ 5.24. is the asymptotic solution as per expressions (44) and (48) and f (ζ) is the numerical resolution of Equation (36). For Re = 10 (left), the global error between the analytical and numerical solution is ≤10 −3 for ζ ≥ 1.13, while for Re = 100 (right), the global error is ≤10 −3 for ζ ≥ 2.34.

Conclusions
The analysis discussed in this paper has provided results on the regularity of solutions to a fluid flow formulated with a degenerate diffusivity Eyring-Powell viscosity term and a Darcy-Forchheimer law in porous medium. Travelling wave profiles have been obtained and asymptotic solutions have been shown to hold based on the Geometric Perturbation Theory. The analytical assessments presented have been validated with a numerical exploration applicable for moderate values in the Reynolds number, i.e., in medium of high porosity. As a main finding to conclude upon, we highlight the existence of an exponential profile of the solution under an asymptotic approximation. This is not a trivial result and reflects that such an exponential profile holds as well for the case of degenerate diffusivity introduced. As a future research topic related with the presented analysis, the authors would like to bring attention to the introduction of further advanced numerical analyses in line with the references cited. In addition, other analytical techniques based in perturbation to solitons, (G'/G)-expansion methods or maximal-minimal profiles may be explored to further precise the exponential solution found.