Stokes steady-state approximation  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 . Many classical problems of high-viscous liquids hydrodynamics were solved in this approximation: the destruction of bubbles , the formation of the cusp on the liquid surface [4,5], the coalescence of liquid particles  and some problems of nanofluidics .
The idea of using inertia-less approximation to calculate wave relaxation on the surface of a fluid with high viscosity was suggested by . The relaxation of harmonic perturbations in this approximation was considered by [9,10]. This approach has also been used in studying thin film flow , 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 located on the unperturbed surface and the axis directed along the external normal to it. Then the perturbation of the surface can be defined as , where is the vector with components and is the time. We express the equations of motion in the Stokes approximation , the continuity equation, the equation for the renormalized pressure (where 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 , have the form
The following conditions in case of infinite half-space must also be met:
Here and are the components of velocity, and are the derivatives along the coordinates and accordingly, 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 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
Integrating (5) with boundary conditions (6), we obtain
It is apparent that the Fourier component has an effect in the layer of depth . 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 (). In this case, the solution of equation (9) with initial condition becomes
where is the characteristic relaxation time of the corresponding harmonics. For and the relaxation time tends to zero assuming non-zero gravity and surface tension, respectively. At the relaxation time has its maximal value . 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 , therefore at . 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 .
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 and the multiplier . The maximum value of this factor is equal to , and corresponds to . Let be the values of at which the multiplier takes a value in times less than the maximum-see Fig. 1.
By solving the equation , we obtain
At . On the other hand, at .
To estimate perturbation evolution we replace function with the step function which is equal to zero for and for . In the interval this function is equal to -see Fig.1:
where 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 . At the same time the amplitude of the Fourier transform of the perturbation decreases as . 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 , then slowly-due to symmetric increasing and decreasing . At large times, the ring is localized in the neighborhood of a circle of radius and has a characteristic width . 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 . It leads to the subtraction of additional perturbation localized in a circle with radius from the initial perturbation. When . 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 . 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 leads to an expansion of its external borders, in the coordinate representation, to a circle of radius . In dimensional variables , 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 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 , then the Fourier transform of this disturbance is concentrated in the region with a typical size . 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 . In the coordinate representation this additive is concentrated in a circle of radius and has an amplitude
At the initial stage of relaxation process , respectively, ,
Through this long-wavelength addition, the initial perturbation goes down (in case of positive volume) on , practically without distortion of the shape. Outside the region of the initial perturbation arises an additional perturbation with the amplitude and attenuation radius . 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 , when , 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 (), the region of localization of the Fourier transform has the inverse aspect ratio . Removal from the spectrum of the initial perturbation the circle with radius leads initially to the formation of a symmetric long-wave perturbation of the radius and the amplitude . 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 .
From the time moment , when , the area of influence of long-range perturbations loses circular symmetry. In the direction of the axis the scale of the area becomes, and subsequently maintained, equal to , and along the axis this area continues to shrink as . The amplitude of the perturbation begins to decrease linearly .
At the moment , when 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.
2.5. Relaxation of short-wave perturbations due to capillary forces
The evolution of perturbation with the scale is different. The Fourier transform of this disturbance is concentrated in the region with a typical size of . 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 , when . The size of the region occupied by the perturbation begins to grow as , and its amplitude decreases at the initial stage, as . Characteristic decay time of the perturbation is of order unity. Fig. 5 shows the relaxation of short-wave axisymmetric perturbation.
Finally let’s describe the relaxation of small-scale perturbations with significantly different scales along the axes and : . The evolution of this perturbation begins at the moment . The scale begins to increase like in the axis direction and in axis direction it varies slightly. The amplitude of the perturbation decreases according to .
After the moment 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: .
The relaxation process ends substantially at the moment . Fig.5 shows relaxation of short-wave non-axisymmetric perturbations. One can see that over time the perturbation becomes axisymmetric.
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 is the characteristic relaxation time of corresponding harmonic.
At , coincides with (10), i.e. with the case of a fluid of infinite depth. Therefore, for the relaxation time due to the surface tension. On the other hand, for , unlike the case of a fluid of infinite depth. Then, according to (20),
or , 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 is the dimensionless thickness of layer.
3.2. Finite layer over a half-space
Consider a plane layer with thickness , filled with liquid with viscosity . Layer is located over the half-space which is filled with liquid viscosity . 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 . As expected, the characteristic relaxation time at or transferred to the corresponding expression (13) obtained for the half-space, at , (22) transferred to (21)-expression for finite layer over solid bottom. At , . 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 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.
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.
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).
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.
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) , where is the thickness of the diffusion boundary layer, 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 . At correspondent harmonic is unstable, and at 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 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.
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.