Direct Numerical Simulation of Acoustic Turbulence: Zakharov–Sagdeev Spectrum

We present the results of direct numerical simulation of three-dimensional acoustic turbulence in medium with weak positive dispersion. It is shown that at the beginning of the long-wavelength region in the turbulence energy distribution in the k-space, there are formed jets in the form of narrow cones. At higher wavenumbers, the cones broaden, and the distribution accordingly tends to isotropic. In this region of wavenumbers, the angle-averaged turbulence spectrum acquires a power-law character, $E(k) \propto {{k}^{{ - \alpha }}}$, with the exponent close to 3/2, which corresponds to the Zakharov–Sagdeev weak acoustic turbulence spectrum.

INTRODUCTION

As well known (see, e.g., [1–3]), the weak turbulence theory is based on the assumption that the nonlinear interaction of waves is weak in comparison with the linear wave dispersion, which is determined by the dependence of the second derivative of the eigen frequency ω of small amplitude waves relative to the wave vector k. As well known, for both deep water gravitational and capillary waves, the second derivative of $\omega ({\mathbf{k}})$ is nonzero in the entire range of wavenumbers. Thereby, initially Gaussian-distributed linear waves with different k almost remain this property when weak nonlinearity is taken into account. Each wave moving with its own frequency and wave vector experiences the influence of other waves at distances L greater its wavelength, $\propto {\kern 1pt} {{k}^{{ - 1}}}$. This makes it possible to describe the system of waves statistically in the random phase approximation using kinetic equations for the number of waves (quasiparticles) ${{n}_{k}}$ [1]. In the leading approximation in nonlinearity, the kinetic equations describe either decay processes ($1 \to 2$), when

$${{\omega }_{k}} = {{\omega }_{{{{k}_{1}}}}} + {{\omega }_{{{{k}_{2}}}}},\quad {\mathbf{k}} = {{{\mathbf{k}}}_{1}} + {{{\mathbf{k}}}_{2}},$$

(1)

or in the case of four-wave interaction, the processes of scattering ($2 \to 2$) with the resonance condition

$${{\omega }_{k}} + {{\omega }_{{{{k}_{1}}}}} = {{\omega }_{{{{k}_{2}}}}} + {{\omega }_{{{{k}_{3}}}}},\quad {\mathbf{k}} + {{{\mathbf{k}}}_{1}} = {{{\mathbf{k}}}_{2}} + {{{\mathbf{k}}}_{3}}.$$

(2)

From these resonance conditions immediately f-ollows that for the acoustic type waves, when ${{\omega }_{k}}$ depends linearly on the wave vector k, the resonance conditions (1) and (2) are satisfied automatically if the wave vectors ${{{\mathbf{k}}}_{i}}$ of all interacting waves are collinear. This means that all resonant waves propagate in the same direction at the same speed. Strictly speaking, the kinetic equation cannot be applied in such a situation. It is necessary to take into account the dispersion of waves, which should ensure the applicability of the weak turbulence theory. In the dispersion law ${{\omega }_{k}}$, one needs to account the next cubic term in the long-wavelength limit. In isotropic media, this expansion is w-ritten as

$${{\omega }_{k}} = k{{c}_{{\text{s}}}}(1 \pm {{a}^{2}}{{k}^{2}} + ...),$$

(3)

where c_s is the speed of sound and $a$ is the scale parameter characterizing the wave dispersion. In the case of the sign + in (3), one says about positive dispersion, the different sign corresponds to negative dispersion. Resonance condition (1) is satisfied only for positive dispersion, and accordingly is not satisfied for a different sign of the dispersion (the latter means that with a weak nonlinearity, a four-wave process becomes main, when condition (2) is fulfilled).

In this work, we consider the case of positive dispersion for three-dimensional ($D = 3$) acoustic turbulence. As clarified in [4, 5], the kinetic equation for $D = 3$ has a scale-invariant solution of the Ko-lmogorov type in the long-wavelength limit, independent of the dispersion parameter. This isotropic solution is the Zakharov–Sagdeev turbulence spectrum $E(k) \propto {{k}^{{ - 3/2}}}$, which describes a direct cascade with a constant energy flux ε from the long-wave region to the region of high wavenumbers. The fact that the spectrum is independent of the dispersion parameter is related to the integrable singularity in the kinetic equation due to the presence of two δ-functions, i.e., the resonances (1). In the two-dimensional case, in the region of small k, the singularity is non-integrable. As recently found out in [6], the spectrum in this case is a power-law one, $\propto {\kern 1pt} {{k}^{{ - 1}}}$, but depends explicitly on the dispersion parameter a.

