Analysis and profiles of travelling wave solutions to a Darcy-Forchheimer fluid formulated with a non-linear diffusion

1.   Introduction

The proposed analysis discussed along this paper is based on a Darcy-Forchheimer fluid equation extended with a non-linear diffusion.

The kind of materials presenting pores, that may be filled by a fluid, are typically referred as porous materials (or simply porous media). In 1856, Henry Darcy set the foundation of flow passing through such porous media, when he was working on flow of water through sand beds. His theory was known as Darcy Law. Later on, it was observed that the Darcy law experienced limitations to correctly model the inertial and boundary effects at high flow rates. To avoid this phenomena, in 1901 Forchheimer ^[1] extended the results of Henry Darcy concluding on a type of fluid known as Darcy-Forchheimer. The study of porous media has noticeable implications for engineers, geologists, architectures, and has important implications in mathematical modelling. As a set of representative examples discussing porosity and applications, it is possible to discuss about water in reservoirs, oil production, heat exchanger and catalytic reactors. In this direction, some important studies about Darcy-Forchheimer models can be consulted in ^{[2,3,4,5,6,7,8,9,10,11,12,13,14]}.

The Darcy-Forchheimer model is a focus of research currently. In ^[3], the heat transfer principles of a Darcy-Forchheimer flow in nanofluids with radiation are studied. Further, in ^[4], the authors provide a description of the convective flow of an MHD nanoliquid in an odd-shaped cavity filled with carbon nanotube-iron. In both cases, the modelling equations were provided in accordance with the well known Forchheimer model. In ^[5], a macro-scale momentum equation is presented departing from the Navier-Stokes equations to derive a free from turbulence model in porous medium. The cited references have been considered based on the introduced novelty at dealing with porous medium. In all cases, the driving equation were based on the classical order two diffusion (with gaussian kernel). The Darcy-Forchheimer model can be described as well with non-linear diffusion principles. In ^[29], a Darcy law is employed, supported by a non-linear diffusive formulation, to study the behaviour of nanofluids. In other cases, the non-linear diffusion has been introduced to account for specific properties to model, for instance this is the case of coagulation effects in blood vessels ^[31] or the peristaltic movement in a kind of Jeffrey fluid ^[30].

Alternatively, the main novelty introduced in this article is based on the formulation of a Darcy-Forchheimer fluid with a non-linear diffusion of the Porous Medium type. This kind of diffusion is also referred as degenerate due to the impossibility to show a maximum principle in a general sound. Some remarkable assessments related with non-linear diffusion can be consulted in ^[19,22] where important results, involving non-linear diffusion, are introduced (especially topics related with finite propagating speed and support). Note that the use of non-linear diffusion is well extended in the applied sciences and responds to the need of selecting an appropriate diffusion to correctly model the phenomena under study. As an example, non-linear diffusion appears in transverse problems related with the study of cells chemotasis in biology supported by the Keller and Segel model ^[25]. More recent studies can be consulted in ^{[24,26,27,28]}.

One the main objectives pursuit are related with the exploration of exact solution profiles. The proposed solutions dealt along this study are based on a travelling waves formulation. This kind of solutions have been explored in numerous studies in different applied sciences, as they permit to model the involved dynamic along a fixed shape profile moving with constant velocity. Such profile can be determined by the resolution of an Ordinary Differential Equation. Some interesting applications can be consulted in studies ^{[15,16,17,18]}.

As discussed, along the presented analysis we provide a non-linear diffusion formulation of a Darcy-Forchheimer flow so as to account for a diffusion fitting the modelling principles in a Porous Medium. Once the non-linear diffusion is introduced, solutions are shown to exist, to be regular and unique. Afterwards, solutions are explored in the travelling wave domain supported by the Geometric Perturbation Theory. Precise profiles of travelling wave solutions are shown. Such profiles are described asymptotically exhibiting an exponential behaviour to the travelling wave frame and in the proximity of the equilibrium solutions. The existence of an exponential profile is not trivial in the non-linear diffusion formulation and provide remnants of the classical gaussian order two diffusion. Finally, a numerical validation exercise is presented for particular values in the involved model parameters.

2.   Model formulation

We start from the classical 2D Darcy-Forchheimer model driven by the following equation:

where $\rho$ the density, $P$ is the pressure field, $\mu$ the dynamic viscosity, $\nu = \frac{\mu }{\rho }$ is the kinematic viscosity, $F$ the nonuniform inertia coefficient of porous medium, $\phi$ the porosity and $K_{1}$ the medium permeability.

