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

Authors: Irina Eglite; Andrei Kolyshkin
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

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. The 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.


---

## Body

## 1. 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. These formulas contain a constant friction coefficient determined from empirical relationships [1] if the Reynolds number of the flow and surface roughness are given. There 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. Thus, it is necessary to analyse factors affecting shallow flow instability and development of perturbations in a weakly nonlinear regime.Three different approaches for the investigation of shallow water flows are suggested in [4]: experimental studies, numerical modelling, and stability analysis. Experimental data [5–7] indicate that bottom friction has a stabilizing influence on the flow. Temporal linear stability analysis of shallow mixing layers [8–12] shows that bed friction reduces the width of a mixing layer and stabilizes the flow. The development of perturbations can also be analysed from a spatial point of view [13–16]. Calculations show that bed friction reduces the spatial growth rate of perturbations.There 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. These aspects of the linear stability problem are analysed in [13–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–22]. The analysis in [17–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. The Froude number is of order 0.2-0.3 in experiments [2] so that rigid-lid assumption is justified. The following cases of shallow flows are considered in [17–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 two-phase flows in [17–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. The 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. The equation is derived from the shallow water equations under the rigid-lid assumption. The 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.


## 2. Linear Stability Problem

Consider the system of shallow water equations under the rigid-lid assumption written in the form [25](1)∂u∂x+∂v∂y=0,(2)∂u∂t+u∂u∂x+v∂u∂y+∂p∂x+cfy2huu2+v2=0,(3)∂v∂t+u∂v∂x+v∂v∂y+∂p∂y+cfy2hvu2+v2=0,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 cfy is the friction coefficient depending on the transverse coordinate y. The function cfy is assumed to be of the form(4)cfy=cf2γy,where cf2 is a constant and γ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). Thus, 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).The boundary conditions are(5)v±∞=0.The 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(6)u=∂ψ∂y,v=−∂ψ∂x,and eliminating the pressure from (1)–(3), we obtain(7)Δψt+ψyΔψx−ψxΔψy+cfy2hΔψψx2+ψy2+cfy2hψx2+ψy2ψy2ψyy+2ψxψyψxy+ψx2ψxx+cfyy2hψyψx2+ψy2=0,where cfyy=cf2γ′y is the derivative of cfy with respect to y and Δ=∂2/∂x2+∂2/∂y2.Consider a perturbed solution to (7) of the form(8)ψx,y,t=ψ0y+εψ1x,y,t+ε2ψ2x,y,t+ε3ψ3x,y,t+⋯,where Uy=ψ0yy is the base flow velocity. The 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=Uy, we obtain(9)L1ψ1=0,where(10)L1ψ≡ψxxt+ψyyt+ψ0yψxxx+ψoyψyyx−ψ0yyyψx+cf2hψ0yψxx+2ψoyyψy+2ψ0yψyy+cfyhψ0yψy.Using the method of normal modes, we represent the functionψ1 in the form(11)ψ1x,y,t=ϕ1yeikx−ct,where ϕ1y is the amplitude of the normal perturbation and k is the wavenumber.Substituting (11) into (9), we obtain(12)Uk−ck−iγSUϕ1″−iγSUy+γySUϕ1′+k3c−k3U−kUyy+ik2γS2Uϕ1=0,where S=cf2b/h is the bed friction number and b is the half-width of the mixing layer.The boundary conditions follow from (5) and have the form(13)ϕ1±∞=0.Complex eigenvaluesc=cr+ici determine the linear stability of the base flow. The flow is linearly stable if all ci<0 and unstable if at least one ci>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=Uy for shallow water equations (1)–(3) due to the empirical character of the Chézy formula in (2) and (3). The paper [28] gives an example of a longitudinal velocity distribution in the form of a generalized power law. This 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]). The obtained formula for the longitudinal dispersion coefficient in [28] is based on the assumption of straight symmetrical channel and uniform flow. The 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. The mixing layer case (the case considered in our paper) has a completely different longitudinal velocity distribution [2]. The 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). The model profile has an inflection point. Thus, instabilities can occur. The hyperbolic tangent profile is widely used in stability analyses of shallow flows [8–12]. The 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 isU˜1 and U˜2 at y˜=−∞ and y˜=+∞, respectively, the function U˜=U˜y is chosen in the form(14)U˜y˜=U˜1+U˜22+U˜2−U˜12tanhy˜,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(15)Uy=C+tanhy,where C=U˜2+U˜1/U˜2−U˜1 is the velocity ratio. All calculations in the paper are performed for the case C=2. The base flow velocity distribution is given in Figure 1. The 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 U1=1 and U2=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).Figure 1
Base flow velocity distribution.The shape functionγy is selected in the form(16)γy=β+12+β−12tanhαy,where β=cf1/cf2 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.The following arguments are used to selectγy in (16). The friction coefficients for the undisturbed flow (away from the shear layer) in the main channel and floodplain are denoted by cf1 and cf2, 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 α. The 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. The flow becomes more unstable as α increases. This 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). The critical values of the parameters S and k, namely, Sc and kc, respectively, correspond to the maximum of the marginal stability curves. Figure 3 plots the critical bed friction numbers Sc 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). Thus, the increase in the ratio of the friction coefficients stabilises the flow.Figure 2
Marginal stability curves for the casesβ=0.3,α=1 (bottom curve) and α=20 (top curve).Figure 3
Critical bed friction numbersSc versus β for the case α=1.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).


