On an equation arising in natural convection boundary layer flow in a porous medium

A generalization of a boundary value problem treated previously by Magyari et al. (ZAMP 53:782–793, 2002a; Transp Porous Medium 46:91–102, 2002b), Zhang (Math Anal Appl 417:361–375, 2014) and Paullet (Appl Math E Notes 14:123–126, 2014) for a prescribed wall temperature and by Zhang (Appl Math Comput 339:367–373, 2018) for a prescribed wall heat flux is considered. The problem involves a parameter λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} and an exponent p. Numerical solutions are obtained for representative values of λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} and p, a feature of which is the existence of critical values λc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _c$$\end{document} of λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} with two solution branches in λ>λc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda >\lambda _c$$\end{document}, with λc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _c$$\end{document} dependent on p. Asymptotic solutions for large λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} and large p are derived for both types of boundary condition. For large λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}, the nature of the solution is essentially different on upper and on the lower branches with similar feature being seen in the behaviour for p large.


Introduction
The boundary value problem θ + λ θ + θ 2 = 0, θ(0) = 1, θ → 0 as y → ∞, (1) where primes denote differentiation with respect to the independent variable y has been treated by Magyari et al. [1,2]. This problem arose in looking for similarity solutions in the free convection boundary layer flow on a vertical surface in a porous material. In particular, a wall temperature proportional to x m was assumed, where the coordinate x measures distance along the surface, with Eq. (1) appearing as the special case of m = −1. The analysis for this particular case was shown to be different to the other values of m and so is a limiting case for m. Equations of the type given in (1) have a relatively long history and have also arisen in flows on stretching surfaces, see for example the original work of Sakiadis [3], Crane [4] and Banks [5] with a full review being provided by Wang [6]. The problem described by Eq. (1) has also been considered in a more theoretical way by Zhang [7]. Several properties of the solution were established including the existence of a critical value λ c with two solution branches for λ > λ c , seen also in the numerical integrations of [1,2]. However, existence or nonexistence of solutions for 0 < λ < λ c remains an open question. Paullet [8] further proved that for some range less than λ c , no solution to Eq. (1) exists. Bounds on the critical value λ c were derived consistent with that computed by [1,2]. Zhang [9] also considered the 'prescribed wall flux condition', i.e. the same equation given in (1) but now subject to the boundary conditions As for the 'prescribed wall temperature condition' in equation (1), properties of the solution were established theoretically again including bounds on the critical value leading to the existence of two solution branches.
Our aim here is to extend the problems considered by Zhang [7,9] to the more general case for an exponent p > 1 subject to prescribed wall temperature condition given in (1) or to the prescribed wall flux condition (2). We assume throughout that λ ≥ 0 and, following Zhang [7,9], we limit our attention to those solutions for which the integral ∞ 0 θ(y) dy is bounded. As a consequence, here we consider only those solutions which are exponentially small, of O(e −λ y ), for y large. We note that there are also other solutions which are algebraic, of O(y −1/(p−1) ) for y large. Equations which involve a term of the type θ p , as in (3), have arisen in the somewhat different context of convective flow in a porous medium with local heat generation being modelled by θ p (p > 1), see [10][11][12][13][14] for example. We start with the case of a prescribed wall heat flux.

Prescribed wall heat flux
Here we examine the problem given by restricting our attention to those solutions which are exponentially small for y large. In Fig. 1a, we plot θ(0) against λ for p = 2, the case treated in [9], obtained from the numerical solution of Eq. (4). There is a critical value λ c of λ, where λ c 1.41, with two solution branches for λ > λ c . On the lower solution branch, θ(0) decreases slowly to zero as the value of λ is increased, whereas on the upper branch θ(0) increases very rapidly as λ is increased. In Fig. 1b On integrating Eq. (4) and applying the boundary conditions, we obtain so that θ(0) > λ −1 . This lower bound on θ(0) is seen for the lower branch in Fig. 2 and provides a reasonable estimate for the lower branch solutions seen in Fig. 1. We can calculate the critical values numerically following the approach described in [15,16], for example, and in Fig. 3a we plot the values of λ c and in Fig. 3b the values of θ c (0), the value of θ(0) at λ = λ c , against p. Our direct calculation of the critical values gives a better estimate of λ c 1.40920 for p = 2, with θ c (0) 1.18710, towards the higher end of the range 1.37824 < λ c < 1.41665 given [9]. Both the values of λ c and θ c (0) decrease as p is increased, as perhaps could be expected from Fig. 1, with λ c → 2 as p → 1.
We now consider the nature of the solution when λ is large as seen in Fig. 1.

