The effects of wall transpiration on the unsteady free convection boundary-layer flow near a stagnation point in a heat generating porous medium

The natural convection boundary-layer flow near a stagnation point on a permeable surface embedded in a porous medium is considered when there is local heat generation within the boundary layer at a rate proportional to (T-T∞)p,p≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(T-T_\infty )^p,\ p \ge 1$$\end{document}, where T is the fluid temperature and T∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_\infty$$\end{document} the ambient temperature. There is mass transfer through the surface characterized by the dimensionless parameter γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document}, with γ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma >0$$\end{document} for fluid injection and γ<0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma <0$$\end{document} for fluid withdrawal. The steady states are considered where it is found that, for p>1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p >1$$\end{document}, there is a critical value γc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma _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}$$\gamma$$\end{document} with solutions existing for γ≥γc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \ge \gamma _c$$\end{document} if 1<p<2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1<p<2$$\end{document} and for γ≤γc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \le \gamma _c$$\end{document} if p>2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p>2$$\end{document}. The initial-value problem reveals that, for 1≤p<2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 \le p<2$$\end{document}, the nontrivial steady states are stable and the solution evolves to this state at large times. However, for p>2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p>2$$\end{document} these steady states are unstable and the solution either approaches the trivial state with the local heating dying out or a finite-time singularity develops for sufficiently large initial inputs.


