Extended two dimensional equation for the description of nonlinear waves in gas-liquid mixture

TL;DR

This paper derives a new two-dimensional nonlinear PDE for gas-liquid mixture waves, investigates its integrability, constructs exact solutions, and performs numerical analysis to understand nonlinear wave behavior.

Contribution

It introduces a novel extended 2D nonlinear equation for gas-liquid mixtures, including high-order terms, and analyzes its integrability and solutions.

Findings

The derived equation is integrable under certain parameter conditions.

Exact solutions of the equation are constructed.

Numerical simulations reveal wave dynamics described by the equation.

Abstract

We consider a system of equations for the description of nonlinear waves in a liquid with gas bubbles. Taking into account high order terms with respect to a small parameter, we derive a new nonlinear partial differential equation for the description of density perturbations of mixture in the two-dimensional case. We investigate integrability of this equation using the Painlev\'e approach. We show that travelling wave reduction of the equation is integrable under some conditions on parameters. Some exact solutions of the equation derived are constructed. We also perform numerical investigation of the nonlinear waves described by the derived equation.

Equations

$$\begin{array}{c}\frac{\partial \tilde{\rho }}{\partial \tau }+\nabla \tilde{\text{u}}+\nabla (\tilde{\rho }\tilde{\text{u}})=0,\\ (1+\tilde{\rho })\left(\frac{\partial \tilde{\text{u}}}{\partial \tau }+\tilde{\text{u}}\nabla \tilde{\text{u}}\right)+\frac{1}{\alpha }\nabla \tilde{p}=0,\\ p=\alpha \tilde{\rho }+{\alpha }_1{\tilde{\rho }}^2+{\alpha }_2{\tilde{\rho }}^3+\beta {\tilde{\rho }}_{\tau \tau }-\left({\beta }_1+{\beta }_2\right)\tilde{\rho }{\tilde{\rho }}_{\tau \tau }-\left({\beta }_1+\frac{3}{2}{\beta }_2\right){\tilde{\rho }}_{\tau }^2+\varkappa {\tilde{\rho }}_{\tau }+{\varkappa }_1\tilde{\rho }{\tilde{\rho }}_{\tau }.\end{array}$$

$$x=ϵ(\xi -\tau ),y={ϵ}^{\frac{3}{2}}\delta \eta ,t={ϵ}^2\tau .$$

$$\begin{array}{c}{\tilde{u}}^{(1)}=ϵu_1^{(1)}+{ϵ}^2{u}_2^{(1)}+\cdots ,{\tilde{u}}^{(2)}=ϵu_1^{(2)}+{ϵ}^2{u}_2^{(2)}+\cdots ,\\ \tilde{\rho }=ϵ{\rho }_1+{ϵ}^2{\rho }_2+\cdots ,\tilde{p}=ϵp_1+{ϵ}^2{p}_2+\cdots .\end{array}$$

$$u_1^{(1)}={\rho }_1,p_1=\alpha {\rho }_1.$$

$$\begin{array}{c}{\rho }_{1t}-{\rho }_{2x}+u_{2x}^{(1)}+({\rho }_1{u}_1^{(1)}{)}_x+ϵ\delta u_{1y}^{(2)}=0,\\ u_{1t}^{(1)}-u_{2x}^{(1)}+u_1^{(1)}{u}_{1x}^{(1)}+{\rho }_{2x}+\frac{2{\alpha }_1}{\alpha }{\rho }_1{\rho }_{1x}-\frac{\varkappa }{\alpha }{\rho }_{1xx}-{\rho }_1{u}_{1x}^{(1)}+\\ +ϵ\left({\rho }_1{u}_{1t}^{(1)}+{\rho }_1{u}_1^{(1)}{u}_{1x}^{(1)}+\frac{{\alpha }_2}{\alpha }({\rho }_1^3{)}_x+\frac{\beta }{\alpha }{\rho }_{1xxx}+\frac{\varkappa }{\alpha }{\rho }_{1tx}-\frac{{\varkappa }_1}{\alpha }({\rho }_1{\rho }_{1x}{)}_x\right)=0,\end{array}$$

