Transient Numerical Solutions of an Extended Korteweg-De Vries Equation Describing Solitary Waves in Open-Channel Flow

A solitary wave in two-dimensional, incompressible, turbulent free-surface flow over a plane bottom with small, constant slope is considered. The flow is assumed to be slightly supercritical with Froude numbers close to 1. If the flow far upstream and far downstream is fully developed, a simple argument based on the law of momentum shows that for a solitary wave to exist, the bottom friction cannot be constant all along the channel bed. In [1] the situation was considered where the bottom roughness of the channel is constant over some distance and slightly higher than in the rest of the channel bed, giving rise to a higher bottom friction coefficient. In an asymptotic analysis in [1] an extended Korteweg-de Vries (KdV) equation was derived to describe the surface elevation of the fluid. Adopting this equation, we solved it numerically by posing a coupled boundary-value eigenvalue problem and obtained results for stationary and transient wave solutions as well as for the eigenvalue, which corresponds to distinct values of the bottom friction coefficient. While the numerical solutions as compared to the asymptotic solutions agree qualitatively in the stationary case, there were major differences found in case of the transient solutions. Preliminary studies to this work were reported in [2].


Introduction
A solitary wave in two-dimensional, incompressible, turbulent free-surface flow over a plane bottom with small, constant slope α is considered, see Fig.If the flow far upstream and far downstream is fully developed, a simple argument based on the law of momentum shows that for a solitary wave to exist, the bottom friction cannot be constant all along the channel bed.Since the momentum flow rate far upsteam and far downstream is the same, the forces acting on a control volume of fluid must balance.Forces arising from the hydrostatic pressure cancel, too.If surface tension and surface shear forces are neglected, the increased weight of the fluid under the solitary wave must therefore be balanced by increased bottom friction forces.
Following that argument, a situation as in Fig. 1 is considered, where the bottom roughness of the channel is assumed to be constant over some distance  and slightly higher than in the rest of the channel bed, giving rise to a higher bottom friction coefficient.The flow is assumed to be slightly supercritical with Froude numbers close to 1.
In an analysis given by Schneider [1], a coupled asymptotic expansion of the Reynolds-averaged Navier-Stokes equations was performed for differences of the Froude number to 1, 1 Institute of Mechanics, Montanuniversitaet Leoben, Franz-Josef-Strasse 18, A-8700 Leoben, Austria.
(Fr −1) ≡ ε → 0 , the Reynolds number Re → ∞ , and the slope of the channel bottom α → 0 , with the remarkable feature that any recourse to turbulence modelling was avoided.r h and r u were defined to denote the reference quantities of Reynolds-averaged fluid elevation and volumetric mean velocity, respectively, in fully developed flow.The main result of the asymptotic analysis was an extended KdV equation describing the surface elevation h H h r = + ( ) and reads as follows: The term extending the classical KdV equation includes a parameter β , which is small within the assumptions of the asymptotic expansion and which characterizes the effect of turbulent dissipation.The function Γ(X) describes the increased bottom friction due to the increased roughness over the bottom length  and was assumed to be given as where L h r = 3 ε  was taken to be 1 in [1].The eigenvalue λ determines the amount of increase in bottom roughness as necessary for the fluid to fulfill the conditions of the law of momentum.Equation (1), together with the obvious boundary conditions for a solitary wave was solved in [1] again by asymptotic methods for small values of β .