We present the results of direct numerical simulation of three- dimensional acoustic turbulence in a medium with weak positive dispersion. It is shown that jets in the form of narrow cones appear in the energy distribution in the k-space in the region close to pumping. At high wavenumbers, the cones broaden, and the distribution tends to isotropic. In this region of wavenumbers, the angle-averaged turbulence spectrum acquires a power-law character, $E(k) \propto {{k}^{{ - \alpha }}}$, with an exponent close to 3/2, which corresponds to the Zakharov–Sagdeev spectrum of weak acoustic turbulence [4].

BASIC EQUATIONS

Direct numerical simulation of acoustic turbulence was carried out in the framework of the nonlinear string equation for a scalar function u depending on three spatial coordinates ${\mathbf{r}} = \{ x,y,z\}$ and time t:

$${{u}_{{tt}}} = \Delta u - 2{{a}^{2}}{{\Delta }^{2}}u + \Delta ({{u}^{2}}),$$

(4)

where a is the dispersion parameter introduced above and Δ is the Laplace operator. Note that this equation in the one-dimensional case refers to equations integrable by the inverse scattering transform [7]. In 3D geometry, this model was first used by Zakharov to study weak acoustic turbulence in [5]. In the linear approximation, Eq. (4) has the dispersion law

$${{\omega }^{2}} = {{k}^{2}} + 2{{a}^{2}}{{k}^{4}},\quad k = {\text{|}}{\mathbf{k}}{\text{|}},$$

(5)

coinciding with the Bogoliubov spectrum for the condensate oscillations of a weakly non-ideal Bose gas. The speed of sound c_s in this expression is equal to 1. In the case of weak dispersion $ka \ll 1$, this law is transformed into Eq. (3).

Equation (4) refers to Hamiltonian systems, it can be represented as a system of two equations:

$${{u}_{t}} = \frac{{\delta H}}{{\delta \phi }},\quad {{\phi }_{t}} = - \frac{{\delta H}}{{\delta u}},$$

(6)

where ϕ has the meaning of the hydrodynamic potential, and the Hamiltonian H is written as

$$\begin{gathered} H = \frac{1}{2}\int \left[ {{{{\left( {\nabla \phi } \right)}}^{2}} + {{u}^{2}}} \right]d{\mathbf{r}} + \int {{a}^{2}}{{(\nabla u)}^{2}}d{\mathbf{r}} + \frac{1}{3}\int {{u}^{3}}d{\mathbf{r}} \\ \equiv {{H}_{1}} + {{H}_{2}} + {{H}_{3}}. \\ \end{gathered} $$

(7)

The Hamiltonian (7) has three terms. The first ${{H}_{1}}$ is the sum of the kinetic and potential energies of linear dispersionless waves. The second term ${{H}_{2}}$ is responsible for the dispersive part of the wave energy, and ${{H}_{3}}$ describes the nonlinear interaction of waves.

Making the spatial Fourier transform and introducing the normal variables ${{a}_{k}}$ and $a_{k}^{ * }$,

$${{u}_{k}} = {{\left( {\frac{{{{k}^{2}}}}{{2{{\omega }_{k}}}}} \right)}^{{1/2}}}({{a}_{k}} + a_{{ - k}}^{ * }),$$

$${{\phi }_{k}} = - i{{\left( {\frac{{{{\omega }_{k}}}}{{2{{k}^{2}}}}} \right)}^{{1/2}}}({{a}_{k}} - a_{{ - k}}^{ * }){\kern 1pt} {\kern 1pt} .$$

Equations (6) take the standard form [8]

$$\frac{{\partial {{a}_{k}}}}{{\partial t}} = - i\frac{{\delta H}}{{\delta a_{k}^{ * }}},$$

(8)

where

