Linear and Weakly Nonlinear Instability of Shallow Mixing Layers with Variable Friction

Linear and weakly nonlinear instability of shallow mixing layers is analysed in the present paper. It is assumed that the resistance force varies in the transverse direction. Linear stability problem is solved numerically using collocation method. It is shown that the increase in the ratio of the friction coefficients in the main channel to that in the floodplain has a stabilizing influence on the flow.*e amplitude evolution equation for the most unstable mode (the complex Ginzburg–Landau equation) is derived from the shallow water equations under the rigid-lid assumption. Results of numerical calculations are presented.


Introduction
Understanding the interaction between fast and slow fluid streams at river junctions and in compound channels is important for a proper description of mass and momentum transfer in shallow flows.In order to analyse the problem, hydraulic engineers apply depth-averaged shallow water equations [1] where either Chézy or Manning formulas are used to take into account the effect of bottom friction.ese formulas contain a constant friction coefficient determined from empirical relationships [1] if the Reynolds number of the flow and surface roughness are given.ere are cases, however, where the resistance force changes considerably in the transverse direction [2,3].One example of such a situation is flow in compound channels or rivers in case of floods.Hence, the analysis of instability characteristics of shallow mixing layers with variable friction is important from environmental point of view.In particular, contaminants and residues can accumulate in some parts of the flow due to instability.us, it is necessary to analyse factors affecting shallow flow instability and development of perturbations in a weakly nonlinear regime.ree different approaches for the investigation of shallow water flows are suggested in [4]: experimental studies, numerical modelling, and stability analysis.Experimental data [5][6][7] indicate that bottom friction has a stabilizing influence on the flow.Temporal linear stability analysis of shallow mixing layers [8][9][10][11][12] shows that bed friction reduces the width of a mixing layer and stabilizes the flow.e development of perturbations can also be analysed from a spatial point of view [13][14][15][16].Calculations show that bed friction reduces the spatial growth rate of perturbations.
ere are other factors affecting stability characteristics of shallow mixing layers: (a) flow curvature, (b) presence of solid particles in a fluid stream, and (c) variable friction force in the transverse direction.ese aspects of the linear stability problem are analysed in [13][14][15][16].In particular, a stably curved mixing layer has a stabilizing influence on the flow, while an unstably curved mixing layer destabilises the flow.In addition, the presence of solid particles in the stream reduces the growth rate of perturbations.Linear stability analysis is a useful tool for calculating critical values of the parameters characterizing the flow.However, linear theory cannot be used to analyse the development of perturbations in unstable regime.Assuming that the bed friction number is slightly smaller than the critical value (i.e., the flow is linearly unstable with a small growth rate), weakly nonlinear theories can be used to obtain an amplitude evolution equation for the most unstable mode [17][18][19][20][21][22]. e analysis in [17][18][19][20][21] shows that the amplitude evolution equation is the complex Ginzburg-Landau equation (GLE) under the rigid-lid assumption (free surface is treated as a rigid lid).As it is shown in [10], rigid-lid assumption works well for small Froude numbers.e Froude number is of order 0.2-0.3 in experiments [2] so that rigid-lid assumption is justified.e following cases of shallow flows are considered in [17][18][19][20][21]: (a) Rayleigh friction [19], (b) bottom friction modelled by the Chézy formula [21], and (c) generalization of the analysis in [21] for the case of slightly curved mixing layers and twophase flows in [17][18][19].Recent application of weakly nonlinear theory to shallow water flows without rigid-lid assumption [22] shows that the amplitude evolution equation for the most unstable mode is the complex Schrödinger equation.Since the solutions of the GLE vary from truly deterministic to almost chaotic [23] (depending on the values of the coefficients), in many cases, the GLE is used as a phenomenological equation for the analysis of spatiotemporal dynamics of complex flows.
e coefficients of the equation are estimated using experimental data.Complex phenomena in fluid mechanics (see, e.g., [24]) can be modelled using the obtained equation.
In the present paper, we perform linear and weakly nonlinear analysis of shallow mixing layers under the assumption that the friction force changes in the transverse direction.Weakly nonlinear analysis shows that the amplitude evolution equation for the most unstable mode is the complex GLE. e equation is derived from the shallow water equations under the rigid-lid assumption.e coefficients of the equation are calculated in closed form using linear stability characteristics of the flow.Results of numerical calculations are presented.Practical applications of the proposed model are also discussed.