$$u_{1x}^{(2)}=\delta {\rho }_{1y}+ϵ\left(u_{1t}^{(2)}+u_1^{(2)}{u}_{1x}^{(1)}\frac{\delta {\alpha }_1}{\alpha }{\rho }_{1y}^2-\frac{\delta \varkappa }{\alpha }{\rho }_{1xy}-{\rho }_1{u}_{1x}^{(2)}\right).$$

$$\begin{array}{c}{\left({\rho }_{1t}+\left(1+\frac{{\alpha }_1}{\alpha }\right){\rho }_1{\rho }_{1x}-\frac{\varkappa }{2\alpha }{\rho }_{1xx}\right)}_x+\\ +ϵ\frac{1}{2}{\left(\frac{3{\alpha }_2-{\alpha }_1}{\alpha }{\rho }_1^2{\rho }_{1x}+\left(\frac{\beta }{\alpha }+\frac{{\varkappa }^2}{2{\alpha }^2}\right){\rho }_{1xxx}+\frac{\varkappa }{2\alpha }{\rho }_1{\rho }_{1xx}-\frac{\varkappa \left(2\alpha +{\alpha }_1\right)}{{\alpha }^2}({\rho }_1{\rho }_{1x}{)}_x\right)}_x+\\ +ϵ\frac{{\delta }^2}{2}{\rho }_{1yy}=0.\end{array}$$

$${\rho }_1=\rho +ϵ\left({\lambda }_1{\rho }^2+{\lambda }_2{\rho }_x{\partial }_x^{-1}\rho \right),$$

$${\left({\rho }_t+a_1\rho {\rho }_x+a_2{\rho }_{xx}+a_3(\rho {\rho }_x{)}_x+a_4{\rho }^2{\rho }_x+a_5{\rho }_{xxx}\right)}_x+b{\rho }_{yy}=0,$$

$$\begin{array}{c}{a}_1=\left(1+\frac{{\alpha }_1}{\alpha }\right),a_2=-\frac{\varkappa }{2\alpha },a_3=ϵ\frac{\varkappa (1-2\alpha -{\alpha }_1)-4\alpha ({\lambda }_1+2{\lambda }_2)}{4{\alpha }^2},\\ a_4=ϵ\left(\frac{3{\alpha }_2-{\alpha }_1}{2\alpha }+\left(1+\frac{{\alpha }_1}{\alpha }\right)\left({\lambda }_1+{\lambda }_2\right)\right),a_5=ϵ\frac{2\beta \alpha +{\varkappa }^2}{4{\alpha }^2},b=ϵ\frac{{\delta }^2}{2}.\end{array}$$

$${\lambda }_1=-{\lambda }_2+\frac{\varkappa \left(1-6\alpha -3{\alpha }_1\right)}{\left(8\varkappa +4\right)\alpha }.$$

$$\begin{array}{c}{x}^{\prime }=Ax+Bt,t^{\prime }=Ct,{\rho }^{\prime }=D\rho +E,y^{\prime }=Fy,\end{array}$$

$$\begin{array}{c}B=\frac{a_1^{2}A}{4a_4},C=a_5{A}^3,D=-\frac{a_{5}A}{a_3},E=-\frac{a_1}{2a_4},F=\sqrt{\frac{a_5}{b}}{A}^2,\mu =\frac{a_4{a}_5}{a_3^2},\end{array}$$

$${\left({\rho }_t+\mu {\rho }^2{\rho }_x+{\rho }_{xxx}-(\rho {\rho }_x{)}_x\right)}_x+{\rho }_{yy}=0.$$

$$\rho (x,y,t)={\Phi }^p\sum _{j=0}^{\infty }{u}_j{\Phi }^j,$$

$$p=-1,u_0=\frac{\left(-3\pm \sqrt{9-24\mu }\right){\Phi }_x}{2\mu }.$$

$$\rho (x,y,t)=\frac{u_0(x,y,t)}{\Phi (x,y,t)}+u_{0j}\Phi (x,y,t{)}^{j-1},$$

