1. Introduction and Objective
The Navier–Stokes equations constitute the foundations of fluid mechanics. Published by Claude Navier in two articles in 1821 and 1823, they were definitively consolidated with the work of G.G. Stokes, who settled the no-slip condition at the fluid–solid boundary in 1845 [
1]. Due to the nonlinearity of the equations, analytical solutions are limited to relatively few cases, typically those in which the nonlinear terms are absent, a condition that arises in specific flow scenarios, for instance, when parallel streamlines are present. This is the case of some flows solved in the seminal articles by Navier and Stokes. Thus, Navier derived the velocity distribution for a steady uniform free-surface flow in a rectangular channel [
2]. Stokes, in his study on the motion of pendulums, solved “the extremely simple case of an oscillating plane” [
3] (p. 20). In a note to the same study, Stokes also provided the solution of the velocity distribution due to the sudden motion of a plate initially at rest [
3] (Note B, p. 100). Problems where the solid boundary is not a plane, but a deformed surface are much more challenging to solve analytically. In this case, analytical solutions are available only when the nonlinear terms can be addressed through linearization. Focusing on the velocity distribution in free-surface flows, perhaps the classic problem corresponds to a wavy bed with undulations transverse to the flow, as solved by Wang [
4], who addressed the problem for a low Reynolds number film flow over a wavy bed with amplitude much higher than the wavelength. The propagation of water waves over rigid rippled beds was studied numerically by Huang and Dong [
5].
In the case of unsteady fluid motion, aside from instances where the equations can be linearized, there is limited knowledge about the structure of solutions as the physical and geometric parameters of the problem undergo variations [
6]. Wierschem et al. [
7,
8] determined the velocity distribution of film flows down an inclined wavy channel, considering the effect of surface tension. The flow over a moving wavy bottom was approached by Benjamin [
9]. Another flow configuration for which an analytical result was obtained is that studied by Lyne [
10], who, in addition to the velocity field, determined the steady stream generated by an oscillatory viscous flow over a wavy wall. Lyne’s solution, obtained by the method of perturbations was extended by Kaneko and Honji [
11], and subsequently by Vittori [
12].
The free surface behavior, including the stability, of the flow over a wavy bed has been studied by Pozrikidis [
13], Bontozoglou and Papapolymerou [
14], Scholle et al. [
15], Trifonov [
16,
17,
18], and Tseluiko et al. [
19], among many others. Pozrikidis [
13] and Scholle et al. [
15] made the analysis considering the lubrication approximation. Trifonov [
18] also computed the velocity field numerically. Dávalos-Orozco [
20] studied the nonlinear instability of a thin film flowing down a smoothly deformed surface. Although his analysis is general as the surface can vary as much as in the direction of the main flow as in the direction normal to it, for simplicity, the numerical analysis is carried out considering a sinusoidal wall that deforms in the flow direction.
The oscillating flow over a wavy bottom problem has been approached by many authors due to its implication in sediment transport and ripple formation at the bottom of sea waves. In addition to the analytical approaches developed by Lyne [
10], Kaneko and Honji [
11], and Vittori [
12], the works by Hara and Mei [
21], Nishimura et al. [
22], and Shen and Chan [
23] can be mentioned. In all the above studies, the generatrix of the bed undulations is transverse to the direction of the main flow, as sketched in
Figure 1.
The configuration in which the main flow is parallel to the bed deformation has also been studied, although all the literature found by the authors refers to the case in which the fluid is confined between two plates, and the flow is driven by a pressure gradient. Thus, Szumbarski [
24] computed semi-analytically the velocity distribution and performed a linear stability analysis of the flow. Results of linear stability analysis have also been presented by Kowalewski et al. [
25], Szumbarski and Błoński [
26], and Moradi and Floryan [
27]. Szumbarski et al. [
28] studied numerically the effect of the wall corrugation on the hydraulic resistance of the flow. Among the results presented in their article, they show the circulation pattern induced by the velocity’s normal to the main flow direction, generating closed circulation cells. Gepner et al. [
29] focused on the nonlinear solutions of the velocity field generated by the flow limited by two corrugated plates and the effect of the induced secondary currents enhancing diffusive transport. Yadav and Gepner [
30] studied the effect on the flow instabilities for a geometry consisting in a fixed corrugated bottom plate and an upper moving plane plate. Later, Gepner and Floryan [
31] approached numerically the flow driven by a pressure gradient of a fluid bounded by two sinusoidal walls which are in phase in flows at low and moderate sub-turbulent Reynolds numbers. As the fluid moves transverse to the plate deformations, it is prone to separation and generation of a recirculating region in the valleys of the walls. Previously, Smith [
32] had studied the viscous Poiseuille 2-D flow, the walls of which had symmetrical and non-symmetrical cavities. His approach was theoretical, and the final equations were solved numerically. Bordner [
33] had solved numerically the laminar boundary layer over a periodic stationary wavy surface, noticing the separation effect near the valleys of the wavy surface. A similar problem was solved by Sobey [
6], who modeled and computed numerically the flow in an extracorporeal membrane oxygenator for both the steady and oscillatory flow conditions. Caponi et al. [
34], also numerically, found that separation is obtained from the linearized Navier–Stokes equations for the laminar flow near a moving wavy surface. They also determined the range of validity of the linear approximation of the Navier–Stokes equations. Tsangaris and Leiter [
35], using a perturbation method, and solving numerically some resulting equations, studied the separation of viscous laminar flows through wavy-walled channels. The separation characteristics of fluid flow inside two parallel wavy plates were also studied numerically by Mahmud et al. [
36], who found that separation occurs for a critical Reynolds number which depends on dimensionless values of the amplitude and wavelength of the bed undulation.
Figure 1.
Configuration analyzed in previous studies with free-surface or unconfined flows: (
a) free-surface fluid flow over a fixed bottom [
1]; (
b) flow induced by the sudden motion of a plane plate [
3]; (
c) flow induced by the oscillatory motion of a plane plate [
3]; (
d) flow over a fixed wavy bed [
4,
5,
9,
13,
14,
15,
16,
17,
18,
19,
20]; (
e) flow resulting from the motion of a wavy bottom [
9,
34]; and (
f) oscillatory flow over a wavy bottom [
10,
11,
12,
21,
22,
23].
Figure 1.
Configuration analyzed in previous studies with free-surface or unconfined flows: (
a) free-surface fluid flow over a fixed bottom [
1]; (
b) flow induced by the sudden motion of a plane plate [
3]; (
c) flow induced by the oscillatory motion of a plane plate [
3]; (
d) flow over a fixed wavy bed [
4,
5,
9,
13,
14,
15,
16,
17,
18,
19,
20]; (
e) flow resulting from the motion of a wavy bottom [
9,
34]; and (
f) oscillatory flow over a wavy bottom [
10,
11,
12,
21,
22,
23].
So far, the cited references on flows over deformed walls correspond to pre-established geometries that do not change over time. However, there are studies in which the upper and lower walls confining the fluid fluctuate temporarily. Thus, Carlsson et al. [
37,
38] investigated a scenario involving the confinement of fluid within flexible walls, exhibiting transverse motion in the form of small amplitude standing waves. The velocity field was ascertained through a perturbation approach, revealing that streaming occurs at the second order, leading to the formation of cellular flow patterns in the channel. Lebbal et al. [
39] analyzed the linear instabilities of the flow between flexible walls, considering that the base condition of the flow and the displacement of the walls are sinusoidal in both time and the direction of the main flow. The resulting equations were solved numerically. Vladimirov [
40] considered that the flow domain is the half-space, with a boundary subject to tangential vibrations. The solution of the problem contains both oscillating and non-oscillating components and the streaming flows (either steady or unsteady) are obtained.
The goal of this article is to determine the velocity field of a free-surface flow when the bottom undergoes harmonic deformation in the direction transverse to the main flow (i.e., the undulations’ generatrix is parallel to the main flow), asymptotically approaching a maximum amplitude from a flat initial condition. The problem addressed here arose from the analysis of the generation of rills in their incipient stage due to overland flow over erodible soil [
41]. The original problem is highly complex due to the interaction between the flow and the soil. As the soil erodes, it alters the geometry, subsequently changing the flow characteristics, and so on, until eventually reaching a final equilibrium condition. The primary objective of the article is not to solve the physical problem, as the equations accounting for soil erosion are not considered. The problem presented is merely an idealization of the bed deformation process, or a mimicry of it. Nevertheless, as far as the authors are aware, this has not been accomplished previously. Thus, the aim is to determine the flow characteristics of the rill formation process considering the depth of erosion to be small compared to the width of the rill. The main flow is assumed to be uniform (i.e., the flow depth does not change along the main flow direction), with the primary velocity component oriented in the direction transverse to the bed undulation (the main flow is parallel to the undulation generatrix), as illustrated in
Figure 2. The specific goal of this paper is to present analytical solutions for the three velocity components and the distortion of the free surface caused by the deformation of the channel bottom. The bottom starts as initially flat and evolves to reach an equilibrium wavy shape.
4. Analysis and Visualization of the Solution
4.1. Dimensionless Parameters
To generalize the results,
,
, and
are used as scales to make dimensionless
, and
, respectively:
The minus sign in
appears because
(Equation (38)). The dimensionless bed perturbation (Equations (37) and (38)) is written as
.
is also defined, such that
. Thus, the following dimensionless parameters arise:
is always positive and it appears naturally from the scaling of the equations for the perturbed components of the transverse velocities
or
(Equations (16) or (17)). Let
be a scale for the velocity
such that the dimensionless velocity is given by
. In terms of the dimensionless variables, Equation (16) becomes:
Calling
the velocity scale to make dimensionless the perturbation of the velocity in the z direction,
. Thus, the dimensionless version of Equation (17) is:
As the momentum equations for disturbances in the and directions are linear and depend on only one velocity component, the final result is independent of the chosen velocity scale, depending on two parameters, and . From the above equations, clearly is a kind of Reynolds number that represents the ratio between the fluid’s inertia due to the vertical deformation of the bed and the fluid’s viscous resistance. Note that from the continuity equation for the perturbed velocity components (Equation (14)), the following relationship between the velocity scales is obtained: .
For ,
is a real number, with , and the solutions for the stream function, velocity field, flow depth, and free surface deformation involve hyperbolic functions (Equations (47)–(49), (52), (53), and (59)). When ,
is an imaginary number and the hyperbolic functions become trigonometric ones (Equations (60)–(65)), changing the behavior of the solutions, as will be seen below. Defining , and .
4.1.1. Behavior of the Solutions When ( Real)
Taking the real part of Equation (44) and rendering it dimensionless results in:
Similarly, we can define
and
as new dimensionless velocity components, such that
The transverse velocities have been scaled with and not with the scales and , solely to facilitate the graphical presentation of the results. Indeed, by having the product appear as a coefficient in the transverse velocities, generating the new variables and reduces the number of parameters to be considered.
Denoting
as the ratio between the deformation of the free surface and the bed, Equation (53) in terms of the dimensionless parameters is
For (i.e., ), is in the range for any value of . This means that the amplitude of the free surface deformation is always smaller than the bed deformation. Also, the term in brackets in Equations (73) and (74) is always less than 1, for any .
Finally, the dimensionless perturbation of the longitudinal component of the velocity is as follows:
It is easier to realize the analytical results presented above if they are plotted as in
Figure 4,
Figure 5,
Figure 6,
Figure 7 and
Figure 8. It must be noted that the general behavior of the streamlines, velocity components, and deformation of the free surface essentially does not change with the values of
and
, and the general shapes are those presented in the figures. Arbitrary values of these parameters were used in the computations. The spanwise dimension (
-axis) of the plots correspond to a wavelength
of the bed deformation. The obtained results require small values of
to maintain the validity of the linear approximation. However, using
produces results with variations that are practically imperceptible in the graphs. Thus, to facilitate visualizing the variations introduced by the first-order term in the base solution, the amplitudes of the perturbations were greatly amplified in the figures depicting flow properties.
4.1.2. Behavior of the Solutions When ( Imaginary)
By employing relationships that transform hyperbolic functions of a complex variable into trigonometric ones and taking their real parts, we obtain the dimensionless versions of the equations for the streamline, velocity components, perturbed deformation of the free surface, and the ratio between the deformation of the free surface and the bed. They are presented below:
It is worth discussing the behavior of
before analyzing the stream function, velocities, and free surface deformation.
becomes indeterminate when
, fluctuating between positive and negative values, as it is observed in
Figure 9. Although the figure corresponds to the particular case
, the general behavior is the same for any
. The change in sign of
indicates that the free surface deformation is not always in phase with the bed deformation. Furthermore, due to viscous dissipation, it is not possible to have values of the amplitude of the free surface greater than the amplitude of the bed, limiting the range of
to the interval
. This last condition is important because it imposes a relation between
and
. For the upper limit
,
with
must be satisfied. For the lower limit
, the relationship between
and
is given by
. For a given interval of
, the larger the value of
, the greater the number of branches of
.
The behavior of
when
is imaginary is completely different from the case when it is real. The presence of an oscillatory function like cosine could lead to a change in the sign of
along
for a given transversal location
. In contrast to the case of real
, many configurations of streamlines arise depending on the values of
and
, as shown in the examples in
Figure 10. This sometimes generates cells to make compatible the velocities induced simultaneously by both bottom and free surface deformations.
The dimensionless transverse velocities
and
are presented in
Figure 11 and
Figure 12, respectively, for the three values of
corresponding to
shown in
Figure 10. The change in sign of the velocities is a result of the presence of circulation cells. The deformations of the free surface and bottom over time for
are shown in
Figure 13. As it was indicated, both deformations can be in phase or out of phase by
, depending on the value of
(
can be positive or negative). Finally, the dimensionless main flow velocity
is presented in
Figure 14, where the effect of the perturbation on the base velocity has been exaggerated for easier visualization.
4.2. Flow Separation
The presence of a flow over protrusions on the channel bed may lead to flow separation. Currently, there is a lack of studies determining the separation conditions for the specific configuration explored in this article—convergent flows over a sinusoidally undulated bed with a free surface. Nevertheless, as an initial approximation, extrapolating results from studies on pressure-driven flows in two-dimensional pipes with sinusoidally varying walls along the main flow direction is feasible. To estimate the separation condition for transverse flow with a free surface, we rely on a previous result obtained for flows constrained by wavy walls. This estimation is made despite differences in flow configurations, such as the constant net flow existing in the -direction in pressure-driven cases (where separation occurs in the plane), which contrasts with the variable cross-sectional flow ( plane) of the problem addressed in this article.
In this regard, the analytical results of Tsangaris and Leiter [
35] based on the separation of flow from a laminar boundary layer over a sinusoidally varying surface is applied. These authors determined that the separation is conditioned by a critical Reynolds number (
) that is related to the dimensionless amplitude of bed deformation (
). In general, this relationship is not simple and must be determined by numerically integrating the vorticity equation in terms of the stream function, obtaining a relationship
, which they graphically present in their article. The critical Reynolds number at which separation occurs is defined as
, where
is half of the mean separation between walls, and
is the average velocity.
, with
being the amplitude of the wall deformation. For large Reynolds numbers, the authors provide the approximation
.
To apply the separation point condition of Tsangaris and Leiter to the cases generated in our study, the velocity
is taken as the maximum transverse velocity
, and
is set equal to the flow height coordinate,
. Thus, the critical Reynolds number is given by
, which, expressed in terms of the dimensionless variables given by Equations (66) and (67), is written as
. For
and considering the most unfavorable scenario (
an upper bound for
is obtained:
Similarly, for the case
, the flow separation condition is bounded by
As an example, for the case
, if
, and
, the result is
. Similarly, when
, considering
, and
, the result is
< 5.5. In both cases, the values of
are much lower than
, obtained from the graphical relationship presented in the article by Tsangaris and Leiter [
35].
Choosing represents merely an arbitrary value for the dimensionless amplitude of the bed deformation (). While the preceding results confirm the absence of flow separation, they do not guarantee the extent to which the linear solution deviates from the one obtained by accounting for the non-linear terms in the momentum equations. In strict terms, linearizing the problem requires . Therefore, it is reasonable to expect that the results obtained from the linear approach are valid for . However, this a priori assumption needs validation through a nonlinear analytical solution of the problem or via numerical modeling. Both approaches are beyond the scope of this work, and the authors are not aware of this problem being addressed previously.
4.3. Conditions for the Validity of the Results Obtained from the Linearized Equations
In order to determine the limitations of the applicability of the first-order results obtained in this article, the conditions that the parameters must satisfy when non-dimensionalizing the linearized equations are established below. The analysis can be carried out interchangeably for the momentum equations of perturbations in the or directions, taking into account that the velocity scales for both components are related as .
The momentum equation for the z-component of the disturbance for a uniform flow is given by:
Considering the scales mentioned before,
and
, the scalings for the different terms of Equation (79) are obtained:
To neglect the nonlinear terms means:
The velocity in the
-direction is induced by the bed deformation (Equation (36)), so it is permissible to consider
. Then, inequality (88) is reduced to:
The previous relationship corresponds to
and is nothing more than the condition of small amplitude of bed deformation that has been imposed from the beginning and is the basis for the linearization of the equations. The nonlinear terms of the momentum flow are also non-negligible compared to the forces of viscous origin, that is:
The scaling of relation (90) results in:
To simplify the writing, we define . Given condition (89), inequality (92) can be satisfied in various ways:
- i.
and . The second condition implies or even , a condition that corresponds to deep waters.
- ii.
and , such that . The condition means that , an approximation of long waves. The condition includes , which is interpreted as a transverse flow with a large Reynolds number but within the laminar regime.
- iii.
and , such that . This case cannot occur because, to have , its denominator must be less than one.
- iv.
and . This case implies and that viscous effects dominate in the transverse flow.
Note that the conditions
and
correspond to the cases analyzed in
Section 4.1.1 and
Section 4.1.2, respectively. The condition for deep waters
, when applied to shallow depths means very short wavelengths and the effect of surface tension should be included in the dynamic boundary condition at the free surface.
5. Discussion
The steady and unsteady flow in ducts with wavy walls driven by pressure gradients as well as the free-surface flow with a bed deformed in the downstream direction have been widely studied. In this study, we have examined free-surface flow, determining the characteristics of the velocity field and the surface when the channel bed deforms over time and in the direction transverse to the flow. To the best of the authors’ knowledge, this has not been studied before.
The pattern of the perturbation component of the velocity field, streamlines, and deformation of the free surface depends on the discriminant of
, given by
, where
. The parameter
is a kind of Reynolds number that represents the fluid inertia induced by bed deformation relative to viscous resistance to motion. Thus, when the viscous effect dominates,
, and the root of
is real. This results in a monotonic decay from the source of motion (bed deformation) to the free surface, generating simple and generalizable patterns of streamlines, velocity components, and surface deformation, as observed in
Figure 4,
Figure 5,
Figure 6 and
Figure 7. In this case, the deformation of the free surface is always in phase with the bed deformation.
More interesting and cumbersome to analyze is the case when fluid inertia dominates over viscous resistance, as indicated by
and the root of
being imaginary. This incorporates trigonometric functions into the solutions, leading to an oscillatory behavior of the perturbation velocity components. Consequently, circulation cells of flow transverse to the main flow direction can be observed, as illustrated in
Figure 10,
Figure 11,
Figure 12 and
Figure 13. With an imaginary
, the solution for the ratio of the amplitudes of free surface and bed deformations (
) is unbounded. However, due to viscous dissipation, the free surface deformation cannot exceed that of the bed, necessitating that
. This limitation constrains the range of
and
to ensure physically meaningful solutions. Depending on the relative importance of inertia to viscosity, the deformation of the free surface may be in phase with the bed deformation or out of phase by
. In the latter case, circulation cells are generated, serving as a way to reconcile the different signs imposed on the vertical velocity component (
) by the deformation of both the free surface and the bed.
In this work, the bed deformation is described by the perturbation with (Equations (37) and (38)). However, the perturbation can be generalized for any function whose variation in the -direction can be expressed in terms of Fourier series, with being a bounded function, while maintaining the restriction that is small.
Finally, it is emphasized that the velocity fields and streamlines obtained here are first-order approximations of the nonlinear equations of motion. Therefore, the results are valid to the extent that the disturbances are small. The removal of this restriction can now be done without major difficulties if the problem is approached through proper numerical modeling, which is free from the constraints imposed by the linearization of equations, enabling the calculation of the velocity field and deformation of the free surface (including the effect of surface tension) and the identification of conditions leading to flow separation in real cases. Furthermore, numerical modeling is not limited to laminar flow, as is the case with most analytical approaches, and can also address turbulent flows. However, in the authors’ judgment, the approximate and linearized analytical solution can offer a clearer and more conceptual initial understanding of the system’s behavior. This is evident in this work through the role played by the parameter , interpreted as a Reynolds number associated with the flow inertia induced by bed deformation. This parameter arises easily and naturally during the nondimensionalization of the (linear) equations governing first-order perturbations. Depending on whether is greater or less than 1, the behavior of the solutions can be very different.