Solution for λ large
We first consider the solution on the lower branch by writing Equation (4) then becomes where primes now denote differentiation with respect to ζ. The leading-order solution is H 0 = e −ζ . For the term H 1 at O(λ −(p+1) ), we obtain giving so that, on the lower branch consistent with the results shown in Fig. 1 and expression (5). We note that this result is, to leading order, independent of the exponent p and that the asymptotic limit is approached more rapidly as the value of p is increased.
For the upper branch solution, we put with ζ given in (6). Equation (4) becomes The leading-order problem is This is a homogeneous problem, and for a nontrivial solution, we suppose that h(0) = a 0 for some constant a 0 > 0 which we expect to depend on p. We then put so that Eq. (13) is now Equation (15) is, in effect, an eigenvalue problem for a 0 . In Fig. 4, we plot a 0 against p obtained from the numerical solution of Eq. (15). We see that a 0 reaches a maximum value of approximately 1.2418 at p 4.9 with a 0 → 1 slowly for large values of p and going almost to zero close to p = 1. For p = 2, a 0 0.8587 so that θ(0) ∼ 0.8587 λ 2 + · · · , for p = 5, a 0 1.1821 giving θ(0) ∼ 1.1821 λ 1/2 and for p = 10, a 0 1.19265 with θ(0) ∼ 1.19265 λ 2/9 + · · · as λ → ∞. The difference in the powers of λ in these cases accounts for the way the upper branch solutions increase with λ, rapidly for p = 2 (Fig. 1a) and much more slowly for p = 10 (Fig. 1c).

Solution for p large
If θ < 1, then the term θ p in Eq. (4) is negligible for p large and hence has the solution, to leading order, of so that θ(0) ∼ λ −1 + · · · as p → ∞ on the lower solution branch, as can be seen in Fig. 2 for λ = 1.5.
The nature of the solution on the upper branch for p large can be seen in Fig. 5 where we plot θ and θ against y for p = 61.8 and λ = 1.5. For these parameter values, θ(0) 1.06346. We see that θ decreases monotonically from its wall value. However, there is a thin region next to the wall where θ decreases rapidly before increasing to its outer value. To describe this solution, we start in the outer region by putting y = Y 0 + y, where Y 0 depends on p. This leaves equation (4) unaltered and to the leading-order solution θ = e −λ y in the outer region.
We can calculate Y 0 from showing that the thickness of the inner region decreases as the value of p is increased. The asymptotic expression given in (21) gives Y 0 0.0445 for p = 61.8 and p = 1.5 in good agreement with the value 0.0446 obtained from the numerical solution of (4). The above analysis indicates why λ c → 1 as p → ∞, as indicated in Fig. 3a. Expression also (20) shows that θ c (0) → 1 in this limit, as can seen in Fig. 3b.

Prescribed wall temperature
Here we consider the problem for p > 1, where primes again denote differentiation with respect to y and again limiting our attention to those solutions which are exponentially small for y large. In Fig. 6, we plot θ (0) against λ for p = 2, Fig. 6a, p = 5, Fig. 6b and for p = 10, Fig. 6c, obtained from the numerical solution of Eq. (22). In each case, we see the existence of a critical value λ c of λ with λ c 1.079 for p = 2, in agreement with the value given in [1], λ c 0.648 for p = 5 and λ c 0.404 for p = 10. The critical points appear to arise where θ (0) = 0 and gives rise to two solution branches in λ > λ c . On the lower branch, the values of θ (0) are negative and appear to decrease linearly with λ. On the upper branch θ (0) > 0 with the values of θ (0) increasing with λ, much more rapidly for p = 2 than for p = 5 or for p = 10. In Fig. 7, we plot θ (0) against p for λ = 1.5. Again we see the existence of a critical value p c , here at p c 1.238 with two solution branches in p > p c . The lower solution branch has θ (0) negative and approaches a constant value which appears to be −λ as p is increased. The upper branch has θ (0) > 0, reaches a maximum value of approximately 3.144 at p 2.059 before decreasing as p is increased and appears to approach the value of λ very slowly for large p, seen in numerical solutions to much larger values of p than used to plot Fig. 7.
We can integrate Eq. (22) in this case to obtain giving a lower bound θ (0) > −λ and providing a reasonable estimate for the lower branch solutions seen in Figs. 6 and 7.
The critical values λ c are determined by making a linear perturbation θ 1 to Eq. (22). This leads to where θ c is the solution to (22)  to Eq. (24), we can construct the full solution as which can be seen in Fig. 3b.
Equation (22) then gives The leading-order problem satisfies homogeneous boundary conditions and, if we assume that g (0) = b 0 (b 0 > 0) and Equation (33) is then an eigenvalue problem for b 0 . In Fig. 9,

Solution for p large
For the lower branch solutions θ < 1 on 0 < y < ∞ and so the term θ p becomes negligible for p large. Solving the resulting equation gives θ = e −λ y so that as can be seen in Fig. 7. The structure of the solution for p on the upper branch can be seen in Fig. 10 where we plot θ and θ against y for p = 96.94. There is an outer region y ≥ Y 1 (p) where θ ≤ 1 and the term θ p in Eq. (22) is negligible. If we put y = Y 1 + y, Eq. (22) becomes, to leading order, θ + λ θ = 0, θ → 0 as y → ∞, where primes now denote differentiation with respect to y. Equation (36) is solved subject to the condition that θ = 1 at y = 0, i.e. defining the position of Y 1 as where θ = 1. This gives θ = e −λ y . For the inner region, θ ≥ 1 and the term λ θ in Eq. (22) is negligible. Hence we again have equation (17) at leading order now to be solved subject to, on matching with the outer region, This gives From (38) θ (0) ∼ λ + · · · as p → ∞.