SIAM Journal on Applied Mathematics, 57(1997), 1-14.

IMPULSIVE STRETCHING OF A SURFACE IN A VISCOUS FLUID

C.Y. WANG^{∗}, Q. DU^{†}, M. MIKLAVČIČ^{‡}, AND C.C. CHANG^{§}

Abstract. The induced transient viscous flow due to a suddenly stretched surface is studied. After a similarity transform, the unsteady Navier-Stokes equation is solved by several methods, including perturbation for small times, numerical integration and asymptotic analysis for large times. It is found that the validity of the small-time series can be greatly extended and the approach to steady state is exponential.

Key words. Navier-Stokes equations, perturbations, asymptotics, stability

AMS subject classifications. 76D05, 41A60, 41A58,76E99

1. Introduction. The study of the flow resulting from a stretching boundary is important in processes such as the extrusion of sheet material into coolant fluids. The tangential velocity imparted by the sheet induces motion to the surrounding fluid which alters the convection cooling of the sheet. Due to the much higher viscosity of the extrusate, one can assume the fluid is affected by the sheet but not vice versa. Thus the fluid dynamic problem can be idealized to the case of a fluid disturbed by a tangentially moving boundary. Experiments show that the velocity of the boundary is approximately proportional to the distance from the orifice [11].

The steady state solution to this problem belongs to an important class of exact solutions of the Navier-Stokes equations [12]. Since the tangential velocities are proportional to the distance from the orifice, the Navier-Stokes equations reduce to an ordinary differential equation through a similarity transformation. The two dimensional stretching of a sheet is particularly interesting since a closed form solution has been obtained by Crane [3].

The aim of the present paper is to study the impulsive starting from rest of a stretching sheet in a viscous fluid and its approach towards equilibrium. As time approaches infinity the transient solution should approach Crane’s steady state solution. Although strictly impulsive motion does not occur in practice (infinite stresses at time zero), it is nevertheless an excellent approximation of sudden changes in the boundary conditions. The study of bodies impulsively started in a viscous fluid is fundamental to the understanding of unsteady viscous flows [9, 10]. For example, the impulsive starting of a translating circular cylinder is now a classical problem whereby experimental, analytical and numerical results can be compared. In general, at the moment an impulsive start is initiated, the flow field is a potential. Then vorticity begins to diffuse from the boundary. A steady state is attained when convection of vorticity balances diffusion of vorticity. Due to nonlinear convection, there are no analytic solutions for impulsive start, except for parallel or concentric flows where the nonlinear terms are identically zero. Most solutions are found by a variety of numerical schemes.

For a very short time interval after impulsive start, nonlinear effects are secondary to diffusion and a perturbation method may be used to obtain an analytic solution. Following [2, 4], numerous authors used small-time expansion to solve impulsive starting problems. The solution, as the method suggests, is valid only for small times. As time t approaches infinity, the solution, which is a power series in t, becomes infinite.

∗ Department of Mathematics, Michigan State University, E. Lansing, MI 48824 USA(cywang@math.msu.edu) † Department of Mathematics, Michigan State University, E. Lansing, MI 48824 USA(du@math.msu.edu). The

research was supported in part by DOE under Grant DE-FG02-93ER25172.

‡ Department of Mathematics, Michigan State University, E. Lansing, MI 48824 USA(milan@math.msu.edu) § Institute of Applied Mechanics, Taiwan University, Taipei, Taiwan, ROC.

An interesting work by Benton[1] studied the impulsive rotation from rest of a disc in an infinite viscous fluid. Assuming the time scale is small, Benton expanded the unknowns in power series in time. The unsteady Navier-Stokes equations then become decoupled for each successive order of the time variable. Benton solved for the first two orders analytically, then numerically computed the next two orders. Although such a method was used to treat impulsive start problems before, what is significant is the way Benton presented the small-time results. He plotted, on the same graph, the solutions of various different orders in time. The range of validity of the solutions up to order N was determined by observing the location where the solution of order N + 1 bifurcated. This is equivalent to the quantitative error estimation of an asymptotic series. Benton showed that the solution up to the fourth order was valid for “small” normalized times as large as 1.5, at which time the stresses and velocities were very close to the steady state values of the classical Von Karman rotating disc problem. Thus in certain transient problems the solution for small times may infer the steady state.

