Open access peer-reviewed chapter

Evolution of Small Perturbations of the Free Surface of Viscous Fluid in the Stokes Approximation

By Tamara Pogosian, Ivan Melikhov, Anastasia Shutova, Alexey Amosov and Sergey Chivilikhin

Submitted: April 30th 2014Reviewed: September 14th 2014Published: April 1st 2015

DOI: 10.5772/59226

Downloaded: 995

1. Introduction

Stokes steady-state approximation [1] is widely used for calculation of viscous flows with low Reynolds numbers. This approach was first applied to describe the flow with free boundaries in [2]. Many classical problems of high-viscous liquids hydrodynamics were solved in this approximation: the destruction of bubbles [3], the formation of the cusp on the liquid surface [4,5], the coalescence of liquid particles [6] and some problems of nanofluidics [7].

The idea of using inertia-less approximation to calculate wave relaxation on the surface of a fluid with high viscosity was suggested by [8]. The relaxation of harmonic perturbations in this approximation was considered by [9,10]. This approach has also been used in studying thin film flow [11], two-dimensional Stokes flow with free boundary [12-13].

In this work we study an aperiodic relaxation of a small localized perturbation of the planar surface of an incompressible fluid under the influence of gravity and surface tension. The description is given in the Stokes approximation. It is shown that this imposes limitations on the three-dimensional spectrum of the considered perturbations. An equation describing perturbation damping in Fourier space was obtained. It is shown that the volume of the perturbation becomes zero in the time that is small by comparison to the time of perturbation damping. In the short wave limit, the law of evolution of the perturbation permits a simple geometric interpretation.

The second part of the work is focused on theoretical investigation of instability of the surface shape under the influence of diffusion mass flux through the surface. A local eminence on the surface can lead to perturbation growth. Stability and instability of the perturbation will depend on which mechanism prevails. Diffusion causes distortion of surface shape, gravitational and capillary forces have the opposite effect.

2. Evolution of half-space free boundary perturbation

2.1. The general equations

We consider a half-space filled with a highly viscous incompressible Newton liquid of constant viscosity. At the initial moment of time the planar surface of half-space obtains a small local perturbation. We will investigate further the evolution of this perturbation due to gravitation, capillary forces and external mass flux through the free surface. We introduce the Cartesian coordinate system with axes (x1,x2)located on the unperturbed surface and the axis zdirected along the external normal to it. Then the perturbation of the surface can be defined as z=h(x,t), where xis the vector with components x1,x2and tis the time. We express the equations of motion in the Stokes approximation [1], the continuity equation, the equation for the renormalized pressure Π=P+ρgz(where Pis the true pressure) as


The indices α,βcan be 1 or 2; summation over repeated indices is assumed. The boundary conditions which express the continuity of the normal and tangential stresses on the free surface, in the approximation linear by h, have the form


The following conditions in case of infinite half-space must also be met:


Here Vzand Vαare the components of velocity, zand αare the derivatives along the coordinates zand xαaccordingly, gis the acceleration due to gravity; μ, ρand σare the viscosity, density and surface tension coefficient of the liquid.

Equation of free surface evolution and its initial condition have the form


where Wdescribes the influence of external mass fluxes on the velocity of free boundary.

Applying to (1)-(4) the two-dimensional Fourier transform with respect to the longitudinal coordinates


we obtain


Integrating (5) with boundary conditions (6), we obtain


It is apparent that the Fourier component hkhas an effect in the layer of depth k1. Based on (7), (8), the equation of the free surface motion can be written as


2.2. Relaxation of the initial perturbation of the liquid surface, exact solution

Consider the relaxation of the initial perturbation of the surface due to capillary and gravitation forces in the absence of mass flux at the boundary (Wk=0). In this case, the solution of equation (9) with initial condition hk|t=0=hk0becomes


where T(k)=2μkσk2+ρgis the characteristic relaxation time of the corresponding harmonics. For k0and kthe relaxation time tends to zero assuming non-zero gravity and surface tension, respectively. At k=k*=(ρg/σ)1/2the relaxation time has its maximal value T*=μ/(ρgσ)1/2. Small relaxation time for high-frequency harmonics means that the surface tension forces lead to a rapid smoothing of sharp irregularities. On the other hand T|k=0=0, therefore hk|k=0=0at t>0. Noting that