$$\begin{gathered} H = \int {{\omega }_{k}}{\text{|}}{{a}_{k}}{{{\text{|}}}^{2}}d{\mathbf{k}} + \frac{1}{2}\int {{V}_{{{{k}_{1}}{{k}_{2}}{{k}_{3}}}}}\left( {a_{{{{k}_{1}}}}^{ * }{{a}_{{{{k}_{2}}}}}{{a}_{{{{k}_{3}}}}}} \right. \hfill \\ \left. { + \;{{a}_{{{{k}_{1}}}}}a_{{{{k}_{2}}}}^{ * }a_{{{{k}_{3}}}}^{ * }} \right)\delta \left( {{{{\mathbf{k}}}_{1}} - {{{\mathbf{k}}}_{2}} - {{{\mathbf{k}}}_{3}}} \right)d{{{\mathbf{k}}}_{1}}d{{{\mathbf{k}}}_{2}}d{{{\mathbf{k}}}_{3}}. \hfill \\ \end{gathered} $$

In this Hamiltonian, we left only one resonance term corresponding to decay processes (1) in the nonlinear term.

In the weak dispersion approximation, we will take into account the dispersion in the quadratic Hamiltonian H,

$${{\omega }_{k}} = k(1 + {{a}^{2}}{{k}^{2}}),$$

but in the matrix element ${{V}_{{{{k}_{1}}{{k}_{2}}{{k}_{3}}}}}$, it will be neglected:

$${{V}_{{{{k}_{1}}{{k}_{2}}{{k}_{3}}}}} = \frac{3}{{4{{\pi }^{{3/2}}}}}{{\left( {{{k}_{1}}{{k}_{2}}{{k}_{3}}} \right)}^{{1/2}}}.$$

Hence, the kinetic equation for the pair correlator ${{n}_{k}}$ ($\langle a_{k}^{ * }{{a}_{{{{k}_{1}}}}}\rangle = {{n}_{k}}\delta ({\mathbf{k}} - {{{\mathbf{k}}}_{1}})$) in the weak turbulence approximation is written as

$$\frac{{\partial {{n}_{k}}}}{{\partial t}} = 2\pi \int d{{{\mathbf{k}}}_{1}}d{{{\mathbf{k}}}_{2}}\left( {{{T}_{{k{{k}_{1}}{{k}_{2}}}}} - {{T}_{{{{k}_{1}}k{{k}_{2}}}}} - {{T}_{{{{k}_{2}}k{{k}_{1}}}}}} \right),$$

(9)

where:

$$\begin{gathered} {{T}_{{k{{k}_{1}}{{k}_{2}}}}} = {\text{|}}{{V}_{{k{{k}_{1}}{{k}_{2}}}}}{{{\text{|}}}^{2}}({{n}_{{{{k}_{1}}}}}{{n}_{{{{k}_{2}}}}} - {{n}_{k}}{{n}_{{{{k}_{2}}}}} \\ - \;{{n}_{k}}{{n}_{{{{k}_{1}}}}})\delta \left( {{\mathbf{k}} - {{{\mathbf{k}}}_{1}} - {{{\mathbf{k}}}_{2}}} \right)\delta \left( {{{\omega }_{k}} - {{\omega }_{{{{k}_{1}}}}} - {{\omega }_{{{{k}_{2}}}}}} \right). \\ \end{gathered} $$

(10)

The turbulence spectrum $E(k)$, i.e., the dependence of the energy on the absolute value of k, is found from the solution of this equation after averaging the quantity ${{\omega }_{k}}{{n}_{k}}{{k}^{2}}$ over the entire solid angle Ω:

$$E(k) = {{k}^{2}}{{\omega }_{k}}\int {{n}_{k}}d\Omega .$$

For isotropic distributions, obviously $E(k) = 4\pi {{k}^{2}}{{\omega }_{k}}{{n}_{k}}$.

As was first noted by Zakharov [5], in the kinetic Eq. (9) in the case of isotropic distributions, the dispersion contribution in ${{\omega }_{k}}$ can be neglected, despite the presence of the product of two delta functions with respect to frequencies and wave vectors in the collision term giving a singularity in the kinetic equation. This singularity in the kinetic equation turns out to be integrable after averaging over the angles. As a result, the kinetic equation admits a stationary power-law solution: ${{n}_{k}} \propto$ ${{k}^{\alpha }}$. The exponent α for the Kolmogorov-type spectrum is found using the Zakharov transformations (see [1]): $\alpha = - 11{\text{/}}2$, which corresponds to the Zakharov–Sagdeev spectrum [4]:

$$E(k) = C{{\varepsilon }^{{1/2}}}{{k}^{{ - 3/2}}}.$$

(11)

Here C is the Kolmogorov–Zakharov constant, and ε is the energy dissipation rate per unit volume, which is the energy flux over the spectrum. The power dependence on ε in the spectrum ((11)) with the exponent 1/2 corresponds to the resonant three-wave interaction.

The existence of the Zakharov–Sagdeev spectrum in the inertial interval as a spectrum of the Kolmogorov type was confirmed in a number of papers [9–11] by numerically solving the kinetic Eq. (9) in the presence of long-wavelength pumping and high-frequency damping. It is important to note that the numerical study of weak acoustic turbulence in the framework of the kinetic Eq. (9) in the three-dimensional case was carried out only for isotropic distributions. In this paper, we will show that in direct numerical simulation of three-dimensional acoustic turbulence described by Eq. (4) supplemented by damping at high k and pumping in the region of long waves, the structure of the spectra is not isotropic, especially in the region of small k.

With account of both pumping and damping terms, Eqs. (6) are written in the form:

$${{u}_{t}} = - \Delta \phi + \mathcal{F}({\mathbf{k}},t) - {{\gamma }_{k}}u,$$

(12)

$${{\phi }_{t}} = - u + 2{{a}^{2}}\Delta u - {{u}^{2}},$$

(13)

where the operator ${{\gamma }_{k}}$ responsible for dissipation and the forcing term $\mathcal{F}({\mathbf{k}},t)$ are given in Fourier space as:

$${{\gamma }_{k}} = 0,\quad k \leqslant {{k}_{d}},$$

$${{\gamma }_{k}} = {{\gamma }_{0}},\quad k > {{k}_{d}},$$

$$\mathcal{F}({\mathbf{k}},t) = F(k)\exp [iR({\mathbf{k}},t)],$$

$$F(k) = {{F}_{0}}\exp [ - {{(k - {{k}_{1}})}^{4}}{\text{/}}k_{2}^{4}],\quad k \leqslant {{k}_{2}},$$

$$F(k) = 0,\quad k > {{k}_{2}}.$$

Here, $R({\mathbf{k}}{\mathbf{,}}t)$ are random numbers uniformly distributed in the interval $[0,2\pi ]$, ${{\gamma }_{0}}$ and ${{F}_{0}}$ are constants. The quantity ${{k}_{1}}$ determines the wavenumber on which the maximum pumping amplitude is reached, ${{k}_{2}}$ sets its width, and ${{k}_{3}}$ corresponds to the scale at which dissipation occurs.

NUMERICAL SCHEME AND PARAMETERS

Numerical integration of the system of Eqs. (12) and (13) was carried out in a cubic domain ${{(2\pi )}^{3}}$ with periodic boundary conditions over all three coordinates. Time integration was performed using explicit scheme by means of the Runge–Kutta method of the fourth order accuracy with step $\delta t = 2.5 \times {{10}^{{ - 3}}}$. The equations were integrated over spatial coordinates using spectral methods with the total number of harmonics ${{N}^{3}} = {{512}^{3}}$. To suppress the aliasing effect, we used a filter that nulls higher harmonics with a wavenumber above ${{k}_{a}} \geqslant N{\text{/}}3$. Below we present results of numerical simulation for the following parameters: ${{k}_{d}} = 125$, ${{k}_{1}} = 3$, $k_{2}^{4} = 6$, ${{\gamma }_{0}} = 100$, $a = 2.5 \times {{10}^{{ - 3}}}$, ${{F}_{0}} = 5 \times {{10}^{5}}$. With this choice of parameters, the inertial interval was more than one decade. The maximum dispersion addition at the end of the inertial interval, $k = {{k}_{d}}$, was ${{({{k}_{d}}a)}^{2}} \approx 0.1$.

SIMULATION RESULTS