The present paper studies the flow due to an impulsively stretched boundary. In Section 2 the mathematical problem is formulated. In Section 3 the solution in terms of small-time expansion is obtained to order 6. Its range of validity will be ascertained by both Benton’s method and direct numerical integration (Section 4) of the partial differential equations. Analytical behavior of the solution for large time , its approach to the steady state Crane’s solution, and the stability of Crane’s solution are presented in Sections 5 and 6. Thus the complete transient history from impulsive start to steady state is analysed.

2. Formulation. Figure 1 shows the two-dimensional impulsive stretching of an extensible surface in a viscous fluid. Let the y axis denote the symmetry plane. Velocity is zero for time t < 0. For t > 0 the surface at y = 0 is suddenly stretched such that the local tangential velocity is bx where b is a constant of dimension [T^{−1}]. The Navier-Stokes equations are

q_{t} + (q · ∇)q = −\frac{1}{ρ}∇p + ν∇^{2}q, ∇ · q = 0,

where q is the velocity, p is the pressure, ρ is the density and ν is the kinematic viscosity. In two dimensions, the stream function ψ is related to the velocity by

q = k × (∇ψ).

After taking the curl, the Navier-Stokes equations become

(∇^{2}ψ)_{t} + \frac{∂(∇^{2}ψ, ψ)}{∂(x, y)} = ν∇^{4}ψ. (1)

The initial and boundary conditions are

(2)  t < 0,  ψ ≡ 0
(3)  t > 0,  y = 0, ψ_{y} = bx, ψ = 0

y → ∞, ψ_{y} → 0, ψ bounded. (4)

Note that (4) indicates zero lateral velocity (and pressure gradient) at infinity. This is quite different from the potential towards a plate where there exists, even far from the plate, a lateral velocity caused by a lateral pressure gradient. Since a natural length scale is absent, we introduce the similarity transformation

ψ = \sqrt{bν}xf(η, ξ), (5)

Impulsive stretching of a surface 3

η = \sqrt{b/ν}y, ξ = bt. (6) The x dependency is due to the boundary condition in (3). If there were a natural length scale, such a similarity transformation would not be applicable. The governing equation, after integrating once, becomes f_{ηηη} − f_{ηξ} − (f_{η})^{2} + ff_{ηη} = 0. (7) The initial and boundary conditions are f(ξ < 0, η) = 0 (8) f(ξ > 0, 0) = 0 (9) f_{η}(ξ > 0, ∞) = 0 (10) f_{η}(ξ > 0, 0) = 1 (11) As ξ → ∞, we expect the solution to approach Crane’s steady state solution [3] f(ξ → ∞, η) = 1 − e^{−η}. (12)

Fig. 1. The physical problem

3. Perturbation solution for small times. For small times, after impulsive start, we set ξ = ε^{2}τ, η = εζ (13) where τ, ζ are of order unity and ε << 1. This form is used because the unsteady Navier-Stokes equation includes ζ/√τ = η/√ξ as an invariant group. We perturb from the quiescent state in a power series in time

f = ∑_{n=1}^{∞} ξ^{n− \frac{1}{2}} G_{n}(z) (14)

4 WANG, DU, MIKLAVČIČ AND CHANG

where the power n − 1/2 is obtained by applying the method of stretching used in (7) and

z ≡ \frac{ζ}{2\sqrt{τ}} (15)

is the Stokes’ variable. (7) reduces to the successive ordinary differential equations G′′′_{1} + 2zG′′_{1} = 0 (16) G′′′_{2} + 2zG′′_{2} − 4G′_{2} = 2[(G′_{1})^{2} − G_{1}G′′_{1}] (17) G′′′_{3} + 2zG′′_{3} − 8G′_{3} = 2[2G′_{1}G′_{2} − G_{1}G′′_{2} − G_{2}G′′_{1}] (18) or in general

G′′′_{n} + 2zG′′_{n} − 4(n − 1)G′_{n} = 2 ∑_{i=1}^{n−1}[G′_{i}G′_{n−i} − G_{i}G′′_{n−i}]. (19)

The boundary conditions are G_{1}(0) = 0, G′_{1}(0) = 2, G′_{1}(∞) = 0 (20) G_{n}(0) = G′_{n}(0) = G′_{n}(∞) = 0, n ≥ 2. (21) The solutions are obtained analytically whenever possible. The first order solution is equivalent to Stokes’ first problem of the impulsively translated plate: ( Note that the impulsively translated plate solution has been erroneously credited to Rayleigh, but Stokes published the same solution 60 years earlier.)

