On an equation arising in the boundary-layer flow of stretching/shrinking permeable surfaces

In a recent paper, Al-Housseiny and Stone (J Fluid Mech 706:597–606, 2012) considered the dynamics of a stretching surface and how this interacts with the boundary-layer flow it generates. These authors discussed the cases $$c= -3$$c=-3 for an elastic sheet and $$c= -1$$c=-1 for the viscous fluid, c being representative for the stretching velocity of the sheet. The aim of the present paper is to extend the analysis of Al-Housseiny and Stone (2012) to the general values of c, to allow for both a stretching and a shrinking sheet and for the surface to be permeable through the parameter S, where $$S>0$$S>0 for the fluid withdrawal and $$S < 0$$S<0 for fluid injection. Both the cases $$S=0$$S=0 (impermeable surface) and $$S \ne 0$$S≠0 (permeable surface) are considered for both stretching surfaces and shrinking surfaces. In all these cases, asymptotic solutions are presented for large values c and S (both withdrawal and injection).


Introduction
The boundary-layer flow induced by a surface moving in a direction parallel to its length has been much studied since the original contributions by Sakiadis [1], Crane [2] and Banks [3]. This basic problem has been extended in several different ways, for example for a moving cylindrical surface [4], to include the effect of a permeable surface, Magyari and Keller [5] and Fang et al. [6], shearing surfaces, Weidman [7], and cross-flow, Magyari and Weidman [8] and Weidman [9]. The effects of MHD on stretching surfaces has been treated by Ishak et al. [10] and Fang al. [11]. This latter paper also involved the effects of velocity slip on the surface, a topic also addressed by Mukhopadhyay and Andersson [12], Mukhopadhyay [13] and Wang and Ng [14] with convective boundary conditions being taken for the heat transfer by Yao et al. [15]. Much of the work on flows on moving surfaces is reviewed by Wang [16].
In a recent paper, Al-Housseiny and Stone [17] considered the dynamics of the stretching surface itself and how this interacts with the boundary-layer flow it generates, in previous studies, the velocity of the moving surface was taken as given. They considered two specific cases, namely an elastic sheet and where the surface is a highly viscous Newtonian fluid. This resulted in the equations for the boundary-layer flow as where primes denote differentiation with respect to the independent variable y, as defined in [17], and c is representative of the stretching velocity of the sheet. In [17] they discussed the cases c = −3 for an elastic sheet and c = −1 for the viscous fluid. Here, our aim is to extend their analysis for general values of c in an attempt to put their original problem into a much wider context. We allow for both a stretching and a shrinking sheet which we characterize by the parameter λ = ±1, and for the surface to be permeable through the parameter S, S > 0 for fluid withdrawal, S < 0 for fluid injection. This then leads to the boundary conditions where we take λ = 1 for a stretching surface and λ = −1 for a shrinking surface. We note in passing that, when λ = 0, the solution to (1) is simply f ≡ S. We start by considering an impermeable surface, S = 0.

