Accepted: December 09, 2017 | Published Online: December 11, 2017

Citation: Sarychev V, Nevskii S, Alexey G, Chinakhov D (2017) Simulation of Nanoparticles Formation by Mechanism of Kelvin-Helmholtz Instability. Int J Nanoparticles Nanotech 3:012.

Copyright: © 2017 Sarychev V, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

The Kelvin-Helmholtz instability is analyzed at scales of micro-and nano-structures. We have studied the stability of a plane motionless layer of a two-layer incompressible viscous fluid using the Navier-Stokes equations for linear and nonlinear analyses. The effects of viscosity were assumed to occur at an interface with a flow inside the layers to be irrotational. A dispersion equation was developed for small perturbations, which is similar in form to the dispersion equation appeared in earlier papers, provided the short-wave approximation is applied. We first determined: 1) The dependence of the perturbation decrement for a viscous two-layer fluid that has two maxima, with the first maximum being within the wavelengths ranging from 100 to 300 nm and the second from 1 to 3 µm; 2) The approximate analytic dependence of wavenumber that contains the maximum for the perturbation decrement on input parameters for a problem (densities of both layers, their thicknesses, viscosity, surface tension, velocity of layer motion, a coefficient of resistance). It allows us to size up the emerging vortex structures under different conditions. The range of input parameters for the problem was determined, where two maxima relating to the dependence of disturbance decrement were observed. To verify the results after performing linear analysis, the Level Set Method was used for analyzing nonlinear equations. The results of calculations prove that linear analysis adequately describes how vortex structures in various sizes are formed.

Keywords

Kelvin-Helmholtz instability, Navier-Stokes equations, Level set method

Introduction

Hydrodynamic concept to determine the mechanism for forming nanostructured elements under external actions was first used when studying effects of heterogeneous plasma flows, created by electric explosion of conductors, on metals [1]. This work introduces a hypothesis that nanostructures in a metal are formed because of the Kelvin-Helmholtz instability. An ideal fluid models a layer of plasma and a viscous fluid models a molten layer. Using the Navier-Stokes equations for linear analysis in hydrodynamics relating to small perturbations allowed us to obtain a dispersion equation, with its numerical solution making possible to obtain the dependence of increment on wavelength. The maximum in this dependence is achieved at a wavelength λ_m. With certain values of parameters for the problem (density, layer thickness, viscosity, tension coefficients and relative velocity), the wavelength λ_m is in the nanometer range. This means that the waves within this wavelength region will progress, and after the molten layer solidifies these waves can be observed using an electron microscope. Photographs representing the characters of metal surfaces are given in [1]. Thus, it can be said that the proposed model for forming nanostructures allows us to obtain, in principle, the dependence of critical wavelength on the parameters. However, the dispersion equation is complicated by allowing for viscosity, even while considering one layer and, therefore, the parameterization can be carried out only in terms of numbers in the case of semi-infinite layers. The finiteness of layers taken into account leads to a more complicated dispersion equation [2]. Under the assumption of a viscous potential flow, a simplified dispersion equation was obtained in [3], but its parametric analysis in these and subsequent works of these authors was not carried out. The work in [2] gives a numerical analysis of the complete equation and an analytic-graphical analysis of the dispersion equation under the assumption of the viscous potential flow. The dependence of the increment on the wavelength having two maxima was found for certain parameters. This suggests that even in the first linear approximation there can exist two scales. With a purpose of obtaining a simplified dispersion equation, the case of short waves was considered in [4,5]. In this approximation, the dispersion equation was derived that coincides with the equation for the visco-potential model [3]. It is shown that from the obtained equation it is possible to explicit the increment as a function of wavenumber and parameters that facilitates to find the conditions required for forming two maxima, and to give approximate expressions for the two values of λ_m. In our opinion, this allows us to carry out parameterization in various physical conditions with the purpose of identifying conditions for the formation of nanostructures. Therefore, an important role in our study plays a visco-potential flow.