$$\begin{array}{c}{j}_1^{(1,2)}=-1,j_2^{(1,2)}=3,j_3^{(1,2)}=4,j_4^{(1,2)}=\frac{\pm \sqrt{9-24\mu }+8\mu -3}{2\mu }.\end{array}$$

$$\begin{array}{c}{j}_1^{(1)}=-1,j_2^{(1)}=0,j_3^{(1)}=3,j_4^{(1)}=4,\\ j_1^{(2)}=-1,j_2^{(2)}=0,j_3^{(2)}=3,j_4^{(2)}=4.\end{array}$$

$$\begin{array}{c}{j}_1^{(1)}=-1,j_2^{(1)}=3,j_3^{(1)}=3,j_4^{(1)}=4;\\ j_1^{(2)}=-1,j_2^{(2)}=3,j_3^{(2)}=4,j_4^{(2)}=6.\end{array}$$

$$\begin{array}{c}{j}_1^{(1)}=-1,j_2^{(1)}=1,j_3^{(1)}=3,j_4^{(1)}=4,\\ j_1^{(2)}=-1,j_2^{(2)}=-2,j_3^{(2)}=3,j_4^{(2)}=4.\end{array}$$

$$\begin{array}{c}{u}_0=-3{\Phi }_x,u_2=\frac{u_1^2}{3{\Phi }_x}-\frac{u_{1x}}{{\Phi }_x}-\frac{{\Phi }_{xx}{u}_1}{{\Phi }_x^2}+\frac{{\Phi }_y^2}{{\Phi }_x^3}+\frac{{\Phi }_{xxx}}{{\Phi }_x^2}+\frac{{\Phi }_t}{{\Phi }_x^2},\end{array}$$

$$3\frac{({\Phi }_y{)}^2{\Phi }_{xx}}{{\Phi }_x}+3{\Phi }_{yy}{\Phi }_x-6{\Phi }_y{\Phi }_{xy}=0.$$

$$\rho (x,y,t)=\frac{A_0(x,y,t)}{\Phi (x,y,t)}+A_1(x,y,t).$$

$$\omega =\frac{1}{2k_1}\left(k_x^4+2k_y^2\right),$$

$$\begin{array}{c}{A}_0=-4k_x{e}^{k_{x}x+k_{y}y-\omega t+{\varphi }_0},A_1=2k_x.\end{array}$$

$$\rho =-2k_x\frac{\left(e^{k_{x}x+k_{y}y-\omega t+{\varphi }_0}-1\right)}{1+e^{k_{x}x+k_{y}y-\omega t+{\varphi }_0}}.$$

$$\begin{array}{c}\omega =-\frac{1}{4k_x}\left(k_x^4-4k_y^2\right),A_0=-k_x{e}^{k_{x}x+k_{y}y-\omega t+{\varphi }_0},A_1=\frac{k_x}{2},\end{array}$$

$$\rho =-\frac{k_x}{2}\frac{\left(e^{k_{x}x+k_{y}y-\omega t+{\varphi }_0}-1\right)}{1+e^{k_{x}x+k_{y}y-\omega t+{\varphi }_0}}.$$

$$\begin{array}{c}\omega =\frac{k_x^4+k_y^2}{k_x},A_0=-6k_x{e}^{k_{x}x+k_{y}y-\omega t+{\varphi }_0},A_1=3k_x,\end{array}$$

Full text

Keywords: Nonlinear equation; Nonlinear wave; Liquid with gas bubbles; Reductive perturbation method; Painlevé test; Exact solutions

1 Introduction

A liquid with gas bubbles is a complex dissipative and dispersive nonlinear media. Nonlinear character of waves in such medium brings essential difficulties for investigation, although there are some interesting properties of wave processes in gas-liquid mixture and mathematical models for the description of such systems occurs widely in different sciences: chemistry, biology, physic and etc. (see [1, 2, 3]).