Stretching surface
For this case we take the boundary conditions When c = 0, Eqs. (1) and (3) have the solution f = 1 − e −y and then f (0) = −1. The solution for c = −1 is given in [17] as We plot f (0) against c in Fig. 1 obtained from the numerical solution of Eqs. (1) subject to boundary conditions (3). We see from the figure that f (0) < 0 is decreasing as c is increased with the solution continuing to large positive values of c. To obtain a solution valid for c large and positive we put f = c −1/2 g, η = c 1/2 y to give where primes now denote differentiation with respect to η. Equation (5) suggests expanding in inverse powers of c The leading-order term g 0 satisfies Problem (7) has arisen previously, see [18,19] for example, and has g 0 (0) = −0.62756. We solve the equations arising at O(c −1 ) numerically finding that g 1 (0) = −0.76263. Thus we have Asymptotic expression (8) is shown in Fig. 1 by the broken line showing good agreement with the numerically determined values even at relatively small values of c.
For c > −1 we have 'exponential solutions' in that f ∼ A 0 e −ay as y → ∞, for some a = a(c) > 0 and a constant A 0 . When c = −1 the solution is given by (4) and is the start of the 'algebraic solutions' for c ≤ −1 in which, for c < −1, This means that, in obtaining numerical solutions to Eqs. (1) and (3), we have to increase the value of y ∞ , the value of y where the outer boundary condition is applied, when c ≤ −1. For Fig. 1 we required at least y ∞ = 200.0, whereas we found y ∞ = 15.0 sufficient for c > −1. We terminated our numerical integrations at c −1.5 as it became increasingly more difficult to satisfy the outer boundary condition with sufficient accuracy. The reason for this can be seen in (9) where f approaches zero more slowly as the value of c is decreased. To describe the transition from 'exponential solutions' to 'algebraic solutions' at c = −1 we start by assuming that c > −1 putting c = −1 + δ, where 0 < δ 1. We look for a solution in an inner region by expanding where f 0 is given by (4). At O(δ) we then have A consideration of f as given by expansion (10) and expression (11) indicates that the solution breaks down when y is of O(δ −1 ). This suggests an outer region where we put where primes now denote differentiation with respect to ξ . Equation (13) indicates an expansion in powers of δ, with the leading-order term G 0 satisfying The boundary conditions are, on matching with the solution in the inner region, log ξ + · · · as ξ → 0, G 0 → 0 as ξ → ∞.
Boundary conditions (15) are consistent with this equation and we note that G 0 ∼ e − √ 6 ξ as ξ → ∞ or that f ∼ δ 2 e − √ 6 δ y . The above approach only applies when δ > 0 and recovers the exponential nature of the outer condition for c > −1.
When c < −1 we now put c = −1 − δ, where again 0 < δ 1. This alters the sign of f 1 in expression (11) and, perhaps more significantly, it changes the sign of the second term in Eq. (14) which now becomes The approach of G 0 to its outer condition can no longer be exponential and is now algebraic, of O(ξ −1 ). It is worth commenting that the solutions for c ≤ −1, including the cases given by Al-Housseiny and Stone [17], which have algebraic decay at large distances from the surface are not complete. As pointed out by Brown and Stewartson [20] and later by Merkin [21] these should be regarded as inner solutions and an outer region is also required in which the decay is exponential so as to match smoothly with the ambient conditions. This structure can be seen in the above discussion for (c + 1) small.

Shrinking surface
For this case, we apply the boundary conditions We are unable to get a numerical solution when c ≥ −1. We can show directly that there are no solutions when c = −1 and c = 0. To see that the solutions must terminate at a finite negative value of c we multiply Eq. (1) by f and integrate. On applying boundary conditions (17) we obtain Expression (18) leads to a contradiction when c = −3/2. For c = −2 the solution is In Fig. 2 we plot f (0) against c for this case, starting our numerical integrations with the solution given by (19). We see that the solution continues to large negative values of c and becomes singular as c → −1.5 from below with f (0) becoming large and negative, as could perhaps be expected from (18). For c = −2, f is exponentially small for y large from (19). This follows for the other solutions for c < −3/2 with these also being exponentially small for y large.
We can obtain the solution for large |c|, c < 0, as above, now writing f = −|c| −1/2 g, η = |c| 1/2 y. This results in Eq. (7) at leading order so that We plot asymptotic expression (20) in Fig. 2 by a broken line, showing good agreement with the numerically determined values even to relatively small values of |c|.
To determine how the solution behaves as c → −3/2 we put c = −3/2 − , where 0 < 1, writing We substitute transformation (21) into Eq. (1) and then look for a solution by expanding At leading order we have  1) and (17). Asymptotic expression (20) for |c| large is shown by a broken line primes now denoting differentiation with respect to ζ . To obtain a nontrivial solution to (23) we assume that subject to homogeneous boundary conditions. We construct a solution to Eq. (24) in the form  (18), puts an upper bound on c, requiring that c < −3/2 for solutions for an impermeable shrinking wall.

Permeable surface, S = 0
We can get some insight into this case by considering the solution for c = 0 for which Eqs. (1) and (2) have the solution For a stretching surface, λ = 1 and only a + is applicable giving only one solution with f (0) = −a + , valid for all S, with a + ∼ S as S → ∞ and a + ∼ |S| −1 as S → −∞. For a shrinking wall, λ = −1 and both solutions for a ± in (28) can hold though now require S 2 ≥ 4. When S > 2 (fluid withdrawal), both a + and a − are positive with then two solutions emerging from the critical point at S = 2 and continuing to large S with f (0) ∼ S on one branch and f (0) ∼ S −1 on the other branch. When −2 < S < 2, a ± is complex and there are no real solutions and when S ≤ −2 both roots are negative.