G_{1} = 2[zE + \frac{1}{\sqrt{π}} (1 − E)], (22)

where

E ≡ erfc(z) ≡ \frac{2}{\sqrt{π}} ∫_{z}^{∞} e^{−s^{2}}ds, E ≡ exp(−z^{2}) (23)

and

G′′_{1}(0) = −\frac{4}{\sqrt{π}} . (24)

(17) then becomes

G′′′_{2} + 2zG′′_{2} − 4G′_{2} = 8[E^{2} + \frac{2}{\sqrt{π}} zEE + \frac{2}{π} (E − E^{2})]. (25)

After some work we obtain that the solution satisfying the boundary conditions is

G_{2} = 2[\frac{1}{3}z^{3}E^{2} − \frac{2}{3\sqrt{π}} z^{2}EE + \frac{11}{6\sqrt{π}} EE + \frac{1}{3π} zE^{2} − \frac{1}{2}zE^{2} + \frac{2}{3\sqrt{π}} E

− \frac{4\sqrt{2}}{3\sqrt{π}} erfc(\sqrt{2}z) + (\frac{1}{2} − \frac{2}{3π})(zE − \frac{2}{3\sqrt{π}} E + \frac{2}{3}z^{3}E − \frac{2}{3\sqrt{π}} z^{3}E) (26)

− \frac{13}{6\sqrt{π}} + \frac{4\sqrt{2}}{2\sqrt{π}} − \frac{4}{9π\sqrt{π}}].

Impulsive stretching of a surface 5

Thus,

G′′_{2}(0) = \frac{2(8 − 9π)}{3π\sqrt{π}} = −2.427340. (27)

The solution of G_{3} is more complicated. We substitute (22) and (26) into (18). Partial integration yields

G′_{3} = 2\{Φ(z)+ \frac{2}{15π^{2}} (16 + 3π−8\sqrt{2}π)E + \frac{2}{3π\sqrt{π}} zE^{3}+ \frac{8}{9π\sqrt{π}} zE − \frac{14}{3π} E^{2}E

+ \frac{4}{9π\sqrt{π}} (−4 + 3π)zE− \frac{5}{3π} z^{2}E^{2}E+ \frac{4}{\sqrt{π}} zEE^{2}+ \frac{4}{3\sqrt{π}} z^{3}EE^{2}+ \frac{1}{4}E^{3}−z^{2}E^{3}

− \frac{1}{3}z^{4}E^{3}−(−\frac{1}{8} + \frac{16}{45π^{2}} − \frac{31}{10π})(−\frac{10}{3\sqrt{π}} zE − \frac{4}{3\sqrt{π}} z^{3}EE+E+4z^{2}E+ \frac{4}{3}z^{4}E) (28)

− \frac{4}{3π\sqrt{π}} zE^{2}− \frac{2}{3π} EE+ \frac{4}{3π} z^{2}EE+ \frac{4}{3\sqrt{π}} zE^{2}+(1 − \frac{4}{3π})(\frac{4}{3π} E^{2}+ \frac{1}{3π} z^{2}E^{2}

+ \frac{23}{12\sqrt{π}} zEE− \frac{1}{2\sqrt{π}} z^{3}EE− \frac{3}{8}E^{2}+ \frac{1}{2}z^{2}E^{2}+ \frac{1}{6}z^{4}E^{2})+ \frac{16\sqrt{2}}{15π}E erfc(\sqrt{2}z)\},

where Φ(z) is governed by

Φ′′ + 2zΦ′ − 8Φ = \frac{72}{5π\sqrt{π}} zE^{3}, Φ(0) = Φ(∞) = 0. (29)

Unfortunately no closed form solution of Φ(z) can be found. The two homogeneous solutions are

Φ_{1} = z^{4} + 3z^{2} + \frac{3}{4} (30)

Φ_{2} = −\frac{10}{\sqrt{π}} zE − \frac{4}{\sqrt{π}} z^{3}E + 3E + 12z^{2}E + 4z^{4}E. (31)

By variation of parameters the solution in integral form is