The fluid mechanics theory of potential flow goes back to Euler in 1761. The concept of viscosity was not known in Euler's time. The effects of viscous stresses were introduced by Navier (1822) and Stokes (1845). Stokes (1851) considered potential flow of a viscous fluid in an approximate sense, but later authors restrict their attention to "potential flow of an inviscid fluid". All the books on fluid mechanics and all courses in fluid mechanics have chapters on "potential flow of inviscid fluids" and none on the "potential flow of a viscous fluid" or "Viscous Potential Flow" (VPF). An authoritative and readable exposition of irrotational flow theory and its applications can be found in chapter 6 of the book on fluid dynamics [6]. He speaks of the role of the theory of flow of an inviscid fluid: "various aspects of the flow of a fluid regarded as entirely inviscid (and incompressible) will be considered. The results presented are significant only in as much as they represent an approximation to the flow of a real fluid at large Reynolds number and the limitations of each result must be regarded as important as the result itself". In book [7] has considered irrotational flows of a viscous fluid. That is of the opinion that when one is considering irrotational solutions of the Navier-Stokes equations it is never necessary and typically not useful for one to put the viscosity to zero. This observation runs counter to the idea frequently expressed that potential flow is a topic that is useful only for inviscid fluids; many people think that the notion of a viscous potential flow is an oxymoron. Incorrect statements like "irrotational flow implies inviscid flow but not the other way around" can be found in popular textbooks. The [3] argues that the Navier-Stokes equations are satisfied by potential flow; the viscous term is identically zero when the vorticity is zero but the viscous stresses are not zero. It is not possible to satisfy the no-slip condition at a solid boundary or the continuity of the tangential component of velocity and shear stress at a fluid-fluid boundary when the velocity is given by a potential. KH instability is included by a discontinuity of the velocity at a two-fluid interface. This discontinuity is inconsistent with the no-slip condition for Navier-Stokes studies of viscous fluids, but is consistent with the theory of potential flow of a viscous fluid. Viscous potential flow analysis gives good approximations to fully viscous flows in cases where the shears from the gas flow are negligible. In addition, in [4], it is shown that in the case of short waves the dispersion equation for the viscous layer coincides with the dispersion equation of the VPF. Then analysis is based on the viscous potential flow theory in this paper.

We can observe nanostructured elements in materials not only under the influence of concentrated energy flows, but also intense plastic deformation. An example is equal-channel angular pressing [8], shift under pressure [9], friction and wear processes [10,11] and long-term service life of rails [12,13]. We pay special attention to the papers in [10,11]. In these works, the shear instability of subsurface layers in friction was studied. With the help of scanning electron microscopy [10] it was established that, as moving from the depth of metal to the surface of friction, several characteristic zones can be distinguished: The zone of plastic deformation and texturing (I), the zone of intense fragmentation (II), the zone of turbulent flow (III), and the zone of laminar flotation (VI). The authors [10] designate Zones I and II as the zones of ordinary plastic deformation, whereas Zones III and IV, in their opinion, relate to the zones of shear instability similar to the Kelvin-Helmholtz instability at the shear boundary. To evaluate the availability of such a phenomenon, the authors [10] used the approach proposed yet by Heisenberg [14] and Lin [15]. The use of these results in works [10] and [11] helps to conclude that the pattern of plastic flow is similar to the flow of a fluid. The work [11] describes a material with the layer of a nanoscale grain structure, produced by the experimental process. However, these articles only mention the Kelvin-Helmholtz instability, but its analytical results are not used.

There is no single viewpoint on the applicability of the hydrodynamic approach to the formation of nanostructures in a deformed solid because there are no direct observations of the hydrodynamic instability of the material at this level. Electron microscopy enables the scanning of nanoscale structural elements that have already formed in a liquid metal near the shearing boundary of tangential velocity. Note that this method helps to establish that electron diffraction patterns for the nanostructures are of a quasi-ring structure [16-18] specific to liquid and quasi-liquid media that allows us to conclude that the hydrodynamic vision is applicable. The possibility of using the hydrodynamic approach to describe the intense plastic deformation is also implied by the phenomenon of deformation-induced amorphization investigated in [17,18]. The linear analysis of the Kelvin-Helmholtz instability developed in our works [1-3] to various technological problems is presented in [19]. A mechanism and a mathematical model for the formation of shear bands from the position of the Kelvin-Helmholtz instability are proposed in [20]. Within the framework of this model, a neutral curve is obtained, which allows us to determine only the conditions for the onset of instability. Dependence of the decrement on the wave number was not given. Therefore, it is impossible to determine on what scale the development of instability will occur. In addition, the developed approach allows us to propose technologies for the formation of nanostructured states, for example, nanodroplets produced with employing plasma-arc processes [21]. In this case, λ_m corresponds to the size of the nanodrops.