constitutes the "volume of perturbation". According to (10) this value becomes equal to zero due to the gravity force. Immediacy of this process is determined using a quasi-stationary Stokes approximation. Detailed discussion of this topic see in [14].

Switching to the case of arbitrary wave numbers, we introduce the dimensionless variables


Then (10) takes the form


Applying inverse Fourier transform to (11) we have


The expression (12) allows calculation of the surface perturbation in any time for an arbitrary initial perturbation.

2.3. Qualitative description of the perturbation relaxation

For a qualitative description of the relaxation of the initial perturbation is convenient to use the original expression (11). We have


From (13) it is clear that the configuration at the moment τis determined by the Fourier transform of the initial perturbation hκ0and the multiplier eτ/θ(κ). The maximum value of this factor is equal to eτ, and corresponds to τ=1. Let κ1,2(τ)be the values of κat which the multiplier takes a value in etimes less than the maximum-see Fig. 1.

Figure 1.

Graph of the function e−τ/θ(κ) at τ=1 and its approximation with a step function.

By solving the equation τθ(κ)=τ+1, we obtain


At τ0,κ1=τ2,κ2=2τ. On the other hand, at τ,κ1,2=12τ.

To estimate perturbation evolution we replace function eτ/θ(κ)with the step function f(τ,κ)which is equal to zero for κ<κ1(τ)and for κ>κ2(τ). In the interval κ1(τ)κκ2(τ)this function is equal to eτ-see Fig.1:


where ϑ(x)is the Heaviside step function.

Returning to the expression (13) we see that the evolution of perturbations can be considered as a result the initial perturbation “cut-off” in Fourier space. The system under consideration behaves as filter which "cuts" in the spectrum of the initial perturbation wave numbers the annular region with characteristic dimensions κ1(τ)|κ|κ2(τ). At the same time the amplitude of the Fourier transform of the perturbation decreases as eτ. It can be seen that the width of the annular region decreases with increasing time, rapidly at the beginning-mainly due to the reduction of the outer radius of the ring κ2, then slowly-due to symmetric increasing κ1and decreasing κ2. At large times, the ring is localized in the neighborhood of a circle of radius κ=1and has a characteristic width Δκ=22τ. Such filtering process determines the qualitative features of the relaxation of the initial perturbation.

Spatial filtering of the Fourier transform of the initial perturbation in the region of small wave numbers is equivalent to subtracting its part localized in a circle of radius κ1. It leads to the subtraction of additional perturbation localized in a circle with radius ξ1=2πk1from the initial perturbation. When τ0,ξ1=4πτ. This long-range effect of the initial perturbation is determined by the force of gravity. This is clearly seen in the transition to dimensional movements: the radius of the long-range perturbations is r1=ξ1k*=4πμρgt. Figuratively speaking, the perturbation is drowning in the ambient fluid, distorting the shape of the surface at long distance.

Localization of the Fourier transform of the perturbation in a circle of radius κ2leads to an expansion of its external borders, in the coordinate representation, to a circle of radius ξ2=2πκ2. In dimensional variables r2=ξ2k*=πσtμ, i.e. this effect is due to the action of surface tension forces.

To evaluate the value of the perturbation at the origin of coordinates we use of the fact that




The above statements are based on the approximate replacement of real eτ/θ(κ)functions with the step function (15). Therefore, we can only expect an approximate correspondence between the description, based on these statements, and the actual process of the perturbation relaxation.

Let us consider the relaxation of perturbations with different spatial scale and shape. To illustrate the qualitative dependencies we will use graphics produced on the basis of the exact solution (12).

2.4. Relaxation long-wave perturbations by gravity

If the characteristic dimensionless scale of the region occupied by the perturbation a>>1, then the Fourier transform of this disturbance is concentrated in the region with a typical size 2πa<<1. At the initial stage the relaxation of such perturbation is mainly determined by filtering out small wave numbers of its Fourier representation, i.e., by gravity. In this case, as shown above, to the initial perturbation is added, taken with the opposite sign, the Fourier inverse transform of the initial perturbation, lying in a circle of radius κ1(τ). In the coordinate representation this additive is concentrated in a circle of radius ξ1(τ)=2π/κ1(τ)and has an amplitude