## 3. Weakly Nonlinear Analysis

Suppose that the bed friction numberS is slightly smaller than the critical value Sc:(17)S=1−ε2Sc.Thus, the small parameterε in asymptotic expansion (8) can be defined by the formula(18)ε=1−SSc.This relationship shows how close the bed friction number is to the corresponding critical value. Since the value ofS can be calculated for any set of experimental conditions, the corresponding value of ε also can be computed.Using the method of multiple scales and following [29], we introduce the slow time and longitudinal variable by the relations(19)τ˜=ε2t,ξ˜=εx−cgt,where cg is the group velocity. Substituting (17) and (19) into (8) and collecting the terms containing ε,ε2,ε3,…, we obtain the system of equations(20)L1ψ1=0,(21)L1ψ2=f1,(22)L1ψ3=f2,where the operator L1 is defined in (10). The function f1 in (21) has the form(23)f1=cgψ1xxξ˜+ψ1yyξ˜−2ψ1xξ˜t−3Uψ1xxξ˜−ψ1yψ1xxx−ψ1yψ1yyx−Uψ1ξ˜yy+ψ1xψ1xxy+ψ1xψ1yyy+Uyyψ1ξ˜−cf2hψ1xxψ1y+2Uψ1xξ˜+2ψ1yyψ1y−2UUy+2ψ1xψ1xy−cfy2h12ψ1x2+ψ1y2−U2.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(24)ψ1x,y,t=Aξ˜,τ˜ϕ1yeikx−ct+c.c.,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. The 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. The Ginzburg–Landau equation has the form(25)∂A∂τ˜=σA+δ∂2A∂ξ˜2−μA2A.The sign of the real part ofμ (known as the Landau constant in the literature) plays an important role. The Landau constant has the “wrong sign” in [26]. Thus, 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–21] (in contrast to [29]) that, for shallow water flows, the Landau constant in (25) has the “right sign” so that (25) can be used (and is successfully used in [17–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(26)τ=τ˜σr,ξ=ξ˜σrδr,A˜=Aμrσre−ic0σrτ˜,c0=σiσr,c1=δiδr,c2=μiμr.Thus,(27)∂A˜∂τ=A˜+1+c1i∂2A˜∂ξ2−1+c2iA˜2A˜.Some closed-form solutions of (27) are known in the literature [23]. One of the simplest solutions is the solution of the form(28)A˜=a0eiqξ+iωτ,where(29)a0=1−q2,ω=−c2−c1−c2q2.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(30)A˜=a0+a^eikx+λt+a^∗e−ikx+λteiqx+iϖt,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(31)1+c1c2>0,provided that q satisfies the inequality(32)q2<1+c1c22c22.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).


## 4. Numerical Results

The 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]. The 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.3andS=0.7,0.9,and0.11.The results of calculations are shown in Figure4, where S < Sc and the dots represent the calculated growth rates λr=kci for three values of S. The 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. Thus, we conclude that the GLE can be successfully used to analyse the dynamics of the flow above the threshold.The results of the numerical computations of the linear stability characteristics and the coefficients of the GLE are shown in Table1.Table 1
Linear and weakly nonlinear calculations for four values ofα and β=0.3.
α
0.25
0.5
1.0
1.5
k
0.442
0.437
0.438
0.437
S
0.198
0.205
0.211
0.214
C
1.972
1.985
2.004
2.018
σ
0.195 − 0.487i
0.195 − 0.080i
0.183 − 0.133i
0.174 − 0.173i
μ
2.374 + 0.690i
2.151 + 0.687i
2.090 + 0.516i
2.092 + 0.262i
c
g
1.956
1.981
2.007
2.018
δ
7.077 + 13.243i
7.330 + 12.434i
7.403 + 10.752i
7.255 + 9.645iTable2 shows numerical values of the coefficients c1 and c2 (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. The Landau constant (the real part of μ) is positive. Thus, finite amplitude saturation is possible.Table 2
Numerical values of the coefficientsc1 and c2 of the GLE (27) for four values of α and β=0.3.
β
0.3000
0.3000
0.3000
0.3000
α
0.25
0.5
1.0
1.5
c
1
=
δ
i
/
δ
r
1.8713
1.6963
1.4524
1.3294
c
2
=
μ
i
/
μ
r
0.2906
0.3194
0.2469
0.1252
1
+
c
1
c
2
2.1619
2.0157
1.6993
1.4547Numerical solutions of the GLE (27) are presented below for different values of the parameters c1 and c2 and different initial conditions. The problem is formulated as follows: find the solution of (27) for the given boundary conditions(33)A˜ξ=0=0,A˜ξ=L=0,and the initial condition(34)A˜τ=0=fξ.Method of lines implemented in Mathematica is used for the numerical solution of problems (27), (33), and (34). The value of L=100 is used for all simulations below.The first computation is performed for the casec1 = 1.3294 and c2 = 0.1252. The values of these parameters are taken from Table 2. The function fξ in (34) is assumed to be a small random noise of order 0.01. The results are shown in Figure 5. Since we are interested in the analysis of spatiotemporal development of a perturbation, the modulus of the amplitude A˜ is plotted versus the longitudinal coordinate ξ and time τ. The parameters of the problem satisfy (31) and (32) (in other words, are in the region of stability). Thus, the modulus of the amplitude reaches a constant value.Figure 5
Plot of theA˜ for the case c1 = 1.3294 and c2 = 0.1252 (the initial condition corresponds to small random noise).The second set of computations is performed for the casefξ=1−q2eiqξ where q = 0.5. The results are shown in Figure 6. The 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.Figure 6
Plot of theA˜ for the case c1 = 1.3293 and c2 = 0.1251 (the initial condition is of the form fξ=1−q2eiqξ where q = 0.5).Finally, we consider the case where the Benjamin–Feir stability condition (31) is not satisfied. The values of the parameters are taken from [21]: c1 = −0.799564 and c2 = 2.189654 (these parameters correspond to weakly nonlinear analysis of wake flows). Random noise of order 0.01 is used as the initial condition. The 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). This 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.Figure 7
Plot of theA˜ for the case c1 = −0.799564 and c2 = 2.189654.These examples illustrate the well-known fact that the Ginzburg–Landau model is quite rich in terms of different solutions. Illustrative computations in Figures5–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). The 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). This 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. Thus, periodic solutions of the GLE are stable (this fact is supported by experimental observations in [2, 3]).


## 5. 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. The coefficients of the equation are obtained in closed form available for calculations. Numerical results show that plane wave solutions of the GLE are stable.



---

*Source: 8079647-2018-03-19.xml*