Therefore, our work is aimed at developing a mathematical model that enables parameterization in an approximate model for linear analysis of the Kelvin-Helmholtz instability, and carry out numerical calculations of nonlinear equations in the context of dependences of wavenumbers on parameters derived at maximum decrement values and on the basis of the identified parameters.

This paper is structured as follows: Section 2 sets up the problem. In Section 3 there is a linear analysis of equations of motion with regard to cases relating to viscous and viscoelastic fluids, as well as at the interface with a porous medium. Section 4 gives the results after solving nonlinear equations using numerical analysis and their comparison with the results obtained by linear analysis.

Problem Formulation

Let us consider the stability of a plane motionless layer of a two-layer incompressible fluid (Figure 1). We choose the direction of the x-axis along the interface between the layers, and the y-axis is perpendicular to x and directed toward the second layer. The first layer $(-\infty <x<\infty ,-h<y<\eta (x,t))$ is a viscous, viscoelastic or porous motionless medium. The second layer $(-\infty <x<\infty ,\eta (x,t)<y<H)$ moves with a velocity u₀.

The equation of motion for every layer is expressed as:

$${\rho }_n\frac{d{\overrightarrow{u}}_n}{dt}\text{ = }\frac{\partial {\sigma }_{ik,n}}{\partial x_{k,n}}\text{ (1)}$$

Where σ_ik is the stress components, n = 1,2 is a number of layers. Kinematic and dynamic boundary conditions are explicated as:

$$\begin{array}{l}y\text{ = }H:\text{}{\text{v}}_2\text{ = }0;\\ y\text{ = }-h:{\text{v}}_1\text{ = }0;\\ y\text{ = }0:\frac{\partial \eta }{\partial t}+u_2\frac{\partial \eta }{\partial x}{\text{ = v}}_2,\frac{\partial \eta }{\partial t}+u_1\frac{\partial \eta }{\partial x}{\text{ = v}}_1\text{ (2)}\\ {\sigma }_{1n}-{\sigma }_{2n}-\sigma {\rm K}=0\end{array}$$

Where $K\text{ = }\frac{{\eta }_{xx}}{\sqrt{{\left(1+{({\eta }_x)}^2\right)}^3}}$ is the surface curvature. The condition for shear stresses is not stated. These stresses can be calculated if necessary. For an isotropic viscous fluid, the components of stress tensor have the form:

$${\sigma }_{ik}\text{ = }-p{\delta }_{ik}+2\mu {\dot{\epsilon }}_{ik}\text{ (3)}$$

Where ${\dot{\epsilon }}_{ik}\text{ = }\frac{1}{2}\left(\frac{\partial u_i}{\partial x_k}+\frac{\partial u_k}{\partial x_i}\right)$ is the rate of strain tensor, p is the pressure, µ is the fluid viscosity. For cases when the medium is viscoelastic, the stress tensor has the form:

$${\sigma }_{ik}\text{ = }-p{\delta }_{ik}+2\mu {\dot{\epsilon }}_{ik}+\lambda \theta {\delta }_{ik}+2G{\epsilon }_{ik}\text{ (4)}$$

Where λ and G are constant, θ is the volumetric strain.

According to the results [5] obtained in an approximation of short waves, the viscosity has an effect occurring only at the interface between the layers. Therefore, the Navier-Stokes equations can be replaced with the Euler equations, with the flow in the layers assumed to be free of vortex. Such flows are called viscous-potential. The analysis [5] proves that in this approximation the dispersion equation for waves at the interface of a viscous and an ideal fluid agrees with the dispersion equation in a viscous-potential approximation. Therefore, we intend apply a visco-potential model to viscous, viscoelastic and porous media. According to this model, one can formulate an expression for pressure. To write the equations to describe the viscous-potential model, we introduce the potential of velocity:

$${\overrightarrow{u}}_n\text{ = }\nabla {\Phi }_n\text{ (5)}$$

So that the normal components of tensor (3) take the form:

$${\sigma }_n\text{ = }-p_n+2{\mu }_n\frac{\partial v_n}{\partial y}\text{ (6)}$$

For pressure in viscosity potential flow conditions: $p_n\text{ = }-{\rho }_n\left(\frac{\partial {\Phi }_n}{\partial t}+\frac{1}{2}{\left|\nabla \Phi {}_n\right|}^2+C(t)\right)$. The proof of this formula is similar to the conclusion of the Cauchy-Lagrange and for the stationary case, the Bernoulli integral. Arbitrary function C(t) can be included in Φ.