In a numerical experiment with the above parameters, we observed a transition to weak turbulence regime. Figure 1 shows how the total energy of the system (7) evolves. It can be seen the rather quick transition to the regime of quasi-stationary chaotic motion when the effect of pumping in the region of small k is completely compensated by dissipative effects. The inset to Fig. 1 shows the time dependences of the dispersive part of the energy ${{H}_{2}}$ and the nonlinear interaction energy ${{H}_{3}}$. Both contributions ${{H}_{2}}$ and ${{H}_{3}}$ turn out to be small compared to the total energy of the system (respectively, ${{H}_{1}}$). The dispersive part of the energy exceeds the energy of the nonlinear interaction by almost an order of magnitude, which indicates on the realization of a weakly nonlinear regime. Thus, the total energy in the inertial interval is approximately equal to ${{H}_{1}} \approx \int {{\epsilon }_{k}}d{\mathbf{k}}$, where ${{\epsilon }_{k}} = k{\text{|}}{{a}_{k}}{{{\text{|}}}^{2}}$ is the wave energy density in k-space.

Fig. 1.

(Color online) Time dependence of the total energy of the system (7) for $a = 2.5 \times {{10}^{{ - 3}}}$. The inset shows the time dependences of the dispersive part of the energy ${{H}_{2}}$ and the nonlinear interaction energy ${{H}_{3}}$.

The behavior of the spectrum of space-time Fourier transform of the function $u({\mathbf{r}},t)$ shown in Fig. 2 also testifies to the weakly nonlinear character of wave propagation. The figure shows that the wave energy is concentrated along the linear dispersion relation (5). Line broadening is due to nonlinearity. For almost the entire inertial interval, this broadening does not exceed the linear dispersion. For small k, the broadening is comparable to the dispersion. For higher k, the dispersion exceeds the nonlinear broadening, which agrees with the ratio of the corresponding contributions ${{H}_{2}}$ and ${{H}_{3}}$.

Fig. 2.

(Color online) Space–time Fourier transform ${\text{|}}u({\mathbf{k}},\omega ){{{\text{|}}}^{2}}$ in a logarithmic scale. The black dotted line corresponds to the exact value of the dispersion curve (5), the white dotted line corresponds to the non-dispersive wave propagation, $\omega = {\text{|}}{\mathbf{k}}{\text{|}}$.

The numerical experiment shows that after the system enters the quasi- stationary state, the behavior of $u({\mathbf{r}})$ acquires a complex (turbulent) character. In Fig. 3 this behavior demonstrates the dependence of the function $u({\mathbf{r}})$ in the $z = 0$ plane for the quasi-stationary state at the time $t = 2500$. At the same time, the distribution of the energy density ${{\epsilon }_{k}}$ of turbulent fluctuations in the $k$-space is not isotropic. The anisotropy is especially pronounced in the region of small wavenumbers near the pumping.

Fig. 3.

(Color online) Section of the function $u({\mathbf{r}})$ by the $z = 0$ plane is shown at the time $t = 2500$ corresponding to the quasi-stationary state.

In Fig. 4 we present three isosurfaces of the function ${\text{|}}{{u}_{{\mathbf{k}}}}{\text{|}}{\kern 1pt} ( = \epsilon _{k}^{{1/2}})$. As seen, in the region of small wavenumbers, structures with a large number of jets in the form of narrow cones appear in the distribution of turbulent fluctuations. The onset of such structures is the result of resonant wave interactions (1) at very small k close to the pumping region when dispersion can be neglected. As k increases, the cones broaden and the distribution tends to be isotropic, see Fig. 5. In this figure, the blue color (at $k \geqslant 30$) shows a tendency to spectrum isotropization, which is associated with an increase in dispersion with growing $k$ and accordingly with an angular broadening of the resonant surface (1) by an angle of the order of $ka$.

Fig. 4.

(Color online) Isosurfaces of the Fourier spectrum of ${\text{|}}{{u}_{{\mathbf{k}}}}{\text{|}} = {{u}_{0}}$, (a), (b) and (c) correspond to the values ${{u}_{0}} = 5 \times {{10}^{{ - 5}}}$, 0.5× 10^–5, and 0.25× 10^–5, respectively, $t = 2500$.

Fig. 5.

(Color online) Fourier spectrum of ${\text{|}}{{u}_{k}}{\text{|}} \equiv \epsilon _{k}^{{1/2}}$ in section ${{k}_{z}} = 0$ (logarithmic scale), $t = 2500$.