For the first time nonlinear evolution equations like Burgers, Korteweg–de Vries and Burgers–Korteweg–de Vries were obtained for the description of long weakly nonlinear waves in a gas–liquid mixture in works [4, 5, 6] for the one-dimensional case. The three–dimensional case was considered in work [7], but only first–order terms in an asymptotic series have been taken into account. On the other hand, considering high–order corrections in asymptotic expansions, we are able to obtain more complicated nonlinear equations. It allows us to describe wave processes more accurate then in [7]. Besides, we can discover some new physical effects. In work [8] models for non–linear waves in a gas–liquid mixture were generalized, taking into account higher order terms with respect to small parameters. Models of work [8] take into consideration an interphase heat transfer, surface tension and weak liquid compressibility, although only one–dimensional case is considered. Thus, it is interesting to study long weakly nonlinear waves in a liquid with a gas bubbles in two–dimensional case, taking into consideration both high order terms in the asymptotic expansions and physical properties mentioned above.

Here we derive a new nonlinear partial differential equation for the description of long weakly nonlinear two–dimensional waves in a bubbly gas–liquid mixture. We consider waves propagating in a certain direction. We assume that perturbations in perpendicular directions are less essential but we take them into account. We take into account high order terms in the asymptotic expansions, interphase heat transfer, surface tension and weak liquid compressibility. We also investigate equation derived analytically and numerically. To the best of our knowledge, this equation has not been obtained and investigated before.

In order to investigate integrability of the nonlinear equation we apply the Painlevé approach. It is shown that the equation does not have the Painlevé property in general case. However solitary wave solutions are constructed by means of the truncated expansion method. Using travelling wave variables it is shown that the equation passes the Painlevé test under some conditions on parameters. With the Hopf–Cole transformations the equation is linearised and its solutions are obtained in different forms. Nonlinear waves described by the equation are also investigated numerically using the spectral method. It is shown, that this method have good accuracy and stability.

The rest of this work is organized as follows. In section 2 we derive a nonlinear partial differential equation for the description of waves in gas–liquid mixture, taking into consideration second order terms with respect to the small parameters. In section 3 we apply the Painlev? approach to investigate integrability of the equation. In section 4 the new nonlinear equation is investigated using travelling wave variables. It is shown that the equation is integrable under some conditions on parameters. In section 5 we present the results of the numerical simulation of waves, described by the equation. In the last section we briefly discuss our results.

2 Extended equation for the description of waves in a liquid with gas bubbles in two–dimensional case.

In this section we obtain a two–dimensional nonlinear equation for the description of waves in a liquid with gas bubbles. We use the system of equations for the description of waves in bubbly liquid, presented in [7]. We suppose that the gas–liquid mixture is a homogeneous medium with an average pressure and temperature. We assume that the liquid is incompressible and gas bubbles are spherical. We do not consider destruction, formation, interaction and coalescence of bubbles. We suppose that total amount of gas in a bubble and the amount of gas bubbles in unit of mass of liquid are constants. Gas in bubble is an ideal and a pressure in bubble is described by the politropic law. Liquid viscosity is considered only on the interphase boundary. Taking into account assumptions mentioned above, the following system of equations for the description of waves in liquid with gas bubbles is used (see [7] )

[TABLE]

Here $p,\,\tilde{\rho},\,\tilde{\textbf{u}}$ are the non–dimensional pressure, density and velocity of the mixture correspondingly, $\xi,\,\eta$ are Cartesian coordinates and $\tau$ is the time; $\alpha,\,\alpha_{1},\,\alpha_{2},\,\beta,\,\beta_{1},\,\beta_{2},\,\varkappa,\,\varkappa_{1}$ are non–dimensional parameters [7].

For the derivation of an equation for the description of nonlinear waves, we use the reductive perturbation method (see e.g. [9, 10, 11, 12, 13]). Let us introduce ’slow’ variables

[TABLE]

We suppose that perturbations in $x$ direction are more essential then in $y$. We chose power of $\epsilon$ in ’slow’ variables in order to obtain equations for the case of dissipation main influence. We search for a solution of system (1) in the form of asymptotic series:

[TABLE]

Substituting (2) and (3) into (1) and collecting coefficients at $\epsilon^{0}$ we obtain

[TABLE]

Collecting coefficients at the same powers of $\epsilon$ and using (4) we have the following equations:

[TABLE]

[TABLE]