Linear analysis of the equations of motion

We first consider the two viscous fluids. In this case, it is necessary to put (3) in (1) and the dynamic boundary condition (2), then after their modification the equations of motion will have the form:

$$\begin{array}{l}\frac{\partial u_n}{\partial t}+u_n\frac{\partial u_n}{\partial x}+v_n\frac{\partial u_n}{\partial y}\text{ = }-\frac{1}{{\rho }_n}\frac{\partial p_n}{\partial x}+{\nu }_n\left(\frac{{\partial }^2{u}_n}{\partial x^2}+\frac{{\partial }^2{u}_n}{\partial y^2}\right),\\ \frac{\partial v_n}{\partial t}+u_n\frac{\partial v_n}{\partial x}+v_n\frac{\partial v_n}{\partial y}\text{ = }-\frac{1}{{\rho }_n}\frac{\partial p_n}{\partial y}+{\nu }_n\left(\frac{{\partial }^2{v}_n}{\partial x^2}+\frac{{\partial }^2{v}_n}{\partial y^2}\right),\frac{\partial u_n}{\partial x}+\frac{\partial v_n}{\partial y}\text{ = }0\text{ (7)}\end{array}$$

Where ν_n is the kinematic viscosity of the n-th layer. We linearize (7) and the boundary conditions (2) relating to small perturbations. In the result, we get the following:

$$\frac{\partial U_1}{\partial t}\text{ = }-\frac{1}{{\rho }_1}\frac{\partial P_1}{\partial x}+{\nu }_1\left(\frac{{\partial }^2{U}_1}{\partial x^2}+\frac{{\partial }^2{U}_1}{\partial y^2}\right),\frac{\partial V_1}{\partial t}\text{ = }-\frac{1}{{\rho }_1}\frac{\partial P_1}{\partial y}+{\nu }_1\left(\frac{{\partial }^2{V}_1}{\partial x^2}+\frac{{\partial }^2{V}_1}{\partial y^2}\right),\frac{\partial U_1}{\partial x}+\frac{\partial V_1}{\partial y}\text{ = }0$$

$$\begin{array}{l}\frac{\partial U_2}{\partial t}+u_0\frac{\partial U_2}{\partial x}\text{ = }-\frac{1}{{\rho }_2}\frac{\partial P_2}{\partial x}+{\nu }_2\left(\frac{{\partial }^2{U}_2}{\partial x^2}+\frac{{\partial }^2{U}_2}{\partial y^2}\right),\\ \frac{\partial V_2}{\partial t}+u_0\frac{\partial V_2}{\partial x}\text{ = }-\frac{1}{{\rho }_2}\frac{\partial P_2}{\partial y}+{\nu }_2\left(\frac{{\partial }^2{V}_2}{\partial x^2}+\frac{{\partial }^2{V}_2}{\partial y^2}\right),\frac{\partial U_2}{\partial x}+\frac{\partial V_2}{\partial y}\text{ = }0\text{ (8)}\end{array}$$

The boundary conditions (2) take the following form:

$$\begin{array}{l}y\text{ = }H:V_2\text{ = }0;\\ y\text{ = }-h:V_1\text{ = }0\\ y\text{ = }0:\frac{\partial \eta }{\partial t}+u_0\frac{\partial \eta }{\partial x}\text{ = }{V}_2,\frac{\partial \eta }{\partial t}\text{ = }{V}_1;\\ -P_1+2{\mu }_1\frac{\partial V_1}{\partial y}+P_2-2{\mu }_2\frac{\partial V_2}{\partial y}\text{ = }\sigma {\eta }_{xx}\text{ (9)}\end{array}$$

As mentioned above, we use the model of a viscous-potential fluid. To this end, we linearize equation (5). As a result, we get:

$$U_1\text{ = }\frac{\partial {\Phi }_1}{\partial x},V\text{ = }\frac{\partial {\Phi }_1}{\partial y};U_2\text{ = }\frac{\partial {\Phi }_2}{\partial x},V_2\text{ = }\frac{\partial {\Phi }_2}{\partial y}\text{ (10)}$$

After that, the continuity equations give the following:

$$\frac{{\partial }^2{\Phi }_1}{\partial x^2}+\frac{{\partial }^2{\Phi }_1}{\partial y^2}\text{ = }0,\frac{{\partial }^2{\Phi }_2}{\partial x^2}+\frac{{\partial }^2{\Phi }_2}{\partial y^2}\text{ = }0\text{ (11)}$$

From the momentum conservation equation, it follows

$$P_1\text{ = }-{\rho }_1\left(\frac{\partial {\Phi }_1}{\partial t}\right);P_2\text{ = }-{\rho }_2\left(\frac{\partial {\Phi }_2}{\partial t}+u_0\frac{\partial {\Phi }_2}{\partial x}\right)\text{ (12)}$$

Allowing for (10) and (12), the boundary conditions (9) are written as

$$y\text{ = }H:\frac{\partial {\Phi }_2}{\partial y}\text{ = }0,y\text{ = }-h:\frac{\partial {\Phi }_1}{\partial y}\text{ = }0$$

$$y\text{ = }0:\frac{\partial \eta }{\partial t}+u_0\frac{\partial \eta }{\partial x}\text{ = }\frac{\partial {\Phi }_2}{\partial y},\frac{\partial \eta }{\partial t}\text{ = }\frac{\partial {\Phi }_1}{\partial y}\text{ (13)}$$

$${\rho }_1\left(\frac{\partial {\Phi }_1}{\partial t}\right)+2{\rho }_1{\nu }_1\frac{{\partial }^2{\Phi }_1}{\partial y^2}-{\rho }_2\left(\frac{\partial {\Phi }_2}{\partial t}+u_0\frac{\partial {\Phi }_2}{\partial x}\right)-2{\rho }_2{\nu }_2\frac{{\partial }^2{\Phi }_2}{\partial y^2}\text{ = }\sigma \frac{{\partial }^2\eta }{\partial x^2}$$

Thus, the formulated mathematical problem reduces to solving the Laplace's equations (12) with the boundary conditions (13). The solution to the Laplace's equations satisfying (13) can be written in the form:

$$\begin{array}{l}{\Phi }_1\text{ = }{A}_1\exp (\omega t-ikx)\cosh (k(y+h)),{\Phi }_2\text{ = }{A}_2\exp (\omega t-ikx)\cosh (k(y-H)),\text{ (14)}\\ \eta \text{ = }B\exp (\omega t-ikx),k\text{ = }2\pi /\lambda \end{array}$$

Substitution of (14) into Equations (11) and (13) gives a system of algebraic equations for A₁, A₂, B. If this determinant set to zero, we obtain the dispersion equation:

$\left(R_1+R_2\right){\omega }^2+\omega \left(2i{R}_2{u}_{0}k+(2R_1{\nu }_1+2R_2{\nu }_2)k^2\right)-u_0(2i{R}_2{\nu }_2{k}^3)-R_2{k}^2{u}_0^2+\sigma k^3\text{ (15)}$

Where $R_1\text{ = }{\rho }_1\coth (kh),$ $R_2\text{ = }{\rho }_2\coth (kH).$ We rewrite this equation in the form:

$${\omega }^2+2\omega (a+ib)+c+i{c}_1\text{�\#x00A0;= }0\text{ (16)}$$

$$a\text{ = }\frac{1}{1+A}\left(({\nu }_1+A{\nu }_2)k^2\right),b\text{ = }\frac{Ak{u}_0}{1+A},с\text{ = }\frac{{\omega }_0^2-k^2{u}_0^{2}A}{1+A},{с}_1\text{ = }\frac{2k^3{\nu }_2{u}_{0}A}{1+A}$$

Where $A\text{ = }{R}_2/R_1.$ We note that the dispersion equation (16) follows the dispersion equation derived in [4] with the use of the short-wave approximation. Thus, we assume the flow of material is irrotational in the layers and the effect of viscosity at the interface is significant and these assumptions make possible to obtain the dispersion equation for short waves without additional conditions. From (16), by substituting $\omega \text{ = }\alpha +i\Omega ,$ with subsequent modification, we express the dependence of decrement α on wavenumber k:

$$\alpha \text{ = }\left(\sqrt{\frac{\delta +\sqrt{{\delta }^2+4a^2{b}^2{({\delta }_1-1)}^2}}{2}}-a\right)\text{ (17)}$$

Where $\delta \text{ = }{a}^2-b^2-c,{\delta }_1\text{ = }\frac{c_1}{2ab}.$ We now search for values of wavenumber at which (18) takes on the maximum value. If α > 0, then is instability regions, and contra: if α < 0, then stability. Achieving the maximum value by α = α(k) dependence means that the waves with a wavenumber of ~k_max will progress. These waves will generate particles with dimensions ~λ_max [22]. Equation (17) is suitable for numerical analysis. However, obtaining the dependence of wavelength that includes the maximum α(k) on the problem parameters suggests a cluttered appearance in equations, the physical meaning of which is difficult to comprehend. Therefore, it is necessary to solve the problem of finding an approximate analytical dependence λ_max on the problem parameters. To this end, we transform (17) to the following form:

$$\alpha \text{ = }a\left(\sqrt{\frac{f+\sqrt{f^2+4f_1^2}}{2}}-1\right)\text{ (18)}$$

Where $f\text{ = }\frac{a^2-b^2-c}{a^2},f_1\text{ = }\frac{ab-c_1/2}{a^2}$. Thus, the behavior of decrement with regard to wavenumber depends on the behavior of functions f and f₁. Let us consider two cases. The first case is $\left|f\right|>2\left|f_1\right|$. In this case, (18) will have the form:

$$\alpha \approx a\left(\sqrt{\frac{f+\left|f\right|\left(1+2\epsilon \right)}{2}}-1\right)\text{ = }a\left\{\begin{array}{c}\sqrt{f}-1,f>0;\\ \left|f_1\right|/\sqrt{-f}-1,f<0.\end{array}\right\}\text{ (19)}$$

$$\epsilon \text{ = }{\left(f_1/f\right)}^2$$

The second case is $\left|f\right|<2\left|f_1\right|$:

$$\begin{array}{l}\alpha \approx a\left(\sqrt{\frac{f+2\left|f_1\right|\left(1+{(f/2f_1)}^2/2\right)}{2}}-1\right)\text{ = }a\left(\sqrt{\left|f_1\right|}\sqrt{z^2/2+z+1}-1\right),\text{ (20)}\\ z\text{ = }f/2f_1\end{array}$$

The substitution of f and f₁ in (19) and subsequent transformations result in

$$\begin{array}{l}{\alpha }_1\text{ = }\frac{1}{1+A}\left(k\sqrt{-(\frac{k\sigma }{R_1}(1+A)-A{u}_0^2)}-({\nu }_1+A{\nu }_2)k^2\right)\text{ (21)}\\ {\alpha }_2\text{ = }\frac{k^2}{1+A}\left(\frac{A{u}_0({\nu }_1-{\nu }_2)}{\sqrt{\frac{\sigma k}{R_1}(1+A)-A{u}_0^2}}-({\nu }_1+A{\nu }_2)\right)\end{array}$$

Figure 2 shows the results of numerical calculations according to (17) at ν₁ = 10^-6 m²/s, ν₂ = 10^-8 m²/s, H = 10^-6 m, h = 6 × 10_-6 m, σ = 1.91 N/m, ρ₁ = ρ₂ = 7800 kg/m³, u₀ = 30 m/s. It is seen that the dependence α(k) has two maxima (Curve 1).

The first maximum falls on the wave number km = 4.22 × 106 m^-1, λ_m = 1.487 µm, and the second within km = 1.64 × 10⁷ m^-1 and λ_m = 383 nm. We note that at ν₁ = ν₂ = 0 one maximum of the decrement is observed that is at λ_m = 1.633 µm and km = 3.85 × 10⁶ m^-1. Based on the above, we can conclude (Figure 3) that the viscosity present in one of the layers contributes to generating waves at the submicro and nanoscales. Comparing the diagram (Figure 2, Curves 2 and 3) for the approximations (21) with the dependence (17) shows that these approximations are adequate. To find the maximum values of (21) we reduce them to the dimensionless form:

$$\begin{array}{l}{\overline{\alpha }}_1\text{ = }\left(\overline{k}\sqrt{-(\overline{k}-1)}-{\epsilon }_1{\overline{k}}^2\right)\text{ (22)}\\ {\overline{\alpha }}_2\text{ = }{\overline{k}}^2\left({\epsilon }_2\sqrt{\frac{A}{\overline{k}-1}}-{\epsilon }_1\right)\end{array}$$

$$k_0\text{ = }\frac{A{u}_0^2{R}_1}{(1+A)\sigma },{\alpha }_0\text{ = }\frac{A\sqrt{A}{u}_0^3{R}_1}{{(1+A)}^2\sigma },{\epsilon }_1\text{ = }\frac{u_0{R}_1({\nu }_1+A{\nu }_2)\sqrt{A}}{\sigma (1+A)},{\epsilon }_2\text{ = }\frac{u_0{R}_1({\nu }_1-{\nu }_2)\sqrt{A}}{\sigma (1+A)}$$

Approximate values of the dimensionless critical wavenumbers including the maximum (22) have the form:

$${\overline{k}}_{m1}\text{ = }\frac{2}{3}-\frac{16}{15}{\epsilon }_1,$$

$${\overline{k}}_{m2}\text{ = }1+{\left(\frac{8A^{5/2}{\epsilon }_1-\sqrt{A}{\epsilon }_2(3-2A+3A^2)}{4A(1+3A){\epsilon }_1-2{\epsilon }_2(3A^2+1)}\right)}^2\text{ (23)}$$

Substituting the problem parameters in (23) gives k_m1 = 1.05 × 10⁶ m^-1, k_m1 = 7.55 × 10⁶ m^-1, and this result differ from those obtained using methods of numerical calculations. On the whole, however, one can conclude that the order of critical wavenumber and wavelength can be defined with using (23), which has the maximum instability decrement, and to evaluate further the dimensions of structural elements.

Let us compare the results of calculation using (17) with those of obtained using numerical calculations [1] for the interface between layers of an ideal and a viscous fluid. The dispersion equation obtained in this work has the form:

$${(z^2+1)}^2+\mu {\left((z^2-1)+i{p}_0\right)}^2+{\omega }_{01}-4z\text{ = }0\text{ (24)}$$

where $p_0\text{ = }\frac{u_0}{{\nu }_{1}k},{\omega }_{01}\text{ = }\frac{\sigma }{{\nu }_1^{2}k{R}_1},z\text{ = }\sqrt{1+\frac{\omega }{{\nu }_1{k}^2}}$

Figure 4а represents the dependencies expressed using (17) and (24) when ν₁ = 10^-6 m²/s, ν₂ = 0, H = 10^-6 m, h = 10^-2 m, σ = 1.91 N/m, ρ₁ = ρ₂ = 7800 kg/m³, u₀ = 30 m/s. In both cases (17) and (24), two maxima are observed in the dependency of the decrement on the wavelength. But the values of the critical wavelengths are slightly different. If according to (17) λ_m1 = 318 nm and λ_m2 = 2.76 µm, then according to (24) λ_m1 = 340 nm and λ_m2 = 2.34 µm. Increasing a velocity of up to 50 m/s results in λ_m1 = 114 nm and λ_m2 = 612 nm according to (17), and λ_m1 = 120 nm according to (24), without the second maximum observed (Figure 4b).

On the basis of the above stated, it can be concluded that (17) adequately describes the Kelvin-Helmholtz instability at interfaces between layers.

We now proceed to the case relating to viscoelastic media. For such a medium, the equation for conservation of momentum is obtained by putting (4) into the equations of motion (1). As shown in [23], the viscoelastic medium differs from the viscous one because viscosity depends on relaxation time according to the Maxwell's equation:

$${\nu }_1^{*}\text{ = }\frac{{\nu }_1}{1+\omega t_{1*}},{\nu }_2^{*}\text{ = }\frac{{\nu }_2}{1+\omega t_{2*}}\text{ (25)}$$

Here, $t_{n*}$ is the relaxation time. The dispersion equation can be derived in the same way, and, in the short wave approximation, the dispersion equation itself coincides with equation (16) in form, but with different coefficients:

$$\begin{array}{l}{\omega }^2+2\omega (a_{*}+i{b}_{*})+c_{*}+i{c}_{1*}\text{ = }0,\\ a_{*}\text{ = }\frac{a}{1-\omega t_{*}},b_{*}=\frac{b}{1-2\omega t_{*}},c_{*}\text{ = }\frac{c}{1-2\omega t_{*}},c_{1*}=\frac{c}{1-2\omega t_{*}}\text{ (26)}\end{array}$$

For simplicity, we assume that the time of relaxation in Layers 1 and 2 are the same. Then the dependence of the decrement α on the wavenumber is written as:

$${\alpha }_1\text{ = }\frac{1}{1-\omega t_{*}}\left(\sqrt{\frac{\delta -\omega t_{*}+\sqrt{{(\delta -\omega t_{*})}^2+4a^2{b}^2{\left({\delta }_1-\omega t_{*}-1\right)}^2}}{2}}-a\right)\text{ (27)}$$

Thus, it follows from the linear analysis of the equations of motion that the two-mode Kelvin-Helmholtz instability is observed in the case of a viscous, viscoelastic and porous medium. This means that the waves with a wavelength of λ_m1 and λ_m2 will progress. As is known, linear analysis is valid only for small perturbations. In the case where the boundary is displaced from its equilibrium position, nonlinear analysis needs to be done. The next section is devoted to nonlinear numerical analysis using the finite element method.

Nonlinear numerical analysis

We use the finite element method as a base for numerical analysis to solve nonlinear equations. In recent times, the Level-Set Method (LSM) has been used to solve the hydrodynamic equations of multilayer immiscible fluids [24]. This method is based on the idea that the motion of interfaces is described by means of the so-called level-set functions that determine a distance from a given point to an interface. The function is assumed to take positive values for points inside the fluid and negative values for points in the empty space. Due to a smooth change in this function when crossing the interface, the effect of diffusion is insignificant at the interface, but for large times the values of the function already cease to have their meaning as distances to the interface, and from time to time the distance to the interface must be redefined in the calculation [25-27]. In this article, calculations are made in the COMSOL Multiphysics package under the following initial conditions:

$$u_1(0)\text{ = }0,v_1(0)\text{ = }0,u_2(0)\text{ = }{u}_0,v_2(0)\text{ = }{A}_m\sin \left(\frac{2\pi x}{\lambda }\right)\text{ (28)}$$

Where A_m is the amplitude of perturbation on the velocity. Figure 5 and Table 1 shows the boundary conditions. Table 2 shows the materials characteristics and input parameters of problem.

Figure 6 and Figure 7 demonstrate wave evolution at different times in different modes. It is seen that the perturbations at the interface gradually increase first in Regime 1 (Figure 6a and Figure 6b). Then, at t = 10 ns (Figure 6c) the wave front breaks down. The process of forming a vortex starts at a time of 22 ns (Figure 6d). At t > 22 ns, the vortex structure is fragmented and droplets are formed. The dimensions of these droplets range from ~90 nm to ~242 nm. In addition to the vortex breaking down into droplets, the process of merging small droplets into larger drops is observed in Regime 1. Thus, when we allow for the nonlinearity while investigating the Kelvin-Helmholtz instability, we can observe the two processes of forming nanosized droplets due to the decay of the vortex and their merging into submicrosized particles. According to linear analysis, in theses modes, the single-mode dependences of decrement should be observed, with their maxima falling on a wavelength of 194 nm (Regime 1) In mode 2, the situation is radically different at t = 10 ns (Figure 7a), the perturbations first increase. Then, at a time of 30 ns (Figure 7b), a large drop begins to form, which at t = 50 ns (Figure 7c) breaks off. The size of these drops is ~1.53 µm. At t > 70 ns small droplets begin to form (Figure 7d) whose size is ~150 nm. According to the linear analysis (see the previous section), in this case there will be a two-mode instability. Thus, a nonlinear analysis of the Kelvin-Helmholtz instability fully confirms the results of the linear analysis.

Conclusion

1. A simplified dispersion equation is obtained for short-wave perturbations at the boundary of two viscous liquids. The dependence of the decrement on the wavelength for the total and the reduced dispersion equations for short waves coincide.

2. A simplified equation makes it possible to obtain an analytical dependence of the decrement on the wavenumber. An analysis of this dependence showed that there exist one- and two-mode Kelvin-Helmholtz instability.

3. Parameters are determined under which one of the maxima of this dependence occurs in the submicro and nanorange. For viscoelastic, the same two-mode dependences are obtained.

4. The modeling of the Kelvin-Helmholtz instability in the nonlinear case for the parameters corresponding to single- and double-mode instabilities has shown a fundamental difference: For a single-mode regime, the droplet sizes lie in a narrow range of parameters, and for the two-mode regime the maximum and minimum dimensions differ by a factor of 10.

5. The obtained results will allow not only to explain, but also to predict the conditions for the formation of nanomaterials and to select technological parameters in the formation of nanostructured states.

This research is supported by the Russain Science Foundation (№ 15-12-00010).

---