At the initial stage of relaxation process κ1<<2π/a, respectively, ξ1>>a,


Through this long-wavelength addition, the initial perturbation goes down (in case of positive volume) on Δh(τ), practically without distortion of the shape. Outside the region of the initial perturbation arises an additional perturbation with the amplitude Δh(τ)and attenuation radius ξ1(τ). At large times, the shape of additional perturbation is somewhat distorted, its radius decreases and becomes comparable to the size of the initial perturbation.

By the time τ1=4π/a, when κ1=2π/a, the process of relaxation of the perturbation is mainly completed. To illustrate these effects on the Fig.2 presented the relaxation of axisymmetric long-wave perturbation shape. On the Fig.3 we can see the increasing of this perturbation amplitude.

If the characteristic dimensions of the initial perturbation in two mutually perpendicular directions are very different, the relaxation process is more complicated. Note that for the perturbation, elongated along an axis ξ1(a1>>a2>>1), the region of localization of the Fourier transform has the inverse aspect ratio 2π/a1<<2π/a2. Removal from the spectrum of the initial perturbation the circle with radius κ1(τ)=τ/2leads initially to the formation of a symmetric long-wave perturbation of the radius 4π/τand the amplitude Δh=τ216πh0(ξ)d2ξ. At this stage, the difference of the transverse scales the perturbation does not affect practically the nature of its evolution. Perturbation settles as a whole, and its amplitude decreases according to h(0,τ)=h0(0)Δh(τ).

Figure 2.

The relaxation of axisymmetric long-wave perturbations. The initial shape of perturbation is h0=Ae−ξ2/a2, a=20.

Figure 3.

Increasing of amplitude of long-wave axisymmetric perturbation. The initial shape of perturbation is h0=Ae−ξ2/a2, a=20.

From the time moment τ=4π/a1, when k1(τ)=4π/a1, the area of influence of long-range perturbations loses circular symmetry. In the direction of the axis ξ1the scale of the area ξ1*becomes, and subsequently maintained, equal to a1, and along the axis ξ2this area continues to shrink as ξ2*=4π/τ. The amplitude of the perturbation begins to decrease linearly h(0,τ)=h0(0)Δh(τ),Δh=τ4πL1h0(ξ)d2ξ.

At the moment τ=4π/a2, when κ1(τ)=2π/a2the perturbation is practically relaxed. Figure 4 shows the characteristic configuration of the long-wavelength perturbation, the extent of which the two axes are substantially different.

Figure 4.

The configuration of the long-wave non-axisymmetric perturbations in the time τ=0.04. The initial shape of perturbation is h0=Aexp(−ξ12a12−ξ22a22), a1=100, a2=20.

2.5. Relaxation of short-wave perturbations due to capillary forces

The evolution of perturbation with the scale a<<1is different. The Fourier transform of this disturbance is concentrated in the region with a typical size of 2π/a>>1. In this case occurs the evolution of the perturbation in the initial stage due to filtration of its Fourier transform through small wave numbers. However, the contribution of this region to the total amplitude of the perturbation is small, and therefore, the effects of subsidence of the perturbation and the formation of the region of the far influence, in this case, is irrelevant.

Relaxation of short-wave perturbations is mainly determined by the action of surface tension forces, i.e. filtration its Fourier transform by the large wave numbers. The relaxation process begins at the moment τ=a/π, when κ2=2π/a. The size of the region occupied by the perturbation begins to grow as ξ(τ)=2π/κ2=πτ, and its amplitude decreases at the initial stage, as 1πτ2h0(ξ)d2ξ. Characteristic decay time of the perturbation is of order unity. Fig. 5 shows the relaxation of short-wave axisymmetric perturbation.

Figure 5.

Relaxation of short-wave axisymmetric perturbation. The initial shape of perturbation h0=Ae−ξ2/a2, a=0.1.