Differentiating (5) and (6) with respect to $x$ and $y$ correspondingly and using obtained relations to avoid velocity, we get:

[TABLE]

Now, using the near-identity transformations [14, 15]

[TABLE]

we obtain equation

[TABLE]

where

[TABLE]

and

[TABLE]

Let us use shift and scaling transformations in the form:

[TABLE]

where $A$ is an arbitrary parameter, $B,\;C,\;D,\;E,\;F$ are described by the following relations:

[TABLE]

constant $\lambda_{2}$ is found from the equation $a_{1}a_{3}=2a_{2}a_{4}$. Taking into account transformations (12) we can rewrite equation (7) in the form

[TABLE]

Below we study equation (14).

3 Painlevé test to equation (14)

To investigate integrability of (14) let us apply the Weiss–Tabor–Carnevale (WTC) test [16, 17]. We look for a solution of (14) in the form of series

[TABLE]

where $\Phi=\Phi(x,y,t)$ is a new function and $u_{j}=u_{j}(x,y,t)$ are coefficients in expansion (15). The leading terms of equation (14) are $\mu\left(\rho^{2}\rho_{x}\right)_{x},\,\rho_{xxxx}$ and $-\left(\rho\rho_{x}\right)_{xx}$. Substituting $\rho=u_{0}\Phi^{p}$ into leading terms and equating coefficients at the lowest order of $\Phi(x,y,t)$ we obtain

[TABLE]

Then, substituting into leading terms the expression

[TABLE]

and equating the coefficient at $u_{j}$ to zero, the following Fuchs indexes are found for different branches of solution

[TABLE]

Here upper index is a number of a branch of the solution and lower index is a number of an index for definite branch. Let us recall that an equation pass the Painlevé test only if Fuchs indexes are integer. This holds only for the limited number of values of $\mu$. In case of $\mu=\frac{3}{8}$ we have the flowing Fuchs indexes on each branch of the solution:

[TABLE]

We see that equation (14) does not pass Painlevé test in case of $\mu=\frac{3}{8}$, because if the Fuchs index equal to zero, $u_{0}$ must be an arbitrary function, but it is determined, thus we can not take three arbitrary constants in expression (15).

In case of $\mu=-3$ we have two sets of Fuchs indexes:

[TABLE]

Thus equation (14) does not pass Painlevé test when $\mu=-3$ because $j_{4}^{(1)}=j_{2}^{(1)}=3$ and expression (15) has to have logarithmic terms.

Let us consider case $\mu=\frac{1}{3}$. We have Fuchs indexes

[TABLE]

According to the idea of Painlevé test we recall that (14) have the Painlevé property if functions $\Phi,\,u_{1},\,u_{3},\,u_{4}$ are arbitrary in (15) for one branch of solution and for the other branch we need to use special modification of Painlevé test. Substituting (15) into (14) and consistently equating terms at various powers of $\Phi$ to zero, we find

[TABLE]

and $u_{1}$ is an arbitrary function. Coefficient at $\Phi^{-2}$ has to be zero as $u_{3}$ has to be an arbitrary function. However we obtain that $u_{3}$ can be arbitrary only if the following equation is satisfied:

[TABLE]

It is obvious from (23) that $u_{3}$ can not be taken as an arbitrary function. Thus equation (14) does not have Painlevé property and seems not to be integrable. However one is able to find a solution for equation (14) applying truncated expansion method [19, 18]. We search for a solution of equation (14) in the form

[TABLE]

Substitute expression (24) into equation (14) and equating coefficients at different powers of $\Phi$ to zero we find coefficients $A_{0},A_{1}$. As general case leads to cumbersome formulas we fix parameter $\mu=\frac{3}{8}$, $\mu=-3$ or $\mu=\frac{1}{3}$. These cases correspond to integer values of Fuchs indexes. We also search for the solution in the form of solitary wave so we set $\Phi=e^{k_{x}x+k_{y}y-\omega t+\phi_{0}}$. After all operations in case of $\mu=\frac{3}{8}$ we get the dispersion relation

[TABLE]

and functions

[TABLE]