Linear Stability Problem
Consider the system of shallow water equations under the rigid-lid assumption written in the form [25] zu zx where u and v are the depth-averaged velocity components in the x and y directions, respectively, h is water depth, p is the pressure, and c f (y) is the friction coefficient depending on the transverse coordinate y. e function c f (y) is assumed to be of the form where c f2 is a constant and c(y) is an arbitrary differentiable "shape" function.Since friction terms in (2)-(3) do not contain derivatives, the boundary conditions for the longitudinal velocity component are not needed (otherwise the problem would be overdetermined).us, the boundary conditions for the velocity components are essentially the same as in the case of inviscid flow (we need only zero normal velocity at the boundaries).
e boundary conditions are e assumption of an unbounded layer in the transverse direction is widely used in the stability analysis of mixing layers and wakes (see, e.g., [26,27]) and is adopted here.
Introducing the stream function ψ by the relations and eliminating the pressure from ( 1)-( 3), we obtain where c f y (y) � c f2 c ′ (y) is the derivative of c f (y) with respect to y and Δ � (z 2 /zx 2 ) + (z 2 /zy 2 ).Consider a perturbed solution to (7) of the form ψ(x, y, t) � ψ 0 (y) + εψ 1 (x, y, t) + ε 2 ψ 2 (x, y, t) where U(y) � ψ 0y (y) is the base flow velocity.e role of the small parameter ε will be clarified later.Substituting ( 8) into (7) and linearizing the resulting equation in the neighbourhood of the base flow U � U(y), we obtain where Using the method of normal modes, we represent the function ψ 1 in the form where ϕ 1 (y) is the amplitude of the normal perturbation and k is the wavenumber.Substituting (11) into (9), we obtain where S � c f2 b/h is the bed friction number and b is the halfwidth of the mixing layer.
2 Advances in Civil Engineering e boundary conditions follow from (5) and have the form Complex eigenvalues c � c r + ic i determine the linear stability of the base flow.e flow is linearly stable if all c i < 0 and unstable if at least one c i > 0. Problems ( 12) and ( 13) are solved numerically by means of the collocation method based on the Chebyshev polynomials [15].In the classical hydrodynamic stability theory (see, e.g., [26,27]), the base flow is obtained as a simple one-dimensional steady solution of the equations of motion (e.g., a plane Poiseuille flow is obtained as a steady one-dimensional solution of the Navier-Stokes equations).It is impossible to find a simple analytical solution of the form U � U(y) for shallow water equations ( 1)-( 3) due to the empirical character of the Chézy formula in ( 2) and ( 3). e paper [28] gives an example of a longitudinal velocity distribution in the form of a generalized power law. is profile does not have an inflection point (recall that the presence of an inflection point is the necessary condition for instability of inviscid flows [26]).
e obtained formula for the longitudinal dispersion coefficient in [28] is based on the assumption of straight symmetrical channel and uniform flow.e resulting formula is multiplied by a "catch-all" revision coefficient taking into account many kinds of nonuniformities in natural rivers.However, this adjustment does not change the shape of the longitudinal velocity distribution.e mixing layer case (the case considered in our paper) has a completely different longitudinal velocity distribution [2].e presence of a porous layer with much higher resistance than in the main channel results in a highly nonuniform velocity distribution resembling a hyperbolic tangent function with some asymmetry (however, the asymmetric profile is obtained after all the instabilities have been taken into account).e model profile has an inflection point.us, instabilities can occur.e hyperbolic tangent profile is widely used in stability analyses of shallow flows [8][9][10][11][12].e goal of our investigation is to analyse instabilities of the flow, obtain the linear stability characteristics, and consider the development of perturbations in a weakly nonlinear regime.
Assuming that the velocity of undisturbed flow is  U 1 and where all the dimensional quantities are represented with tildes.Using (  U 2 −  U 1 )/2 as the velocity scale and b as the length scale, we rewrite (14) in the dimensionless form where C � ( is the velocity ratio.All calculations in the paper are performed for the case C � 2. e base flow velocity distribution is given in Figure 1. e domain of the flow is defined as follows: the transverse coordinate y varies from −∞ to +∞.A porous layer occupies the region y < 0, while the region y > 0 corresponds to the main channel.Away from the shear layer, the base flow velocity approaches the undisturbed values U 1 � 1 and U 2 � 3, respectively.We assume that the volume fraction is small (in such a case, as it is shown in [2] shallow water equations of the form (1)-( 3) can be used to analyse the problem).
e shape function c(y) is selected in the form where β � c f1 /c f2 is the ratio of the friction coefficients in the main channel and floodplain (β < 1) and α is the factor reflecting the shape of the transitional zone from floodplain to the main channel.e following arguments are used to select c(y) in ( 16). e friction coefficients for the undisturbed flow (away from the shear layer) in the main channel and floodplain are denoted by c f1 and c f2 , respectively.A stronger resistance force results in a smaller base flow velocity so that ( 15) and ( 16) are consistent.In addition, we would like to remove discontinuity in the friction force used in [2] and consider a more realistic case of a continuous resistance changing with respect to the transverse coordinate.Examples of marginal stability curves are shown in Figure 2 for two different values of α. e flow is stable above the curves and unstable below the curves.Small circles in Figures 2-4 indicate the calculated points, while the solid lines represent the interpolating function.e flow becomes more unstable as α increases.is is related to the fact that larger α result in a steeper transition from the value of the friction coefficient in the main channel to the one in floodplain (16).
e critical values of the parameters S and k, namely, S c and k c , respectively, correspond to the maximum of the marginal stability curves.Figure 3 plots the critical bed friction numbers S c versus β.It is seen from the graph that the flow becomes more unstable as the parameter β decreases (i.e., if the ratio of the friction coefficients in the main channel to that in floodplain becomes smaller).
us, the increase in the ratio of the friction coefficients stabilises the flow.

Weakly Nonlinear Analysis
Suppose that the bed friction number S is slightly smaller than the critical value S c : us, the small parameter ε in asymptotic expansion (8) can be defined by the formula is relationship shows how close the bed friction number is to the corresponding critical value.Since the value of S can be calculated for any set of experimental conditions, the corresponding value of ε also can be computed.

Advances in Civil Engineering
Using the method of multiple scales and following [29], we introduce the slow time and longitudinal variable by the relations where c g is the group velocity.Substituting ( 17) and ( 19) into (8) and collecting the terms containing ε, ε 2 , ε 3 , . .., we obtain the system of equations where the operator L 1 is defined in (10).e function f 1 in (21) has the form Equation ( 20) coincides with (9) and corresponds to the linearized problem where the function ψ 1 is assumed to be of the form (11). Suppose that where k and c are, respectively, the critical wavenumber and speed of the most unstable perturbation in accordance with the linear theory, A(  ξ,  τ) is the unknown amplitude function, and the abbreviation c.c. means complex conjugate.e goal of the weakly nonlinear analysis is to derive the amplitude evolution equation for the function A(  ξ,  τ).Derivation of the amplitude equation can be found in the Appendix.e Ginzburg-Landau equation has the form e sign of the real part of μ (known as the Landau constant in the literature) plays an important role.e Landau constant has the "wrong sign" in [26].us, finite amplitude saturation is not possible, and higher order terms (with respect to A) quickly become important so that (25) can be used for a very short time (in other words, practical application of ( 25) is very limited).It is shown in [17][18][19][20][21] (in contrast to [29]) that, for shallow water flows, 4 Advances in Civil Engineering the Landau constant in (25) has the "right sign" so that (25) can be used (and is successfully used in [17][18][19][20][21]) in order to describe some important features of shallow wake flows.Equation ( 25) can be written in a more convenient form using the substitutions us, Some closed-form solutions of ( 27) are known in the literature [23].One of the simplest solutions is the solution of the form where Since ( 28) represents a periodic solution in both space and time, it is important to investigate the linear stability of such a plane wave.If solution ( 28) is stable, then periodic wave-like solutions can also be observable in experiments.
Assuming that linear stability of ( 28) can be analysed [23].Substituting (30) into (27), we obtain the equation for λ.For the case of small k, the stability condition has the form provided that q satisfies the inequality Condition ( 31) is known as the Benjamin-Feir stability condition.If (31) is not satisfied, then plane wave solutions of the form (28) are unstable (and, therefore, are not observable in experiments).

Numerical Results
e complex Ginzburg-Landau equation as the amplitude evolution equation for the most unstable mode is analysed in the previous section using asymptotic expansion (8).Since it is not known in general whether ( 8) is convergent, it would be nice to have a criterion which can be used to convince us that the Ginzburg-Landau model can adequately represent the dynamics of a fully nonlinear model at least at the initial stage of transition period when the base flow becomes linearly unstable.One such criterion is given in [31].e criterion is as follows: if the growth rates of the unstable perturbation can be well approximated by a parabola in the whole range of unstable wavenumbers, then the GLE can be used to analyse the dynamics of the flow (at least in the beginning of the nonlinear regime).In order to test this assertion, we computed the growth rates for the range of unstable wavenumbers for the following values of the parameters of the problem: β � 0.3 and S � 0.7, 0.9, and 0.11.
e results of calculations are shown in Figure 4, where S < S c and the dots represent the calculated growth rates λ r � kc i for three values of S. e solid lines represent the best-fit parabolas.As can be seen from the graph, the calculated growth rates are well approximated by parabolas in all regions of unstable wavenumbers.us, we conclude that the GLE can be successfully used to analyse the dynamics of the flow above the threshold.e results of the numerical computations of the linear stability characteristics and the coefficients of the GLE are shown in Table 1.
Table 2 shows numerical values of the coefficients c 1 and c 2 (26) for different values of α.As can be seen from Table 2, condition (31) is satisfied for all cases considered.Hence, periodic solutions of the GLE are stable and can be observed in experiments.e Landau constant (the real and the initial condition Method of lines implemented in Mathematica is used for the numerical solution of problems ( 27), (33), and (34).e value of L � 100 is used for all simulations below.e first computation is performed for the case c 1 �1.3294 and c 2 � 0.1252.e values of these parameters are taken from Table 2. e function f(ξ) in ( 34) is assumed to be a small random noise of order 0.01.e results are shown in Figure 5. Since we are interested in the analysis  us, the modulus of the amplitude reaches a constant value.
e second set of computations is performed for the case f(ξ) � ����� 1 − q 2  e iqξ where q � 0.5.e results are shown in Figure 6. e value of q in this case corresponds to the stable range of wavenumbers (both conditions (31) and (32) are satisfied in this case).As a result, perturbations do not propagate throughout the domain and finite amplitude saturation occurs.
Finally, we consider the case where the Benjamin-Feir stability condition (31) is not satisfied.e values of the parameters are taken from [21]: c 1 � −0.799564 and c 2 � 2.189654 (these parameters correspond to weakly nonlinear analysis of wake flows).Random noise of order 0.01 is used as the initial condition.e results are shown in Figure 7.As can be seen from the figure, stabilization of the amplitude does not occur in this case.Moreover, instability is spreading throughout the whole domain (see the parts of the graph for large times).is fact is confirmed by experiments (see, e.g., discussion of wake flow instability in [21] and references therein)-periodical structures do not appear in experiments in the weakly nonlinear regime.
ese examples illustrate the well-known fact that the Ginzburg-Landau model is quite rich in terms of different solutions.Illustrative computations in Figures 5-7 show that both initial conditions and the values of the coefficients are responsible for spatiotemporal dynamics of the amplitude.
Direct comparison of linear and weakly nonlinear instability results with field data is not an easy task.In order to compare directly the results of weakly nonlinear theory, one has to use carefully designed experiments (or numerical experiments).Available experimental data (see, e.g., [2,3]) give the velocity distribution in the longitudinal direction which is averaged over sufficiently large period of time (in this case, nonlinear effects have been already taken into account).e goal of weakly nonlinear theory is to describe the initial stage of the development of instability.Using weakly nonlinear calculations, we can answer the following questions: Should we expect finite amplitude saturation?Would a periodic solution be stable (and, therefore, observable in experiments)?Illustrative calculations for the solution of the boundary value problem containing GLE are consistent with experimental observations.For example, if the Benjamin-Feir instability condition (31) is not satisfied for wake flows, then periodic solutions of the GLE are unstable (and are not observable in experiments).
is fact is supported by experimental observations for wake flows (see the references in [21]).On the other hand, our calculations show that condition (31) is satisfied for shallow mixing layers.us, periodic solutions of the GLE are stable (this fact is supported by experimental observations in [2,3]).

Conclusions
Linear stability of shallow mixing layers is investigated under the assumption that the resistance force varies in the transverse direction.Stability of the base flow depends on the ratio of the friction coefficients in the main channel to that in the floodplain.Weakly nonlinear theory is used to derive the amplitude evolution equation for the most unstable mode.It is shown that the amplitude evolution equation is the Ginzburg-Landau equation with complex coefficients.e coefficients of the equation are obtained in closed form available for calculations.Numerical results show that plane wave solutions of the GLE are stable.

Figure 4 :
Figure 4: Quadratic approximations of the growth rates for the following values of the parameters: S � 0.07, 0.09, and 0.11 (from top to bottom).

Figure 5 :Figure 6 :Figure 7 :
Figure 5: Plot of the |  A| for the case c 1 � 1.3294 and c 2 � 0.1252 (the initial condition corresponds to small random noise).

Table 1 :
Linear and weakly nonlinear calculations for four values of α and β � 0.3.

Table 2 :
Numerical values of the coefficients c 1 and c 2 of the GLE (27) for four values of α and β � 0.3.