Stretching surface
We now have In Fig. 3, we plot f (0) against S for c = −0.5, 0.0, 0.5, 1.0, 2.0. In each case, the solution continues to both large-positive and large-negative values of S with the curves becoming linear for large S. The curves for these values of c look as though they pass through the same point, at S −0.29, however, this is not really the case and is just a feature of the plot. We plot f (0) against c in Fig. 4 for S = −0.5, 0.5, 1.0, 2.5. All curves pass through f (0) = −0.8165 at c = −1 and continue to large positive values of c. For S > 0, | f (0)| increases as c is increased, linearly for larger c. However, for S < 0, the values of f (0) decrease slowly to zero as c is increased. For c < −1 there is a change from 'exponential' to 'algebraic' behaviour for y large and, as noted above, it becomes increasingly difficult to obtain a solution that satisfies the outer boundary condition with sufficient accuracy. To obtain a solution valid for S large we put Eqs. (1) and (29) become subject to where primes now denote differentiation with respect to ζ . Equation (31) leads to an expansion in powers of S −2 , the leading-order term F 0 being given by Expression (33) holds only when Expression (34) shows the linear behaviour with S seen in Fig. 3. We note that, for a shrinking wall F (0) = −1, giving a change in sign for F 0 and then f (0) ∼ (1 + c) S + · · · for large S, again for (c + 1) > 0.
When c = −1, we have the general solution For c < −1 we now put At leading order we obtain Equation (37) has the solution Then from (36) The different nature of the solution for S large between c > −1 and c < −1 can be seen in Fig. 5 where we plot f (0) against S for c = −0.8 and c = −1.1, the constant value of f (0) −0.8165 for c = −1 is shown by a broken line. In the former case, | f (0)| increases linearly with S whereas in the latter case f (0) decreases slowly as S is increased. The transition from 'exponential solutions' to 'algebraic solutions' at c = −1 can be seen in expressions (33) and (38) follows, with only a very slight modification, from that described previously for an impermeable surface.
Strong fluid injection, S < 0, |S| 0 To obtain a solution for |S| large we put with Eqs. (1) and (29) becoming the outer boundary condition is relaxed at this stage, and where primes now denote differentiation with respect to Y . Equation (41) suggests looking for a solution by expanding in powers of |S| −2 , with leading-order problem being When c = 0, Eq. (42) has the solution H 0 = e −Y which also satisfies the required outer condition and gives f (0) ∼ −|S| −1 + · · · as S → −∞. Otherwise Eq. (42) has the solution We recover the previous expression for H 0 in the limit as c → 0. For −1 < c < 0, c/(c + 1) < 0 and expression (43) satisfies the outer condition that H 0 → 0 as Y → ∞ through algebraically small terms. For c outside this range, expression (43) is zero at Y = (1 + c)/c and an outer region is required in this case to satisfy the outer boundary condition. We put and then This leaves Eq. (1) essentially unaltered, though primes now denote differentiation with respect to η, subject to, on matching with the inner region, From (40) and (43) we have, for c = −1, consistent with the results shown in Fig. 3.