Our aim is to provide a new formulation on the involved diffusion by the introduction of a non-linear term of the Porous Medium Type. This is supported by the fact that such non-linear diffusion has important properties such us finite propagation (see ^[19,22]) that can further fit the modelling exercise with the physical reality to explore in porous media. Hence, the driving equation is given by:

where $m > 1$. Making the differentiation of (2.3) with regards to $x$, the following holds:

Using this last obtained value for $-\frac{1}{\rho }\frac{dP}{dx}$ in (2.3):

with initial condition:

where $v_{0}(y)$ refers to the initial velocity. For $m = 1$, the expression (2.4) represents the classical Darcy-Forchheimer flow model formulated with a linear gaussian diffusion as provided in (2.1).

3.   Preliminaries

The degenerate diffusion requires the definition of a weak formulation for (2.4) to support the analysis of existence and uniqueness in the solutions:

Definition 1. Consider a test function $\phi _{1}\in C^{\infty }\left(R\right)$ such that for $0 < \tau < t < T$, the following weak formulation to (2.4) holds:

Based on (3.1), the following definition holds:

Definition 2. Given a finite spatial location $r_0$, admit a ball $B_{r}$ centered in $r_0$ and with radium $r > > r_{0}$. In the proximity of the borders $\partial B_r$ the following equation is defined:

with the following boundary and initial conditions:

4.   Regularity, existence and uniqueness of the solutions

Firstly, the following theorem enunciates the bounded properties of any solution to (3.1), i.e. solutions do not show global blow-up behaviour.

Theorem 4.1. Given $v_{0}\left(y\right) \in L^{1}\left(R\right) \cap L^{\infty}\left(R\right)$, then any solution to (3.1) is bounded for all $\left(y, t\right) \in B_{r}\times (0, T]$ with $r > > 1.$

Proof. Consider an arbitrary $\zeta \in R^{+}$, so that the following cut-off function is defined:

so that

being $A_a > 0$ an arbitrary constant. Now, Multiplying (3.2) by $\psi _{\zeta }$ and integrating in $B_{r}\times \left[\tau, T\right]$:

Note that for some large $r > > r_{0} > 1$ the following holds (^[19,22]):

where $C_1$ is a suitable constant depending on the initial distribution. For $m = 2$:

Now, it is the intention to assess the integral involving the non-linear diffusion term along the domain $B_r$. To this end:

Admit $r > > 1$ and that $\phi_1$ satisfies $\frac{\partial \phi_1}{\partial y} \psi_\zeta < < 1$ over $\partial B_r.$ For $r < \zeta,$ we have

Then

Using the expression (4.4) into (4.1), the following holds:

Now,

Next, consider a test function $\phi_1$ of the form:

Note that the exponent $\alpha$ is chosen in a way that (4.6) is convergent, therefore:

For $\alpha = \frac{2m}{m-1}$ and $r\rightarrow \infty,$ then we have

Considering the results in (4.9) and introducing into (4.1), the following holds:

which shows that the integrals at infinity are small given the decreasing behaviour of the associated test function. Similarly, it is possible to state that integrals are finite for finite value values of $r$ in $\tau < s < t < T$. Consequently, we conclude on the theorem postulation about the boundness of solutions in $B_r \times (0, T]$ for any $T > 0$.

To show the uniqueness in the solutions, the following theorem holds:

Theorem 4.2. Admit $v_{11}\left(y, t\right)$ and $v_{12}\left(y, t\right)$ are two solutions to Eq (2.4), with the same initial data i.e. $v_{11}\left(y, 0\right) = v_{12}\left(y, 0\right) = v_{0}\left(y\right) > 0.$ Then, solutions are unique, i.e. $v_{11}\left(y, t\right) = v_{12}\left(y, t\right)$.

Proof. To prove the uniqueness of solutions, admit a test function $\psi \left(y, t\right) \in C^{\infty }\left(Q_{T}\right)$ where $Q_{T} =$ $R\times \left(0, T\right)$. Firstly, we assume that $v_{11}\left(y, t\right) > v_{12}\left(y, t\right)$ and a weak solution to (2.4) defined as:

and

Subtracting Eqs (4.10) and (4.11):

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

where $K_{3} = \underset{\left(y, t\right) \in Q_{T}}{\max } \lbrace v_{11} \rbrace.$

Admit the following choice for the test function $\psi \left(y, t\right) \in C^{\infty }\left(Q_{T}\right)$:

where $\alpha$ is such that

For simplicity, we can choose a test function under a normalization condition of the form