Thus, the solution of (14) in the form of solitary wave have the form

[TABLE]

On figure 1 the exact solution for the equation (14) in the form of solitary wave (27) is illustrated.

Analogously to case of $\mu=\frac{3}{8}$, setting $\mu=-3$ dispersion relation and parameters take the form

[TABLE]

Substituting (28) into (24) analytical solution of the equation (14) is found in the form

[TABLE]

And in case of $\mu=\frac{1}{3}$ one can obtain relations

[TABLE]

and the solution of equation (14)

[TABLE]

It is obvious that solutions (29) and (31) has the same type as (27).

4 Travelling wave solutions of equation (14).

Let us investigate equation (14) using the travelling wave variables. Assuming that $\rho=\rho(z),\;z=x+y-C_{0}t$ and integrating equation (14) with respect to $z$, we get

[TABLE]

where $C_{1}$ is a constant of integrating. Substituting $\rho(z)=a_{0}z^{p}+a_{j}z^{p+j}$ into the leading terms of equation (32), we find constants $p=1$, $a_{0}=-\frac{3\pm\sqrt{9-24\,\mu}}{2\mu}$ and three Fuchs indexes which are the same as (18) except index $j=4$. Therefore equation (32) may pass the Painlevé test only if $\mu=\frac{1}{3}$. In this case the Fuchs indexes are

[TABLE]

Here upper index stands for the number of branch of solution of (32) and lower index is the number of index on the certain branch. On the first branch, according to Painlevé approach solution of equation (32) is being searched in the form of expansion

[TABLE]

Substituting (34) into (32) and equating coefficients at various powers of $z$ to zero, we obtain that $a_{1},\;a_{3}$ are arbitrary constants and

[TABLE]

The other constants can be defined too. On the second branch of solution there is the Fuchs index below zero. One can confirm that equation (14) pass the Painlevé test in travelling wave variable using Conte–Fordy–Pickering algorithm (see, [20, 21]).

Let us apply Hopf–Cole transformations $\rho(z)=-3\frac{\psi^{\prime}(z)}{\psi(z)}$ for linearization of equation (32) in case of $\mu=\frac{1}{3}$. Using this transformations we obtain linear equation for $\psi(z)$

[TABLE]

Taking into account the general solution of equation (36) we obtain general solution of equation (32)

[TABLE]

Under some conditions on parameters, for example $c_{0}=1,c_{1}=1,c_{2}=1,C_{0}=1.25$ we are able to obtain solution in the form of solitary wave. The result is demonstrated on figure 2.

When $c-1$ is below zero, we are able to obtain the periodical solution. In case of $t=0,\;C_{0}=-3,\;c_{0}=1,\;c_{1}=1,\;c_{2}=1$ solution have form

[TABLE]

and is given on figure 3. It is important to obtain solution of equation (14) in the form of periodical wave (38) because it can be used for the testing of program for numerical solution of equation (14) with the spectral method and periodical boundary value conditions.

Solution (37) is obtained under strong conditions on parameters. Let us apply method of the logistic function [22, 23] to get solution of equation (32) in case of arbitrary value of parameter $\mu$. We search for the solution of equation (32) in the following form:

[TABLE]

where $B_{0},\;B_{1}$ are constants and $\Theta\left(z\right)$ is the so–called logistic function. Following [22], we expend (39) into Laurent series, substitute obtained expansion into equation (32) and equate coefficients at different powers of $z$ to zero. We find that there are two solutions of equation (32) in the form (39) whith the following constants:

[TABLE]

It is worth noting that since $B_{1}^{(1)}+B_{1}^{(2)}\neq 0$, there are no elliptic solutions of equation (32). We can also assume that $kz\rightarrow kz+i\pi$ and obtain solution in the form (39) without poles on the real line.

5 Numerical periodical solutions of the equation (14).

In this section we study numerically nonlinear waves governed by equation (14). To this end we consider boundary value problem with periodical boundary conditions. Let us rewrite eq. (14) in the form

[TABLE]

