- Email: [email protected]
An averaging method for the weakly unstable shallow water equations in a flat inclined channel
Fluid Dynamics Research 38 (2006) 469 – 488
An averaging method for the weakly unstable shallow water equations in a flat inclined channel Richard Spindler∗ Department of Mathematics, Bemidji State University, Bemidji, MN 56601, USA Received 21 May 2005; received in revised form 31 January 2006; accepted 27 February 2006 Communicated by E. Knobloch
Abstract The asymptotic behavior of small disturbances as they evolve in space from boundary conditions in a flat inclined channel under weakly unstable conditions is determined using a method of averaging. The shallow water equations are first transformed to a normal form using a variation of parameters technique, and then a new method of averaging is defined. The method of averaging developed differs from averaging methods for ordinary differential equations and other partial differential equations by averaging along the characteristics of the linearized equations rather than along the time variable. It is found that the solution is dominated by the evolution of the solution along one characteristic. Both asymptotic and numerical results for periodic disturbances are presented. In addition, connections to dynamical systems concepts are described. © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. PACS: 47.20. − k; 47.35.+i; 47.60.+i MSC: 34E10 Keywords: Averaging; Instability; Shallow water
1. Introduction and linear instability Normalization and asymptotic averaging techniques are developed to analyze the semi-infinite boundary value problem of shallow water waves in an inclined constant-slope channel, where the evolution is in ∗ Tel.: +1 802 865 5245.
E-mail address: [email protected]. 0169-5983/$30.00 © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. doi:10.1016/j.fluiddyn.2006.02.006
470
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
the spatial variable. The system is assumed to be in a “weakly unstable” state, meaning the Froude number F is slightly greater than the neutrally stable value (Yu and Kevorkian, 1992). The methods developed here could also work in the initial value problem with minor modifications. The averaging technique created in this study differs from previous averaging methods in several respects. First, and most importantly, the averaging is performed over two variables, namely the characteristics of the reduced (linear) problem. It appears that in previous averaging applications, the averaging was performed on all equations along one variable, namely the “time” variable. An important exception should be noted. Chikwendu and Kevorkian (1972) developed a method of multiple scales for a class of nonlinear wave equations which in the final result is an average along characteristics. Hence, there is some precedent for this type of averaging. Second, the normalization process (where the equation(s) are transformed to Lagrange standard form, to be defined in Section 2) in many, though not all, averaging methods are analogous to the variation of parameters used for first order linear ordinary differential equations with exponential solutions. A variation of parameters idea is used in this study, but it is not exactly analogous to first order linear ordinary differential equations. In the shallow water equations, shallow means that the depth is small compared to some characteristic length, such as the wavelength. In addition, it is assumed the depth is small compared to the breadth of the conduit. Using dimensionless variables, the shallow water equations are then ht + uhx + hux = 0, ut + uux + hx = 1 −
(1) u2 . hF 2
(2)
See Kevorkian (1990, Section 5.1.1) or Spindler (2005, Chapter 1). Here h is the height of the free surface normal to the channel bottom and is normalized by the uniform height. Also, u is the flow speed parallel to the channel bottom averaged over the cross-section of the channel, and is normalized by the uniform flow. See Fig. 1. (The uniform solutions are the constant solutions of the dimensional shallow water equations essentially corresponding to equilibrium, smooth, stratified flow.) Note that this is an inviscid flow and yet has a quadratic drag law. This is due to the channel frictional effects modeled by a Chezý law. In the derivation of these equations, the Froude Number F is defined to be tan F := , C where is the angle of inclination and C is the coefficient of friction. It can also be written as Us2 , F= gH s cos Thus, the Froude number is of the where Hs and Us are the uniform depth and flow speed, respectively. order of the ratio of the uniform fluid velocity to the velocity gH s of the gravity wave. It represents a balance between inertia and gravity. In general, the one-dimensional shallow water equations have been analytically insolvable. The uniform solutions of (1) and (2) are hs := 1 and us := F . In dynamical systems terms, this is a rest, or fixed, point. The next obvious question to ask is: under what conditions is the system locally stable or unstable near this rest point? There is only one parameter, the Froude number, in this model and hence
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
471
1
Free Surface h (x,t)
u (x,t)
Bottom
θ X
Fig. 1. Flat inclined channel.
the question of stability focuses on the Froude number. Linear instability occurs, by definition, at values of the parameter(s) of a given system where the solution to the linearized problem near the rest point grows without bound. Whitham (1974, pp. 85–86) showed that the flow is linearly unstable for F > 2 and Spindler (2005) confirms this for the boundary value problem. Recall that the usual procedure is to substitute in traveling wave solutions ei(kx−t) with k real. Then the regions of instability are determined by conditions under which the imaginary part of is positive. Here, the same procedure was executed, but because it is a boundary value problem, is assumed real. Then the regions of instability are determined by conditions when the imaginary part of k is negative. The instability results say in effect that, for the linearized system, a transcritical bifurcation occurs at F = 2. Here the system’s (linear) stability at this rest point changes from stable to unstable as F increases through 2. In addition, the original system is assumed stable for F < 2 and unstable for F > 2 in some neighborhood of the rest point and in some neighborhood of F = 2. As described in slightly more detail in Section 3.1, this is similar to the statement in ordinary differential equation dynamical systems theory that the 1-jet of the system is sufficient to determine local stability near the rest point for F = 2. For the nonlocal situation, the nonlinear terms of the model come into play by modifying this instability in the far-field (large x). Weak instability, defined by 0 < F − 2>1, is the region of concern in this study. Assuming a small boundary disturbance to uniform flow (the inlet conditions) (1)
h(0, t) = 1 + hb (t), (1)
u(0, t) = F + ub (t),
(3) (4)
472
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488 x
Characteristic of slope one Q
x=t
P
Ω
Characteristics of slope one andt hree
t Boundary Values
Fig. 2. Domain of interest: .
the method of averaging is developed and applied to this weakly unstable situation. In dynamical system terms, we “start” at a point near the rest point where the above linear instability result applies and develop a method of averaging to observe how the system evolves. Thus, inlet conditions perturbed from the rest state are being specified. To model the weak nonlinearity, set F = 2 + where > 0 is “small”. Note that this ties together the physical description of the system (the Froude number) with the boundary condition. Before proceeding, the domain of interest is described. Recall that in a hyperbolic system the solution at a point is determined by the backward characteristics to the Cauchy data, which in√this case is√the boundary data. The slopes of the characteristics of the shallow water equations are u + h and u − h. To a first approximation, these are F + 1 and F − 1, since u = F and h = 1 are the uniform solutions. In the case here, F ≈ 2, so the characteristic slopes are given approximately by 3 and 1. The relevant domain for the boundary value problem is found by observing that the backward characteristics from the points in the domain := {(t, x) ∈ R2 : 0 < x < t} intersect the positive t axis on which the inlet conditions are specified. In Fig. 2, P is in the domain because the solution there is determined solely by the evolution from the boundary data. On the other hand, the backward characteristic of slope 1 from the points in the region {(t, x) ∈ R2 : 0 < t < x} intersect the positive x axis, on which data is not specified. Thus, to determine the solution at point Q in Fig. 2, initial data must also be specified. In this study, only boundary data is considered, and hence only the region is considered. 2. Introduction to averaging The method of averaging has been used to derive asymptotic evolution equations and to prove uniform convergence of asymptotic expansions over expanding intervals of order 1/ for ordinary differential
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
473
equations in Lagrange standard form (to be defined shortly). Sanders and Verhulst (1985) provide an excellent account of this theory, as does Murdock (1991). In this study, a method of averaging that applies to the shallow water equations is created. Krol (1991) developed a method for an advection–diffusion equation by averaging along the time variable. Buitelaar (1993) defines a time averaging method for systems of equations in abstract Banach spaces of the form dw + Aw = f (w, t, ), dt where A is the generator of a C0 group of unitary operators on a Hilbert Space. By Stone’s theorem, this implies iA is self-adjoint. In all of the applications considered by Buitelaar (e.g. nonlinear wave equations), this self-adjointness requires homogeneous boundary conditions. This study’s approach differs from that of Krol or Buitelaar in a couple of ways. First, the average is applied along the characteristics of the linearized equations, and second, the averaging is applied to a hyperbolic system. As noted above, Chikwendu and Kevorkian (1972) developed a method of multiple scales for a class of nonlinear wave equations which in the final result is an average along characteristics. The idea behind averaging methods is as follows. Assume an ordinary differential equation can be transformed into Lagrange standard form, which is defined to be an ordinary differential equation of the form dx = f (t, x, ), dt where these could be vector quantities and where f is real analytic in . In addition, initial conditions are specified at t = 0. In effect, this form implies that t is a slow variable. Now suppose we expand f in a Taylor’s series in to obtain dx = f (1) (t, x) + 2 f (2) (t, x, ). dt
(5)
The system is now in a normal form of order 1 (Murdock 2003, pp. 12–13). Usually f is required to be Lipschitz continuous in x. Also, often the existence of the average quantity 1 T (1) (1) f (x) := lim f (t, x) dt (6) T →∞ T 0 is imposed. In that case, Sanders and Verhulst (1985) call f a KBM-vectorfield (after Krylov, Bogoliubov, and Mitrolosky). An alternative assumption used is that f is periodic in time. To understand the mathematical origin of averaging, perform what is sometimes called a near-identity transformation x = x˜ + g(t, x) ˜
(7)
to transform Eq. (5) to the “autonomous” (in t) problem dx˜ = p(1) (x) ˜ + 2 p(2) (t, x, ˜ ). dt
(8)
474
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
What makes this work is for p(1) to be the average of f (1) holding x constant, namely Eq. (6). To see why, differentiate Eq. (7) with respect to t dx dx˜ jg = + x˜ g · x˜t + dt dt jt jg = p(1) (x) + O(2 ) ˜ + (9) jt using Eq. (8) and the fact that x˜t is O(). Now express xt in terms of x˜ instead of x by expanding the right-hand side of Eq. (5) around x˜ in a Taylor’s series. (Recall that x and x˜ differ by an order of .) Equate the order terms of Eqs. (5) and (9) to obtain f (1) (t, x) ˜ = p(1) (x) ˜ +
jg jt
.
Now average over t, 1 1 (g(t, x) ˜ − g(0, x)) ˜ = lim lim T →∞ T T →∞ T
T
f (1) (t, x) ˜ dt − p(1) (x). ˜
0
If g is uniformly bounded, the left-hand side is zero, and p(1) is found to be equal to the average of f (1) . Thus, dx˜ = f (1) (x) ˜ + 2 p(2) (t, x, ˜ ) dt is the new differential equation. In effect, the t dependence has been “pushed” up to an order of . It can be shown that the solution to the averaged system dz = f (1) (z) dt is an o(1) approximation to the solution to the original problem on expanding intervals of 1/, i.e. on intervals 0 t L/ for some constant L (Sanders and Verhulst, 1985). Murdock (2003) calls the method of using a near-identity transformation to simplify, or normalize, the right-hand side of Eq. (5) format 1a, or the direct/iterative format. To see how normal form theory works in this context, append a zeroth component to x and f (i) as follows: x0 = t, (i)
f0 = 0,
i 1
and define a zeroth order column vector f (0) = [1 0 · · · 0]T so now Eq. (5) can be written as dx = f (0) + f (1) (t, x) + 2 f (2) (t, x, ). dt Define the homological operator Lu := Duf (0) − Df (0) u,
(10)
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
475
where u is a vector field as a function of the vector x and D indicates the total derivative with respect to x. (Thus, these derivatives are matrices.) Then, for format 1a, Eq. (10) transforms to dx˜ = f (0) + [f (1) (t, x) ˜ − (Lg)(t, x)] ˜ + 2 f (2) (t, x, ˜ ). dt From the perspective of normal form theory, the idea then is to choose p(1) so that f (1) − p(1) is in the range of L (which then determines g) and in such a way as to make the problem simpler. The way in which this is done is called the style. In the case here, Lg = jg/jx0 = jg/jt. For the style, p(1) is chosen to be the projection of f (1) into the kernel of L, i.e. p (1) is the time average of f (1) as found above. The approach in this study is similar to the procedure outline above. Note that the progression through the above derivation showed how the method of averaging comes about. If one knows how to average already, then just define t g(t, z) = [f (1) (s, z) − f (1) (z)] ds 0
in the near-identity transformation and work with that (Verhulst, 1996). This also arises naturally from the normal form theory discussed above. In addition, note that if := t (slow time), then the averaged system is dz = f (1) (z), d which indicates the connection to multiple scales. Finally, if f (1) is periodic with period T, the same method is used except the integration is over a period. (In that case g turns out to be periodic.) 3. The method of averaging applied to the weakly unstable shallow water equations A method of averaging for the weakly unstable boundary value shallow water system is now developed, applying similar reasoning as used for ordinary differential equations described above. In dynamical system problems for ordinary differential equations, the systems are often assumed to be in Lagrange standard form. Transforming a system of partial differential equations into this form is not trivial, and that is true in this case too. One common method of attack, mentioned earlier and applied in the next section, is to use a variation of parameters procedure. In that approach, the reduced (linear) problem is solved, which results in some “constants” of integration. Then the full problem is transformed by letting these “constants” of integration depend on all independent variables. Once this is finished, the method of averaging can be defined. Note that all of the calculations related to this work can be found in Spindler (2005, Chapter 3), and also are available by contacting the author. 3.1. The reduced problem and Lagrange standard form The first task is to transform the shallow water equations (1) and (2) into Lagrange standard form. To do this, the reduced problem is solved first. Center the shallow water equations (1) and (2) around the
476
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
uniform values by the trivial change of variables h = 1 + h(1) , u = F + u(1) ,
(11)
which define h(1) and u(1) for = 0, and substitute into the shallow water equations. (This is not a perturbation expansion.) Divide by , multiply the second equation through by 1 + h(1) , and simplify the right-hand side of the second equation to arrive at (1)
+ F hx(1) + ux(1) = −{u(1) hx(1) + h(1) ux(1) }, 1 2 (1) (1) ut + hx(1) + F ux(1) − h(1) − u(1) = − 2 u(1)2 − [u(1) ux(1) + h(1) ut + h(1) hx(1) F F (1) (1) (1) (1) (1) + F h ux ] − h u ux .
ht
(12)
Of course, this change of variables centers the rest point at the origin, but it also in effect “lowers” the order by 1 compared to the original shallow water equations. Note the time derivative on the right-hand side of the second equation. Before inserting F = 2 + , the reduced problem is determined and solved. Set = 0, so that Eqs. (12) reduce to the linearized shallow water equations with F = 2. Define hˆ (1) and uˆ (1) to be the solutions of these linearized equations. Add and subtract the two equations (12) (with = 0) and define Sˆ := hˆ (1) + uˆ (1)
and
Rˆ := hˆ (1) − uˆ (1)
to obtain Sˆt + 3Sˆx − Rˆ = 0, Rˆ t + Rˆ x + Rˆ = 0. The characteristic slopes of this system are 1 and 3 and so the characteristics of this system are x − t= constant and x − 3t= constant. Change variables along these characteristics, that is (x, t) := x − t
and
(x, t) := x − 3t,
(13)
which transforms the first order to ˆ , ), Sˆ (, ) = 21 R(
(14)
ˆ , ). Rˆ (, ) = 21 R(
(15)
In effect, the linearized system has been transformed to normal form in the sense of hyperbolic partial differential equations (Forsythe and Wasow, 1960). The solution of Eq. (15) is found to be ˆ , ) = rˆ ()e/2 , R(
(16)
where rˆ () is in effect a “constant” of integration. Using this in the Eq. (14) above gives ˆ , ) = R ˆ ()e/2 + sˆ (), S(
(17)
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
477
Ω
= 3 Fig. 3. Domain of interest: .
where sˆ () is in effect another “constant” of integration and 1 ˆ () := R rˆ ( ) d . 2
(18)
Note that the primes here do not designate differentiation. The variable is just a dummy integration variable. The domain will be needed in the new (, ) coordinate system. When x = 0, Eqs. (13) are = −t and = −3t. Eliminating t determines the boundary = 3. Similarly, x = t implies = 0. Also, x 0 maps to the region 3, and x < t maps to the region < 0 under the transformation (13). (The boundary t 0 maps to the region which is less restrictive than 3.) Thus, = {(, ) : 3 0}. See Fig. 3. Note that using these transformations, the concept of jet sufficiency of ordinary differential equations can be related to the situation here. Applying the change of variables S := h(1) + u(1) , R := h(1) − u(1) , := x − (F − 1)t, and := x − (F + 1)t, Eqs. (12) are then transformed to a normal form
1 1 1 1 1 − S+ + R + m1 , S = 2 2 F 2 F
1 1 1 1 1 R = − S+ + R + m2 . (19) 2 2 F 2 F The remainder terms m1 and m2 terms are quadratic polynomials in S, R, and their derivatives. Because only linear transformations have been applied to the system, its stability properties are not changed. In ordinary differential equation dynamical systems theory, a system in the form of (19) where the right-hand side is expanded into a Taylor’s series in the dependent variables with F = 2, is 1-jet sufficient since the linear polynomial terms determine the stability locally around the rest point. If F = 2, the quadratic terms
478
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
would need to be examined to determine if the system is 2-jet sufficient, and so on. See Murdock (2003, p. 8) for a definition and discussion of sufficiency of jets in the context of ordinary differential equations. The reason for determining the solution to the reduced problem is to transform the original problem (12) using this reduced solution by the method of variation of parameters. First, to determine the full dependence of the full equations on , substitute F = 2 + into the second equation of (12) and multiply both sides by F 2 to obtain (1)
4(ut
(1)
+ hx(1) ) + 8ux(1) − 4(h(1) − u(1) ) + 4(ut (1)
= {−u(1)2 − 4[u(1) ux(1) + h(1) ut
+ hx(1) ) + (12ux(1) − 2h(1) ) − 2(h(1) − u(1) )
+ h(1) hx(1) + 2h(1) ux(1) ]} + O(2 ).
All 2 terms are summarized in O(2 ). (Note that here the O(2 ) refers to the terms in this equation, of course, which would imply O(3 ) in the original shallow water equations.) Moving terms to the right-hand side, the equations become (1)
+ 2hx(1) + ux(1) = −(u(1) hx(1) + h(1) ux(1) + hx(1) ),
(1)
+ hx(1) + 2ux(1) − (h(1) − u(1) ) = {−(1 + h(1) )(ut + hx(1) + 2ux(1) ) − (1 + u(1) )ux(1) + h(1) − 21 u(1) − 41 u(1)2 } + O(2 ).
ht ut
(1)
(20)
Observe that (20) is recursive in the sense that part of the left-hand side is in the right-hand side, namely (1) (1) (1) ut + hx + 2ux . Substitute this in for itself on the right-hand side, moving higher order terms into O(2 ) to obtain (1)
+ 2hx(1) + ux(1) = −(u(1) hx(1) + h(1) ux(1) + hx(1) ),
(1)
+ hx(1) + 2ux(1) − (h(1) − u(1) ) = {−(1 + u(1) )ux(1) + 21 u(1) − (h(1) − 21 u(1) )2 }
ht ut
+ O(2 ).
(21)
(22)
As in the reduced problem, transform the left-hand side of (21) and (22) into normal form (in the sense of hyperbolic partial differential equations), so let S := h(1) + u(1)
and
R := h(1) − u(1)
and eliminate h(1) and u(1) in the equations. (Add and subtract the equations). Again, transform along the characteristics of the linearized equations, (x, t) = x − t and (x, t) = x − 3t, to obtain S − 21 R = 21 {−(1 + 43 S − 41 R)(S + S ) + 41 (S + R)(R + R ) + 41 (S − R) −
1 16 (S
+ 3R)2 } + O(2 ),
(23)
R − 21 R = 21 {(1 + 41 S − 43 R)(R + R ) + 41 (S + R)(S + S ) + 41 (S − R) −
1 16 (S
+ 3R)2 } + O(2 ).
(24)
Next, change variables using the reduced solutions, Eqs. (16) and (17), by the method of variation of parameters. Recall that sˆ and rˆ are “constants” of integration for the reduced equations with the first
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
479
depending on and the second on . Now let s depend also on and r also on . Therefore, define r(, ) := R(, )e−/2 , s(, ) := S(, ) − R(, )e/2 , where R is similar to Eq. (18): 1 R(, ) = r( , ) d . 2
(25)
Determine the derivatives of R and S and substitute these into the equations for S and R , (23) and (24), to obtain s = 21 {−[1 + 43 (s + Re/2 ) − 41 re/2 ](s + 21 re/2 + s + R e/2 + 21 Re/2 )
+ 41 (s + Re/2 + re/2 )(r e/2 + r e/2 + 21 re/2 ) + 41 (s + Re/2 − re/2 )
−
1 16 (s
+ Re/2 + 3re/2 )2 } + O(2 )
and r = 21 {(1 + 41 s + Re/2 − 43 re/2 )(r e/2 + r e/2 + 21 re/2 )
+ 41 (s + Re/2 + re/2 )(s + 21 re/2 + s + R e/2 + 21 Re/2 ) + 41 (s + Re/2 − re/2 )
−
1 16 (s
+ Re/2 + 3re/2 )2 }e−/2 + e−/2 O(2 ).
Re-arrange and combine terms. Also, because s and r are equal to times terms on the right, move s , r , and R terms on the right up into O(2 ) and e−/2 O(2 ) to arrive at
1 3 1 2 s = − 1 + s s + s − s 2 4 4 16
1 1 1 5 1 3 + − (3R − r)s + − r − R + r s + − R − r e/2 4 8 2 4 4 4
5 2 1 7 2 1 + − r − r R − R + r (r + R) e + O(2 ) (26) 16 2 16 4 and
1 1 2 −/2 1 1 1 1 1 1 ss + s − s e r − r s + (r + R)s r = + r + r + R + 2 4 4 16 4 4 4 8 4
13 2 1 2 1 3 + − r + R + R − r r e/2 + e−/2 O(2 ). (27) 16 16 4 4
Basically these independent variables introduced into s and r correspond to slowly varying variables for their respective function, meaning that corresponds to a slowly varying independent variable for s and corresponds to a slowly varying independent variable for r. Observe that everything has been exact so far. (Although the higher order terms are hidden in the O(2 ) term.) Various transformations and algebraic operations have just been applied to the problem. Note that for = 0 the reduced solutions are recovered. Also observe that, at least up to order O(), s = 0 and r = 0 are rest points as they should be, corresponding to h(1) = 0 and u(1) = 0.
480
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
3.2. Near-identity transformation The above system of equations, (26) and (27), is now expressed in normal form in the sense of partial differential equations. This is similar to the normal form as used in dynamical systems for ordinary differential equations (Murdock, 2003), but there are differences. For example, on the left-hand side, the derivatives are in terms of two independent variables instead of one, of course, since these are partial differential equations. Also, derivatives appear on the right-hand side. Despite these differences, an averaging theory can be developed along the lines as described in Section 2. This will also be analogous to the transformation of normal forms (Murdock, 2003). Now that the equations are in normal form, the method of averaging along characteristics arises formally by applying the idea of a near-identity transformation defined by
g1 (, s˜ , s˜ , r˜ , r˜ ) s˜ s = + . (28) r r˜ 2 e−/2 g2 (, s˜ , s˜ , r˜ , r˜ ) The goal is to define this transformation such that the Lagrange standard form equations (26) and (27) are “autonomous” and uncoupled from each other
s˜ = p1 (˜s , s˜ ) + O(2 ), 2
r˜ = p2 (˜r , r˜ ) + e−/2 O(2 ) 2
(29) (30)
solving the same boundary conditions, namely s˜ (, 3) = s(, 3) and r˜ (, 3) = r(, 3). This is similar to the example in Section 2 for ordinary differential equations where the independent variable t was eliminated. The difference here is that dependent variables are being eliminated from the right-hand side functions instead of an independent variable. For convenience, a factor e−/2 of g2 is inserted in the second equation. Assume that gi functions and their derivatives are uniformly bounded. Also, by definition of change of variables, the near-identity transformation must be invertible in a neighborhood of =0 which is assumed. To proceed, differentiate (28) to obtain
s s˜ = r r˜
g1 + g1˜s s˜ + g1˜s s˜ + g1˜r r˜ + g1˜r r˜ + . (31) 2 e−/2 g2 − 21 e−/2 g2 + e−/2 (g2˜s s˜ + g2˜s s˜ + g2˜r r˜ + g2˜r r˜ + g2 ) Recall that s˜ = O() and r˜ = e−/2 O(). Thus, those terms are moved up into higher orders. Also assume the same is true of higher derivatives of these terms. Then
g1 + g1˜r r˜ + g1˜r r˜ s O(2 ) s˜ + = + e−/2 O(2 ) r r˜ 2 e−/2 g2 − 21 e−/2 g2 + e−/2 (g2˜s s˜ + g2˜s s˜ + g2 ) ⎤ ⎡ d g1
2) ⎥ d s˜,˜s constant ⎢ s˜ O ( ⎥ + −/2 2 . = + ⎢ (32) ⎦ e r˜ O ( ) 2 ⎣ d −/2 (e g2 ) d r˜ ,˜r constant
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
481
Observe that by the near-identity transformation (28), s and s˜ differ by O() and r and r˜ differ by e−/2 O(). Expand the right-hand side terms of Eqs. (26) and (27) (as functions of s and r and their derivatives) in a Taylor’s series around s˜ and r˜ and their derivatives, to simply replace s with s˜ and r with r˜ , with the remainder terms of the Taylor’s Series moved up into the O(2 ) and e−/2 O(2 ) terms. (Note that to apply this Taylor Series, g1 and g2 and their derivatives must be uniformly bounded.) Next substitute Eqs. (26) and (27) (with s˜ and r˜ , etc., instead due to the discussion in the last paragraph) into the left-hand side of Eq. (32), and Eqs. (29) and (30) into the right-hand side of Eq. (32). After cancelling an from both sides and re-arranging, d p1 (˜s , s˜ ) + g1 d s˜,˜s constant
1 3 1 = − 1 + s˜ s˜ + s˜ − s˜ 2 4 4 16
1 1 5 1 3 1 ˜ ˜ ˜ + − (3R − r˜ )˜s + − r˜ − R + r˜ s˜ + − R − r˜ e/2 4 8 2 4 4 4
5 1 7 2 1 ˜ ) e + O(), ˜ − R ˜ + r˜ (˜r + R + − r˜ 2 − r˜ R (33) 16 2 16 4 d −/2 p2 (˜r , r˜ ) + (e g2 ) d r˜ ,˜r constant
1 2 −/2 1 1 s˜ s˜ + s˜ − s˜ e = 4 4 16
1 1 1 1 1 ˜ + ˜ )˜s r˜ − r˜ s˜ + (˜r + R + r˜ + r˜ + R 4 4 4 8 4
1 2 3 1 13 2 /2 ˜ ˜ R − r˜ r˜ e + − r˜ + R + 16 16 4 4 + e−/2 O().
(34)
˜ (Note that the cancellation of lowers the order 1 more relative to the shallow water equations.) Here, R is defined similar to R except using r˜ instead of r in the integral defining R. Examine Eq. (33) first. Move p1 to the other side of Eq. (33) and integrate as follows:
1 2 3 1 g1 = − 1 + s˜ s˜ + s˜ − s˜ 4 4 16
1 1 1 5 1 3 /2 ˜ ˜ ˜ − (3R − r˜ )˜s + − r˜ − R + r˜ s˜ + − R − r˜ +e 4 8 2 4 4 4
5 1 7 2 1 ˜ ) − p1 (˜s , s˜ ) d ˜ − R ˜ + r˜ (˜r + R +e − r˜ 2 − r˜ R 16 2 16 4 s˜ ,˜s constant
+ O()
482
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
or
1 2 1 3 1 /2 ˜ − r˜ ) d − s˜ (3R g1 = − 1 + s˜ s˜ + s˜ − s˜ + e 4 4 16 4
1 5 1 1 3 ˜ + r˜ d + ˜ − r˜ d − r˜ − R − R +˜s 8 2 4 4 4
5 1 7 2 1 ˜ − R ˜ + r˜ (˜r + R ˜ ) d − r˜ 2 − r˜ R + e 16 2 16 4 − p1 (˜s , s˜ ) + c1 (, s˜ , s˜ ) + O(). Next divide both sides by and let → −∞ and continue to treat s˜ and s˜ as constant. This limit could conceivably be taken in a number of ways, as long as t > x 0, that is 3 < 0 (to remain in the domain as discussed in the first section). However, another requirement for choosing the curve along which to take the limit is that singularities (shocks) must not be encountered in the limit. The boundary data curve x = 0 certainly satisfies this requirement, and hence the limit is defined along this curve. Since g1 is assumed uniformly bounded, the left-hand side tends to zero. Suppose also there exists a positive constant C(x) 0 for each fixed x such that ˜ |, |R ˜ | C(x)e−3/2 |˜r |, |˜r |, |R
(35)
uniformly in t. (It will be seen in Section 4 that the solution is consistent with this assumption.) By this assumption, each of the terms above, except the first and the p1 and c1 terms, will be of the form e1/2(−3)
=
e−x
or
e(−3)
=
e−2x
,
which approaches zero uniformly as → −∞ since x is fixed. Assuming the integration “constant” c1 is uniformly bounded, all that remains is p1 (˜s , s˜ ) = −(1 + 43 s˜ )˜s + 41 s˜ −
1 2 16 s˜
+ O().
Substitute this into (29) to obtain one evolution equation. Process Eq. (34) in a similar manner. First, integrate (34) holding r˜ and r˜ constant,
1 2 − /2 1 1 1 1 ˜ g2 = e s˜ s˜ + s˜ − s˜ e + r˜ + r˜ + R 4 4 16 4 4
1 1 1 ˜ )˜s r˜ − r˜ s˜ + (˜r + R + 4 8 4
1 2 1 3 13 2 /2 ˜ ˜ R − r˜ r˜ e − p2 (˜r , r˜ ) d + − r˜ + R + 16 4 4 16 r˜ ,˜r constant /2
+ e/2 c2 (, r˜ , r˜ ) + O().
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
483
Then re-write the equation as
1 2 − /2 1 /2 1 s˜ s˜ + s˜ − s˜ e e g2 = 4 4 16
1 ˜ 1 ˜ 1 ˜ 1 1 ˜ ˜ ˜ ˜ R − R s˜ + (R + RR )˜s + R + R + RR + 4 4 4 8 4
−)/2 1 2 1 3 ˜ ˜ 13 ˜ 2 ( ˜ R˜ ) d ˜R + ˜ R − R R e R − p2 (R, + − R + R ˜ ˜ 16 16 4 4 R,R constant + e/2 c2 (, r˜ , r˜ ) + O() ˜ . (In the first equation above, r˜ and r˜ have been held constant ˜ R := e/2 R using R˜ = e/2 r˜ and where R in the integral, which is equivalent to holding R˜ = e/2 r˜ and R˜ = e/2 r˜ constant.) Also, it has been assumed that p2 is homogeneous in its arguments (meaning yp 2 (x1 , x2 ) = p2 (yx 1 , yx 2 )). Factor the R˜ terms out of the integrals to obtain
1 2 − /2 1 1 1 ˜ 1 /2 ˜ ˜ g2 = e s˜ s˜ + s˜ − s˜ e d + R + R + RR 4 4 16 4 4 1 ˜ 1 1 ˜ R )˜s R − R˜ + s˜ d + (R˜ + R 4 8 4
1 2 3 ˜ ˜ 1 13 ˜ 2 ˜ R˜ ) + e/2 c2 (, r, r ) + O(). ˜ ˜ RR − R R − p2 (R, + 2 − R + RR + 16 16 4 4 Divide both sides by , treat R˜ and R˜ as constant, and let → −∞. Again, choose the curve x = 0 along which to take this limit. Since g2 is assumed uniformly bounded, the left-hand side tends to zero. Also suppose that the various integrals involving s and s in this last equation are uniformly bounded as well as the integration “constant” e/2 c2 (, r˜ , r˜ ). In effect the solutions have been restricted to the class of functions satisfying these specifications. Then ˜ R˜ ) = [R˜ + 1 R˜ + 1 R ˜ p2 (R, 4 4 R ] + O() or ˜ ] + e−/2 O(). p2 (˜r , r˜ ) = [˜r + 41 r˜ + 41 R This result provides a second evolution equation by substituting into Eq. (30). Also observe that p2 , at least up to order , is homogeneous and so that assumption is consistent. From the work above, the averaging operation that appears to transform the equations as desired is to average Eq. (26) by 1 lim [· · ·] d →−∞ s,s constant
484
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
along the line x = 0, where s and s are treated as constants in the limit. Average Eq. (27) by e/2 −/2 e lim [· · ·] d →−∞ R,R constant
along the line x = 0, where R, R , and RR are treated as constants in the limit. The near-identity transformations gi are 1 1 5 1 /2 ˜ + r˜ d + ˜ − r˜ ) d + s˜ − s˜ − r˜ − R (3R g1 = e 4 8 2 4
5 1 7 1 ˜ − R ˜ 2 + r˜ (˜r + R ˜ ) d + e − r˜ 2 − r˜ R 16 2 16 4 and
1 1 1 ˜ 1 ˜ 1 2 − /2 d + s˜ s˜ + s˜ − s˜ e R − R g2 = e 4 4 16 4 8
13 1 2 3 1 ˜R + ˜ R − R˜ R˜ . + 2 − R˜ 2 + R R 16 16 4 4 /2
1 3 ˜ − r˜ d − R 4 4
1 ˜ R )˜s s˜ d + (R˜ + R 4
These are uniformly bounded if the assumptions on s˜ and r˜ hold, as required. To summarize, the equations in Lagrange standard form (26) and (27) transform to the evolution equations
3 1 2 1 − 1 + s˜ s˜ + s˜ − s˜ + O(2 ) (36) s˜ = 2 4 4 16 and
1 1 ˜ r˜ + r˜ + R + e−/2 O(2 ) r˜ = 2 4 4
(37)
under the near-identity transformation. In ordinary differential equations, the system of equations retaining the higher order O(2 ) terms is sometimes called the transformed system. Of course, the O(2 ) terms are dropped to approximate the actual solution. When these higher order terms are dropped,
3 1 2 1 − 1 + s˜ s˜ + s˜ − s˜ (38) s˜ = 2 4 4 16 and r˜ =
1 1 ˜ r˜ + r˜ + R 2 4 4
(39)
are sometimes called the averaged system. If ¯ = is defined, then the averaged Eq. (38) becomes s˜¯ = 21 [−(1 + 43 s˜ )˜s + 41 s˜ −
1 2 16 s˜ ]
and, if ¯ = is defined, then the second averaged (39) becomes ˜ ]. r˜¯ = 21 [r˜ + 41 r˜ + 41 R
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
485
In effect, the new variables define far-field scales for each equation. In ordinary differential equations, this system of equations is sometimes called the guiding system (Murdock, 1991). The development of the averaging method can be placed in the context of normal forms theory of dynamical systems. Here, the homological operator is defined to be ⎡ ⎤ d g1 ⎢ ⎥ d s˜,˜s constant ⎥. Lg := ⎢ ⎣ d −/2 ⎦ (e g2 ) d r˜ ,˜r constant Let f (1) be the order term on the right-hand side of Lagrange standard form equations (26) and (27). Eqs. (33) and (34) result from transforming the normal form system (26) and (27) into format 1a by applying the near-identity transformation. The abstract form of Eqs. (33) and (34) is Lg = f (1) − p (1) .
The average operation and the final evolution equations then arise by defining p (1) to be the projection of f (1) into the kernel of L. 4. A comparison The solutions to the averaged evolution equations (38) and (39) are determined, assuming a sine wave boundary condition, namely h(0, t) = 1 + sin(t)
and
u(0, t) = F .
It is hoped that the solutions to the averaged equations are a good approximation to the solution of the transformed Eqs. (36) and (37). This means that the solutions are within o(1) over expanding intervals of x of order 1/ as in the ordinary differential equations case. Also, recall that s= s˜ + O() and r = r˜ +e−/2 O(). Hence R = R˜ + O() (where R˜ := e/2 r˜ ) and S = S˜ + O(). Thus, up to order , R and S, and hence h(1) ˜ In this section, the solution of Eq. (39) is solved explicitly for r˜ , and u(1) , are determined by R˜ and S. while Eq. (38) is solved numerically. Finally, the shallow water equations are numerically integrated and compared to the solution of the evolution equations, at a perturbation of = 0.1. The averaged r-evolution equation (39) is linear in r˜ with constant coefficients, so it can be solved analytically. The boundary condition for r˜ can be written r˜ (, 3) = − sin()e−3/2 ˜ Assume a solution of the form using the definitions of r˜ and R. r˜ (, ) = −e−(3/2)−(/2)a(3−) sin + b(3 − ) , 2
(40)
where the parameters a and b are to be determined. Substitute the proposed solution into the r˜ -evolution equation (39) to solve for a and b to obtain 5 3/4 a=− − 4 9 + 4 2
and
b=
/2
9 + 4 2
− .
(41)
486
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488 Direct Average
0.08 0.06 0.04
h(1)
0.02 0 -0.02 -0.04 -0.06 -0.08 37
38
39
40 Time t
41
42
43
Fig. 4. Height h1 of wave at x = 30.
Note that a is less than zero. Also observe that r˜ is consistent with assumption (35). Transforming back to the R˜ variable ˜ , ) = − e−(1/2)(1+a)(3−) sin + b(3 − ) R( 2 −(1+a)x = −e sin + b(3 − ) , 2 where of course R˜ = e/2 r˜ . This implies that the part of the solution from the characteristic decays exponentially in the x direction as long as < 1/|a|. Thus, with this restriction on , the contribution of R˜ to h(1) and u(1) is negligible as x becomes large. The far-field evolution arises from the contribution ˜ contribution to S˜ decays exponentially in the x of S˜ which is essentially determined by s˜ (since the R (1) (1) direction also). This also implies that h and u become nearly identical rather quickly. A basic forward–backward (FB) method is used to numerically integrate the s-evolution equation. In addition, the differences are designed to satisfy the CFL condition. The boundary condition in terms of s˜ is /2 ˜ s˜ , = − sin −e R , . 3 3 3 From Eq. (25) for R and the solution for r˜ 1 s˜ 3 sin , = − sin − + 2 cos . 2 3 3 9 + 4 3 3
(42)
Eq. (42) is the boundary condition for the s-evolution equation (38) used in the numerical integration. The shallow water equations are numerically integrated by applying the artificial viscosity method (Lapidus and Pinder, 1999, Section 6.2.3).Again, the differences are designed to satisfy the CFL condition.
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
487
Direct Average
0.08 0.06 0.04
h(1)
0.02 0 -0.02 -0.04 -0.06 -0.08 0
5
10
15 20 Distance x
25
30
35
Fig. 5. Height h1 of wave at t = 30.
Finally, all of the programs were written in MatLab 6.5. Fig. 4 is an overlay of the two solutions at x = 30 and = 0.1. Observe that the waves are periodic in t as expected, and the formation of a bore occurs. The waves are within about 2 of each other up to (and past) x ≈ 1/ as expected from perturbation theory. Fig. 5 is an overlay of the two solutions at t = 30 and = 0.1. Notice the gradual exponential growth of the wave in the x direction, as the linear stability theory predicted. Also observe that the waves steepen to a bore as x increases. Finally, the solutions from the two methods at t = 30 are still within about 2 of each other up to (and past) x ≈ 1/. 5. Conclusion The averaging method introduced in this study for the shallow water equations indicates a method of averaging for more general two-dimensional systems of hyperbolic partial differential equations. By averaging along the characteristics of the linearized equations, an approximate solution may be found in the far-field. It is hoped that this method can be extended to that more general situation. Note that Yu and Kevorkian (1992) analyzed the initial value problem for a flat inclined constant-slope channel using a multiple scales method, while Spindler and Yu (2006) analyzed the boundary value problem using multiple scales in conjunction with the Fredholm alternative theorem. The s-evolution equations found for mutliple scales there and for averaging here differ by one term, but the numerical solutions of these equations are both within order 2 of each other. These results suggest that the two solutions may be asymptotically equivalent, although this question of equivalence merits further study. For ordinary differential equations, Murdock (1991) shows that the first order two-scales method (using one far-field scale) is identical to averaging. In fact, Perko (1969) shows that the two methods for ordinary differential equations are asymptotically equivalent, though not necessarily identical, up to all orders of .
488
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
Although deriving the averaging technique is more involved than applying multiple scales, once the averaging operation is known, it is no more difficult to apply than determining the secular terms in multiple scales. However, the asymptotic validity of multiple scales for partial differential equations is an open mathematical question as well as its relationship to averaging. The advantage of the normalization and averaging technique is that it may provide a more tractable framework for deriving rigorous error estimates for both techniques applied to partial differential equations or systems of partial differential equations, just as with ordinary differential equations. Acknowledgments I would like to thank Ferdinand Verhulst for his helpful comments on this project. I am also grateful for the many suggestions the reviewers made to improve the manuscript. References Buitelaar, R.P., 1993. The method of averaging in Banach spaces: theory and applications. Ph.D. Thesis, Rijksuniversiteit Utrecht, Netherlands. Chikwendu, S.C., Kevorkian, J., 1972. A perturbation method for hyperbolic equations with small nonlinearies. SIAM J. Appl. Math. 22 (2), 235–258. Forsythe, G.E., Wasow, W.R., 1960. Finite-Difference Methods for Partial Differential Equations. Wiley, New York. Kevorkian, J., 1990. Partial Differential Equations: Analytical Solution Techniques. Wadsworth and Brooks/Cole, Pacific Grove, CA. Krol, M.S., 1991. Nearly time-periodic advection–diffusion problems. SIAM J. Appl. Math. 51 (6), 1622–1637. Lapidus, L., Pinder, G., 1999. Numerical Solution of Partial Differential Equations in Science and Engineering. WileyInterscience, New York. Murdock, J.A., 1991. Perturbation Theory and Methods. Wiley-Interscience, New York. Murdock, J., 2003. Normal Forms and Unfoldings for Local Dynamical Systems. Springer, New York. Perko, L.M., 1969. Higher order averaging and related methods for perturbed periodic and quasi-periodic systems. SIAM J. Appl. Math. 17, 698–724. Sanders, J.A., Verhulst, F., 1985. Averaging Methods in Nonlinear Dynamical Systems. Springer, New York. Spindler, R., 2005. The one-dimensional shallow water boundary value problem and instability. Ph.D. Thesis, University of Vermont, Burlington, VT. Spindler, R., Yu, J., 2006. Nonlinear spatial evolution of small disturbances from boundary conditions in flat inclined-channel flow. Stud. Appl. Math. 116, 173–200. Verhulst, F., 1996. Nonlinear Differential Equations and Dynamical Systems, second ed. Springer, Berlin. Whitham, G.B., 1974. Linear and Nonlinear Waves. Wiley, New York. Yu, J., Kevorkian, J., 1992. Nonlinear evolution of small disturbances into roll waves in an inclined open channel. J. Fluid Mech. 243, 575–594.
An averaging method for the weakly unstable shallow water equations in a flat inclined channel Richard Spindler∗ Department of Mathematics, Bemidji State University, Bemidji, MN 56601, USA Received 21 May 2005; received in revised form 31 January 2006; accepted 27 February 2006 Communicated by E. Knobloch
Abstract The asymptotic behavior of small disturbances as they evolve in space from boundary conditions in a flat inclined channel under weakly unstable conditions is determined using a method of averaging. The shallow water equations are first transformed to a normal form using a variation of parameters technique, and then a new method of averaging is defined. The method of averaging developed differs from averaging methods for ordinary differential equations and other partial differential equations by averaging along the characteristics of the linearized equations rather than along the time variable. It is found that the solution is dominated by the evolution of the solution along one characteristic. Both asymptotic and numerical results for periodic disturbances are presented. In addition, connections to dynamical systems concepts are described. © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. PACS: 47.20. − k; 47.35.+i; 47.60.+i MSC: 34E10 Keywords: Averaging; Instability; Shallow water
1. Introduction and linear instability Normalization and asymptotic averaging techniques are developed to analyze the semi-infinite boundary value problem of shallow water waves in an inclined constant-slope channel, where the evolution is in ∗ Tel.: +1 802 865 5245.
E-mail address: [email protected]. 0169-5983/$30.00 © 2006 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. doi:10.1016/j.fluiddyn.2006.02.006
470
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
the spatial variable. The system is assumed to be in a “weakly unstable” state, meaning the Froude number F is slightly greater than the neutrally stable value (Yu and Kevorkian, 1992). The methods developed here could also work in the initial value problem with minor modifications. The averaging technique created in this study differs from previous averaging methods in several respects. First, and most importantly, the averaging is performed over two variables, namely the characteristics of the reduced (linear) problem. It appears that in previous averaging applications, the averaging was performed on all equations along one variable, namely the “time” variable. An important exception should be noted. Chikwendu and Kevorkian (1972) developed a method of multiple scales for a class of nonlinear wave equations which in the final result is an average along characteristics. Hence, there is some precedent for this type of averaging. Second, the normalization process (where the equation(s) are transformed to Lagrange standard form, to be defined in Section 2) in many, though not all, averaging methods are analogous to the variation of parameters used for first order linear ordinary differential equations with exponential solutions. A variation of parameters idea is used in this study, but it is not exactly analogous to first order linear ordinary differential equations. In the shallow water equations, shallow means that the depth is small compared to some characteristic length, such as the wavelength. In addition, it is assumed the depth is small compared to the breadth of the conduit. Using dimensionless variables, the shallow water equations are then ht + uhx + hux = 0, ut + uux + hx = 1 −
(1) u2 . hF 2
(2)
See Kevorkian (1990, Section 5.1.1) or Spindler (2005, Chapter 1). Here h is the height of the free surface normal to the channel bottom and is normalized by the uniform height. Also, u is the flow speed parallel to the channel bottom averaged over the cross-section of the channel, and is normalized by the uniform flow. See Fig. 1. (The uniform solutions are the constant solutions of the dimensional shallow water equations essentially corresponding to equilibrium, smooth, stratified flow.) Note that this is an inviscid flow and yet has a quadratic drag law. This is due to the channel frictional effects modeled by a Chezý law. In the derivation of these equations, the Froude Number F is defined to be tan F := , C where is the angle of inclination and C is the coefficient of friction. It can also be written as Us2 , F= gH s cos Thus, the Froude number is of the where Hs and Us are the uniform depth and flow speed, respectively. order of the ratio of the uniform fluid velocity to the velocity gH s of the gravity wave. It represents a balance between inertia and gravity. In general, the one-dimensional shallow water equations have been analytically insolvable. The uniform solutions of (1) and (2) are hs := 1 and us := F . In dynamical systems terms, this is a rest, or fixed, point. The next obvious question to ask is: under what conditions is the system locally stable or unstable near this rest point? There is only one parameter, the Froude number, in this model and hence
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
471
1
Free Surface h (x,t)
u (x,t)
Bottom
θ X
Fig. 1. Flat inclined channel.
the question of stability focuses on the Froude number. Linear instability occurs, by definition, at values of the parameter(s) of a given system where the solution to the linearized problem near the rest point grows without bound. Whitham (1974, pp. 85–86) showed that the flow is linearly unstable for F > 2 and Spindler (2005) confirms this for the boundary value problem. Recall that the usual procedure is to substitute in traveling wave solutions ei(kx−t) with k real. Then the regions of instability are determined by conditions under which the imaginary part of is positive. Here, the same procedure was executed, but because it is a boundary value problem, is assumed real. Then the regions of instability are determined by conditions when the imaginary part of k is negative. The instability results say in effect that, for the linearized system, a transcritical bifurcation occurs at F = 2. Here the system’s (linear) stability at this rest point changes from stable to unstable as F increases through 2. In addition, the original system is assumed stable for F < 2 and unstable for F > 2 in some neighborhood of the rest point and in some neighborhood of F = 2. As described in slightly more detail in Section 3.1, this is similar to the statement in ordinary differential equation dynamical systems theory that the 1-jet of the system is sufficient to determine local stability near the rest point for F = 2. For the nonlocal situation, the nonlinear terms of the model come into play by modifying this instability in the far-field (large x). Weak instability, defined by 0 < F − 2>1, is the region of concern in this study. Assuming a small boundary disturbance to uniform flow (the inlet conditions) (1)
h(0, t) = 1 + hb (t), (1)
u(0, t) = F + ub (t),
(3) (4)
472
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488 x
Characteristic of slope one Q
x=t
P
Ω
Characteristics of slope one andt hree
t Boundary Values
Fig. 2. Domain of interest: .
the method of averaging is developed and applied to this weakly unstable situation. In dynamical system terms, we “start” at a point near the rest point where the above linear instability result applies and develop a method of averaging to observe how the system evolves. Thus, inlet conditions perturbed from the rest state are being specified. To model the weak nonlinearity, set F = 2 + where > 0 is “small”. Note that this ties together the physical description of the system (the Froude number) with the boundary condition. Before proceeding, the domain of interest is described. Recall that in a hyperbolic system the solution at a point is determined by the backward characteristics to the Cauchy data, which in√this case is√the boundary data. The slopes of the characteristics of the shallow water equations are u + h and u − h. To a first approximation, these are F + 1 and F − 1, since u = F and h = 1 are the uniform solutions. In the case here, F ≈ 2, so the characteristic slopes are given approximately by 3 and 1. The relevant domain for the boundary value problem is found by observing that the backward characteristics from the points in the domain := {(t, x) ∈ R2 : 0 < x < t} intersect the positive t axis on which the inlet conditions are specified. In Fig. 2, P is in the domain because the solution there is determined solely by the evolution from the boundary data. On the other hand, the backward characteristic of slope 1 from the points in the region {(t, x) ∈ R2 : 0 < t < x} intersect the positive x axis, on which data is not specified. Thus, to determine the solution at point Q in Fig. 2, initial data must also be specified. In this study, only boundary data is considered, and hence only the region is considered. 2. Introduction to averaging The method of averaging has been used to derive asymptotic evolution equations and to prove uniform convergence of asymptotic expansions over expanding intervals of order 1/ for ordinary differential
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
473
equations in Lagrange standard form (to be defined shortly). Sanders and Verhulst (1985) provide an excellent account of this theory, as does Murdock (1991). In this study, a method of averaging that applies to the shallow water equations is created. Krol (1991) developed a method for an advection–diffusion equation by averaging along the time variable. Buitelaar (1993) defines a time averaging method for systems of equations in abstract Banach spaces of the form dw + Aw = f (w, t, ), dt where A is the generator of a C0 group of unitary operators on a Hilbert Space. By Stone’s theorem, this implies iA is self-adjoint. In all of the applications considered by Buitelaar (e.g. nonlinear wave equations), this self-adjointness requires homogeneous boundary conditions. This study’s approach differs from that of Krol or Buitelaar in a couple of ways. First, the average is applied along the characteristics of the linearized equations, and second, the averaging is applied to a hyperbolic system. As noted above, Chikwendu and Kevorkian (1972) developed a method of multiple scales for a class of nonlinear wave equations which in the final result is an average along characteristics. The idea behind averaging methods is as follows. Assume an ordinary differential equation can be transformed into Lagrange standard form, which is defined to be an ordinary differential equation of the form dx = f (t, x, ), dt where these could be vector quantities and where f is real analytic in . In addition, initial conditions are specified at t = 0. In effect, this form implies that t is a slow variable. Now suppose we expand f in a Taylor’s series in to obtain dx = f (1) (t, x) + 2 f (2) (t, x, ). dt
(5)
The system is now in a normal form of order 1 (Murdock 2003, pp. 12–13). Usually f is required to be Lipschitz continuous in x. Also, often the existence of the average quantity 1 T (1) (1) f (x) := lim f (t, x) dt (6) T →∞ T 0 is imposed. In that case, Sanders and Verhulst (1985) call f a KBM-vectorfield (after Krylov, Bogoliubov, and Mitrolosky). An alternative assumption used is that f is periodic in time. To understand the mathematical origin of averaging, perform what is sometimes called a near-identity transformation x = x˜ + g(t, x) ˜
(7)
to transform Eq. (5) to the “autonomous” (in t) problem dx˜ = p(1) (x) ˜ + 2 p(2) (t, x, ˜ ). dt
(8)
474
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
What makes this work is for p(1) to be the average of f (1) holding x constant, namely Eq. (6). To see why, differentiate Eq. (7) with respect to t dx dx˜ jg = + x˜ g · x˜t + dt dt jt jg = p(1) (x) + O(2 ) ˜ + (9) jt using Eq. (8) and the fact that x˜t is O(). Now express xt in terms of x˜ instead of x by expanding the right-hand side of Eq. (5) around x˜ in a Taylor’s series. (Recall that x and x˜ differ by an order of .) Equate the order terms of Eqs. (5) and (9) to obtain f (1) (t, x) ˜ = p(1) (x) ˜ +
jg jt
.
Now average over t, 1 1 (g(t, x) ˜ − g(0, x)) ˜ = lim lim T →∞ T T →∞ T
T
f (1) (t, x) ˜ dt − p(1) (x). ˜
0
If g is uniformly bounded, the left-hand side is zero, and p(1) is found to be equal to the average of f (1) . Thus, dx˜ = f (1) (x) ˜ + 2 p(2) (t, x, ˜ ) dt is the new differential equation. In effect, the t dependence has been “pushed” up to an order of . It can be shown that the solution to the averaged system dz = f (1) (z) dt is an o(1) approximation to the solution to the original problem on expanding intervals of 1/, i.e. on intervals 0 t L/ for some constant L (Sanders and Verhulst, 1985). Murdock (2003) calls the method of using a near-identity transformation to simplify, or normalize, the right-hand side of Eq. (5) format 1a, or the direct/iterative format. To see how normal form theory works in this context, append a zeroth component to x and f (i) as follows: x0 = t, (i)
f0 = 0,
i 1
and define a zeroth order column vector f (0) = [1 0 · · · 0]T so now Eq. (5) can be written as dx = f (0) + f (1) (t, x) + 2 f (2) (t, x, ). dt Define the homological operator Lu := Duf (0) − Df (0) u,
(10)
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
475
where u is a vector field as a function of the vector x and D indicates the total derivative with respect to x. (Thus, these derivatives are matrices.) Then, for format 1a, Eq. (10) transforms to dx˜ = f (0) + [f (1) (t, x) ˜ − (Lg)(t, x)] ˜ + 2 f (2) (t, x, ˜ ). dt From the perspective of normal form theory, the idea then is to choose p(1) so that f (1) − p(1) is in the range of L (which then determines g) and in such a way as to make the problem simpler. The way in which this is done is called the style. In the case here, Lg = jg/jx0 = jg/jt. For the style, p(1) is chosen to be the projection of f (1) into the kernel of L, i.e. p (1) is the time average of f (1) as found above. The approach in this study is similar to the procedure outline above. Note that the progression through the above derivation showed how the method of averaging comes about. If one knows how to average already, then just define t g(t, z) = [f (1) (s, z) − f (1) (z)] ds 0
in the near-identity transformation and work with that (Verhulst, 1996). This also arises naturally from the normal form theory discussed above. In addition, note that if := t (slow time), then the averaged system is dz = f (1) (z), d which indicates the connection to multiple scales. Finally, if f (1) is periodic with period T, the same method is used except the integration is over a period. (In that case g turns out to be periodic.) 3. The method of averaging applied to the weakly unstable shallow water equations A method of averaging for the weakly unstable boundary value shallow water system is now developed, applying similar reasoning as used for ordinary differential equations described above. In dynamical system problems for ordinary differential equations, the systems are often assumed to be in Lagrange standard form. Transforming a system of partial differential equations into this form is not trivial, and that is true in this case too. One common method of attack, mentioned earlier and applied in the next section, is to use a variation of parameters procedure. In that approach, the reduced (linear) problem is solved, which results in some “constants” of integration. Then the full problem is transformed by letting these “constants” of integration depend on all independent variables. Once this is finished, the method of averaging can be defined. Note that all of the calculations related to this work can be found in Spindler (2005, Chapter 3), and also are available by contacting the author. 3.1. The reduced problem and Lagrange standard form The first task is to transform the shallow water equations (1) and (2) into Lagrange standard form. To do this, the reduced problem is solved first. Center the shallow water equations (1) and (2) around the
476
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
uniform values by the trivial change of variables h = 1 + h(1) , u = F + u(1) ,
(11)
which define h(1) and u(1) for = 0, and substitute into the shallow water equations. (This is not a perturbation expansion.) Divide by , multiply the second equation through by 1 + h(1) , and simplify the right-hand side of the second equation to arrive at (1)
+ F hx(1) + ux(1) = −{u(1) hx(1) + h(1) ux(1) }, 1 2 (1) (1) ut + hx(1) + F ux(1) − h(1) − u(1) = − 2 u(1)2 − [u(1) ux(1) + h(1) ut + h(1) hx(1) F F (1) (1) (1) (1) (1) + F h ux ] − h u ux .
ht
(12)
Of course, this change of variables centers the rest point at the origin, but it also in effect “lowers” the order by 1 compared to the original shallow water equations. Note the time derivative on the right-hand side of the second equation. Before inserting F = 2 + , the reduced problem is determined and solved. Set = 0, so that Eqs. (12) reduce to the linearized shallow water equations with F = 2. Define hˆ (1) and uˆ (1) to be the solutions of these linearized equations. Add and subtract the two equations (12) (with = 0) and define Sˆ := hˆ (1) + uˆ (1)
and
Rˆ := hˆ (1) − uˆ (1)
to obtain Sˆt + 3Sˆx − Rˆ = 0, Rˆ t + Rˆ x + Rˆ = 0. The characteristic slopes of this system are 1 and 3 and so the characteristics of this system are x − t= constant and x − 3t= constant. Change variables along these characteristics, that is (x, t) := x − t
and
(x, t) := x − 3t,
(13)
which transforms the first order to ˆ , ), Sˆ (, ) = 21 R(
(14)
ˆ , ). Rˆ (, ) = 21 R(
(15)
In effect, the linearized system has been transformed to normal form in the sense of hyperbolic partial differential equations (Forsythe and Wasow, 1960). The solution of Eq. (15) is found to be ˆ , ) = rˆ ()e/2 , R(
(16)
where rˆ () is in effect a “constant” of integration. Using this in the Eq. (14) above gives ˆ , ) = R ˆ ()e/2 + sˆ (), S(
(17)
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
477
Ω
= 3 Fig. 3. Domain of interest: .
where sˆ () is in effect another “constant” of integration and 1 ˆ () := R rˆ ( ) d . 2
(18)
Note that the primes here do not designate differentiation. The variable is just a dummy integration variable. The domain will be needed in the new (, ) coordinate system. When x = 0, Eqs. (13) are = −t and = −3t. Eliminating t determines the boundary = 3. Similarly, x = t implies = 0. Also, x 0 maps to the region 3, and x < t maps to the region < 0 under the transformation (13). (The boundary t 0 maps to the region which is less restrictive than 3.) Thus, = {(, ) : 3 0}. See Fig. 3. Note that using these transformations, the concept of jet sufficiency of ordinary differential equations can be related to the situation here. Applying the change of variables S := h(1) + u(1) , R := h(1) − u(1) , := x − (F − 1)t, and := x − (F + 1)t, Eqs. (12) are then transformed to a normal form
1 1 1 1 1 − S+ + R + m1 , S = 2 2 F 2 F
1 1 1 1 1 R = − S+ + R + m2 . (19) 2 2 F 2 F The remainder terms m1 and m2 terms are quadratic polynomials in S, R, and their derivatives. Because only linear transformations have been applied to the system, its stability properties are not changed. In ordinary differential equation dynamical systems theory, a system in the form of (19) where the right-hand side is expanded into a Taylor’s series in the dependent variables with F = 2, is 1-jet sufficient since the linear polynomial terms determine the stability locally around the rest point. If F = 2, the quadratic terms
478
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
would need to be examined to determine if the system is 2-jet sufficient, and so on. See Murdock (2003, p. 8) for a definition and discussion of sufficiency of jets in the context of ordinary differential equations. The reason for determining the solution to the reduced problem is to transform the original problem (12) using this reduced solution by the method of variation of parameters. First, to determine the full dependence of the full equations on , substitute F = 2 + into the second equation of (12) and multiply both sides by F 2 to obtain (1)
4(ut
(1)
+ hx(1) ) + 8ux(1) − 4(h(1) − u(1) ) + 4(ut (1)
= {−u(1)2 − 4[u(1) ux(1) + h(1) ut
+ hx(1) ) + (12ux(1) − 2h(1) ) − 2(h(1) − u(1) )
+ h(1) hx(1) + 2h(1) ux(1) ]} + O(2 ).
All 2 terms are summarized in O(2 ). (Note that here the O(2 ) refers to the terms in this equation, of course, which would imply O(3 ) in the original shallow water equations.) Moving terms to the right-hand side, the equations become (1)
+ 2hx(1) + ux(1) = −(u(1) hx(1) + h(1) ux(1) + hx(1) ),
(1)
+ hx(1) + 2ux(1) − (h(1) − u(1) ) = {−(1 + h(1) )(ut + hx(1) + 2ux(1) ) − (1 + u(1) )ux(1) + h(1) − 21 u(1) − 41 u(1)2 } + O(2 ).
ht ut
(1)
(20)
Observe that (20) is recursive in the sense that part of the left-hand side is in the right-hand side, namely (1) (1) (1) ut + hx + 2ux . Substitute this in for itself on the right-hand side, moving higher order terms into O(2 ) to obtain (1)
+ 2hx(1) + ux(1) = −(u(1) hx(1) + h(1) ux(1) + hx(1) ),
(1)
+ hx(1) + 2ux(1) − (h(1) − u(1) ) = {−(1 + u(1) )ux(1) + 21 u(1) − (h(1) − 21 u(1) )2 }
ht ut
+ O(2 ).
(21)
(22)
As in the reduced problem, transform the left-hand side of (21) and (22) into normal form (in the sense of hyperbolic partial differential equations), so let S := h(1) + u(1)
and
R := h(1) − u(1)
and eliminate h(1) and u(1) in the equations. (Add and subtract the equations). Again, transform along the characteristics of the linearized equations, (x, t) = x − t and (x, t) = x − 3t, to obtain S − 21 R = 21 {−(1 + 43 S − 41 R)(S + S ) + 41 (S + R)(R + R ) + 41 (S − R) −
1 16 (S
+ 3R)2 } + O(2 ),
(23)
R − 21 R = 21 {(1 + 41 S − 43 R)(R + R ) + 41 (S + R)(S + S ) + 41 (S − R) −
1 16 (S
+ 3R)2 } + O(2 ).
(24)
Next, change variables using the reduced solutions, Eqs. (16) and (17), by the method of variation of parameters. Recall that sˆ and rˆ are “constants” of integration for the reduced equations with the first
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
479
depending on and the second on . Now let s depend also on and r also on . Therefore, define r(, ) := R(, )e−/2 , s(, ) := S(, ) − R(, )e/2 , where R is similar to Eq. (18): 1 R(, ) = r( , ) d . 2
(25)
Determine the derivatives of R and S and substitute these into the equations for S and R , (23) and (24), to obtain s = 21 {−[1 + 43 (s + Re/2 ) − 41 re/2 ](s + 21 re/2 + s + R e/2 + 21 Re/2 )
+ 41 (s + Re/2 + re/2 )(r e/2 + r e/2 + 21 re/2 ) + 41 (s + Re/2 − re/2 )
−
1 16 (s
+ Re/2 + 3re/2 )2 } + O(2 )
and r = 21 {(1 + 41 s + Re/2 − 43 re/2 )(r e/2 + r e/2 + 21 re/2 )
+ 41 (s + Re/2 + re/2 )(s + 21 re/2 + s + R e/2 + 21 Re/2 ) + 41 (s + Re/2 − re/2 )
−
1 16 (s
+ Re/2 + 3re/2 )2 }e−/2 + e−/2 O(2 ).
Re-arrange and combine terms. Also, because s and r are equal to times terms on the right, move s , r , and R terms on the right up into O(2 ) and e−/2 O(2 ) to arrive at
1 3 1 2 s = − 1 + s s + s − s 2 4 4 16
1 1 1 5 1 3 + − (3R − r)s + − r − R + r s + − R − r e/2 4 8 2 4 4 4
5 2 1 7 2 1 + − r − r R − R + r (r + R) e + O(2 ) (26) 16 2 16 4 and
1 1 2 −/2 1 1 1 1 1 1 ss + s − s e r − r s + (r + R)s r = + r + r + R + 2 4 4 16 4 4 4 8 4
13 2 1 2 1 3 + − r + R + R − r r e/2 + e−/2 O(2 ). (27) 16 16 4 4
Basically these independent variables introduced into s and r correspond to slowly varying variables for their respective function, meaning that corresponds to a slowly varying independent variable for s and corresponds to a slowly varying independent variable for r. Observe that everything has been exact so far. (Although the higher order terms are hidden in the O(2 ) term.) Various transformations and algebraic operations have just been applied to the problem. Note that for = 0 the reduced solutions are recovered. Also observe that, at least up to order O(), s = 0 and r = 0 are rest points as they should be, corresponding to h(1) = 0 and u(1) = 0.
480
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
3.2. Near-identity transformation The above system of equations, (26) and (27), is now expressed in normal form in the sense of partial differential equations. This is similar to the normal form as used in dynamical systems for ordinary differential equations (Murdock, 2003), but there are differences. For example, on the left-hand side, the derivatives are in terms of two independent variables instead of one, of course, since these are partial differential equations. Also, derivatives appear on the right-hand side. Despite these differences, an averaging theory can be developed along the lines as described in Section 2. This will also be analogous to the transformation of normal forms (Murdock, 2003). Now that the equations are in normal form, the method of averaging along characteristics arises formally by applying the idea of a near-identity transformation defined by
g1 (, s˜ , s˜ , r˜ , r˜ ) s˜ s = + . (28) r r˜ 2 e−/2 g2 (, s˜ , s˜ , r˜ , r˜ ) The goal is to define this transformation such that the Lagrange standard form equations (26) and (27) are “autonomous” and uncoupled from each other
s˜ = p1 (˜s , s˜ ) + O(2 ), 2
r˜ = p2 (˜r , r˜ ) + e−/2 O(2 ) 2
(29) (30)
solving the same boundary conditions, namely s˜ (, 3) = s(, 3) and r˜ (, 3) = r(, 3). This is similar to the example in Section 2 for ordinary differential equations where the independent variable t was eliminated. The difference here is that dependent variables are being eliminated from the right-hand side functions instead of an independent variable. For convenience, a factor e−/2 of g2 is inserted in the second equation. Assume that gi functions and their derivatives are uniformly bounded. Also, by definition of change of variables, the near-identity transformation must be invertible in a neighborhood of =0 which is assumed. To proceed, differentiate (28) to obtain
s s˜ = r r˜
g1 + g1˜s s˜ + g1˜s s˜ + g1˜r r˜ + g1˜r r˜ + . (31) 2 e−/2 g2 − 21 e−/2 g2 + e−/2 (g2˜s s˜ + g2˜s s˜ + g2˜r r˜ + g2˜r r˜ + g2 ) Recall that s˜ = O() and r˜ = e−/2 O(). Thus, those terms are moved up into higher orders. Also assume the same is true of higher derivatives of these terms. Then
g1 + g1˜r r˜ + g1˜r r˜ s O(2 ) s˜ + = + e−/2 O(2 ) r r˜ 2 e−/2 g2 − 21 e−/2 g2 + e−/2 (g2˜s s˜ + g2˜s s˜ + g2 ) ⎤ ⎡ d g1
2) ⎥ d s˜,˜s constant ⎢ s˜ O ( ⎥ + −/2 2 . = + ⎢ (32) ⎦ e r˜ O ( ) 2 ⎣ d −/2 (e g2 ) d r˜ ,˜r constant
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
481
Observe that by the near-identity transformation (28), s and s˜ differ by O() and r and r˜ differ by e−/2 O(). Expand the right-hand side terms of Eqs. (26) and (27) (as functions of s and r and their derivatives) in a Taylor’s series around s˜ and r˜ and their derivatives, to simply replace s with s˜ and r with r˜ , with the remainder terms of the Taylor’s Series moved up into the O(2 ) and e−/2 O(2 ) terms. (Note that to apply this Taylor Series, g1 and g2 and their derivatives must be uniformly bounded.) Next substitute Eqs. (26) and (27) (with s˜ and r˜ , etc., instead due to the discussion in the last paragraph) into the left-hand side of Eq. (32), and Eqs. (29) and (30) into the right-hand side of Eq. (32). After cancelling an from both sides and re-arranging, d p1 (˜s , s˜ ) + g1 d s˜,˜s constant
1 3 1 = − 1 + s˜ s˜ + s˜ − s˜ 2 4 4 16
1 1 5 1 3 1 ˜ ˜ ˜ + − (3R − r˜ )˜s + − r˜ − R + r˜ s˜ + − R − r˜ e/2 4 8 2 4 4 4
5 1 7 2 1 ˜ ) e + O(), ˜ − R ˜ + r˜ (˜r + R + − r˜ 2 − r˜ R (33) 16 2 16 4 d −/2 p2 (˜r , r˜ ) + (e g2 ) d r˜ ,˜r constant
1 2 −/2 1 1 s˜ s˜ + s˜ − s˜ e = 4 4 16
1 1 1 1 1 ˜ + ˜ )˜s r˜ − r˜ s˜ + (˜r + R + r˜ + r˜ + R 4 4 4 8 4
1 2 3 1 13 2 /2 ˜ ˜ R − r˜ r˜ e + − r˜ + R + 16 16 4 4 + e−/2 O().
(34)
˜ (Note that the cancellation of lowers the order 1 more relative to the shallow water equations.) Here, R is defined similar to R except using r˜ instead of r in the integral defining R. Examine Eq. (33) first. Move p1 to the other side of Eq. (33) and integrate as follows:
1 2 3 1 g1 = − 1 + s˜ s˜ + s˜ − s˜ 4 4 16
1 1 1 5 1 3 /2 ˜ ˜ ˜ − (3R − r˜ )˜s + − r˜ − R + r˜ s˜ + − R − r˜ +e 4 8 2 4 4 4
5 1 7 2 1 ˜ ) − p1 (˜s , s˜ ) d ˜ − R ˜ + r˜ (˜r + R +e − r˜ 2 − r˜ R 16 2 16 4 s˜ ,˜s constant
+ O()
482
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
or
1 2 1 3 1 /2 ˜ − r˜ ) d − s˜ (3R g1 = − 1 + s˜ s˜ + s˜ − s˜ + e 4 4 16 4
1 5 1 1 3 ˜ + r˜ d + ˜ − r˜ d − r˜ − R − R +˜s 8 2 4 4 4
5 1 7 2 1 ˜ − R ˜ + r˜ (˜r + R ˜ ) d − r˜ 2 − r˜ R + e 16 2 16 4 − p1 (˜s , s˜ ) + c1 (, s˜ , s˜ ) + O(). Next divide both sides by and let → −∞ and continue to treat s˜ and s˜ as constant. This limit could conceivably be taken in a number of ways, as long as t > x 0, that is 3 < 0 (to remain in the domain as discussed in the first section). However, another requirement for choosing the curve along which to take the limit is that singularities (shocks) must not be encountered in the limit. The boundary data curve x = 0 certainly satisfies this requirement, and hence the limit is defined along this curve. Since g1 is assumed uniformly bounded, the left-hand side tends to zero. Suppose also there exists a positive constant C(x) 0 for each fixed x such that ˜ |, |R ˜ | C(x)e−3/2 |˜r |, |˜r |, |R
(35)
uniformly in t. (It will be seen in Section 4 that the solution is consistent with this assumption.) By this assumption, each of the terms above, except the first and the p1 and c1 terms, will be of the form e1/2(−3)
=
e−x
or
e(−3)
=
e−2x
,
which approaches zero uniformly as → −∞ since x is fixed. Assuming the integration “constant” c1 is uniformly bounded, all that remains is p1 (˜s , s˜ ) = −(1 + 43 s˜ )˜s + 41 s˜ −
1 2 16 s˜
+ O().
Substitute this into (29) to obtain one evolution equation. Process Eq. (34) in a similar manner. First, integrate (34) holding r˜ and r˜ constant,
1 2 − /2 1 1 1 1 ˜ g2 = e s˜ s˜ + s˜ − s˜ e + r˜ + r˜ + R 4 4 16 4 4
1 1 1 ˜ )˜s r˜ − r˜ s˜ + (˜r + R + 4 8 4
1 2 1 3 13 2 /2 ˜ ˜ R − r˜ r˜ e − p2 (˜r , r˜ ) d + − r˜ + R + 16 4 4 16 r˜ ,˜r constant /2
+ e/2 c2 (, r˜ , r˜ ) + O().
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
483
Then re-write the equation as
1 2 − /2 1 /2 1 s˜ s˜ + s˜ − s˜ e e g2 = 4 4 16
1 ˜ 1 ˜ 1 ˜ 1 1 ˜ ˜ ˜ ˜ R − R s˜ + (R + RR )˜s + R + R + RR + 4 4 4 8 4
−)/2 1 2 1 3 ˜ ˜ 13 ˜ 2 ( ˜ R˜ ) d ˜R + ˜ R − R R e R − p2 (R, + − R + R ˜ ˜ 16 16 4 4 R,R constant + e/2 c2 (, r˜ , r˜ ) + O() ˜ . (In the first equation above, r˜ and r˜ have been held constant ˜ R := e/2 R using R˜ = e/2 r˜ and where R in the integral, which is equivalent to holding R˜ = e/2 r˜ and R˜ = e/2 r˜ constant.) Also, it has been assumed that p2 is homogeneous in its arguments (meaning yp 2 (x1 , x2 ) = p2 (yx 1 , yx 2 )). Factor the R˜ terms out of the integrals to obtain
1 2 − /2 1 1 1 ˜ 1 /2 ˜ ˜ g2 = e s˜ s˜ + s˜ − s˜ e d + R + R + RR 4 4 16 4 4 1 ˜ 1 1 ˜ R )˜s R − R˜ + s˜ d + (R˜ + R 4 8 4
1 2 3 ˜ ˜ 1 13 ˜ 2 ˜ R˜ ) + e/2 c2 (, r, r ) + O(). ˜ ˜ RR − R R − p2 (R, + 2 − R + RR + 16 16 4 4 Divide both sides by , treat R˜ and R˜ as constant, and let → −∞. Again, choose the curve x = 0 along which to take this limit. Since g2 is assumed uniformly bounded, the left-hand side tends to zero. Also suppose that the various integrals involving s and s in this last equation are uniformly bounded as well as the integration “constant” e/2 c2 (, r˜ , r˜ ). In effect the solutions have been restricted to the class of functions satisfying these specifications. Then ˜ R˜ ) = [R˜ + 1 R˜ + 1 R ˜ p2 (R, 4 4 R ] + O() or ˜ ] + e−/2 O(). p2 (˜r , r˜ ) = [˜r + 41 r˜ + 41 R This result provides a second evolution equation by substituting into Eq. (30). Also observe that p2 , at least up to order , is homogeneous and so that assumption is consistent. From the work above, the averaging operation that appears to transform the equations as desired is to average Eq. (26) by 1 lim [· · ·] d →−∞ s,s constant
484
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
along the line x = 0, where s and s are treated as constants in the limit. Average Eq. (27) by e/2 −/2 e lim [· · ·] d →−∞ R,R constant
along the line x = 0, where R, R , and RR are treated as constants in the limit. The near-identity transformations gi are 1 1 5 1 /2 ˜ + r˜ d + ˜ − r˜ ) d + s˜ − s˜ − r˜ − R (3R g1 = e 4 8 2 4
5 1 7 1 ˜ − R ˜ 2 + r˜ (˜r + R ˜ ) d + e − r˜ 2 − r˜ R 16 2 16 4 and
1 1 1 ˜ 1 ˜ 1 2 − /2 d + s˜ s˜ + s˜ − s˜ e R − R g2 = e 4 4 16 4 8
13 1 2 3 1 ˜R + ˜ R − R˜ R˜ . + 2 − R˜ 2 + R R 16 16 4 4 /2
1 3 ˜ − r˜ d − R 4 4
1 ˜ R )˜s s˜ d + (R˜ + R 4
These are uniformly bounded if the assumptions on s˜ and r˜ hold, as required. To summarize, the equations in Lagrange standard form (26) and (27) transform to the evolution equations
3 1 2 1 − 1 + s˜ s˜ + s˜ − s˜ + O(2 ) (36) s˜ = 2 4 4 16 and
1 1 ˜ r˜ + r˜ + R + e−/2 O(2 ) r˜ = 2 4 4
(37)
under the near-identity transformation. In ordinary differential equations, the system of equations retaining the higher order O(2 ) terms is sometimes called the transformed system. Of course, the O(2 ) terms are dropped to approximate the actual solution. When these higher order terms are dropped,
3 1 2 1 − 1 + s˜ s˜ + s˜ − s˜ (38) s˜ = 2 4 4 16 and r˜ =
1 1 ˜ r˜ + r˜ + R 2 4 4
(39)
are sometimes called the averaged system. If ¯ = is defined, then the averaged Eq. (38) becomes s˜¯ = 21 [−(1 + 43 s˜ )˜s + 41 s˜ −
1 2 16 s˜ ]
and, if ¯ = is defined, then the second averaged (39) becomes ˜ ]. r˜¯ = 21 [r˜ + 41 r˜ + 41 R
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
485
In effect, the new variables define far-field scales for each equation. In ordinary differential equations, this system of equations is sometimes called the guiding system (Murdock, 1991). The development of the averaging method can be placed in the context of normal forms theory of dynamical systems. Here, the homological operator is defined to be ⎡ ⎤ d g1 ⎢ ⎥ d s˜,˜s constant ⎥. Lg := ⎢ ⎣ d −/2 ⎦ (e g2 ) d r˜ ,˜r constant Let f (1) be the order term on the right-hand side of Lagrange standard form equations (26) and (27). Eqs. (33) and (34) result from transforming the normal form system (26) and (27) into format 1a by applying the near-identity transformation. The abstract form of Eqs. (33) and (34) is Lg = f (1) − p (1) .
The average operation and the final evolution equations then arise by defining p (1) to be the projection of f (1) into the kernel of L. 4. A comparison The solutions to the averaged evolution equations (38) and (39) are determined, assuming a sine wave boundary condition, namely h(0, t) = 1 + sin(t)
and
u(0, t) = F .
It is hoped that the solutions to the averaged equations are a good approximation to the solution of the transformed Eqs. (36) and (37). This means that the solutions are within o(1) over expanding intervals of x of order 1/ as in the ordinary differential equations case. Also, recall that s= s˜ + O() and r = r˜ +e−/2 O(). Hence R = R˜ + O() (where R˜ := e/2 r˜ ) and S = S˜ + O(). Thus, up to order , R and S, and hence h(1) ˜ In this section, the solution of Eq. (39) is solved explicitly for r˜ , and u(1) , are determined by R˜ and S. while Eq. (38) is solved numerically. Finally, the shallow water equations are numerically integrated and compared to the solution of the evolution equations, at a perturbation of = 0.1. The averaged r-evolution equation (39) is linear in r˜ with constant coefficients, so it can be solved analytically. The boundary condition for r˜ can be written r˜ (, 3) = − sin()e−3/2 ˜ Assume a solution of the form using the definitions of r˜ and R. r˜ (, ) = −e−(3/2)−(/2)a(3−) sin + b(3 − ) , 2
(40)
where the parameters a and b are to be determined. Substitute the proposed solution into the r˜ -evolution equation (39) to solve for a and b to obtain 5 3/4 a=− − 4 9 + 4 2
and
b=
/2
9 + 4 2
− .
(41)
486
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488 Direct Average
0.08 0.06 0.04
h(1)
0.02 0 -0.02 -0.04 -0.06 -0.08 37
38
39
40 Time t
41
42
43
Fig. 4. Height h1 of wave at x = 30.
Note that a is less than zero. Also observe that r˜ is consistent with assumption (35). Transforming back to the R˜ variable ˜ , ) = − e−(1/2)(1+a)(3−) sin + b(3 − ) R( 2 −(1+a)x = −e sin + b(3 − ) , 2 where of course R˜ = e/2 r˜ . This implies that the part of the solution from the characteristic decays exponentially in the x direction as long as < 1/|a|. Thus, with this restriction on , the contribution of R˜ to h(1) and u(1) is negligible as x becomes large. The far-field evolution arises from the contribution ˜ contribution to S˜ decays exponentially in the x of S˜ which is essentially determined by s˜ (since the R (1) (1) direction also). This also implies that h and u become nearly identical rather quickly. A basic forward–backward (FB) method is used to numerically integrate the s-evolution equation. In addition, the differences are designed to satisfy the CFL condition. The boundary condition in terms of s˜ is /2 ˜ s˜ , = − sin −e R , . 3 3 3 From Eq. (25) for R and the solution for r˜ 1 s˜ 3 sin , = − sin − + 2 cos . 2 3 3 9 + 4 3 3
(42)
Eq. (42) is the boundary condition for the s-evolution equation (38) used in the numerical integration. The shallow water equations are numerically integrated by applying the artificial viscosity method (Lapidus and Pinder, 1999, Section 6.2.3).Again, the differences are designed to satisfy the CFL condition.
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
487
Direct Average
0.08 0.06 0.04
h(1)
0.02 0 -0.02 -0.04 -0.06 -0.08 0
5
10
15 20 Distance x
25
30
35
Fig. 5. Height h1 of wave at t = 30.
Finally, all of the programs were written in MatLab 6.5. Fig. 4 is an overlay of the two solutions at x = 30 and = 0.1. Observe that the waves are periodic in t as expected, and the formation of a bore occurs. The waves are within about 2 of each other up to (and past) x ≈ 1/ as expected from perturbation theory. Fig. 5 is an overlay of the two solutions at t = 30 and = 0.1. Notice the gradual exponential growth of the wave in the x direction, as the linear stability theory predicted. Also observe that the waves steepen to a bore as x increases. Finally, the solutions from the two methods at t = 30 are still within about 2 of each other up to (and past) x ≈ 1/. 5. Conclusion The averaging method introduced in this study for the shallow water equations indicates a method of averaging for more general two-dimensional systems of hyperbolic partial differential equations. By averaging along the characteristics of the linearized equations, an approximate solution may be found in the far-field. It is hoped that this method can be extended to that more general situation. Note that Yu and Kevorkian (1992) analyzed the initial value problem for a flat inclined constant-slope channel using a multiple scales method, while Spindler and Yu (2006) analyzed the boundary value problem using multiple scales in conjunction with the Fredholm alternative theorem. The s-evolution equations found for mutliple scales there and for averaging here differ by one term, but the numerical solutions of these equations are both within order 2 of each other. These results suggest that the two solutions may be asymptotically equivalent, although this question of equivalence merits further study. For ordinary differential equations, Murdock (1991) shows that the first order two-scales method (using one far-field scale) is identical to averaging. In fact, Perko (1969) shows that the two methods for ordinary differential equations are asymptotically equivalent, though not necessarily identical, up to all orders of .
488
R. Spindler / Fluid Dynamics Research 38 (2006) 469 – 488
Although deriving the averaging technique is more involved than applying multiple scales, once the averaging operation is known, it is no more difficult to apply than determining the secular terms in multiple scales. However, the asymptotic validity of multiple scales for partial differential equations is an open mathematical question as well as its relationship to averaging. The advantage of the normalization and averaging technique is that it may provide a more tractable framework for deriving rigorous error estimates for both techniques applied to partial differential equations or systems of partial differential equations, just as with ordinary differential equations. Acknowledgments I would like to thank Ferdinand Verhulst for his helpful comments on this project. I am also grateful for the many suggestions the reviewers made to improve the manuscript. References Buitelaar, R.P., 1993. The method of averaging in Banach spaces: theory and applications. Ph.D. Thesis, Rijksuniversiteit Utrecht, Netherlands. Chikwendu, S.C., Kevorkian, J., 1972. A perturbation method for hyperbolic equations with small nonlinearies. SIAM J. Appl. Math. 22 (2), 235–258. Forsythe, G.E., Wasow, W.R., 1960. Finite-Difference Methods for Partial Differential Equations. Wiley, New York. Kevorkian, J., 1990. Partial Differential Equations: Analytical Solution Techniques. Wadsworth and Brooks/Cole, Pacific Grove, CA. Krol, M.S., 1991. Nearly time-periodic advection–diffusion problems. SIAM J. Appl. Math. 51 (6), 1622–1637. Lapidus, L., Pinder, G., 1999. Numerical Solution of Partial Differential Equations in Science and Engineering. WileyInterscience, New York. Murdock, J.A., 1991. Perturbation Theory and Methods. Wiley-Interscience, New York. Murdock, J., 2003. Normal Forms and Unfoldings for Local Dynamical Systems. Springer, New York. Perko, L.M., 1969. Higher order averaging and related methods for perturbed periodic and quasi-periodic systems. SIAM J. Appl. Math. 17, 698–724. Sanders, J.A., Verhulst, F., 1985. Averaging Methods in Nonlinear Dynamical Systems. Springer, New York. Spindler, R., 2005. The one-dimensional shallow water boundary value problem and instability. Ph.D. Thesis, University of Vermont, Burlington, VT. Spindler, R., Yu, J., 2006. Nonlinear spatial evolution of small disturbances from boundary conditions in flat inclined-channel flow. Stud. Appl. Math. 116, 173–200. Verhulst, F., 1996. Nonlinear Differential Equations and Dynamical Systems, second ed. Springer, Berlin. Whitham, G.B., 1974. Linear and Nonlinear Waves. Wiley, New York. Yu, J., Kevorkian, J., 1992. Nonlinear evolution of small disturbances into roll waves in an inclined open channel. J. Fluid Mech. 243, 575–594.