Finally let’s describe the relaxation of small-scale perturbations with significantly different scales along the axes ξ1and ξ2: a1<<a2<<1. The evolution of this perturbation begins at the moment τ=a1/π. The scale begins to increase like πτin the ξ1axis direction and in ξ2axis direction it varies slightly. The amplitude of the perturbation decreases according to h(0,τ)=2πa2τh0(ξ)d2ξ.

After the moment τ=a2/πthe magnitudes of the perturbation in the direction of the two axes become comparable. From this moment the perturbation becomes axially symmetric, and its radius increases as πτ. Amplitude of perturbation thus decreases faster than during the first stage: h(0,τ)=1πτ2h0(ξ)d2ξ.

The relaxation process ends substantially at the moment τ=1. Fig.5 shows relaxation of short-wave non-axisymmetric perturbations. One can see that over time the perturbation becomes axisymmetric.

Figure 6.

Configurations of the short-wave non-axisymmetric perturbation. The initial shape of perturbation is h0=Aexp(−ξ12a12−ξ22a22), a1=0.02, a2=0.2.

3. Finite layer over solid bottom and over a half-space

3.1. Finite layer over solid bottom

In this section we study surface perturbations of a flat liquid layer with depth H resting on a solid bottom. We use system (1) and conditions on the free boundary (2) to calculate perturbation evolution due to gravity and surface tension. The conditions at infinity (3) are replaced by the no-slip condition along the solid surface


Again applying the Fourier transform in the longitudinal coordinates and solving the resulting system of ordinary differential equations, we obtain




The evolution equation of the free surface and the initial condition has the form


The solution of equation (19) is similar to (10)


where T(k)=ch(2kH)+1+2(kH)2sh(2kh)2kH2μkσk2+ρgis the characteristic relaxation time of corresponding harmonic.

At kH>>1, T(k,H)=2μkσk2+ρgcoincides with (10), i.e. with the case of a fluid of infinite depth. Therefore, for kthe relaxation time T(H,k)0due to the surface tension. On the other hand, for k0, T(H,k)unlike the case of a fluid of infinite depth. Then, according to (20),

hk|k=0=hk0|k=0or h(x,t)d2x=h0(x)d2x, i.e. the “volume of perturbation” is kept. Using the same scales as we used for infinite liquid, we present (20) in dimensionless form


where H¯=k*His the dimensionless thickness of layer.

3.2. Finite layer over a half-space

Consider a plane layer with thickness H, filled with liquid with viscosity μ1. Layer is located over the half-space which is filled with liquid viscosity μ2. At the interface between the layer and half-space we require continuity of the velocity components and continuity of the convolution of the stress tensor with the normal vector. Taking into account the gravitational forces and capillary forces on the free surface of the layer, using the technique described above, we obtain the equation of evolution of the free surface of the layer:


where m=μ2μ1. As expected, the characteristic relaxation time θ(κ,H¯,m)at m=1or H¯transferred to the corresponding expression (13) obtained for the half-space, at m, (22) transferred to (21)-expression for finite layer over solid bottom. At κ=0, θ(κ,H¯,m)=0. Consequently, the "volume of perturbation" instantly vanishes, as in the case of the half-space.

3.3. Comparing the cases of finite layer and finite layer over a half-space

Fig.7 shows function eτ/θfor the cases of finite layer and finite layer over a half-space. Qualitative differences in the behavior of this function lead to a difference in the perturbation decay for these cases.

Figure 7.

Graph of the function e−τ/θ at τ=1 for half-space, finite layer and layer with half-space.

Fig. 8 and Fig. 9 show decrease of the amplitude for the perturbation with the same initial characteristics. The figures show that the perturbation on the surface layer over solid bottom decays faster than the perturbation on the surface layer over the half-space.

Figure 8.

Relaxation of axisymmetric perturbation of free surface of the layer with thickness H¯=1.The initial shape of perturbation h0=Ae−ξ2/a2, a=0.1

Figure 9.

Relaxation of axisymmetric perturbation of free surface of the layer with thickness H¯=1 over a half-space. The initial shape of perturbation h0=Ae−ξ2/a2, a=0.1, m=2.

3.4. Comparing the relaxation of finite layer free surface perturbation in linear and non-linear approximation