Φ(z) = \frac{6}{5π} (∫_{∞}^{z} zE^{2}Φ_{2}dz)Φ_{1}(z) + \frac{6}{5π} [\frac{1}{4} ∫_{0}^{∞} zE^{2}Φ_{2}dz − ∫_{0}^{z} zE^{2}Φ_{1}dz]Φ_{2}(z) . (32)

From (32) an exact initial value is obtained:

Φ′(0) = \frac{−24}{5π\sqrt{π}} (\frac{11}{4} − \frac{9}{2\sqrt{3}}) = −0.130961 . (33)

Then

G′′_{3}(0) = \frac{1}{\sqrt{π}} [\frac{512}{135π^{2}} + \frac{(216\sqrt{3} − 730)}{15π} + \frac{49}{6}] = 0.563620 . (34)

The function G_{3}(z) is obtained numerically from (18) and the exact initial values G_{3}(0) = G′_{3}(0) = 0 and (34). Since G_{3}(z) has no explicit solution, all higher order solutions are obtained numerically. Using G_{4}(0) = G′_{4}(0) = 0, we guess G′′_{4}(0) and integrate (19) as an initial value problem by a variable step Runge-Kutta algorithm. At a large value of z (z ∼ 5) we check if |G′_{4}(z)| < 10^{−6}. One dimensional shooting yields

G′′_{4}(0) = −0.192699. (35)

A similar method yields

G′′_{5}(0) = 0.0696358, (36)

G′′_{6}(0) = −0.0253397. (37)

The process may be continued, limited only by computer storage. The functions G′_{2}(z), G′_{3}(z), G′_{4}(z), G′_{5}(z) are shown in Fig 2. It is shown that the functions alternate in sign and decay within a boundary layer of z ∼ 3. All functions approach constants at infinity.

Fig. 2. The similarity function G′_{n}(z)

We also have

\begin{matrix}G_{3}(∞) = 0.239028, & G_{4}(∞) = −0.117979, \\ G_{5}(∞) = 0.057411, & G_{6}(∞) = −0.029042. \end{matrix} (38)

Fig 3 shows a plot of the natural logarithm of |G′′_{n}(0)| and |G_{n}(∞)|. We see these values approach straight lines as n increases. This implies that the amplitudes of G_{n}(z) decay exponentially as n → ∞. The local shear stress is proportional to f_{ηη}(0, ξ). Using (14) we find

f_{ηη}(0, ξ) = \frac{1}{4} ∑_{n=1}^{∞} ξ^{n− \frac{3}{2}} G′′_{n}(0). (39)

4. Numerical solution of the original equation. In order to ascertain the accuracy of the results for small time expansion, (7) through (11) are solved by a fully discrete finite difference method.

Fig. 3. The decay of G′_{n}(z) with n

Since ϕ(ξ, η) = f_{η}(ξ, η), (7) becomes

ϕ_{ξ} − ϕ_{ηη} + ϕ^{2} − ϕ_{η} ∫_{0}^{η} ϕ = 0, for η ≥ 0, ξ ≥ 0 (40)

and its steady state solution is e^{−η}.

If g(ξ, η) = ϕ(ξ, η) − e^{−η} then (40) becomes

g_{ξ} − g_{ηη} − g_{η} + e^{−η} (g_{η} + 2g + ∫_{0}^{η} g) + g^{2} − g_{η} ∫_{0}^{η} g = 0 (41)

with the boundary conditions for ξ > 0:

g(ξ, 0) = 0, g(ξ, ∞) = 0. (42)

We view ξ as the time variable and η as the spatial variable. The numerical scheme involves solving a parabolic type nonlinear equation for g. The original equation is posed on a semiinfinite spatial interval. Numerical computations are made, however, on a truncated interval of finite length L. Standard difference approximations on a uniform grid are used on the spatial derivatives while an explicit Euler method is used in the time integration. For example, let {η_{j}} be the grid points, η_{0} = 0, η_{N} = L and η_{j+1} − η_{j} = h for j = 0, 1, 2, ..., N − 1. At time ξ_{l}, we use

g_{ηη}(ξ_{l}, η_{j}) ≈ \frac{g_{j+1}^{l} − 2g_{j}^{l} + g_{j−1}^{l}}{h^{2}} , ∀j

for the second order derivative terms and

g_{ξ}(ξ_{l}, η_{j}) ≈ \frac{g_{j}^{l+1} − g_{j}^{l}}{Δt} , ∀ l ≥ 0 , 1 ≤ j ≤ N − 1