Solution for c large
When S is of O(1), the previous behaviour for c large no longer applies and we now make the transformation Equations (1) and (29) become We expand P(z; c) = P 0 (z) + c −1 P 1 (z) + · · · (50) Substituting expansion (50) into (49) gives from which it follows that The above asymptotic solution holds only when S > 0 and gives the linear behaviour for c large seen in Fig. 4. Also seen in Fig. 4, the asymptotic behaviour for c large is different when S < 0. In this case, a region next to the wall develops in which f = −|S| + y. There is then a shear layer in which we put y = |S| + y and then f = c −1/2 g, η = c 1/2 y. This results in Eq. (7) at leading order now subject to g ∼ η as η → −∞, g → 0 as η → ∞.
From Eqs. (7) and (53) The asymptotic form given by (48)-(52) for c large is different to that seen above when S = 0. To continue further, we now assume that S is of O(c −1/2 ), putting S = c −1/2 T and proceed as for the case when S = 0. Now, the leading-order term in expansion (6) is given by Eq. (7) subject to the boundary conditions in (7) though now with g 0 (0) = T . In Fig. 6, we plot g 0 (0) against T where we see that the solution continues to large-positive T . To obtain a solution for T large we write g 0 = T + T −1 g 0 , η = T η giving, at leading order, g 0 + g 0 = 0, g 0 (0) = 0, g 0 (0) = 1, g 0 → 0 as η → ∞.
(56) From the above f (0) ∼ −c 1/2 T + · · · = −S c + · · · for c large, (57) agreeing with the leading-order term in expansion (52). Asymptotic expression (57) is shown in Fig. 6 by a broken line and agrees well with the numerically determined values even at relatively small values of T .
The solution continues to large-negative values of T , with g 0 (0) → 0 and g 0 (0) < 0. We can express the solution to Eq. (7) as Thus, the sign of g 0 (η) on 0 < η < ∞ is given by the sign of g 0 (0) and hence a consequence of the boundary condition that g 0 (0) = 1 is that we must have g 0 (0) < 0. As T → −∞ an inner region develops in which g 0 = −|T | + η at leading order. This leads to an outer region where we put η = |T | +η leaving g 0 unscaled which still satisfies Eq. (7) but now subject to g 0 ∼η asη → −∞, g 0 → 0 asη → ∞.
(59) This problem has been treated previously by Chapman [22] from which it follows that g 0 ∼ d 0 + · · · asη → ∞, where d 0 0.8758 and g 0 ∼ B 0 e −η 2 /2 asη → −∞ for some constant B 0 . From this it follows that g 0 (0) ∼ B 0 e −T 2 /2 + · · · as T → −∞ reflecting the rapid fall seen in Fig. 6 in which g 0 (0) deceases as |T | increases. Our numerical integrations suggest that B 0 −0.233. The values of f (0) on the lower branch decrease to zero as S is increased when c ≥ 0 along the lines indicated above for c = 0. However, when c < 0, f (0) becomes zero at a finite value S 0 of S and then decreases rapidly as S is increased further. This is illustrated in Fig. 8 where we plot f (0) against S for c = −0.2, −0.5, −0.8 when f (0) < 0 and in Fig. 9 we give S 0 plotted against c. The values of S 0 increase rapidly as c → −1 as could be expected from (35) and not quite so obvious from the figure as c → 0 suggesting that the behaviour of f (0) with S seen in Fig. 8 is restricted to −1 < c < 0.
We have already noted above how the solution on the upper branch behaves for S large. To obtain how the solution on the lower branch behaves for S large we start by noting that, from expression (4), a − ∼ −S −1 for c = 0 when S is large. This leads us to follow the previous solution for strong injection on a stretching wall by putting At leading order we recover Eq. (42) now subject to H 0 (0) = −1 with solution For c > 0 and an outer region similar to that described above is required to complete the solution. From this we have that It is, for c > 0, through this outer region that the outer boundary condition is approached through exponentially small terms. The asymptotic nature of the solution for S large on the lower branch is different for c < 0, as can be seen in Fig. 8. To obtain an asymptotic solution in this case we now put Boundary conditions (65) indicate looking for a solution by expanding in S −2 , the leading-order term h 0 satisfying To obtain a solution to subject to (66) we assume that h 0 (0) = −d 0 for some constant d 0 > 0. We now put 0 t to obtain the leading-order problem as It is the numerical solution to (67) that determines the value of d 0 for a given value of c. A plot of d 0 against c is shown in Fig. 10 where we see that the solution becomes singular as c → −1 and as c → 0. Thus, limiting this asymptotic solution to the range −1 < c < 0.