Using finite-element software Comsol Multiphysics we developed a 2D axisymmetric model of the viscous flow to calculate transient relaxation of the perturbation on the free surface and compare it with analytical results reported in the previous sections. We used “Laminar Two-Phase Flow, Moving Mesh” module to enable deformation of the computational domain during the solution. Full Navier-Stokes equations are solved on the moving mesh to describe deformation of the domain due to capillary forces acting along the free surface.

Computational mesh was refined near the perturbation to resolve high curvature and velocity gradients (see Fig. 10).

Figure 10.

Initial non-deformed finite-element discretization of the domain. The initial shape of perturbation h0=Ae−ξ2/a2, A=1, a=1. Boundary conditions: 1 – axial symmetry, 2 – free surface, 3 – slip, and 4 – no-slip wall.

Multiple cases were calculated corresponding to small and large non-dimensional perturbation amplitude and width. The evolution of the perturbation amplitude with time was compared for numerical and analytical models (Fig. 11). It can be seen that the analytics and FEM are in excellent agreement for small amplitude A=0.1 for both small and large width (a=0.1, a=1.0). Some discrepancy is observed in the case of large amplitude A=1.0 due to violation of the analytic model assumptions.

Figure 11.

Comparison of the analytical model (solid lines) and FEM (points) results

4. Evolution of small perturbations of half-space free boundary due to external mass flux, gravitation and capillary forces

Let’s describe the influence of external mass flow on stability of the free surface of half-space. A small perturbation on the free surface leads to a perturbation of the diffusion boundary layer above it. The corresponding mass flow perturbations on the free boundary defines the source term in the right hand side of (9) Wk=κVcth(kHD), where HDis the thickness of the diffusion boundary layer, Vis the unperturbed velocity of the liquid level rise due to external mass flow. In this case the equation of perturbation evolution has the following shape


where γ(k)=12μk(σk2+ρg)+kVcth(kHD). At γ>0correspondent harmonic is unstable, and at γ<0it is stable. On the Fig.12 we can see the regions in the space of parameters, correspondent stable and unstable regimes. On the Fig.13 present the graphs of γ(k)for these regions. In the case (a) all the harmonics with the wave number exceeding some critical value are unstable. In the case (b) all harmonics are stable. In the case (c) there is a range of wave numbers corresponding to unstable harmonics. In the case (d) the harmonics in a certain range of wave numbers, and with a wave number exceeding a critical value are unstable.

Figure 12.

Space of parameters with regions, corresponds to regimes a-unstable, b-stable, c-partly unstable, d-mixed.

Figure 13.

Graphs of γ(k) for regimes: a-unstable, b-stable, c-partly unstable, d-mixed

5. Conclusion

The relaxation process of a small perturbation of the free planar surface of a viscous incompressible fluid under the action of capillary and gravitational forces was investigated. The case of half-space, layer over hard bottom and layer over the half-space was analyzed. We propose a method for qualitative analysis of the configuration of the perturbation based on spatial filtering the Fourier transform of the initial perturbation.

It is shown that under the gravity force there is an additional long-wave perturbation arising around the initial perturbation which provides a zero total volume of the perturbation. Capillary forces lead to an increase of the area affected by the perturbation. The influence of amplitude of perturbation on the rate of its relaxation was investigated numerically and compared with the analytics.

The stability of the free surface of the half-space with the mass flux at the boundary was studied.


This work was financially supported by Government of Russian Federation, Grant 074-U01.

© 2015 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tamara Pogosian, Ivan Melikhov, Anastasia Shutova, Alexey Amosov and Sergey Chivilikhin (April 1st 2015). Evolution of Small Perturbations of the Free Surface of Viscous Fluid in the Stokes Approximation, Hydrodynamics - Concepts and Experiments, Harry Edmar Schulz, IntechOpen, DOI: 10.5772/59226. Available from:

chapter statistics

995total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Development of a Coupled Fluid-Structure Model with Application to a Fishing Net in Current

By Chunwei Bi, Yunpeng Zhao and Guohai Dong

Related Book

First chapter

One Dimensional Turbulent Transfer Using Random Square Waves – Scalar/Velocity and Velocity/Velocity Interactions

By H. E. Schulz, G. B. Lopes Júnior, A. L. A. Simões and R. J. Lobosco

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us