for the time derivative. The composite trapezoidal rule is used for the integral

∫_{0}^{η_{j}} gdη ≈ \frac{h}{2}(g_{0} + g_{j}) + h ∑_{k=1}^{j−1} g_{k} .

At each time step, f is updated by summing together the above numerical integration of g in space and the antiderivative of e^{−η}. The three parameters involved in the computations are the truncated interval length L, the spatial mesh size h, and the time-step size Δt. Various values of these parameters were used in our computations to check the accuracy and convergence of the numerical solutions. A stability condition like Δt = ch^{2} was enforced. The value of the constant c was chosen to guarantee that a discrete maximum principle holds for the nonlinear difference equations which ensure the stability of time-integration. This principle can be proved rigorously and has been verified in all the computations we have performed. The initial condition is given by

g_{j}^{0} = −e^{−η_{j}} , 0 ≤ j ≤ N;

and the boundary conditions are imposed as:

g_{0}^{l} = g_{N}^{l} = 0, l > 0 .

Finally, the quantity f_{ηη}(ξ, 0) was computed numerically by taking the difference of g near η = 0. We have also compared the first order difference approximation with higher order differences based on extrapolation. The discrepancy in the values was very insignificant after the first few time-steps. In fact, the error was less than 1% after the initial transient layer. The numerical results are plotted later in Figure 6.

5. Asymptotic behavior for large times. We shall now study the behavior of solutions of (7) that are near Crane’s steady state solution (12). Our physical intuition suggests that Crane’s solution is stable. Thus, we expect that the solution g of (41) decays to 0, as ξ → ∞, for any small initial perturbation. Our numerical experiments confirm this expectation. Moreover, they suggest an exponential decay of g.

In our numerical computations we used the initial condition g(0, η) = −e^{−η} which corresponds to the impulsive stretching of the surface. The domain of the variable η was truncated to [0, ℓ] and the solution was calculated with various grid sizes. In each case the ratio

e(ξ, ℓ) = −\frac{log(|g(ξ + Δξ, η)|/|g(ξ, η)|)}{Δξ} (43)

is calculated as ξ → ∞ for some fixed η, Δξ > 0. Plots of this ratio versus ξ are given in Figure 4 for ℓ = 30 and ℓ = 60. As ξ → ∞ these ratios approach a limit which is asymptotic to the value 0.25 as ℓ increases. In Figure 5 we have several plots of e(ξ, ℓ) − 0.25 versus ξ for ℓ = 30, 45, 60 and 90. These observations suggest that

lim_{ξ→∞} e(ξ, ℓ) ≈ 0.25 + 10ℓ^{−2} (44)

and hence we expect the decay of g(ξ, η) to be like e^{−ξ/4}.

However, just by looking at (41) it is not at all obvious that Crane’s solution is stable. The first three terms of (41) actually suggest that g may decay algebraicly at most when the initial conditions are in L^{p}(0, ∞). This is because the continuous spectrum of the linear operator contains 0. Since the continuous spectrum can be moved by using weighted spaces in place of L^{p}, we find, after some experimentation, that the best weight is e^{−η/2}. This is equivalent to considering perturbations in the form

g(ξ, η) = e^{−η/2}v(ξ, η) (45)

Fig. 4. The ratio e(ξ, ℓ) from (42) as function of ξ. Left l = 30, Right l = 60.

Fig. 5. The value of e(ξ, ℓ) − 0.25 as function of ξ for various values of constant l. Top row: l = 30, l = 45. Bottom row:l = 60, l = 90.

which implies that v has to satisfy the boundary conditions v(ξ, 0) = 0, v(ξ, ∞) = 0 and the equation v_{ξ} + Av = h(v) (46) where Av = −v_{ηη} + e^{−η}v_{η} + ( \frac{1}{4} + \frac{3}{2}e^{−η} )v + e^{−η/2} ∫_{0}^{η} e^{−s/2}v(ξ, s)ds (47)

h(v) = −e^{−η/2}v^{2} + (v_{η} − \frac{1}{2}v) ∫_{0}^{η} e^{−s/2}v(ξ, s)ds. (48)