Along $y\rightarrow \infty,$ the integral mass shall be null. Then, for $r > > 1$:

In the asymptotic $\left \vert y\right \vert \rightarrow \infty,$ we obtain:

then

Now, we define

where

The integral (4.20) is finite under the condition for $\alpha$ given in (4.19):

Integrating in both sides

where $K_{4}\left(\alpha \right) = 4\alpha \left(\alpha +1\right).$ Now, we consider the assessment of each of the integrals involved in (4.12):

Similarly,

and

Combining the different expressions for each of the assessed integrals, the following holds:

Note that $\left \vert \Theta \left(y\right) \right \vert < \infty$. For a sufficiently small $t > 0,$ $\sup \left \vert v_{11}-v_{12}\right \vert \rightarrow 0,$ then

which implies that $v_{11}\left(t\right) \leq v_{12}\left(t\right)$. This contradicts the initial assumption, therefore the only compatible condition is to admit $v_{11}\left(t\right) = v_{12}\left(t\right)$ as intended to prove.

5.   Travelling waves regularity and existence

The travelling wave profiles ($h$) are described as $v_{1}\left(y, t\right) = h\left(\eta \right)$ where $\eta = y-ct\in R,$ $c$ is the travelling wave speed and $h:R\rightarrow \left(0, \infty \right)$ belongs to $L^{\infty }\left(R\right).$

The expression (2.3) is then transformed to the travelling wave domain as:

with $h^{\prime }\left(\eta \right) < 0.$ Admit, now, the new variables (see ^[19] for further justifications):

so that the following system holds

The proposed system in analyzed in the proximity of the critical points. To this end, setting $X^{\prime } = 0$ and $Y^{\prime } = 0,$ yield

The solutions to (5.4) are computed by standards means:

and

Indeed, $\left(X_{1}, 0\right)$ and $\left(X_{2}, 0\right)$ are the system critical points. Now, the objective is to use the Geometric Perturbation Theory to characterize the orbits in the proximity of such critical points.

5.1.   Geometric Perturbation Theory

Along this section, the singular Geometric Perturbation Theory is employed to show the asymptotic behaviour of specific defined manifolds close the system critical points. Firstly, admit the following manifold:

with critical points $\left(X_{1}, 0\right)$ and $\left(X_{2}, 0\right).$ A perturbed manifold $N_{\epsilon }$ close to $N_{0}$ in the critical point $\left(X_{1}, 0\right)$ is defined as:

where $\epsilon$ denotes a perturbation parameter close to equilibrium $\left(X_{1}, 0\right)$ and $C$ is a suitable constant obtained after root factorization. Firstly, define $X_{3} = X-X_{2}$ to apply the Fenichel invariant manifold theorem ^[20] as formulated in ^[21,23]. For this purpose, it is required to show that $N_{0}$ is a normally 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_{\epsilon}:$

The associated eigenvalues are both real $\left(\pm \sqrt{C}\epsilon \right)$. This permits to conclude that $N_{0}$ is a hyperbolic manifold. Now, we shall show that the manifold $N_{\epsilon }$ is locally invariant under the flow (5.3), so that the manifold $N_{0}$ can be described by the asymptotic approach $N_{\epsilon }.$ Based on this, admit the following functions:

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

Since the solutions are bounded and measurable, $\delta = C\left \Vert X_{3}\right \Vert$ is finite. As a consequence, the distance between the manifolds is sufficiently small to hold the normal hyperbolic condition for $\delta \in \left(0, \infty \right)$ and $\epsilon \rightarrow 0$ close to the critical point $\left(X_{1}, 0\right)$.

Now, consider the perturbed manifold $N_{\gamma }$ close to $N_{0}$ in the critical point $\left(X_{2}, 0\right)$:

where $\gamma$ denotes a perturbation parameter close to equilibrium $\left(X_{2}, 0\right)$ and $D$ is a suitable constant which is found after root factorization. Now, admit $X_{4} = X-X_{2}$ to apply the Fenichel invariant manifold theorem in the same way as for the critical point $\left(X_{1}, 0\right)$. The equivalent flow associated to $N_{\gamma}$ is given by:

The associated eigenvalues are both real $\left(\pm \sqrt{D}\gamma \right)$ which show that $N_{0}$ is a hyperbolic manifold close $(X_2, 0)$. Now, it is required to show that the manifold $N_{\gamma }$ is locally invariant under the flow (5.3), so that the manifold $N_{\gamma}$ behaves asymptotically as $N_{0 }.$ To show this last condition, consider the functions:

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