Stationary Solutions
First, the stationary form of the extended KdV, Eq. ( 1) with H 1,T ≡ 0 was considered.The results for this case were already presented in a previous work, see [3].
For the friction-free case β = 0 and Γ ≡ 0 , a well-known analytical solution of Eq. ( 1) is the solitary wave solution with its crest located at the undetermined position X m From the asymptotic theory in [1] it follows that to lowest order in the expansion parameter β the solution of Eq. ( 1) is given by H H (1) starting close to the corresponding stationary solution will stay close to or rather diverge, respectively, from the stable or unstable stationary solution in the limit of large times.Numerical solutions of Eq. ( 1) may be obtained by formulating a coupled boundary-value eigenvalue problem with respect to the space variable X .In addition to the two boundary conditions in Eq. ( 4), Eq. ( 1) needs to be supplemented by two further boundary conditions, which may be posed as to obtain the correct decay behaviour of H 1 for X → ±∞ .Assuming H 1 to decay as H 1 ~ exp(kX) , Eq. ( 1) gives with three roots which are real for β ≤ β max ≈ 0.383 .For small values of β these roots are In the numerical solutions presented in the following, the two additional boundary conditions used were with the corresponding exact solutions of Eq. ( 6).The choice of k 1 instead of k 2 in Eq. ( 9) for X → +∞ gives a "fast" decay behaviour as compared to the choice of k 2 and was motivated by the results of the numerical integration, which was carried out for both choices and revealed rather large gradients close to the right end of the integration domain in case of the choice of k 2 , which was therefore dropped in the following considerations.The numerical solutions of the full Eq.( 1) for the stationary case were obtained with the package MATLAB R2010A using the bvp4c integrator.Typical values chosen were for the stepsize Δ X ≈ 2 • 10 −3 and relative error tolerances of 10 −15 .The discontinuous function φ(X) defined in Eq. ( 3) was replaced by the approximating function δ with δ = 0.01 in order to avoid difficulties in determining numerical derivatives.Figure 2 shows two numerical solutions for H 1 obtained for the case β = 0.1 , L = 1 in comparison to the asymptotic solutions obtained in [1].The stable and unstable asymptotic solution, H 1 0 ( )− and H 1 0 ( )+ , respectively, are shown in grey solid and grey dashed lines, respectively in Fig. 2, and were used as initial guess for the respective numerical solutions, which are shown in black solid and black dashed lines, respectively.Both the stable and unstable numerical solution are slightly shifted upstream.As for the stable solution, its crest is higher compared to the corresponding stable asymptotic solution, and its corresponding eigenvalue is also higher compared to the asymptotic value in the stable case, with the values being indicated in Fig. 2. The opposite is true for the unstable solution. (1) Employing the argument from the preceding section, larger values of the eigenvalue correspond to larger bottom friction forces, which have to be balanced by a larger weight of the fluid.Since the weight of the fluid is proportional to the spacial integral over H 1 , the larger value of the crest in the stabe case complies with this necessity.
It shall be mentioned that the asymptotic solutions depicted Fig. 2 are solutions to zeroth order in the expansion parameter β .From [4], an (so far unpublished) asymptotic solution extended to the order β is available, which is in perfect agreement with the numerical solutions.  .