We can show that the lower bound on the spectrum of the linear operator A is 1/4 which would suggest decay of v like e^{−ξ/4}. However, using the energy method, we shall now demonstrate this decay for all ‘nice’ solutions of the nonlinear equation (46) provided only that v is not too large initially. In the next section this decay will be demonstrated without requiring that the solution be ‘nice’.

Let

W (ξ) = (∫_{0}^{∞} v(ξ, η)^{2}dη)^{1/2} , W_{1}(ξ) = (∫_{0}^{∞} v_{η}(ξ, η)^{2}dη)^{1/2} .

Multiplying (46) by 2v and integrating over η gives

2W W′ + 2 ∫_{0}^{∞} vAv dη = 2 ∫_{0}^{∞} vh(v) dη. (49)

Using (47) and integration by parts we obtain

2 ∫_{0}^{∞} vAv dη = (∫_{0}^{∞} e^{−s/2}v(ξ, s)ds)^{2} + ∫_{0}^{∞} (2v_{η}^{2} + (1/2 + 4e^{−η})v^{2}) dη

hence,

2 ∫_{0}^{∞} vAv dη ≥ 2W_{1}^{2} + W^{2}/2. (50)

Integration by parts gives

2 ∫_{0}^{∞} vv_{η} (∫_{0}^{η} e^{−s/2}v(ξ, s)ds) dη = − ∫_{0}^{∞} v^{3}e^{−η/2}dη

which implies that

2 ∫_{0}^{∞} vh(v) dη = ∫_{0}^{∞} v(6v_{η} − v) (∫_{0}^{η} e^{−s/2}v(ξ, s)ds) dη

and since ∫_{0}^{η} e^{−s/2}|v(ξ, s)|ds ≤ W (ξ), we have

|2 ∫_{0}^{∞} vh(v) dη| ≤ 6W^{2}W_{1} + W^{3}. (51)

If (51) and (50) are used in (49), we obtain

2W W′ ≤ −W^{2}/2 − 2W_{1}^{2} + 6W_{1}W^{2} + W^{3}
= −W^{2}/2 − 2(W_{1} − 3W^{2}/2)^{2} + 9W^{4}/2 + W^{3}
≤ −W^{2}/2 + 9W^{4}/2 + W^{3}
4W′ ≤ −W + 9W^{3} + 2W^{2},

which implies that if W (0) < 0.24 then W (ξ) → 0 as ξ → ∞. Moreover the decay is like e^{−ξ/4} for large time ξ. By analyzing the spectrum of A one can show that we have obtained a sharp estimate on the asymptotic decay rate. This asymptotic decay rate does not depend on nonlinearity h. However, the maximum magnitude of the initial perturbation for which the solutions eventually decay does very much depend on h. Above we have used only one crude bound on h, see (51), so we should alsoexpect decay of W for much larger W (0). In the case of impulsive start, W (0) = 1 and the numerical calculations suggest decay.

To explain ℓ^{−2} term in (44) consider the operator T v = −v_{ηη} + 0.25v with zero boundary conditions on [0, ℓ]. The essential features of A are captured by T. Note that the spectrum of T starts at 0.25 + π^{2}ℓ^{−2} which determines the decay rate of solutions of v_{ξ} + T v = 0.

6. Proof of the decay of perturbations about steady state. Equation (46) will be considered as an abstract evolution equation in the real Hilbert space H = L^{2}(0, ∞), hence, for each ξ ≥ 0 we will consider v(ξ) to be a real valued function of η belonging to H, i.e. v(ξ) ∈ H for ξ ≥ 0. We will need the precise definition of D(A), the domain of A, as well as fractional powers of A. An easy way to do this is through the variational formulation of A.

Let (·, ·) and ∥ · ∥ denote the usual inner product and the norm on H. Define V = W_{0}^{1}(0, ∞) and

(u, v)_{V} = (u′, v′) + (u, v), ∥u∥_{V} = (u, u)_{V}^{1/2} for u, v ∈ V.

Note that v belongs to V if and only if it satisfies the following conditions: v ∈ H, v is absolutely continuous on [0, a] for every a > 0, v(0) = 0 and v′ ∈ H.

Define a bounded sectorial form F : V × V → IR by

F (v, u) = ∫_{0}^{∞} (v′u′ + v′ue^{−η} + vu( \frac{1}{4} + \frac{3}{2}e^{−η} ) + ue^{−η/2} ∫_{0}^{η} e^{−s/2}v(s)ds) dη.

