Open access

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

Written By

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

Submitted: April 30th, 2014 Published: April 1st, 2015

DOI: 10.5772/59226

Chapter metrics overview

1,609 Chapter Downloads

View Full Metrics

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 z directed along the external normal to it. Then the perturbation of the surface can be defined as z=h(x,t), where x is the vector with components x1,x2 and t is 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 P is 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 Vz and Vα are the components of velocity, z and α are the derivatives along the coordinates z and xα accordingly, g is 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 W describes 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 hk has 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=hk0 becomes


where T(k)=2μkσk2+ρg is the characteristic relaxation time of the corresponding harmonics. For k0 and k the relaxation time tends to zero assuming non-zero gravity and surface tension, respectively. At k=k*=(ρg/σ)1/2 the 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=0 at 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κ0 and 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 e times 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 κ1 and decreasing κ2. At large times, the ring is localized in the neighborhood of a circle of radius κ=1 and 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πk1 from 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 κ2 leads 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(τ)=τ/2 leads 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 ξ1 the scale of the area ξ1* becomes, and subsequently maintained, equal to a1, and along the axis ξ2 this 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π/a2 the 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<<1 is 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 ξ1 and ξ2 : a1<<a2<<1. The evolution of this perturbation begins at the moment τ=a1/π. The scale begins to increase like πτ in the ξ1 axis direction and in ξ2 axis 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+ρg is the characteristic relaxation time of corresponding harmonic.

At kH>>1, T(k,H)=2μkσk2+ρg coincides with (10), i.e. with the case of a fluid of infinite depth. Therefore, for k the relaxation time T(H,k)0 due 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=0 or 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*H is 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=1 or 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 HD is the thickness of the diffusion boundary layer, V is 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 γ>0 correspondent harmonic is unstable, and at γ<0 it 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.


  1. 1. Happel S.J., Brenner H. Low Reynolds Number Hydrodynamics. Englewood Cliffs: Prentice-Hall; 1965.
  2. 2. Frenkel, J. Viscous flow of crystalline bodies under the action of surface tension. Journal of Physics 1945;9(5) 385-391.
  3. 3. Tanveer S., Vasconcelos G.L. Bubble Breakup in two-dimensional Stokes flow. Physical Review Letters 1994;73(21) 2845-2848.
  4. 4. Jeong J.-T., Moffatt H.K. Free-surface cusps associated with flow at low Reynolds Number. Journal of Fluid Mechanics 1992;241 1-22.
  5. 5. Pozrikidis C. Numerical studies of singularity formation at free surfaces and fluid interfaces in two-dimensional Stokes flow. Journal of Fluid Mechanics 1997;331 145-167.
  6. 6. Hopper R.W. Coalescence of two equal cylinders: exact results for creeping viscous plane flow driven by capillarity. Journal of the American Ceramic Society 1984;67(12) 262-264.
  7. 7. Rodygina O.A., Chivilikhin S.A., Popov I.Yu., Gusarov V.V. Crystallite model for flow in nanotube caused by wall soliton. Nanosystems: physics, chemistry, mathematics 2014;5(3) 400–404.
  8. 8. Lamb H. Hydrodynamics. Cambridge: University Press: 1932.
  9. 9. Levich V.G. Physicochemical hydrodynamics. Englewood Cliffs: Prentice-Hall; 1962.
  10. 10. Voeght F.D., Joos, P. Damping of a disturbance on a liquid surface. Journal of Colloid and Interface Science 1884;98(1) 2--32.
  11. 11. Voinov O.V. Quasi-steadi-state flow in a liquid film in a gas: comparison of two methods of describing waves. Journal of Applied Mechanics and Technical Physics 2007;48(3) 385-392.
  12. 12. Antonovskii L.K. Interface boundary dynamics under the action of capillary forces. Quasisteady-state plane-parallel motion. Journal of Applied Mechanics and Technical Physics 1988; 29(3) 396-399.
  13. 13. Chivilikhin S.A. Plane capillary flow of a viscous fluid with multiply connected boundary in the Stokes approximation. Fluid Dynamics 1992;27(1) 88-92.
  14. 14. Chivilikhin S.A. Relaxation of a small local perturbation of the surface of a viscous fluid in the Stokes approximation. Fluid Dynamics 1985;20(3) 450-454.

Written By

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

Submitted: April 30th, 2014 Published: April 1st, 2015