Introduction
During the transport and storage of a porous material heat can be generated locally within the medium. This can set up a convective flow within the material which may help to dissipate the heat or alternatively can lead to enhanced heat generation and even thermal runaway with possibly disastrous consequences. Examples include the spontaneous ignition in stock piles of coal [1][2][3][4], in bagasse (the cellulose waste left after the extraction of sugar from sugar cane) [5,6] or in the wetting of cellulosic materials [7][8][9]. Previously we have modelled this problem assuming local heat generation at a rate proportional to ðT À T 1 Þ p , where T 1 is the (constant) ambient temperature and the exponent p ! 1 [10,11]. We have assessed how the effects of the input of heat from the boundary [12] and of an outer flow [13] control the local heat generation, finding that this depended critically on the exponent p and the magnitude of the initial heat input. More recently we have considered modified form of Arrhenius kinetics [14], again seeing that the evolution of the local internal heat generation depended on the initial heat input as well as the local heating rate, finding conditions under which the local heating died away or produced a finite-time blow up in the solution.
The power-law expression for localised heat generation has been used previously in several other contexts as a model for Arrhenius kinetics. These studies have shown at least qualitative agreement with those using the more general Arrhenius form, as seen for a specific case by comparing [12] and [14]. The use of a power-law model is motivated in part by its analytic simplicity as compared to the Arrhenius function and to avoid the 'cold boundary problem'. In this, if an inert outer region is required, as is the case here, the Arrhenius form has to be modified, as in [14]. This can, perhaps, to lead to some artificiality in the kinetic scheme used. An alternative form for the local heating has been suggested, based more on fluid flow, is to model the local heat generation as a source term with a particular spatial variation [15][16][17][18]. However, this approach has the drawback, at least for the present case of spontaneous combustion within a porous material, in that it does not allow the local heat generation to vary with the evolving temperature and velocity fields.
Here we consider how the input from the boundary influences the heat generation within the flow. We again take the above power-law term to model the heat generation within the boundary layer and allow for fluid injection or withdrawal from a thermally insulated surface. Previously, without any flow through the wall [12] and without an outer flow, we found that, when 1 p\2, the local heating died away. However, for p [ 2 a finite-time blow up in the solution could occur for a sufficiently large initial input, otherwise the local heating also died out. Our aim here is to determine how this behaviour is influenced by the fluid input or withdrawal from the wall. We find that, when p ¼ 1, a nontrivial steady state is attained giving a balance between the rate of local heat generation and the wall mass transfer. A similar situation applies when 1\p\2 provided now the rate of mass transfer through the wall is greater than some critical value, dependent on p. Otherwise the heat generation is inhibited. For p [ 2 there is a either finite-time thermal runaway if the initial input is above some threshold value or otherwise the local heating dies away. We start by deriving our model.

Model
We consider the unsteady, two-dimensional natural convection boundary-layer flow that can arise near a forward stagnation point in a fluid-saturated porous material in which there is local heat generation within the boundary layer at a rate proportional to ðT À T 1 Þ p , where T is the fluid temperature and T 1 is the ambient temperature. We assume that the bounding surface is thermally insulated and that there is normal wall velocity v w , where v w can be either positive or negative. We further assume the flow is given by Darcy's law and we make the standard Boussinesq approximation. The basic equations for our model are, see [12,[19][20][21] for example, together with some initial condition and where we assume that p ! 1. Here x measures distance along the body surface and y normal to it, u and v are respectively the velocity components in the x and y directions. Also K is the permeability of the porous media, m the kinematic viscosity of the fluid, b the coefficient of thermal expansion, g the acceleration due to gravity, a m the effective thermal diffusivity and r the heat capacity ratio of the fluid-filled porous medium to that of the fluid. In Eq. (2), S(x) is the sine of the angle between the outward normal and the downward vertical and for a stagnation-point flow To make Eqs. (1-4) dimensionless we write where T s ; U s are respectively temperature and velocity scales given by Equation (2) becomes (on dropping the overbars) u ¼ x h and the continuity Eq. (1) is satisfied by introducing a streamfunction w defined so that u ¼ w y ; v ¼ Àw x . We then write w ¼ x f ðy; tÞ, and consequently h ¼ hðy; tÞ, to obtain subject to the boundary conditions where c ¼ v w R 1=2 =U s . In (8), c [ 0 gives fluid injection from the boundary (blowing) and c\0 gives fluid withdrawal from the boundary (suction). We find different behaviour depending on whether p ¼ 1; 1\p\2 or p [ 2, treating these cases, as well as the transition case when p ¼ 2, see expression (6), separately. We start by considering the steady states of Eq. (7) as these can represent the possible large time behaviour of the full initial-value problem.

Steady states
These represent possible large time limit of the initialvalue problem (7,8) and, for a general value of p, are given by where primes denote differentiation with respect to y. We note that, for p ¼ 2, Eq. (9) can be integrated to give on satisfying the boundary condition as y ! 1.
Applying the boundary conditions on y ¼ 0 then shows that the only possibility when c 6 ¼ 0 is the trivial solution f Àc with a nontrivial solution only for c ¼ 0.
In Fig. 1 we plot f 0 ð0Þ against c for p ¼ 1 obtained from the numerical solution of Eq. (9). The numerical solution terminated at c ¼ À2, where f 0 ð0Þ becomes zero (we were only able to obtain the trivial solution f Àc for c\ À 2), and proceeded to large positive values of c with the values of f 0 ð0Þ increasing as c is increased.
We plot f 0 ð0Þ against c in Fig. 2 for p ¼ 1:5, representative of values of p in the range 1\p\2. In this case we see that there is a critical value c c of c, where c c ' À0:5085 giving two solution branches in c c \c. The upper branch solution continues to large positive c, with f 0 ð0Þ increasing as c is increased, and we find only the trivial solution for c\c c . This leads us to consider the critical values in more detail. To calculate c c numerically we follow the approach given in [22], for example, whereby we make a linear perturbation to Eq. (9) resulting in a linear homogeneous equation subject to homogeneous boundary conditions and it is the solution to this eigenvalue problem that determines c c . We plot c c against p in Fig. 3 where we see that c c ! 0 as p ! 2, as perhaps might be expected, and approaches a finite negative value as p ! 1, with our numerical results suggesting that c c is approaching À2.
We now consider the nature of the solution as p ! 2 from below. We put p ¼ 2 À , with then c ¼ l , and look for a solution valid for small by expanding At leading order we obtain f 0 ¼ 2b 0 tanhðb 0 yÞ for some constant b 0 [ 0 to be found. At OðÞ we have We can integrate Eq. (12) to obtain, on applying the boundary conditions and the expression for f 0 , Meccanica (2018) 53:759-772 761 on using a result given in [12]. From (13) we have Expression (14) gives Hence Asymptotic expression (15) is shown in Fig. 3  We see that c c is positive throughout for this case, approaching zero as p ! 2 from above, and that c c increases as p is increased. The above discussion for the asymptotic behaviour as p ! 2 follows directly for this case on putting p ¼ 2 þ . The result is that This asymptotic result is shown in Fig. 5 by a broken line, again showing good agreement with the numerically determined values. We next consider how the solution behaves with p for a fixed value of c, taking c ¼ 1:0 representative of blowing and c ¼ À1:0 representative of suction. We plot f 0 ð0Þ against p in Fig. 6 for these two cases, in Fig. 6a for p in the range 1 6 p\2 and in Fig. 6b for p [ 2. In Fig. 6a both curves start with their values at p ¼ 1 as shown in Fig. 1 with, for c ¼ 1:0; f 0 ð0Þ increasing and becoming large and the boundary-layer thickness decreasing as p is increased with our numerical integrations suggesting that f 0 ð0Þ ! 1 as p ! 2. For c ¼ À1:0, however, f 0 ð0Þ decreases to the critical point at p ¼ p c ' 1:220 for c c ¼ À1:0 shown in Fig. 3, leading to a lower solution branch in p\p c , not particularly clear in the figure.
For p [ 2, Fig. 6b, there is now a critical value for c ¼ 1:0 at p ¼ p c ' 3:938 giving rise to solutions only in p > p c and also a lower solution branch in p [ p c . For c ¼ À1:0 the values of f 0 ð0Þ increase rapidly as p ! 2, similar to that seen for c ¼ 1:0 in Fig. 6a. In both cases the values of f 0 ð0Þ appear to be approaching the same limit as p is increased.
To see how the solution for c [ 0 in Fig. 6a behaves as p ! 2 from below we put p ¼ 2 À d and look for a solution for d ( 1 by writing f ¼ aðdÞ F; f ¼ aðdÞ y for some aðdÞ to be determined, expecting that a ) 1 for d small. Equation (9) becomes Primes now denote differentiation with respect to f. Since a À2d $ 1 À 2d log a þ Á Á Á, we look for a solution by expanding At leading order we have for some constant c 0 [ 0 to be found. To obtain a solution at next order we require This then gives at Oðd log aÞ, on integrating the equation once and applying the boundary conditions as f ! 1, If we now apply the boundary conditions on f ¼ 0, particularly that F 1 ð0Þ ¼ Àc, and the above expression for F 0 we obtain We can readily adapt this analysis to the case when c\0 and p ! 2 from above, Fig. 6b. We now put p ¼ 2 þ d and make the above transformation. We then follow the above discussion, the only difference being in a sign change on the right-hand side of Eq. (21). This in turn leads to, on applying the boundary conditions on f ¼ 0, with the expression for f 0 ð0Þ being the same as that given in (23).

Strong fluid withdrawal, jcj large
We have seen that, when p [ 2, the solution proceeds to large negative values of c and to obtain a solution valid in this limit we put f ¼ jcj þ jcj ð3ÀpÞ=ðpÀ1Þ G; Y ¼ jcj y. Equations (9) become G 000 þ ð1 þ jcj 2ð2ÀpÞ=ðpÀ1Þ Þ G 00 þ G 0 p ¼ 0; This gives, at leading order, on writing u ¼ G 0 , Now suppose uð0Þ ¼ u 0 . We then put u ¼ u 0 U to obtain the eigenvalue problem for u 0 as A plot of u 0 against p is shown in Fig. 7 obtained from the numerical integration of Eq. (27). The graph increases from its value of u 0 ¼ 0:8587 at p ¼ 2, has a maximum value of u 0 ' 1:2418 at p ' 4:9 before decreasing slowly to u 0 ¼ 1 as p ! 1. We note that the numerical integration of (27) continues smoothly through p ¼ 2 into 1\p\2 though this asymptotic solution is valid only for p [ 2. Hence 4 Initial-value problem Our discussion of the steady states suggests that we should treat the cases p ¼ 1; 1\p\2 and p [ 2 separately. We solved initial-value problem (7,8) numerically using the scheme described in [10,12] for example, taking as our initial condition of oy ¼ a 0 e Ày and prescribing a value for a 0 [ 0.

p ¼ 1
Here we took a 0 ¼ 1:0 and in Fig. 8 we plot values of the wall temperature h w hð0; tÞ against t obtained from the numerical solution of initial-value problem (7,8) for representative values of c. We see that these curves approach at large times the constant, nontrivial value given by the corresponding steady state solution shown in Fig. 1. For values of c À 2 the numerical solutions approached the trivial state f Àc as t increased.

1\p\2
Here we took p ¼ 1:5, as in Fig. 2, and in Fig. 9 we plot the wall temperature h w for representative values of c. In this case, for c\0, there is the possibility of two steady states. Our numerical integrations indicate that it is the solution on the upper solution branch that is approached at large times, indicating that the saddlenode bifurcation at c ¼ c c changes the temporal stability from stable on the upper branch to unstable on the lower branch. For values of c\c c ' À0:5085 we found only the trivial solution f Àc at large times.

p [ 2
In this case we found different behaviour to the previous cases in that now the solution either becomes singular at a finite time or approaches the trivial state at large times. We illustrate this in Fig. 10 for p ¼ 3:0, taking values of c\c c ' 0:6082. In Fig. 10a we plot the wall temperature h w for representative values of c, here taking a 0 ¼ 3:0 in the initial condition. In each case we find that h w increases very rapidly as the singularity at t ¼ t s is approached, finally becoming unbounded. To attain this finite-time singularity require a sufficiently large initial input, i.e. having a 0 greater than some threshold value. We illustrate this in Fig. 10b with plots of h w for c ¼ À1:0, again with p ¼ 3:0. With a 0 ¼ 2:18 we again see the the development of the finite-time singularity. However, with a 0 ¼ 2:17 the solution approaches the trivial state at large times.
In Fig. 10c we plot the times t s at which the finitetime singularity occurred, at the values of c shown by þ. The numerical method involved a procedure for halving the time step Dt to maintain accuracy, with the numerical integration terminating when Dt\5 Â 10 À6 . The values of t s in Fig. 10c are the values of t at this final time step. The values of t s depend on the value of a 0 taken, as well as whether a finite-time singularity occurs, Fig. 10b. For Fig. 10c we took a 0 ¼ 3:0, noting that, for example, with c ¼ À3:0 the solution approached the trivial state, however, with a 0 ¼ 5:0 a finite time singularity was seen. Thus it appears from our numerical investigation that, even though nontrivial steady states exist for p [ 2, these are unstable with the result that the solution evolves to either the trivial state or becomes singular at a finite time. This latter case gives thermal runaway with the wall temperature, as well as the temperature field within the boundary layer, becoming increasingly larger.

Solution near the singularity
The nature of the singularity as t ! t s has been discussed for related cases in [12,13] for a general value of p. Here we restrict attention to p ¼ 3 in line with the results shown in Figs. 10 and 11. We start by putting s ¼ t s À t and looking for a solution for s small. There is an inner region where we put, following [12,13], with Eqs. (7,8) becoming ð30Þ Equation (30) leads to an eigenfunction expansion, following [13], where / 0 ¼ u 0 g; / 1 ¼ Àu 1 ðg 5 À 20g 3 þ 60gÞ; There is then a middle region where we put f ¼ s À1=4 gðn; sÞ; n ¼ y s À1=4 : Matching with the inner region given by (31, 32) suggests an expansion in powers of s 1=4 , the leadingorder term g 0 being given by There is then an outer region in which f and y are left unscaled and in which there is a regular expansion for f, the leading-order term f 0 being indeterminate in the expansion, as could be expected from Stewartson [23].
On matching with the middle region, This outer region is unaffected by the singularity developing near the wall and depends essentially on how the solution evolves from its initial input.
We can see this structure in Fig. 11a where we plot h profiles taken at times close to t s , the outer boundary condition in this numerical integration was taken at The increasingly larger temperatures developing near the wall as t ! t s can clearly be seen with the temperature field essentially the same at relatively small distances from the wall. From (29), h w $ u 0 ðt s À tÞ À1=2 so that a plot of h À2 w against t should be a straight line, at least close to t s . We can see this in Fig. 11b where we estimate the slope to be approximately À1:702 giving u 0 ' 0:767, a little greater than the theoretical value of 0.707. This difference could be accounted for by the difficulty in maintaining accuracy in integrating numerically very close to a singularity.

Conclusions
We have considered the effect that fluid transfer through the wall, both injection and withdrawal, can have on the development the boundary layer near a forward stagnation point when there is local heat generated at a rate proportional to ðT À T 1 Þ p within the boundary layer. Previously we have seen [12,13] that, without the fluid transfer, how the solution evolved depended on whether 1 6 p\2 or p [ 2, respectively approaching a nontrivial steady state, or either dying away or having a finite-time singularity. A similar situation arises here though the injection/ withdrawal of fluid through the wall can have significant effects. For 1 6 p\2, the solution approaches a nontrivial steady state provided the dimensionless parameter c measuring the wall transfer is such that it is greater than some critical value c c , dependent on p. Otherwise only the trivial state is attained at large times. When there is fluid injection, i.e. c [ 0, large temperatures can be achieved within the boundary layer, increasing as c is increased, see Figs. 1, 2, 7 and 8. These critical values depend on the exponent p increasing from À2 for p ¼ 1 to zero as p ! 2, see Fig. 3. Thus if fluid withdrawal is sufficiently large the effect is to inhibit boundary layer development. However fluid injection can greatly increase both the heat transfer and the fluid flow.
The situation is essentially different when the exponent p [ 2. Here the effect of fluid injection is to limit the range of existence of any nontrivial steady state, see Figs. 4 and 5, whereas it is now fluid withdrawal that increases the temperatures seen within the boundary layer. However, these nontrivial steady states are unstable and the time dependent solution either evolves to the trivial state or reaches a finitetime singularity, see Fig. 10a. Which of these states is achieved depends strongly on the initial temperature input, see Fig. 10b.