The linear operator A associated with F is defined as follows: v ∈ D(A) if and only if v ∈ V and there exists g ∈ H such that F (v, u) = (g, u) for all u ∈ V, and in this case Av = g. One can show that A is actually given by (47) and that D(A) = W^{2}(0, ∞) ∩ W_{0}^{1}(0, ∞). Observe that if v ∈ V and g(η) = ∫_{0}^{η} e^{−s/2}v(s)ds, then

F (v, v) = ∫_{0}^{∞} (v^{′2} + v′ve^{−η} + v^{2}( \frac{1}{4} + \frac{3}{2}e^{−η} ) + gg′) dη
= \frac{1}{2}g(∞)^{2} + ∫_{0}^{∞} (v^{′2} + v^{2}( \frac{1}{4} + \frac{3}{2}e^{−η} ) + v′ve^{−η}) dη
= \frac{1}{2}g(∞)^{2} + ∫_{0}^{∞} (v^{′2} + v^{2}( \frac{1}{4} + 2e^{−η} )) dη
≥ \frac{1}{4}∥v∥_{V}^{2} .

This implies, as in [8], that

∥e^{−Aξ}∥ ≤ e^{−ξ/4} for ξ ≥ 0. (52)

It is easy to see that for u, v ∈ V we can express F (v, u) − F (u, v) as

∫_{0}^{∞} (2v′ue^{−η} − vue^{−η} + ue^{−η/2} ∫_{0}^{η} e^{−s/2}v(s)ds − ve^{−η/2} ∫_{0}^{η} e^{−s/2}u(s)ds) dη

which implies that

|F (v, u) − F (u, v)| ≤ 4∥u∥∥v∥_{V}

and therefore D(A^{1/2}) = V, see [8]. Fix now any α ∈ (0, 1/4) and define

∥u∥_{1/2} = ∥(A − α)^{1/2}u∥ for u ∈ V.

Note that the norm ∥ · ∥_{1/2} is equivalent to the norm ∥ · ∥_{V}, see [5, 8].

If v ∈ V, then

v(η)^{2} = 2 ∫_{0}^{η} v′v ≤ 2∥v∥∥v′∥ ≤ ∥v∥_{V}^{2}

hence ∥v∥_{∞} ≤ ∥v∥_{V} which implies that the nonlinearity h : V → H, given by (48), is well defined and that

∥h(v)∥ ≤ 3∥v∥_{V} ∥v∥ ≤ c_{1}∥v∥_{1/2}^{2} for v ∈ V. (53)

Hence

∥h(u) − h(v)∥ ≤ 3∥u − v∥_{V} (∥u∥_{V} + ∥v∥_{V}) for u, v ∈ V,

which implies that h is locally Lipschitz.

Theorem 6.1. For every v_{0} ∈ V there exists τ > 0 and a unique v ∈ C([0, τ), V) such that v(0) = v_{0} and

v′(ξ) + Av(ξ) = h(v(ξ)) for ξ ∈ (0, τ).

Moreover, τ can be chosen so that either τ = ∞ or τ < ∞ and sup_{0<ξ<τ} ∥v(ξ)∥_{V} = ∞. Furthermore, there exist δ > 0 and c < ∞ such that if ∥v_{0}∥_{V} < δ, then τ = ∞ and

∥v(ξ)∥_{V} ≤ c∥v_{0}∥_{V} e^{−αξ} for ξ ≥ 0.

(Recall that α ∈ (0, 1/4) is arbitrary.)

Proof. In view of the above discussion, standard theory implies local existence and the ‘moreover’ part. See for example Appendix 2 of [7] or [5]. We also have

v(t) = e^{−At}v(0) + ∫_{0}^{t} e^{−A(t−s)}h(v(s))ds for t ∈ [0, τ). (54)

Since there exists c_{2} < ∞ such that

∥e^{−At}u∥_{1/2} ≤ c_{2}t^{−1/2}e^{−αt}∥u∥ for u ∈ H, t > 0,

see [5], we have that (52), (53) and (54) imply

∥v(t)∥_{1/2} ≤ e^{−αt}∥v(0)∥_{1/2} + c_{1}c_{2} ∫_{0}^{t}(t − s)^{−1/2}e^{−α(t−s)}∥v(s)∥_{1/2}^{2}ds.

If w(t) = e^{αt}∥v(t)∥_{1/2}, then