where $L(\rho)=-\rho_{xxx}-\partial^{-1}_{x}\rho_{yy}$ and $N(\rho)=\left(\rho\rho_{x}\right)_{x}-\mu\left(\rho^{2}\rho_{x}\right)$ are linear and nonlinear operator of the equation (14) correspondingly. Antiderivative $\partial^{-1}_{x}$ is an integral with respect to $x$ . To solve boundary value problem which includes equation (41), start conditions and periodical boundary conditions, we use the integrating factor with the fourth–order Runge–Kutta approximation method (IFRK4), presented in works [24, 25]. Using the Fourier transformation, we discretize the spatial part of equation (14) and get the system of ordinary differential equations

[TABLE]

where $\widehat{\rho},\;\widehat{L}=ik_{x}^{3}-\frac{ik_{y}^{2}}{k_{x}}$ and $\widehat{N[\rho]}=-\frac{ik_{x}\mu}{3}\widehat{\rho^{3}}-\frac{1}{2}k_{x}^{2}\widehat{\rho^{2}}$ are the Fourier forms of $\rho,\;L,\;N$ correspondingly, $k_{x},\;k_{y}$ are Fourier multipliers. The basic idea of integrating factor (IF) method is to use transformations that allows us to solve the linear part of our problem exactly and then to solve the nonlinear part of our problem numerically. These transformations are the following:

[TABLE]

where $e^{-\widehat{L}t}$ is the IF. As a result, we obtain ordinary differential equation in the following form

[TABLE]

Equation (44) with initial conditions can be solved by the fourth–ordered Runge–Kutta approximation method. It is obvious that equation (44) has to be regularized for $k_{x}=0$ in order to give numerical sense to $\frac{1}{k}_{x}$. Following [26], we add to $k_{x}$ in the denominator a small imaginary part of appropriate sign $i\lambda_{0}$. For $\lambda_{0}$ we use the smallest floating point number that MATLAB can represent $2.2\times 10^{-16}$. In this case the linear operator has form $\widehat{L}=ik_{x}^{3}-\frac{ik_{y}^{2}}{k_{x}+i\lambda}$.

To check our numerical strategy we used the exact solution of the equation (14) in the periodical form (38). On each tome layer we define error for numerical simulation as

[TABLE]

where $i$ and $j$ are discrate coordinates of points in $x$ and $y$ axis correspondingly, $\rho^{i,j}_{num}$ and $\rho^{i,j}_{exact}$ are values of a numerical and exact solutions in the grid point with the coordinate $(i,j)$ correspondingly. Evolution of the error with the time $t$ is given on the figure 4. Taking into account results of numerical simulation of the exact solution we can regard the algorithm effective.

Let us use numerical simulation for the investigation of stability of equation (14) to perturbations of the parameter $\mu$. In case fo $\mu=\frac{1}{3}$, traveling wave reduction of equation (14) is integrable and we use its exact solution (38) as initial conditions for numerical simulation of equation (14) in case of $\mu=\frac{1}{3}+\delta$. The result is presented on figure 5 for $\mu=0.43$. We see that wave hold its shape but acquire different speed.

We also use solitary wave $\frac{5}{\cosh^{2}(x)+\cosh^{2}(y)}$ as an initial condition for numerical calculations. Result of calculations is presented on figure 6. Small perturbations in $y$ direction appears and amplitude of the solitary wave decreases during time evolution.

6 Conclusion

We have obtained new nonlinear equation (14) for the description of waves in a liquid with gas bubbles in two–dimensional case. Using the Painlevé approach, we have shown that this equation is not integrable. We have constructed some analytic solutions of equation (14) with the help of the truncated expansion method. We have shown that traveling wave reduction of eq. (14) can pass the Painlevé test under some conditions on parameters. We have shown that in this case equation (14) can be linearized with the Hopf–Cole transformation. It has been shown that the equation for the description of waves in a gas–liquid mixture admits kink–type and periodical solutions. We have numerically investigated the evolution and stability of nonlinear waves discribed by equation (14).

7 Acknowledgments

This research was partially supported by grant for Scientific Schools 2296.2014.1., by RFBR grant 14–01–00498 and by grant for the state support of young Russian scientists 3694.2014.1.