The generation of jets is associated with two possible causes: linear and non-linear. The first one is the discreteness of the wave vector lattice in the pumping region, $1 \leqslant k \leqslant 6$, which inevitably leads to a small anisotropy in the excitation and differences in the growth of perturbations at the initial stage. The second one is the tendency of the dispersion to zero at small k caused by the equations itself; the three-wave resonance conditions are satisfied for an arbitrary beam. The beams that form the jets have an advantage over other beams. This process, the cooperation of rays into a jet, has a clearly nonlinear character. In particular, this fact follows from the numerical simulation results presented in Figs. 4 and 5, for which the contrasts in intensity in the jets and the regions between them are significant: the difference reaches two orders of magnitude. In our opinion, such a jump in intensity cannot be explained only by a small anisotropy in pumping, but has a nonlinear origin, possibly due to the collapse of sound waves described by the three-dimensional Kadomtsev–Petviashvili equation (see review [12]). It should also be noted that numerical experiments show that similar jets appear at the onset of developed hydrodynamic turbulence [13]. However, this question is beyond the scope of this work: new numerical experiments with a high spatial resolution are required.

To find the turbulence spectrum $E(k)$, as noted above, it is necessary to integrate the expression ${{k}^{2}}{{\epsilon }_{k}}$ over the entire solid angle. Figure 6 shows the result of this averaging over the angle for the time corresponding to the quasi-stationary state. It is clearly seen from Fig. 6 that in the regime of quasi-stationary chaotic motion, the spectrum of $E(k)$ acquires a power law behavior. Recall that the inertial interval was $k \in [6,100]$. As can be seen from Fig. 6, there are two regions with different behavior of the spectrum in the inertial interval. In the region of high wavenumbers, the spectrum of weak acoustic turbulence coincides with the Zakharov–Sagdeev spectrum (11) with high accuracy, and in the long-wave region, deviations from this spectrum are observed, which, in our opinion, arise due to jets, whose role is significant at small $k$. It should be noted that similar large deviations of an oscillatory nature from the Zakharov–Sagdeev spectrum were observed numerically when simulating the condensate in the framework of the Gross–Pitaevskii equation [14]. In our opinion, these deviations can be related to the anisotropy caused by the presence of jets. No jets were found in the experiment [14]. In our numerical experiment, the spectral energy density ${{\epsilon }_{k}}$ approaches an isotropic distribution in the region of high k as it follows from Figs. 4 and 5. The dispersion ${{(ka)}^{2}}$ changes from 10^–3 to 10^–1 in this region, i.e., remains weak. These two factors contribute to the formation of the Zakharov–Sagdeev spectrum. It should be specially emphasized that the found turbulence spectrum is far from the Kadomtsev–Petviashvili spectrum ${{E}_{{KP}}} \propto {{k}^{{ - 2}}}$ [15] related to the dispersionless limit ($a = 0$).

Fig. 6.

(Color online) Turbulence spectrum $E(k)$ measured in the quasi-stationary state, the black dotted line corresponds to the Zakharov–Sagdeev spectrum (11), the red solid line corresponds to the Kadomtsev–Petviashvili spectrum.

CONCLUSIONS

In the present work, direct numerical simulation of three-dimensional acoustic turbulence in a medium with weak positive dispersion is carried out, taking into account the pumping and dissipation of energy. We have established that the system of nonlinear interacting weakly dispersive waves quickly enough passes into the quasi-stationary chaotic state, which is a developed wave turbulence. In the quasi-stationary regime, the wave field acquires a complex, chaotic character. In the long-wavelength region, close to pumping, we observed in the turbulence spectrum the appearance of narrow jets in the form of cones, which expand upon transition to the short-wave region. The latter leads to the fact that the spectral energy density ${{\epsilon }_{k}}$ tends in the region of high k to an isotropic distribution, for which the dispersion remains weak. In this range of scales, the turbulence spectrum calculated in the stationary state agrees with a high accuracy with the analytical Zakharov–Sagdeev spectrum of weak acoustic turbulence. It is numerically demonstrated that the criteria of weak turbulence are fully satisfied for this spectrum. Note that the result obtained is the first reliable observation of the spectrum of weak turbulence of acoustic waves in media with positive dispersion in direct three-dimensional numerical simulations.