Transient solutions
In [1], a method originally due to Scott [5] was applied to obtain an asymptotic solution for transient values of the crest heights H 1m (T ) for L = 1 , and a first-order ordinary differential equation valid for small values of β and slowly varying wave speeds V of the solitary wave was obtained: The transient asymptotic solution thereby is still of the form of a sech 2 function, with its amplitude given by 3(1 − V) and its crest being located at X X VT m ( ) 0 0 = + .For the special case of V ≡ 0 , Eq. (10) gives just an algebraic condition with the asymptotic values for the stationary crest positions X m ( ) 0 − and X m ( ) 0 + as solutions, which were obtained earlier in the asymptotic analysis.
In solving Schneider's full time dependent extended KdV Eq. (1) numerically, transient solitary wave solutions were obtained and compared to the asymptotic transient solutions following from Eq. (10).For the case of the numerical transient solutions of the full Eq.( 1), the stationary solution obtained from the eigenvalue problem of the stationary form of Eq. ( 1) was shifted in X direction by values of −β and + β , respectively, and posed as initial solution of the full time dependent problem.While the time integration was performed by a forward difference scheme, the space integration was performed by fixing the value of λ to the eigenvalue of the respective stationary solution.Since the space integration procedure then reduces to a simple boundary value problem not involving the determination of an eigenvalue, one of the four boundary conditions used in the stationary problem had to be dropped, which was chosen to be the condition (8).The other boundary conditions for the spacewise integration as well as the remaining conditions for the numerical integration were kept the same as in the stationary problem, in particular a fast decay behaviour of the solutions was assumed.The numerical solutions of the full Eq.( 1) were again obtained with the package MATLAB R2010A using the bvp4c integrator in space and a forward difference scheme in time.Typically, a time step size of Δ T ≈ 1 • 10 −1 was used, yielding Courant numbers of about 0.4 or less.
For In accordance between the numerical and the asymptotic analysis of Eq. (1), the transient numerical solutions show that solitary waves starting close to also converge to the stable stationary solution, whereas no convergence to an unstable stationary solution follows.Figure 3 shows this result for β = 0.1 and L = 1 .
For the case of stable solutions, shown as solid black lines, the height of crests initially located at X X with the asymptotic value of the stationary crest height H m 1 0 3 ( )− = , shown with a grey "+" symbol in Fig. 3.In contrast to the numerical solutions, the convergence of the asymptotic solutions is different and shows no spiralling behaviour.This is a consequence of the asymptotic expansion performed to "lowest order" in β , leading to Eq. ( 10) valid for slowly varying wave speeds, which cannot describe the oscillatory fast time behaviour of the wave crest heights seen in the numerical solutions.shown with a grey "x" symbol in Fig. 3.While the unstable asymptotic solution for the solitary wave initially shifted to the right by the value + β decays in the limit of large times, the corresponding numerical solution converges to a different stationary solution, with its values for the crest height and position seen as the lower terminating point of the black dashed line in Fig. 3.An analogon to this novel stationary solution is not present in the asymptotic analysis in [1].The unstable asymptotic solution for the solitary wave initially shifted to the left by the value −β finally converges back to the stable stationary solution.The corresponding numerical solution converges to a stationary solution again different to all the precedingly determined stationary solutions, where the convergence at the left side of the black dashed line takes place in a spiralling manner, with its limiting values for large times missing in Fig. 3.This stationary solution is different from the stable stationary solution determined previously, which is obvious since the value of λ is different, but it is remarkable in the sense that it represents a different eigensolution with a different eigenvalue to the stationary eigenvalue problem subject to the same boundary conditions as compared to the previously determined stable stationary eigensolution H 1 − .Also to this novel stationary solution there is no analogon in the asymptotic analysis considered here.This leads to the conclusion that the eigenvalue problem is degenerate and naturally raises the questions of how many eigensolutions exist and which of them will actually be realized.Figure 4 shows the shapes of the numerical solutions of Eq. ( 1 = + to the right by the amount + β , the crest of which had appeared as the the lower terminating point of the black dashed line in Fig. 3.The line marked T = 68 represents the solution which crest was seen as the the terminating point on the left end of the black dashed line in Fig. 3, resulting from perturbing H 1 + from its original position X X m m = + to the left by the amount −β .This solu- tion may give an idea of how the novel stationary solution it is approaching to will be looking like.At least in the case of the solution marked T = 300 , this makes clear why an analogon to this solution is not accessible via an asymptotic expansion to lowest order at the current stage, since the considerable differences to the solitary wave may not be seen as a small perturbation in terms of the parameter β .

Conclusion
Solitary waves in turbulent open-cannel flow are governed by an extended KdV equation which was derived in [1] without any recourse to turbulence modelling.The problem of finding stationary solitary waves for a piecewise constant enhanced bottom roughness in the cannel may be formulated as a coupled boundary-value eigenvalue problem, with the eigenvalue determining the necessary amount of increase in bottom roughness for the wave to comply with the law of momentum.For given values of the parameters, two eigensolutions were found corresponding to a stable and an unstable stationary wave in the channel.This is in accord with an asymptotic analysis of the stationary problem, given in [1].In determining transient solutions, waves starting close to also converge to the stable stationary solution, whereas no convergence to an unstable solution follows.Instead, the solutions starting close to the unstable stationary solution converge to novel stationary solutions, a feature not being present in the asymptotic analysis.This reveals the problem that multiple solutions to the eigenvalue problem exist, with solutions differing considerably from the classical soliton solution.

Fig. 1 A
Fig. 1 A solitary wave in open-channel flow.The nomenclature was taken from [1], with overbars indicating ensemble-averaged quantities of the velocity u(x, y, t) in direction x parallel to the channel bottom and the surface elevation h(x, t) .The subscript r denotes some reference position very far upstream where the flow is fully developed.The shaded region of length  indicates a region with increased bottom roughness.
termed stable and unstable in[1], corresponding to positions of the wave crests X , refer to the classical notion whether (full time-dependent) solutions of Eq.
the case of the asymptotic transient solutions, transient values of the crest heights H T m 1 0 ( ) ( ) were obtained fromEq.(10) by solving an initial value problem, in which the initial position X m , respectively, were shifted in X direction by −β and + β , respectively.

= − with the stationary crest height H m 1 −
, shown with a black "+" symbol in Fig.3, thereby justifying the identification of these solutions as a stable solutions.Although with different absolute values as compared to the numerical solutions, the same is true for the transient asymptotic solutions, shown as solid grey lines.The height of crests initially located at X

Fig. 3 1 + 1 +
Fig. 3 Crest heights H 1m of the transient numerical solutions of Eq. (1) (in black) in comparison to the asymptotic solutions (in grey) as function of their dimensionless position X m , with β = 0.1 , L = 1 and fixed values of λ .The solid lines represent stable transient solutions in the sense that they converge to the stable stationary solution with crest heights H m 1 − and

1 − 1 + 1 +
) for β = 0.1 , L = 1 .The line marked T = 200 represents the solution resulting from perturbing the stable stationary solution H obtained from the original eigenvalue problem from its original position X X m m = − to the left by the amount −β after a time T = 200 .This solution is almost identical with the stable stationary solution.The line marked H represents the unstable stationary solution obtained from the original eigenvalue problem.The line marked T = 300 represents the novel stationary solution resulting from perturbing H from its original position X X m m