by a broken line and is virtually indistinguishable from the numerically determined values.
There is also a solution following from (19) for c = −2, Expression (69) gives f (0) = −S and we can use this value for the numerical integrations, the results of which are shown in Fig. 11with plots of f (0) against c for S = 0.5, 1.0, 2.0. These results follow closely those seen in Fig. 2 with f (0) becoming singular, with f (0) → −∞, as c → −3/2. For S = 1 and 2, f (0) decreases monotonically to zero as |c| increases, finding the same for other values of S > 1 tried. However, for S = 0.5, we find that f (0) changes sign at c −3.52, is positive for −17.27 S −3.52 before becoming negative again and then approaching zero with f (0) < 0 as |c| increases further. When |c| large and S > 0, there is an inner region in which y is left unscaled and in which we look for a solution by expanding At leading order we have At O(|c| −1 ) we then have (S − y) f 1 = −1 giving From (71) and (72) showing that f (0) → 0 from below as c → −∞. There is also an outer region where we put y = S + y and then f = −|c| −1/2 U and Z = |c| 1/2 y. When this is substituted into Eq. (1) we find that, at leading order, The solution to (74) has already been given by Chapman [22] from which it follows that U → 0.8758 as Z → ∞ so that f (∞) ∼ 0.8758 |c| −1/2 as |c| → ∞.

Critical values
The critical values of S seen in Fig. 7 lead us to consider these in more detail and we give a plot of S c against c in Fig. 12, these values being obtained numerically following the approach discussed in [23,24] for example. The values of S c increase, becoming large, as c approaches −1, as might be expected from (35) which does not give a critical value, and decrease slowly towards zero, with f c (0), the value of f (0) at S c , increasing as the values of c increase. We can apply the above discussion to see how the values of S c vary as c is increased. We now put This results in Eq. (7) at leading order now subject to These equations have a critical value at T c 1.7027 with g c (0) = 1.1751. From this it follows that Expression (77) shows that S c decreases slowly to zero as c is increased, as can be seen in Fig. 12. Asymptotic expression (77) for S c is also plotted in Fig. 12   the dynamics of the moving surface, in an attempt to put the original problem treated in [17] into a much wider framework. A shrinking wall is also considered leading to different flow dynamics to those seen [17] and here for general c for a stretching wall. We have further included the effect of fluid transfer through the wall, characterized by the dimensionless parameter S, S > 0 for fluid withdrawal, S < 0 for fluid injection. We started by considering an impermeable surface for both stretching and shrinking surfaces. For a stretching surface, we found that the solution continued to large values of c, Fig. 1, with an asymptotic solution being derived for large c, expression (8). For negative values of c, there was a change from exponentially decaying solutions for c > −1 to algebraically decaying solutions for c ≤ −1, limiting the range of c where we could get reliable numerical solutions. How the solution behaved as it passed through c = −1 was examined in detail. For a shrinking surface solutions are limited to the range c < −3/2, Fig. 2, though now these solutions are exponential at large distances. The nature of the solution as c → −3/2 from below, expression (26), and as c → −∞ were determined, expression (20). As a consequence, there is no solution for the case c = −1 in [17] for a shrinking surface.
For a stretching permeable surface we found that the solution continued to both large-positive S, strong fluid withdrawal, and large-negative S, strong fluid withdrawal, Fig. 3. When we considered how the solution behaved for large-positive S, we found that this depended on whether c > −1 or c < −1, as noted in Fig. 5. The solution continued to large-positive values of c and again the change from exponential decay solutions to algebraically decaying solutions at c = −1 limited the range of negative values of c where we could obtain reliable numerical solutions. We derived an asymptotic solution for large c, the nature of this being dependent on whether S > 0, expression (52), or S < 0, expression (54).
For a shrinking permeable surface, we found that there was a critical value S c of S, dependent on c, limiting solutions to S ≥ S c with dual solutions for S > S c , Fig. 7. The behaviour of the solution on the upper branch followed that for a stretching surface. However, the nature of the solution on the lower branch for S large depended on whether c ≥ 0 or −1 < c < 0, Fig. 8. There are also solutions for c < −3/2, Fig. 11, which become singular as c → −3/2 from below and continue to large-negative values of c. For a shrinking surface f (0), and consequently the skin friction, can become negative, Figs. 8 and 11, leading to the possibility of the boundary layer separating from its bounding surface. We also discussed the critical values S c , Fig. 12, which exist only for c > −1, an asymptotic form for S c for c large being derived, expression (77).