w(t) ≤ w(0) + c_{1}c_{2} ∫_{0}^{t}(t − s)^{−1/2}e^{−αs}w(s)^{2}ds.

If g(t) = max_{0≤s≤t} w(s) and c_{3} = max_{t≥0} c_{1}c_{2} ∫_{0}^{t}(t − s)^{−1/2}e^{−αs}ds, then

g(t) ≤ g(0) + c_{3}g(t)^{2} for t ∈ [0, τ).

If 0 < g(0) < 1/(4c_{3}), then the quadratic equation c_{3}x^{2} − x + g(0) has roots x_{1}, x_{2} such that

g(0) < x_{1} < 2g(0) < x_{2}.

Since g is continuous we have to have that g(t) ≤ x_{1} for all t ∈ [0, τ) which precludes the case τ < ∞ and therefore

e^{αt}∥v(t)∥_{1/2} ≤ g(t) ≤ 2g(0) = 2∥v(0)∥_{1/2} for t ≥ 0

and the equivalence of norms ∥ · ∥_{1/2}, ∥ · ∥_{V} complete the proof.

7. Results and discussions. Fig. 6 shows our results for the normalized shear stress |f_{nn}(0, ξ)|. The exact numerical solution from Section 4 is compared with the small time series solution. For impulsive start the shear is singular, or infinity, when time or ξ is at 0^{+}. It rapidly decreases, then approaches the exact steady state value of unity as ξ → ∞. Our small time expansion (39), being extremely accurate for small ξ, shows the singularity is of the form ξ^{− \frac{1}{2}} as ξ → 0. Since (39) is divergent for ξ ≥ 1, such expressions are normally useless for ξ ≥ 0(1). However, when various truncations of this small-time series are plotted, we observe that the range of validity can be vastly extended to ξ > 1. As an asymptotic series for ξ → 0, the error is the magnitude of the next order term, which is fairly small since we have shown that the coefficients decay exponentially as the order is increased. We are thus able to extend the validity of (39) to about ξ = 1.5, which is almost steady-state.

Of course the final approach to equilibrium cannot be inferred from the small-time series solution. Our numerical solution showed that the decay to steady state is exponential. This is supported by our analysis of the perturbation about the steady state. As we noted before, the steady state, or Crane’s solution, is a rare closed form nonlinear exact solution of the Navier-Stokes equations which is particularly amenable to analysis.

In conclusion, we have successfully analyzed the complete transient behavior of the unsteady viscous flow caused by a stretching surface. Our methods may also be useful to other practical initial value problems.

Fig. 6. The transient decay of normalized local shear stress, — — — — numerical, − − − − series, − · − · − · − steady state value. Numerals denote the number of terms retained.

REFERENCES

[1] E.R. Benton, On the flow due to a rotating disk, J. Fluid Mech., 24 (1966), pp. 781–800.

[2] H. Blasius, Grenzschichten in Flüssigkeiten mit kleiner Reibung, Zeit. f. Math. Phys., 56 (1908), pp. 1–37.

[3] L.J. Crane, Flow past a stretching plate, ZAMP, 21 (1970), pp. 645–647.

[4] S. Goldstein and L. Rosenhead, Boundary layer growth, Proc. Camb. Philo. Soc.,32 (1936), pp. 392–401.

[5] D. Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, No. 840, Springer, 1981.

[6] M. Katagiri, Flow due to an impulsive rotation of infinite disk, Phys. Fluids, 17 (1974), pp. 1463–1464.

[7] M. Miklavčič, Stability for semilinear parabolic equations with noninvertible linear operator, Pacific J. Math. 118 (1985), pp. 199–214.

[8] M. Miklavčič, Galerkin approximations for weakly nonlinear second order evolution equations, Funkcialaj Ekvacioj, 33 (1990), pp. 291–305.

[9] F. S. Sherman, Viscous Flow, McGraw-Hill, New York, (1990).

[10] D. P. Telionis, Unsteady Viscous Flows, Springer, NewYork, (1981).

[11] J. Vleggaar, Laminar boundary-layer behavior on continuous accelerating surfaces, Chem. Eng. Sci. 32 (1977), pp. 1517–1525.

[12] C. Y. Wang, Exact solutions of the steady-state Navier-Stokes equations, Ann. Rev. Fluid Mech., 23(1991), pp. 159–177.