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

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 $$\left(x_1,x_2\right)$$ 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\left(x,t\right)$$, where $$x$$ is the vector with components $$x_1,x_2$$ 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 $$\Pi =P+\rho gz$$ (where $$P$$ is the true pressure) as

$$\{\begin{array}{c}{\partial }_{zz}^2{V}_z+{\partial }_{\beta \beta }^2{V}_z=\frac{1}{\mu }{\partial }_z\Pi ,\\ {\partial }_{zz}^2{V}_{\alpha }+{\partial }_{\beta \beta }^2{V}_{\alpha }=\frac{1}{\mu }{\partial }_{\alpha }\Pi ,\\ {\partial }_z{V}_z+{\partial }_{\beta }{V}_{\beta }=0,\\ {\partial }_{zz}^2\Pi +{\partial }_{\beta \beta }^2\Pi =0.\end{array}$$E1

The indices $$\alpha ,\beta$$ 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

$${\left(-\Pi +2\mu {\partial }_z{V}_z\right)|}_{z=0}=\sigma {\partial }_{\beta \beta }^{2}h-\rho gh,{\left({\partial }_{\alpha }{V}_z+{\partial }_z{V}_{\alpha }\right)|}_{z=0}=0.$$E2

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

$${\left(\begin{array}{c}{V}_z\\ V_{\alpha }\\ \Pi \end{array}\right)|}_{z\to -\infty }=0.$$E3

Here $$V_z$$ and $$V_{\alpha }$$ are the components of velocity, $${\partial }_z$$ and $${\partial }_{\alpha }$$ are the derivatives along the coordinates $$z$$ and $$x_{\alpha }$$ accordingly, $$g$$ is the acceleration due to gravity; $$\mu$$, $$\rho$$ and $$\sigma$$ are the viscosity, density and surface tension coefficient of the liquid.

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

$${\partial }_{t}h={V_z|}_{z=0}+W,{h|}_{t=0}=h^0\left(x\right),$$E4

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

$$f_k={\displaystyle \int f\left(x\right)}{e}^{-ikx}{d}^{2}x$$

we obtain

$$\{\begin{array}{c}{\partial }_{zz}^2{V}_{zk}-k^2{V}_{zk}=\frac{1}{\mu }{\partial }_z{\Pi }_k\\ {\partial }_{zz}^2{V}_{\alpha k}-k^2{V}_{\alpha k}=\frac{i{k}_{\alpha }}{\mu }{\Pi }_k\\ {\partial }_z{V}_{zk}+i{k}_{\beta }{V}_{\beta k}=0\\ {\partial }_{zz}^2{\Pi }_k-k^2{\Pi }_k=0\end{array}$$E5

$${\left(-{\Pi }_k+2\mu {\partial }_z{V}_{zk}\right)|}_{z=0}=-\left(\sigma k^2+\rho g\right)h_k,{\left(i{k}_{\alpha }{V}_{zk}+{\partial }_z{V}_{\alpha k}\right)|}_{z=0}=0,{\left(\begin{array}{c}{V}_{zk}\\ V_{\alpha k}\\ {\Pi }_k\end{array}\right)|}_{z\to -\infty }=0,$$E6

$${\partial }_t{h}_k={V_{zk}|}_{x_{\perp }=0}+W_k.$$E7

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

$$\left(\begin{array}{c}{V}_{zk}\\ V_{\alpha k}\\ {\Pi }_k\end{array}\right)=\left(\begin{array}{c}kz-1\\ i{k}_{\alpha }z\\ 2\mu k\end{array}\right)\frac{\sigma k^2+\rho g}{2\mu k}{h}_k{e}^{kz}.$$E8

It is apparent that the Fourier component $$h_k$$ has an effect in the layer of depth $$k^{-1}$$. Based on (7), (8), the equation of the free surface motion can be written as

$${\partial }_t{h}_k+\frac{1}{2\mu k}\left(\sigma k^2+\rho g\right)h_k=W_k$$E9

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 ($$W_k=0$$). In this case, the solution of equation (9) with initial condition $${h_k|}_{t=0}=h_k^0$$ becomes

$$h_k=h_k^0{e}^{-t/T\left(k\right)},$$E10

where $$T\left(k\right)=\frac{2\mu k}{\sigma k^2+\rho g}$$ is the characteristic relaxation time of the corresponding harmonics. For $$k\to 0$$ and $$k\to \infty$$ the relaxation time tends to zero assuming non-zero gravity and surface tension, respectively. At $$k=k_{*}={\left(\rho g/\sigma \right)}^{1/2}$$ the relaxation time has its maximal value $$T_{*}=\mu /{\left(\rho g\sigma \right)}^{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 $${h_k|}_{k=0}=0$$ at $$t>0$$. Noting that

$${h_k|}_{k=0}={\displaystyle \int h\left(x,t\right)d^{2}x}$$

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

$$\tau =t/T_{*},\theta =T/T_{*},\xi =k_{*}x,\kappa =k/k_{*}.$$

Then (10) takes the form

$$h_{\kappa }=h_{\kappa }^0{e}^{-\tau /\theta \left(\kappa \right)},\theta \left(\kappa \right)=\frac{2\kappa }{{\kappa }^2+1},h_{\kappa }={\displaystyle \int h\left(\xi ,\tau \right)}{e}^{-i\kappa \xi }.$$E11

Applying inverse Fourier transform to (11) we have

$$h\left(\xi ,\tau \right)={\displaystyle \int g\left(\xi -{\xi }^{\prime },\tau \right)}{h}^0\left({\xi }^{\prime }\right)d^2{\xi }^{\prime },g\left(\xi ,\tau \right)={\displaystyle \int e^{i\kappa \xi -\tau /\theta \left(\kappa \right)}\frac{d^2\kappa }{{\left(2\pi \right)}^2}}.$$E12

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

$$h\left(\xi ,\tau \right)={\displaystyle \int e^{i\kappa \xi }{h}_{\kappa }^0{e}^{-\tau /\theta \left(\kappa \right)}\frac{d^2\kappa }{{\left(2\pi \right)}^2}},\theta \left(\kappa \right)=\frac{2\kappa }{k^2+1}.$$E13

From (13) it is clear that the configuration at the moment $$\tau$$ is determined by the Fourier transform of the initial perturbation $$h_{\kappa }^0$$ and the multiplier $$e^{-\tau /\theta \left(\kappa \right)}$$. The maximum value of this factor is equal to $$e^{-\tau }$$, and corresponds to $$\tau =1$$. Let $${\kappa }_{1,2}\left(\tau \right)$$ be the values of $$\kappa$$ at which the multiplier takes a value in $$e$$ times less than the maximum-see Fig. 1.

Figure 1.

Graph of the function $$e^{-\tau /\theta \left(\kappa \right)}$$ at $$\tau =1$$ and its approximation with a step function.

By solving the equation $$\frac{\tau }{\theta \left(\kappa \right)}=\tau +1$$, we obtain

$${\kappa }_{1,2}=\frac{1}{\tau }\left(\tau +1\pm \sqrt{2\tau +1}\right).$$E14

At $$\tau \to 0,{\kappa }_1=\frac{\tau }{2},{\kappa }_2=\frac{2}{\tau }$$. On the other hand, at $$\tau \to \infty ,{\kappa }_{1,2}=1\mp \sqrt{\frac{2}{\tau }}$$.

To estimate perturbation evolution we replace function $$e^{-\tau /\theta \left(\kappa \right)}$$ with the step function $$f\left(\tau ,\kappa \right)$$ which is equal to zero for $$\kappa <{\kappa }_1\left(\tau \right)$$ and for $$\kappa >{\kappa }_2\left(\tau \right)$$. In the interval $${\kappa }_1\left(\tau \right)\le \kappa \le {\kappa }_2\left(\tau \right)$$ this function is equal to $$e^{-\tau }$$-see Fig.1:

$$f\left(\kappa ,\tau \right)=e^{-\tau }\vartheta \left(\kappa -{\kappa }_1\left(\tau \right)\right)\vartheta \left({\kappa }_2\left(\tau \right)-\kappa \right),$$E15

where $$\vartheta \left(x\right)$$ 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 $${\kappa }_1\left(\tau \right)\le \left|\kappa \right|\le {\kappa }_2\left(\tau \right)$$. At the same time the amplitude of the Fourier transform of the perturbation decreases as $$e^{-\tau }$$. 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 $${\kappa }_2$$, then slowly-due to symmetric increasing $${\kappa }_1$$ and decreasing $${\kappa }_2$$. At large times, the ring is localized in the neighborhood of a circle of radius $$\kappa =1$$ and has a characteristic width $$\Delta \kappa =2\sqrt{\frac{2}{\tau }}$$. 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 $${\kappa }_1$$. It leads to the subtraction of additional perturbation localized in a circle with radius $${\xi }_1=\frac{2\pi }{k_1}$$ from the initial perturbation. When $$\tau \to 0,{\xi }_1=\frac{4\pi }{\tau }\to \infty$$. 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 $$r_1=\frac{{\xi }_1}{k_{*}}=\frac{4\pi \mu }{\rho 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 $${\kappa }_2$$ leads to an expansion of its external borders, in the coordinate representation, to a circle of radius $${\xi }_2=\frac{2\pi }{{\kappa }_2}$$. In dimensional variables $$r_2=\frac{{\xi }_2}{k_{*}}=\frac{\pi \sigma t}{\mu }$$, 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

$$h\left(0,\tau \right)={\displaystyle \int h_{\kappa }}\left(\tau \right)\frac{d^2\kappa }{{\left(2\pi \right)}^2}$$

Then

$$h\left(0,\tau \right)={\displaystyle \underset{{\kappa }_1\left(\tau \right)\le \left|\kappa \right|\le {\kappa }_2\left(\tau \right)}{\int }{h}_{\kappa }^0}\frac{d^2\kappa }{{\left(2\pi \right)}^2}$$

The above statements are based on the approximate replacement of real $$e^{-\tau /\theta \left(\kappa \right)}$$ 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 $$\frac{2\pi }{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 $${\kappa }_1\left(\tau \right)$$. In the coordinate representation this additive is concentrated in a circle of radius $${\xi }_1\left(\tau \right)=2\pi /{\kappa }_1\left(\tau \right)$$ and has an amplitude

$$\Delta h\left(\tau \right)={\displaystyle \underset{\left|\kappa \right|\le {\kappa }_1\left(\tau \right)}{\int }{h}_{\kappa }^0}\frac{d^2\kappa }{{\left(2\pi \right)}^2}$$

At the initial stage of relaxation process $${\kappa }_1<<2\pi /a$$, respectively, $${\xi }_1>>a$$,

$$\Delta h\left(\tau \right)=\frac{1}{4\pi }{\kappa }_1^2\left(\tau \right){h_{\kappa }^0|}_{\kappa =0}=\frac{{\tau }^2}{16\pi }{\displaystyle \int h^0\left(\xi \right)}{d}^2\xi$$

Through this long-wavelength addition, the initial perturbation goes down (in case of positive volume) on $$\Delta h\left(\tau \right)$$, practically without distortion of the shape. Outside the region of the initial perturbation arises an additional perturbation with the amplitude $$\Delta h\left(\tau \right)$$ and attenuation radius $${\xi }_1\left(\tau \right)$$. 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 $${\tau }_1=4\pi /a$$, when $${\kappa }_1=2\pi /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 $${\xi }_1$$ ($$a_1>>a_2>>1$$), the region of localization of the Fourier transform has the inverse aspect ratio $$2\pi /a_1<<2\pi /a_2$$. Removal from the spectrum of the initial perturbation the circle with radius $${\kappa }_1\left(\tau \right)=\tau /2$$ leads initially to the formation of a symmetric long-wave perturbation of the radius $$4\pi /\tau$$ and the amplitude $$\Delta h=\frac{{\tau }^2}{16\pi }{\displaystyle \int h^0\left(\xi \right)}{d}^2\xi$$. 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\left(0,\tau \right)=h^0\left(0\right)-\Delta h\left(\tau \right)$$.

Figure 2.

The relaxation of axisymmetric long-wave perturbations. The initial shape of perturbation is $$h^0=A{e}^{-{\xi }^2/a^2},a=20$$.

Figure 3.

Increasing of amplitude of long-wave axisymmetric perturbation. The initial shape of perturbation is $$h^0=A{e}^{-{\xi }^2/a^2},a=20.$$

From the time moment $$\tau =4\pi /a_1$$, when $$k_1\left(\tau \right)=4\pi /a_1$$, the area of influence of long-range perturbations loses circular symmetry. In the direction of the axis $${\xi }_1$$ the scale of the area $${\xi }_1^{*}$$ becomes, and subsequently maintained, equal to $$a_1$$, and along the axis $${\xi }_2$$ this area continues to shrink as $${\xi }_2^{*}=4\pi /\tau$$. The amplitude of the perturbation begins to decrease linearly $$h\left(0,\tau \right)=h^0\left(0\right)-\Delta h\left(\tau \right),\Delta h=\frac{\tau }{4\pi L_1}{\displaystyle \int h^0\left(\xi \right)}{d}^2\xi$$.

At the moment $$\tau =4\pi /a_2$$, when $${\kappa }_1\left(\tau \right)=2\pi /a_2$$ 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 $$\tau =0.04$$. The initial shape of perturbation is $$h^0=A\exp \left(-\frac{{\xi }_1^2}{a_1^2}-\frac{{\xi }_2^2}{a_2^2}\right),$$ $$a_1=100,a_2=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\pi /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 $$\tau =a/\pi$$, when $${\kappa }_2=2\pi /a$$. The size of the region occupied by the perturbation begins to grow as $$\xi \left(\tau \right)=2\pi /{\kappa }_2=\pi \tau$$, and its amplitude decreases at the initial stage, as $$\frac{1}{\pi {\tau }^2}{\displaystyle \int h^0\left(\xi \right)}{d}^2\xi$$. 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 $$h^0=A{e}^{-{\xi }^2/a^2},a=0.1$$.

Finally let’s describe the relaxation of small-scale perturbations with significantly different scales along the axes $${\xi }_1$$ and $${\xi }_2$$ : $$a_1<<a_2<<1$$. The evolution of this perturbation begins at the moment $$\tau =a_1/\pi$$. The scale begins to increase like $$\pi \tau$$ in the $${\xi }_1$$ axis direction and in $${\xi }_2$$ axis direction it varies slightly. The amplitude of the perturbation decreases according to $$h\left(0,\tau \right)=\frac{2}{\pi a_2\tau }{\displaystyle \int h^0\left(\xi \right)}{d}^2\xi$$.

After the moment $$\tau =a_2/\pi$$ 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 $$\pi \tau$$. Amplitude of perturbation thus decreases faster than during the first stage: $$h\left(0,\tau \right)=\frac{1}{\pi {\tau }^2}{\displaystyle \int h^0\left(\xi \right)}{d}^2\xi$$.

The relaxation process ends substantially at the moment $$\tau =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 $$h^0=A\exp \left(-\frac{{\xi }_1^2}{a_1^2}-\frac{{\xi }_2^2}{a_2^2}\right),a_1=0.02,a_2=\mathrm{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

$${V_z|}_{z=-H}=0,{V_{\alpha }|}_{z=-H}=0.$$E16

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

$$\begin{array}{l}{V}_{zk}=\left(V_{zk0}-\frac{H{\Pi }_{k0}}{2\mu }\right)\frac{sh\left(k\left(z+H\right)\right)}{sh\left(kH\right)}+\frac{\left(z+H\right){\Pi }_k}{2\mu },\\ {\Pi }_k={\Pi }_{k0}\frac{sh\left(k\left(z+H\right)\right)}{sh\left(kH\right)}-{\Pi }_{kH}\frac{sh\left(kz\right)}{sh\left(kH\right)},\end{array}$$E17

where

$$\left(\begin{array}{c}{V}_{zk0}\\ {\Pi }_{k0}\\ {\Pi }_{kH}\end{array}\right)=\left(\begin{array}{c}2kH-sh\left(2kh\right)\\ 2\mu k\left(ch\left(2kH\right)+1\right)\\ 4\mu k\left(kHsh\left(kH\right)+ch\left(kH\right)\right)\end{array}\right)\frac{\sigma k^2+\rho g}{2\mu k\left(ch\left(2kH\right)+1+2{\left(kH\right)}^2\right)}{h}_k.$$E18

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

$${\partial }_t{h}_k+\frac{2kH-sh\left(2kh\right)}{ch\left(2kH\right)+1+2{\left(kH\right)}^2}\frac{\sigma k^2+\rho g}{2\mu k}{h}_k=0,{h_k|}_{t=0}=h_k^0.$$E19

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

$$h_k=h_k^0{e}^{-t/T\left(k\right)},$$E20

where $$T\left(k\right)=\frac{ch\left(2kH\right)+1+2{\left(kH\right)}^2}{sh\left(2kh\right)-2kH}\frac{2\mu k}{\sigma k^2+\rho g}$$ is the characteristic relaxation time of corresponding harmonic.

At $$kH>>1$$, $$T\left(k,H\right)=\frac{2\mu k}{\sigma k^2+\rho g}$$ coincides with (10), i.e. with the case of a fluid of infinite depth. Therefore, for $$k\to \infty$$ the relaxation time $$T\left(H,k\right)\to 0$$ due to the surface tension. On the other hand, for $$k\to 0$$, $$T\left(H,k\right)\to \infty$$ unlike the case of a fluid of infinite depth. Then, according to (20),

$${h_k|}_{k=0}={h_k^0|}_{k=0}$$ or $${\displaystyle \int h\left(x,t\right)d^{2}x={\displaystyle \int h^0}\left(x\right)}{d}^{2}x$$, i.e. the “volume of perturbation” is kept. Using the same scales as we used for infinite liquid, we present (20) in dimensionless form

$$h_{\kappa }=h_{\kappa }^0{e}^{-\tau /\theta \left(\kappa ,\overline{H}\right)},\theta \left(\kappa ,\overline{H}\right)=\frac{ch\left(2\kappa \overline{H}\right)+1+2{\left(\kappa \overline{H}\right)}^2}{sh\left(2\kappa \overline{H}\right)-2\kappa \overline{H}}\frac{2\kappa }{{\kappa }^2+1},h_{\kappa }={\displaystyle \int h\left(\xi ,\tau \right)}{e}^{-i\kappa \xi },$$E21

where $$\overline{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 $${\mu }_1$$. Layer is located over the half-space which is filled with liquid viscosity $${\mu }_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:

$$h_{\kappa }=h_{\kappa }^0{e}^{-\tau /\theta \left(\kappa ,\overline{H}\right)},\theta \left(\kappa ,\overline{H},m\right)=\frac{2m\cdot sh\left(2\kappa \overline{H}\right)+\left(m^2+1\right)ch\left(2\kappa \overline{H}\right)+\left(1+2{\kappa }^2{\overline{H}}^2\right)\left(m^2-1\right)}{\left(m^2+1\right)sh\left(2\kappa \overline{H}\right)+2m\cdot ch\left(2\kappa \overline{H}\right)-2\kappa \overline{H}\left(m^2-1\right)}\frac{2\kappa }{{\kappa }^2+1},$$E22

where $$m=\frac{{\mu }_2}{{\mu }_1}$$. As expected, the characteristic relaxation time $$\theta \left(\kappa ,\overline{H},m\right)$$ at $$m=1$$ or $$\overline{H}\to \infty$$ transferred to the corresponding expression (13) obtained for the half-space, at $$m\to \infty$$, (22) transferred to (21)-expression for finite layer over solid bottom. At $$\kappa =0$$, $$\theta \left(\kappa ,\overline{H},m\right)=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^{-\tau /\theta }$$ 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^{-\tau /\theta }$$ at $$\tau =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 $$\overline{H}=1$$.The initial shape of perturbation $$h^0=A{e}^{-{\xi }^2/a^2},a=0.1$$

Figure 9.

Relaxation of axisymmetric perturbation of free surface of the layer with thickness $$\overline{H}=1$$ over a half-space. The initial shape of perturbation $$h^0=A{e}^{-{\xi }^2/a^2},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 $$h^0=A{e}^{-{\xi }^2/a^2},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) $$W_k=\kappa Vcth(k{H}_D)$$, where $$H_D$$ 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

$${\partial }_t{h}_k+\gamma h_k=0,h_k=h_k^0{e}^{\gamma \left(k\right)t},$$E23

where $$\gamma \left(k\right)=-\frac{1}{2\mu k}\left(\sigma k^2+\rho g\right)+kVcth\left(k{H}_D\right)$$. At $$\gamma >0$$ correspondent harmonic is unstable, and at $$\gamma <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 $$\gamma \left(k\right)$$ 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 $$\gamma \left(k\right)$$ 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.

Acknowledgments

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