Again, as the solutions are bounded and measurable, $\delta = D\left \Vert X_{4}\right \Vert$ is finite. Consequently, the distance between the manifolds keeps the normal hyperbolic condition for $\delta \in \left(0, \infty \right)$ and $\gamma$ sufficiently small close the critical point $\left(X_{2}, 0\right)$.

5.2.   Travelling wave profiles

Based on the normal hyperbolic condition in the manifold $N_{0}$ under the flow (5.3), asymptotic travelling wave profiles can be obtained.

For this purpose, consider firstly the flow (5.3) to obtain the following family of trajectories in the phase plane $\left(X, Y\right)$:

Based on the continuity of $H(X, Y)$ and the changing sign for $X$ sufficiently small and $X$ sufficiently big, it is possible to conclude on the existence of a critical trajectory of the form$:$

which implies that

Solving (5.12) by using standard separation of variables, the following holds:

where

and

Finally and after using (5.2):

which implies that

As it can been seen, the existence of an exponential profile along the travelling wave frame holds. This is not a trivial result for the non-linear diffusion case and reflects the existence of an exponential flow profile as solution to (2.3).

Now, the intention is to show that the defined supporting manifolds $N_{\epsilon}$ and $N_\gamma$ preserve the exponential behaviour close to the critical points. For this purpose, consider the expression (5.8) (idem can be shown for $N_\gamma$ as per (5.9)) so that:

Solving (5.8) by using separation of variables, the following holds:

From (5.8), the above equation becomes

Using (5.8) into (5.15):

Upon resolution:

From (5.2),

This last expression permits to show the conservation of the exponential profile close the critical points defined by the asymptotic manifolds $N_\epsilon$ and $N_\gamma$.

6.   Numerical validation of the analytical solution

The aim along the present chapter is to set a numerical procedure for the expression (5.1) and to compare with the analytical asymptotic solution obtained in expression (5.14) for different values in the travelling wave speed. To this end, an implicit Runge-Kutta algorithm with interpolant extension has been proposed based on the MATLAB software and a scheme proposed in ^[17]. The initial condition has been assumed as a step function, $v_0 (y) = H(-y)$ such that $v_0(y \rightarrow -\infty) = 1$ and $v_0(y \rightarrow -\infty) = 0$. These two conditions permit to study the heteroclinic connection between a positive state at $-\infty$ evolving towards a null state at $+\infty$. Note that to solve the Runge-Kutta algorithm, considering the heteroclinic connection, it is required to add pseudo-boundary conditions in the travelling wave domain given by $h(\eta \rightarrow -\infty) = 1$ and $h(\eta \rightarrow -\infty) = 0$. Hence and to avoid such boundary influence, the spatial travelling wave domain is considered sufficiently large, i.e. $\eta \in (-500,500)$. The domain is split in $10^5$ nodes with an absolute converging error of $10^{-5}$.

The results are presented for different values in the travelling wave speed $c$ (see Figures 1–3) to account for a general picture of solutions. The numerical exercise is performed only for validating the obtained analytical solution. Then the correctness of the analytical profile is checked for dedicated values in the involved model parameters. For the sake of simplicity:

while the travelling wave speed $c$ is a free parameter that can be adjusted in accordance with certain criteria in the involved analytical and exact numerical profiles.

As it can be observed, the structure of solutions follows a similar evolution along the asymptotic travelling wave tail. In addition and in the searching of an optimal travelling wave profile, the following principles have been followed:

Minimum distance between the exponential blends upon arrival to the null condition and, in the asymptotic (beyond $\eta = 1$), a maximum distance of $10^{-4}$ between both profiles.

The last two conditions are met for a travelling wave $c = 16.353$ that can be considered as the most appropriate to account for minimum global errors.

Based on the exposed figures, it is possible to conclude that the there exists an optimal travelling speed for which the analytical solution in (5.14) evolve close to the actual solution.

7.   Conclusions

First of all, assessments on existence, regularity and uniqueness of the solutions have been provided. Afterwards, the problem (2.3) has been studied based on purely analytical techniques applicable to understand the asymptotic behaviour of solutions. The analysis in the travelling wave domain has been complemented with by the Geometric Perturbation Theory to search for asymptotic solutions close the critical points. Even further, specific solution profiles have been shown analytically for the whole manifold evolution $N_0$ (5.7) and in the proximity of the critical points. The main finding is related with the existence of an asymptotic exponential behaviour along the travelling wave frame. Finally, a numerical exercise has been provided to validate the analytical solution (5.14) via exact simulation of the model (5.1).

Conflict of interest

The authors declare no conflict of interests.