\newpage 
\chapter{Слабая волновая турбулентность волн на 
поверхности жидкости}
\label{Weak_Turbulence}
  
В этой главе представлены результаты, опубликованные с работах
\cite{OptTurb, Capil, WeakGrav, WeakTurbGrav2004}

Одним из краеугольных камней в теории слабой турбулентности является резонансное
взаимодействие волн. Как было показано в работе \cite{Kartasheva-1998} и
развито в \cite{Pushkarev-2001}, условия резонанса на характерной
для численных экспериметов дискретной сетке не могут быть выполнены точно. Однако, как
было отмечено выше, в настоящий момент прямое численное моделирование динамических
уравнений является важной задачей. Это один из основных инструментов проверки правильности
предположений и допущений, в рамках которых была получена теория слабой турбулентности.
Поэтому перед нами с очевидностью встают следующие вопросы:
\begin{itemize}
\item Как дискретная сетка волновых чисел в численном моделировании влияет на резонансное
взаимодействие?
\item Может ли нелинейный сдвиг частоты уширить резонансное множество таким образом,
чтобы дискретность перестала быть важной?
\end{itemize}
В данной главе  изучаются эти вопросы для капиллярных и 
гравитационных волн на поверхности
несжимаемой жидкости бесконечной глубины.

Эта глава построена следующим образом. 
В разделе \ref{INTRO4} мы запишем динамические уравнения в виде 
ряда по степеням крутизны волн.
В разделе \ref{RESON} выписаны уравнения на дискретной сетке.
Случай капиллярных волн на поверхности глубокой воды рассмотрен в 
разделе \ref{CapillarRes}.
Произведено численное
моделирование резонансного распада монохроматической волны. Получена картина вторичных
распадов, формирующих Колмогоровский каскад.
В разделе \ref{GravRes} приведены результаты моделирования
распада монохроматической гравитационной волны на поверхности бесконечно глубокой
идеальной жидкости. 
В параграфе \ref{GravTurbu} приводятся результаты численного моделирования слаботурбулентной
системы гравитационных волн, возбуждаемой стохастической внешней силой.
В заключении \ref{GravTurbConcl} приведен анализ полученных 
результатов и выводы,
полезные для построения численной модели слабой турбулентности гравитационных 
поверхностных волн.

\section{Введение}
\label{INTRO4}

Рассмотрим безвихревое потенциальное течение идеальной несжимаемой жидкости
бесконечной глубины. Для простоты будем считать плотность жидкости $\rho = 1$.
Потенциал скорости $\phi$ при этих условиях удовлетворяет уравнению Лапласа:
\begin{equation}
\label{Laplas}
\triangle \phi = 0,
\end{equation}
в области со следующими границами:
\begin{equation}
-\infty < z < \eta (\vec r), \;\;\;\; \vec r = (x,y).
\end{equation}
Граничные условия на потенциал таковы
\begin{equation}
\label{Laplas_boundary}
\begin{array}{c}
\displaystyle
\left. 
\frac {\partial \eta}{\partial t}+
\frac {\partial \phi}{\partial x}\frac {\partial \eta}{\partial x} +
\frac {\partial \phi}{\partial y}\frac {\partial \eta}{\partial y}
= \frac{\partial \phi}{\partial z} \right|_{z= \eta},\\
\displaystyle
\left. \left ( 
\frac {\partial \phi}{\partial t}+
\frac{1}{2}( \nabla \phi )^2
\right ) \right |_{z= \eta} +\sigma (\sqrt{1 + (\nabla \eta)^2} - 1) = 0,
\end{array}
\end{equation}
при $z=\eta$, и:
\begin{equation}
\phi_z|_{z= - \infty} = 0,
\end{equation}
при $z\to -\infty$.
Здесь $\eta = \eta (x,y,t)$ смещение поверхности от состояния покоя.
Энергия системы описывается следующей функцией:
$$
H = T + U,
$$
\begin{equation}
T = \frac{1}{2} \int d^2 r \int \limits_{-\infty}^{\eta} (\nabla \phi)^2 d z,
\end{equation}
\begin{equation}
U = \sigma \int (\sqrt{1 + (\nabla \eta)^2} - 1)d^2 r + \frac{1}{2}g \int \eta^2 d^2 r,
\end{equation}
где $\sigma$ -- коэффициент поверхностного натяжения,
$g$ -- ускорение свободного падения.
В \cite{Z68} было показано, что эта система является гамильтоновой. 
Гамильтоновы переменные в этом случае --- это смещение поверхности
$\eta(x,y,t)$ и потенциал скосроти на поверхности жидкости
$\psi (x,y;t) = \phi (x,y,\eta(x,y;t);t)$. Уравнения Гамильтона имеют следующий вид:
\begin{equation}
\label{Hamiltonian_equations}
\frac{\partial \eta}{\partial t} = \frac{\delta H}{\delta \psi}, \;\;\;\;
\frac{\partial \psi}{\partial t} = - \frac{\delta H}{\delta \eta}.
\end{equation}
Применяя разложение энергии по крутизне, можно получить следующее
выражение для Гамильтониана:
\begin{equation}
\label{Hamiltonian}
\begin{array}{l}
\displaystyle
H = \frac{1}{2}\int\left(  \sigma |\nabla \eta|^2 + g \eta^2 + \psi \hat k  \psi \right) d^2 r +
\frac{1}{2}\int\eta\left[ |\nabla \psi|^2 - (\hat k \psi)^2 \right] d^2 r + \\
\displaystyle
+ \frac{1}{2}\int\eta (\hat k \psi) \left[ \hat k (\eta (\hat k \psi)) + \eta\triangle\psi \right] d^2 r.
\end{array}
\end{equation}
Здесь введен линейный оператор $\hat k$, суть действия которого на функцию
сводится к следующему: надо взять фурье-компоненты функции, умножить их
на модуль $k$ и выполнить обратное преобразование Фурье.
Для правильного описания нелинейного взаимодействия гравитационных волн необходимо использовать
разложение  Гамильтониана до членов четвертого порядка включительно.
Уравнения Гамильтона (\ref{Hamiltonian_equations}) принимают вид:
\begin{equation}
\label{eta_psi_system}
\begin{array}{lcl}
\displaystyle
\dot \eta &=& \hat k  \psi - (\nabla (\eta \nabla \psi)) - \hat k  [\eta \hat k  \psi] +
\hat k (\eta \hat k  [\eta \hat k  \psi]) + \\
\displaystyle
&&+ \frac{1}{2} \triangle [\eta^2 \hat k \psi] + \frac{1}{2} \hat k [\eta^2 \triangle\psi],\\
\displaystyle
\dot \psi &=& \sigma \triangle \eta - g\eta - \frac{1}{2}\left[ (\nabla \psi)^2 - (\hat k \psi)^2 \right] -
[\hat k  \psi] \hat k  [\eta \hat k  \psi] - [\eta \hat k  \psi]\triangle\psi.
\end{array}
\end{equation}
Удобно записать гамильтониан в фурье-гармониках:
\begin{equation}
\label{FourierTrans}
\psi_{\vec k} = \frac{1}{2\pi} \int \psi_{\vec r} e^{i {\vec k} {\vec r}} d^2 r,\;\;
\eta_{\vec k} = \frac{1}{2\pi} \int \eta_{\vec r} e^{i {\vec k} {\vec r}} d^2 r.
\end{equation}
Члены гамильтониана примут следующий вид:
\begin{equation}
\label{Ham_eta_psi_k}
\begin{array}{l}
\displaystyle
H = H_0 + H_1 + H_2 + ...,\\
\displaystyle
H_0 = \frac{1}{2}\int (|k| |\psi_{\vec k}|^2 + \sigma |k|^2|\eta_{\vec k}|^2 + g |\eta_{\vec k}|^2) d\vec k,\\
\displaystyle
H_1 = -\frac{1}{4\pi}\int L_{\vec k_1\vec k_2 }\psi_{\vec k_1}\psi_{\vec k_2}\eta_{\vec k_3}
\delta (\vec k_1 + \vec k_2 + \vec k_3) d\vec k_1 d\vec k_2 d\vec k_3,\\
\displaystyle
H_2 = \frac{1}{16\pi^2}\int M_{\vec k_1\vec k_2\vec k_3\vec k_4}
\eta_{\vec k_1}\psi_{\vec k_2}\eta_{\vec k_3}\psi_{\vec k_4}
\delta (\vec k_1 + \vec k_2 + \vec k_3 + \vec k_4) d\vec k_1 d\vec k_2 d\vec k_3 d\vec k_4,\\
\end{array}
\end{equation}
Здесь введены матричные элементы
\begin{equation}
\begin{array}{l}
\displaystyle
L_{\vec k_1 \vec k_2} = (\vec k_1 \vec k_2) + |\vec k_1||\vec k_2|,\\
\displaystyle
M_{\vec k_1\vec k_2\vec k_3 \vec k_4} = |\vec k_2| |\vec k_4|\left[
\frac{1}{2}(|\vec k_1 + \vec k_2| + |\vec k_1 + \vec k_4| + \right.\\
\displaystyle
\left. +|\vec k_3 + \vec k_2| + |\vec k_3 + \vec k_4|) -
|\vec k_2| - |\vec k_4|\right].\\
\end{array}
\end{equation}
Вследствие вещественности функций $\eta$ и $\psi$ их фурье-образы обладают
симметрией $\eta_{\vec k} = \eta_{-\vec k}^{*}$. Поэтому удобнее ввести канонические
переменные:
\begin{equation}
\label{a_k_substitution}
a_{\vec k} = \sqrt \frac{\omega_k}{2k} \eta_{\vec k} + i \sqrt \frac{k}{2\omega_k} \psi_{\vec k},
\end{equation}
здесь $\omega_k$ --- это закон дисперсии волн:
\begin{equation}
\omega_k = \sqrt {\sigma k^3 + g k}.
\end{equation}
Члены Гамильтониана (\ref{Ham_eta_psi_k}) принимают вид:
\begin{equation}
\label{Ham_a_k}
\begin{array}{lll}
\displaystyle
H_0 = \int \omega_k |a_{\vec k}|^2 d {\vec k},\\
\displaystyle
H_1 = \frac{1}{6}\frac{1}{2\pi}\int E_{\vec k_1 \vec k_2}^{\vec k_0} 
(a_{\vec k_1}a_{\vec k_2}a_{\vec k_0} + a_{\vec k_1}^{*}a_{\vec k_2}^{*}a_{\vec k_0}^{*})
\delta (\vec k_1 + \vec k_2 + \vec k_0) d {\vec k_1}d {\vec k_2}d {\vec k_0} +\\
\displaystyle
+\frac{1}{2}\frac{1}{2\pi}\int C_{\vec k_1 \vec k_2}^{\vec k_0} 
(a_{\vec k_1}a_{\vec k_2}a_{\vec k_0}^{*} + a_{\vec k_1}^{*}a_{\vec k_2}^{*}a_{\vec k_0})
\delta (\vec k_1 + \vec k_2 - \vec k_0) d {\vec k_1} d {\vec k_2}d {\vec k_0},\\
\displaystyle
H_2 = \frac{1}{4}\frac{1}{(2\pi)^2}\int W_{\vec  k_1 \vec k_2 \vec k_3 \vec k_4}
(a_{\vec k_1}a_{\vec k_2}a_{\vec k_3}a_{\vec k_4}
+ a_{\vec k_1}^{*}a_{\vec k_2}^{*}a_{\vec k_3}^{*}a_{\vec k_4}^{*})\times\\
\displaystyle
\times\delta (\vec k_1 + \vec k_2 + \vec k_3 + \vec k_4) d {\vec k_1} d {\vec k_2} d {\vec k_3} d {\vec k_4}+\\
\displaystyle
+\frac{1}{4}\frac{1}{(2\pi)^2}\int F_{\vec  k_1 \vec k_2 \vec k_3 \vec k_4}
(a_{\vec k_1}^{*}a_{\vec k_2}a_{\vec k_3}a_{\vec k_4}
+ a_{\vec k_1}a_{\vec k_2}^{*}a_{\vec k_3}^{*}a_{\vec k_4}^{*})\times\\
\displaystyle
\times\delta (\vec k_1 - \vec k_2 - \vec k_3 - \vec k_4) d {\vec k_1} d {\vec k_2} d {\vec k_3} d {\vec k_4}+\\
\displaystyle
+\frac{1}{4}\frac{1}{(2\pi)^2}\int D_{\vec  k_1 \vec k_2 \vec k_3 \vec k_4}
(a_{\vec k_1}a_{\vec k_2}a_{\vec k_3}^{*}a_{\vec k_4}^{*}
\delta (\vec k_1 + \vec k_2 - \vec k_3 - \vec k_4) d {\vec k_1} d {\vec k_2} d {\vec k_3} d {\vec k_4}.
\end{array}
\end{equation}
Здесь введены такие обозначения матричных элементов:
\begin{equation}
\begin{array}{l}
\displaystyle
E_{\vec k_1 \vec k_2}^{\vec k_0} = V_{\vec k_1 \vec k_2}^{k_0} + V_{\vec k_0 \vec k_2}^{k_1} + V_{\vec k_0 \vec k_1}^{k_2},\\
\displaystyle
C_{\vec k_1 \vec k_2}^{\vec k_0} = V_{\vec k_1 \vec k_2}^{k_0} - V_{-\vec k_0 \vec k_2}^{k_1} - V_{-\vec k_0 \vec k_1}^{k_2},\\
\displaystyle
V_{\vec k_1 \vec k_2}^{k_0} = \sqrt {\frac{\omega_{k_1} \omega_{k_2} k_0}{8 k_1 k_2 \omega_{k_0}}} L_{\vec k_1 \vec k_2},\\
\displaystyle
F_{\vec k_1 \vec k_2 \vec k_3 \vec k_4} = W_{\vec  k_1, -\vec k_2, -\vec k_3, -\vec k_4} 
- W_{\vec  -k_2, \vec k_1, -\vec k_3, -\vec k_4} + W_{-\vec  k_3, -\vec k_2, \vec k_1, -\vec k_4} - \\
\displaystyle
- W_{\vec  -k_4, -\vec k_2, -\vec k_3, \vec k_1},\\
\displaystyle
D_{\vec k_1 \vec k_2 \vec k_3 \vec k_4} = W_{\vec  k_1, -\vec k_3, \vec k_2, -\vec k_4} 
+ W_{\vec  k_3, -\vec k_1, \vec k_4, -\vec k_2} -\\
\displaystyle
-\left(W_{\vec  k_1, \vec k_2, -\vec k_3, -\vec k_4} +
W_{\vec  k_1, \vec k_2, -\vec k_4, -\vec k_3} +
W_{\vec  k_2, \vec k_1, -\vec k_3, -\vec k_4} + W_{\vec  k_2, \vec k_1, -\vec k_4, -\vec k_3}\right),\\
\displaystyle
W_{\vec  k_1 \vec k_2 \vec k_3 \vec k_4} = \sqrt{\frac{k_1 \omega_{k_2} k_3 \omega_{k_4}}
{16 \omega_{k_1} k_2 \omega_{k_3} k_4}}M_{\vec k_1\vec k_2\vec k_3 \vec k_4}.
\end{array}
\end{equation}
Несмотря на кажущуюся громоздкость представления Гамильтониана в таких перемнных,
смысл каждого члена в выражении (~\ref{Ham_a_k}) ясен и понятен. Член $H_0$ отвечает
линейной эволюции начальных условий. Выражение с $H_1$ отвечает трехволновым взаимодействиям,
в том числе: член с $E_{\vec k_1 \vec k_2}^{\vec k_0}$ дает генерацию второй гармоники, а член с
$C_{\vec k_1 \vec k_2}^{\vec k_0}$ отвечает трехволновым процессам распада и слияния.
Выражение с $H_2$ представляет четырехволновые процессы, а именно: член с 
$W_{\vec  k_1 \vec k_2 \vec k_3 \vec k_4}$ дает генерацию третьей гармоники,
$F_{\vec k_1 \vec k_2 \vec k_3 \vec k_4}$ --- распад и слияние одной волны в три,
$D_{\vec k_1 \vec k_2 \vec k_3 \vec k_4}$ --- рассеяние двух волн в две волны.
Уравнения Гамильтона в этих переменных принимают вид:
\begin{equation}
\label{Eq_a_k}
\begin{array}{l}
\displaystyle
\dot a_{\vec k} = -i \frac{\delta H}{\delta a_{\vec k}^{*}} = -i \omega_k a_{\vec k} - \\
\displaystyle
-\frac{i}{2}\frac{1}{2\pi}\int C_{\vec k_1 \vec k_2}^{\vec k} 
a_{\vec k_1}a_{\vec k_2}\delta (\vec k_1 + \vec k_2 - \vec k) d {\vec k_1}d {\vec k_2}-\\
\displaystyle
-\frac{i}{2\pi}\int C_{\vec k \vec k_2}^{\vec k_0} 
a_{\vec k_2}^{*}a_{\vec k_0}\delta (\vec k + \vec k_2 - \vec k_0) d {\vec k_2}d {\vec k_0}-\\
\displaystyle
-\frac{i}{2}\frac{1}{2\pi}\int E_{\vec k_1 \vec k_2}^{\vec k} 
a_{\vec k_1}^{*}a_{\vec k_2}^{*}\delta (\vec k_1 + \vec k_2 + \vec k) d {\vec k_1}d {\vec k_2}-\\
\displaystyle
-\frac{i}{(2\pi)^2}\int W_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
a_{\vec k_1}^{*}a_{\vec k_2}^{*}a_{\vec k_3}^{*}
\delta (\vec k_1 + \vec k_2 + \vec k_3 + \vec k)d {\vec k_1}d {\vec k_2}d {\vec k_3}-\\
\displaystyle
-\frac{i}{4}\frac{1}{(2\pi)^2}\int F_{\vec  k \vec k_2 \vec k_3 \vec k_4}
a_{\vec k_2}a_{\vec k_3}a_{\vec k_4}
\delta (\vec k - \vec k_2 - \vec k_3 - \vec k_4)d {\vec k_2}d {\vec k_3}d {\vec k_4}-\\
\displaystyle
-\frac{i}{4}\frac{3}{(2\pi)^2}\int F_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
a_{\vec k_1}a_{\vec k_2}^{*}a_{\vec k_3}^{*}\delta (\vec k_1 - \vec k_2 - \vec k_3 - \vec k)
d {\vec k_1}d {\vec k_2}d {\vec k_3}-\\
\displaystyle
-\frac{i}{2}\frac{1}{(2\pi)^2}\int D_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
a_{\vec k_1}a_{\vec k_2}a_{\vec k_3}^{*}
\delta (\vec k_1 + \vec k_2 - \vec k_3 - \vec k)d {\vec k_1}d {\vec k_2}d {\vec k_3}.
\end{array}
\end{equation}

Изложенные результаты можно найти во множестве источников (например, 
\cite{Zakharov-1967, ZLF92}),
в том числе и для более общего случая жидкости конечной глубины \cite{Zakharov-1999}.



\section{Резонансные взаимодействия волн на дискретной сетке} 
\label{RESON}

Рассмотрим динамические уравнения (\ref{Eq_a_k}) в периодической области
с линейными размерами $L_x$ и $L_y$. Если мы находимся на дискретной сетке
с количеством точек по направлениям $N_x$ и $N_y$ соответственно, то волновой
вектор и координатный радиус вектор могут быть представлены как:
\begin{equation}
\label{krdiscr}
\vec k = \left( \frac{2\pi}{L_x}k_x,  \frac{2\pi}{L_y}k_y\right),
\;\;\; \vec r = \left( \frac{L_x}{N_x}r_x, \frac{L_y}{N_y}r_y\right),
\end{equation}
где $k_i$ пробегает дискретный ряд значений от $-N_i/2$ до $N_i/2$, а $r_i$ --- от $0$ до $N_i$.
Соответственно, все знаки интегрирования заменятся на знаки суммирования, а $\delta$-функция
Дирака --- на ее дискретный аналог, символ Кронекера $\triangle$. Для сокращения записи будем писать
просто знак суммы $\sum$, подразумевая, что суммирование производится по всем повторяющимся
символам. Тогда можем переписать уравнения (\ref{Eq_a_k}) в следующем виде:
\begin{equation}
\label{Eq_a_k_discr}
\begin{array}{l}
\displaystyle
\dot a_{\vec k} = -i \omega_k a_{\vec k} - \\
\displaystyle
-\frac{i}{2}\frac{2\pi}{L_x L_y}\sum C_{\vec k_1 \vec k_2}^{\vec k} 
a_{\vec k_1}a_{\vec k_2}\triangle_{\vec k_1 + \vec k_2, \vec k}-\\
\displaystyle
-\frac{i2\pi}{L_x L_y}\sum C_{\vec k \vec k_2}^{\vec k_0} 
a_{\vec k_2}^{*}a_{\vec k_0}\triangle_{\vec k + \vec k_2, \vec k_0}-\\
\displaystyle
-\frac{i}{2}\frac{2\pi}{L_x L_y}\sum E_{\vec k_1 \vec k_2}^{\vec k} 
a_{\vec k_1}^{*}a_{\vec k_2}^{*}\triangle_{\vec k_1 + \vec k_2+ \vec k, 0}-\\
\displaystyle
-\frac{i(2\pi)^4}{(L_x L_y)^3}\sum W_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
a_{\vec k_1}^{*}a_{\vec k_2}^{*}a_{\vec k_3}^{*}
\triangle_{\vec k_1 + \vec k_2+ \vec k_3 + \vec k, 0}-\\
\displaystyle
-\frac{i}{4}\frac{(2\pi)^4}{(L_x L_y)^3}\sum F_{\vec  k \vec k_2 \vec k_3 \vec k_4}
a_{\vec k_2}a_{\vec k_3}a_{\vec k_4}
\triangle_{\vec k , \vec k_2 + \vec k_3 + \vec k_4}-\\
\displaystyle
-\frac{i}{4}\frac{3(2\pi)^4}{(L_x L_y)^3}\sum F_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
a_{\vec k_1}a_{\vec k_2}^{*}a_{\vec k_3}^{*}\triangle_{\vec k_1,  \vec k_2 + \vec k_3 + \vec k}-\\
\displaystyle
-\frac{i}{2}\frac{(2\pi)^4}{(L_x L_y)^3}\sum D_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
a_{\vec k_1}a_{\vec k_2}a_{\vec k_3}^{*}
\triangle_{\vec k_1 + \vec k_2, \vec k_3 + \vec k}.
\end{array}
\end{equation}
Удобно избавиться от линейного члена подстановкой:
\begin{equation}
a_{\vec k} = A_{\vec k} e^{i \omega_k t}.
\end{equation}
В этих переменных уравнения принимают вид:
\begin{equation}
\label{Eq_A_k_discr}
\begin{array}{l}
\displaystyle
\dot A_{\vec k} = -\frac{i}{2}\frac{2\pi}{L_x L_y}\sum C_{\vec k_1 \vec k_2}^{\vec k} 
A_{\vec k_1}A_{\vec k_2}\triangle_{\vec k_1 + \vec k_2, \vec k}e^{i(\omega_{k_1}+\omega_{k_2} - \omega_k)t}-\\
\displaystyle
-\frac{i2\pi}{L_x L_y}\sum C_{\vec k \vec k_2}^{\vec k_0} 
A_{\vec k_2}^{*}A_{\vec k_0}\triangle_{\vec k + \vec k_2, \vec k_0}e^{i(-\omega_{k} - \omega_{k_2} + \omega_{k_0})t}-\\
\displaystyle
-\frac{i}{2}\frac{2\pi}{L_x L_y}\sum E_{\vec k_1 \vec k_2}^{\vec k} 
A_{\vec k_1}^{*}A_{\vec k_2}^{*}\triangle_{\vec k_1 + \vec k_2 +\vec k, 0}e^{i(-\omega_{k_1} - \omega_{k_2} - \omega_{k})t}-\\
\displaystyle
-\frac{i(2\pi)^4}{(L_x L_y)^3}\sum W_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
A_{\vec k_1}^{*}A_{\vec k_2}^{*}A_{\vec k_3}^{*}
\triangle_{\vec k_1 + \vec k_2+ \vec k_3 + \vec k, 0}e^{i(-\omega_{k_1} - \omega_{k_2} - \omega_{k_3} - \omega_{k_4})t}-\\
\displaystyle
-\frac{i}{4}\frac{(2\pi)^4}{(L_x L_y)^3}\sum F_{\vec  k \vec k_2 \vec k_3 \vec k_4}
A_{\vec k_2}A_{\vec k_3}A_{\vec k_4}
\triangle_{\vec k , \vec k_2 + \vec k_3 + \vec k_4}e^{i(-\omega_{k} + \omega_{k_2} + \omega_{k_3} + \omega_{k_4})t}-\\
\displaystyle
-\frac{i}{4}\frac{3(2\pi)^4}{(L_x L_y)^3}\sum F_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
A_{\vec k_1}A_{\vec k_2}^{*}A_{\vec k_3}^{*}\triangle_{\vec k_1,  \vec k_2 + \vec k_3 + \vec k}e^{i(\omega_{k_1}-\omega_{k_2} - \omega_{k_3} - \omega_{k_4})t}-\\
\displaystyle
-\frac{i}{2}\frac{(2\pi)^4}{(L_x L_y)^3}\sum D_{\vec  k_1 \vec k_2 \vec k_3 \vec k}
A_{\vec k_1}A_{\vec k_2}A_{\vec k_3}^{*}\triangle_{\vec k_1 + \vec k_2, \vec k_3 + \vec k}e^{i(\omega_{k_1}+\omega_{k_2} - \omega_{k_3} - \omega_{k_4})t}.
\end{array}
\end{equation}
Если проинтегрировать левую и правую части уравнения по времени за некоторый фиксированный
промежуток, то становится очевидно, что наибольший вклад в амплитуду волны дадут те члены
уравнения, которые будут иметь наименьший показатель в экспоненте. В непрерывной среде
мы можем выбрать волны так, чтобы занулить эти показатели точно. Однако на дискретной
сетке волновых векторов сделать это в случае произвольного закона дисперсии, вообще
говоря, нельзя. Заметим, что члены, отвечающие за генерацию
кратных гармоник, в большинстве случаев законов дисперсии дают показатели, существенно
отличающиеся от нуля. Поэтому в дальнейшем мы их рассматривать не будем.

Также заметим, что полученные уравнения удобно решать в области с линейными размерами
$L_x = L_y = 2\pi$. В этом случае все формулы существенно упрощаются, а проекции волнового
вектора $k_x$ и $k_y$, соответственно (\ref{krdiscr}), принимают  только целочисленные значения.
Необходимо заметить, что
выбор линейных размеров нисколько не влияет на универсальность полученных результатов, так
как он приведит лишь к перенормировке времени, при фиксированных остальных константах.

 
\subsection{Капиллярные волны}
\label{CapillarRes}

В случае капиллярных волн на поверхности бесконечно глубокой жидкости закон дисперсии имеет
вид:
\begin{equation}
\omega_k = \sqrt{\sigma k^3},
\end{equation}
здесь $\sigma$ --- коэффициент поверхностного натяжения жидкости.

Рассмотрим трехволновой процесс, отвечающий распаду и слиянию волн. Условия резонанса,
следующие из (\ref{Eq_A_k_discr}), выражаются так:
\begin{equation}
\label{resonant_conditions_cap}
\omega_{k_1} + \omega_{k_2} = \omega_{k_0},\;\;\;
\vec k_1 + \vec k_2 = \vec k_0.
\end{equation}
В данном случае, для определенности, будем говорить о распаде волны $A_{\vec k_0}$ в две волны $A_{\vec k_1}$ и
$A_{\vec k_2}$.
Резонансное многообразие (\ref{resonant_conditions_cap}) для распада:
\begin{equation}
\begin{array}{l}
\displaystyle
\vec k_0 = \left(
\begin{array}{c}
0\\
k_0
\end{array}
\right),\\
\displaystyle
\vec k_1 = \left(
\begin{array}{c}
-k_x\\
k_0 - k_y
\end{array}
\right),\;
\vec k_2 = \left(
\begin{array}{c}
k_x\\
k_0 + k_y
\end{array}
\right).
\end{array}
\end{equation}
приведено на Рисунке~\ref{resonance_curve}.
\begin{figure}[!hbt]
%  \includegraphics[angle=-90,width=\PW]{kadrs_of_energy.eps}
\includegraphics[angle=0,width=\PW]{Fig_4.2/resonance.eps}
\caption{ Резонансное множество для $k_0 = 68$.}
\label{resonance_curve}
\end{figure}

Уравнение (\ref{Eq_A_k_discr}) сводится к системе обыкновенных дифференциальных уравнений:
\begin{equation}
\label{two-waves_decay}
\begin{array}{l}
\displaystyle
\dot A_{\vec k_0} = -\frac{i}{2}\frac{2\pi}{L_x L_y} C_{\vec k_1 \vec k_2}^{\vec k_0} 
A_{\vec k_1}A_{\vec k_2}e^{i\Omega_{k_1 k_2}^{k_0} t},\\
\displaystyle
\dot A_{\vec k_1} = -i\frac{2\pi}{L_x L_y} C_{\vec k_1 \vec k_2}^{\vec k_0} 
A_{\vec k_2}^{*}A_{\vec k_0}e^{-i\Omega_{k_1 k_2}^{k_0} t},\\
\displaystyle
\dot A_{\vec k_2} = -i\frac{2\pi}{L_x L_y} C_{\vec k_1 \vec k_2}^{\vec k_0} 
A_{\vec k_1}^{*}A_{\vec k_0}e^{-i\Omega_{k_1 k_2}^{k_0} t},
\end{array}
\end{equation}
Здесь введена расстройка $\Omega_{k_1 k_2}^{k_0} = \omega_{k_1} + \omega_{k_2} - \omega_{k_0}$.
Пусть $A_{\vec k_1}, A_{\vec k_2}$ малы по сравнению с распадающейся волной
$\left|A_{\vec k_0}\right| \gg \max (|A_{\vec k_1}|, |A_{\vec k_2}|)$ в начальный момент времени $t=0$.
Тогда на начальном этапе роста гармоник уравнения (\ref{two-waves_decay}) можно линеаризовать.
Предполагая ($A_{\vec k_0} \simeq {\mathrm const}$), имеем растущее (при выполнении некоторых условий)
решение:
\begin{equation}
\label{exp-growth}
A_{\vec k_{1,2}}(t) = A_{\vec k_{1,2}}(0) e^{\lambda t},
\end{equation}
здесь:
\begin{equation}
\label{lambda}
\lambda =-\frac{i}{2}\Omega_{k_1 k_2}^{k_0} + \sqrt{\left|\frac{2\pi}{L_x L_y} C_{\vec k_1 \vec k_2}^{\vec k_0} A_{\vec k_0}\right|^2 - \left(\frac{1}{2}\Omega_{k_1 k_2}^{k_0}\right)^2}.
\end{equation}
Отсюда видно, что вблизи от резонансной кривой, при выполнении условий:
\begin{equation}
\label{condition}
\left|\frac{2\pi}{L_x L_y} C_{\vec k_1 \vec k_2}^{\vec k_0} A_{\vec k_0}\right| > \left|\frac{1}{2}\Omega_{k_1 k_2}^{k_0}\right|,
\end{equation}
мы получаем рост гармоник. Выбрав достаточно большую начальную амплитуду мы
можем уширить резонансную кривую так, чтобы в получившуюся область попало 
настолько большое число гармоник, что дискретность сетки волновых векторов перестает
быть важной.

Для описания взаимодействия волн на стадии, когда амплитуды растущих
волн уже нельзя считать малыми по сравнению с амплитудой распадающихся волн,
надо переходить к решению полных динамических уравнений (\ref{eta_psi_system})
или (\ref{Eq_a_k}). Численная схема подробно описана в Приложении \ref{NumScheme}.

Система (\ref{eta_psi_system}) решалась в области с линейными размерами
$L_x = L_y = 2\pi$. Коэффициент поверхностного натяжения $\sigma=1$. Число
точек $512\times512$. В качестве начальных условий была взята
монохроматическая волна с волновым вектором $\vec k_0 = (0;68)$ и амплитудой
$|a_{\vec k_0}|=4 \times 10^{-3}$. Все остальные гармоники имели малую амплитуду
$|a_{\vec k}| \sim 10^{-12}$ и случайную фазу.
Как было отмечено выше, резонансная кривая
практически никогда не проходит точно через точки сетки. При выбранных
линейных размерах области все волновые числа являются целыми. Детальная
картина резонансной кривой на сетке волновых векторов в области, где рядом
с резонансной кривой точки лежат наиболее плотно, показана на
Рисунке~\ref{resonance_curve_local}. Видно, что некоторые точки лежат ближе
к резонансной кривой, чем другие.
\begin{figure}[!ht]
\includegraphics[angle=0,width=\PW]{Fig_4.2/cap_res_local.eps}
\caption{ Фрагмент резонансной кривой. Хорошо видна разная расстройка для разных узлов сетки.}
\label{resonance_curve_local}
\end{figure}
\begin{figure}[!hbt]
\includegraphics[angle=0,width=\PW]{Fig_4.2/flat_grow.eps}
\caption{Эволюция различных гармоник при распадающейся волне $\vec k_0 = (00, 68)$.}
\label{Growing}
\end{figure}
%\clearpage

Как и ожидалось, в начале мы наблюдали экспоненциальный рост ряда гармоник
согласно (\ref{exp-growth}) и (\ref{lambda}).
Эта стадия показана на Рис.\ref{Growing} и Рис.\ref{3d-1}. Видно, что одни гармоники находятся в резонансе, а другие нет.
Через некоторое время практически все гармоники, близкие к резонансному множеству, оказываются вовлечены в процесс распада начальной волны, как показано на Рис.\ref{3d-2}.
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_1.eps}
\caption{Резонансные гармоники начинают расти. Момент времени t=1.4.}
\label{3d-1}
\end{figure}
\clearpage

\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_2.eps}
\caption{Начинаются вторичные процессы распада. Момент времени t=11.}
\label{3d-2}
\end{figure}
\begin{figure}
\includegraphics[angle=0,width=\PW]{Fig_4.2/levels_1.eps}
\caption{Линии уровня поверхности $|a_k|^2$. Хорошо видны вторичные распады. Момент времени t=14.}
\label{levels-1}
\end{figure}
\clearpage
Позже гармоники, являющиеся наиболее близкими к резонансной кривой
(сравните с Рис. \ref{resonance_curve_local}) и находящиеся в той ее части, где
матричный элемент больше, достигают уровня, когда они сами могут служить
источником так называемых <<вторичных распадов>>. 

Амплитуды волн становятся существенно различными. Наибольших амплитуд
достигают волны, для которых инкремент (\ref{lambda}) максимален. На Рис.
\ref{levels-1} видна регулярная структура, порожденная начальной
волной $a_{\vec k_0}$.

Некоторое время спустя вся $k$-плоскость заполнена распадающимися волнами,
как показано на Рис.\ref{3d-3}.
\begin{figure}[hbt]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_3.eps}
\caption{Поверхность $|a_k|^2$ в момент времени t=57.}
\label{3d-3}
\end{figure}


\subsection{Гравитационные волны}
\label{GravRes}
 
В случае гравитационных волн на поверхности бесконечно глубокой жидкости
закон дисперсии принимает вид:
\begin{equation}
\omega_k = \sqrt{g k},
\end{equation}
здесь $g$ --- ускорение свободного падения. Здесь и далее для удобства
положим $g=1$.
В этом случае условия (\ref{resonant_conditions_cap} не могут быть
выполнены, так как данная система не имеет действительных нетривиальных
решений. Следовательно, можно подобрать преобразование, которое исключит
члены третьего порядка, отвечающие этому распаду. Поэтому, если в случае
капиллярных волн нам можно было ограничиться
членами третьего порядка в Гамильтониане, то в случае гравитационных
этого не достаточно. И это именно та причина, по которой мы всюду
рассматривали Гамильтониан с точностью до членов четвертого порядка.

Член, отвечающий генерации третьей гармоники, снова не является резонансным.
Если мы рассмотрим начальные условия, аналогичные случаю капиллярных волн,
то есть одна монохроматическая волна и случайный шум малой амплитуды, то
очевидно, что главными будут процессы, в которые амплитуда начальной волны
входит наибольшее количество раз. Процессы распада одной волны в три и
обратный ему, при условии что начальная волна встречается в уравнениях
два раза, не являются резонансными. Поэтому главным будет процесс, который
подразумевает рассеяние двух волн с одинаковыми волновыми векторами в две
другие волны. Условия резонанса такого процесса, следующие из
(\ref{Eq_A_k_discr}), выражаются так:
\begin{equation}
\label{resonant_conditions_grav}
\omega_{k_1} + \omega_{k_2} = 2\omega_{k_0},\;\;\;
\vec k_1 + \vec k_2 = 2\vec k_0.
\end{equation}
Резонансная кривая для условий (\ref{resonant_conditions_grav}) приведена на
Рис.\ref{res_curv_grav}.
\begin{figure}[!hbtp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/res_curv_grav.eps}
\caption{ Резонансная кривая для $\vec k_0 = (0;30)$.}
\label{res_curv_grav}
\end{figure}
Эта кривая иногда называется кривой Филлипса, по имени ее первооткрывателя. 
Система (\ref{eta_psi_system}) решалась в области с линейными размерами
$L_x = L_y = 2\pi$. Ускорение свободного падения $g=1$. Число
точек $512\times512$. В качестве начальных условий была взята
монохроматическая волна с волновым вектором $\vec k_0 = (0;30)$ и амплитудой
$|a_{\vec k_0}|=2.5 \times 10^{-4}$. Все остальные гармоники имели малую амплитуду
$|a_{\vec k}| \sim 10^{-9}$ и случайную фазу.
В начале процесса ситуация развивается схожим со случаем капиллярных волн образом. Наблюдается
экспоненциальный рост нескольких гармоник, находящихся вблизи от резонансной кривой (Рис.~\ref{grav_res_local}). Данная
ситуация показана на Рис.~\ref{grav_grow}.
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/grav_res_local.eps}
\caption{Фрагмент резонансной кривой. Хорошо заметна разная расстройка для разных узлов сетки.}
\label{grav_res_local}
\end{figure}
\begin{figure}
\includegraphics[angle=0,width=\PW]{Fig_4.2/grav_grow.eps}
\caption{Рост гармоник как функция времени. Видно, что есть резонансные и не резонансные гармоники.}
\label{grav_grow}
\end{figure}
Матричный элемент процесса имеет наибольшие значения в окрестности начальных волн. Этим
обуславливается картина начального роста гармоник: рост наблюдается только около <<перекрестья>>
резонансного множества (Рисунок~\ref{3d_grav_100}).
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_grav_100.eps}
\caption{ Начало роста гармоник при $\vec k_0 = (0;30)$.}
\label{3d_grav_100}
\end{figure}
\begin{figure}
\includegraphics[angle=0,width=\PW]{Fig_4.2/levels_grav_100.eps}
\caption{ Начало роста гармоник. Линии уровня $|a_{\vec k}|^2 = 1.5\times10^{-18}$.}
\label{levels_grav_100}
\end{figure}
Наиболее удобным для анализа картины роста гармоник
является представление поверхности $|a_{\vec k}|^2$ в виде линий уровня. Таким образом,
Рисунок~\ref{3d_grav_100} может быть представлен в ином, более понятном виде (Рисунок~\ref{levels_grav_100}).

Дальнейшая эволюция сводится к росту гармоник (Рис.~\ref{3d_grav_400} и Рисунке~\ref{levels_grav_400})
и началу вторичных распадов, выражающемся
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_grav_400.eps}
\caption{ Продолжение роста гармоник.}
\label{3d_grav_400}
\end{figure}
\begin{figure}
\includegraphics[angle=0,width=\PW]{Fig_4.2/levels_grav_400.eps}
\caption{ Продолжение роста гармоник. Линии уровня $|a_{\vec k}|^2 = 1.5\times10^{-18}$.}
\label{levels_grav_400}
\end{figure}
в <<заплывании>> начальных гармоник окружающими (Рисунок~\ref{3d_grav_830} и Рис.~\ref{levels_grav_830}).
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_grav_830.eps}
\caption{ Начальная волна окружена порожденными
вторичными распадами волн гармониками.}
\label{3d_grav_830}
\end{figure}
\begin{figure}
\includegraphics[angle=0,width=\PW]{Fig_4.2/levels_grav_830.eps}
\caption{ Начальная волна окружена порожденными
вторичными распадами гармониками. Линии уровня 
$|a_{\vec k}|^2 = 1.5\times10^{-18}$.}
\label{levels_grav_830}
\end{figure}

Несколько менее интересным, но не менее поучительным, является случай, когда мы изначально
имее две большие волны $a_{\vec k_0}$ и $a_{-\vec k_0}$. В этом случае резонансной кривой является
окружность с центром в нуле и радиусом $|\vec k_0|$. Причем, нетрудно заметить, что такая же картина
будет наблюдаться и в случае капиллярных волн. Результаты моделирования такого процесса
представленны на Рисунке~\ref{3d_circle_1}. Начальные условия были такими:
$$
\vec k_0 = \left(
\begin{array}{c}
0\\
30
\end{array}
\right); 
|a_{\vec k_0}| = |a_{-\vec k_0}| = 10^{-3};
\vec k \ne \vec k_0, |a_{\vec k}| = 10^{-9}.
$$
\begin{figure}[!hbt]
\includegraphics[angle=0,width=\PW]{Fig_4.2/3d_grav_circle.eps}
\caption{ Рост гармоник в случае начальных условий в виде двух волн с противоположными волновыми векторами.}
\label{3d_circle_1}
\end{figure}

Очевидно, что при моделировании процессов распада и рассеяния скорость роста гармоник напрямую
связана с амплитудой начальной гармоники (по крайней мере на начальной стадии процесса).
Казалось бы, логично взять эту амплитуду как можно больше. Однако в этом предположении
кроется серьезная ошибка. Взяв слишком большую амплитуду начальной волны, мы можем <<открыть>>
нерезонансный процесс. В самом деле, уравнение
$\Omega_{k_1 k_2}^{k_0} = \omega_{k_1} + \omega_{k_2} - \omega_{k_0}$ в случае гравитационных волн
не имеет нетривиальных решений при $\Omega_{k_1 k_2}^{k_0} = 0$, но если положить
$\Omega_{k_1 k_2}^{k_0} > 0$, то решения появляются. Остается только достаточно большой
амплитудой $A_{\vec k_0}$ сделать подкоренное выражение в (\ref{lambda}) положительным, и мы
получим <<квази-резонансное>> трехволновое взаимодействие. Надо понимать, что такая возможность
появляется только на дискретной сетке и в случае непрерывной среды такие процессы проигрывают
действительно резонансным. Таким образом, если выбрать слишком большую амплитуду, то в
моделируемый процесс будет внесен артефакт, искажающий истинную физическую картину.
На Рис.~\ref{Wrong_resonance} приведена <<квази-резонансная>> кривая для трехволнового
взаимодействия при $\Omega_{k_1 k_2}^{k_0} = 2.3$ и волновых векторах.
\begin{equation}
\begin{array}{l}
\displaystyle
\vec k_0 = \left(
\begin{array}{c}
0\\
30
\end{array}
\right),\\
\displaystyle
\vec k_1 = \left(
\begin{array}{c}
-k_x\\
30 - k_y
\end{array}
\right),\;
\vec k_2 = \left(
\begin{array}{c}
k_x\\
30 + k_y
\end{array}
\right).
\end{array}
\end{equation}
\begin{figure}[!hbt]
\includegraphics[angle=0,width=\PW]{Fig_4.2/Wrong_resonance.eps}
\caption{ <<Квази-резонансная>> кривая для трехволнового
взаимодействия гравитационных волн. $\vec k_0 = (0;30)$.}
\label{Wrong_resonance}
\end{figure}


\section{Колмогоровские спектры в турбулентности гравитационных волн}
\label{GravTurbu}

Теория слабой турбулентности многократно проверялась в численных экспериментах 
\cite{OptTurb, Onorato-2002, Dias-2001, Annenkov-2001, Kolm}.
Однако единственным успешным экспериментом по численному моделированию турбулентности
волн на поверхности трехмерной жидкости следует признать \cite{PZ1996}, где моделировалась
турбулентность капиллярных волн. Несмотря на значительные затраченные усилия, в случае
гравитационных волн на поверхности добиться прогресса не удавалось. По нашему мнению,
причина этого заключалась
в трудности правильного выбора параметров численной схемы.

Эта глава построена следующим образом: в разделе \ref{StatDescr} вводится статистическое описание
волнового поля на поверхности жидкости для случая гравитационных волн на глубокой воде;
в разделе \ref{NumSim} приведены параметры численной схемы и результаты расчетов;
в заключении обсуждаются
полученные результаты.


\subsection{Введение. Решения кинетического уравнения.}
\label{StatDescr}


Как было отмечено в параграфе \ref{GravRes}, дисперсионный закон для
гравитационных волн на поверхности глубокой жидкости относится к
так называемому <<нераспадному>> случаю. Уравнения:
\begin{equation}
\omega_{k_1} = \omega_{k_2} + \omega_{k_3},\;\;\; \vec k_1 = \vec k_2 + \vec k_3
\end{equation}
не имеют нетривиальных вещественных решений. Это означает, что кубические
члены в Гамильтониане могут быть исключены соответствующим каноническим
преобразованием $a(\vec k,t) \longrightarrow b(\vec k,t)$ \cite{ZLF92}. Формула
этого преобразования довольно громоздка и хорошо известна
\cite{ZLF92,Zakharov-1999}, поэтому не буудем здесь вдаваться в детали.

Для статистического описания волнового поля  используется
корреляционная функция:
\begin{equation}
<a_{\vec k} a_{\vec k'}^*> = n_{k} \delta (\vec k - \vec k').
\end{equation}
Величина $n_{\vec k}$ непосредственно связана с наблюдаемыми
корреляционными функциями
и может быть получена из экспериментальных данных. Например, из
(\ref{a_k_substitution}) можно получить:
\begin{equation}
\label{I_k_expression}
I_k = <|\eta_{\vec k}|^2> = \frac{1}{2}\frac{\omega_k}{g} (n_k + n_{-k}).
\end{equation}
В случае гравитационных волн удобнее пользоваться другой корреляционной
функцией:
\begin{equation}
<b_{\vec k} b_{\vec k'}^*> = N_{k} \delta (\vec k - \vec k').
\end{equation}
Функция $N_k$ не может быть измерена напрямую. В общем случае связь между
$n_k$ и $N_k$ довольно сложна \cite{Zakharov-1999}, но в случае глубокой воды
все становится проще:
\begin{equation}
\frac{n_k - N_k}{n_k} \simeq \mu,
\end{equation}
здесь $\mu = (ka)^2$, где $a$ --- характерное отклонение поверхности.
В случае слабой турбулентности $\mu << 1$.
Корреляционная функция $N_k$ подчиняется кинетическому уравнению
\cite{Zakharov-1966}:
\begin{equation}
\label{Kinetic_equation}
\frac{\partial N_k}{\partial t} = st(N,N,N) + f_p (k) - f_d (k),
\end{equation}
Здесь:
\begin{equation}
\begin{array}{l}
\displaystyle
st(N,N,N)=4\pi \int \left| T_{\vec k,\vec k_1,\vec k_2,\vec k_3}\right|^2 \times\\
\displaystyle
\times(N_{k_1} N_{k_2} N_{k_3} +N_{k} N_{k_2} N_{k_3} - N_{k} N_{k_1} N_{k_2} -\\
\displaystyle
- N_{k} N_{k_1} N_{k_3})\delta (\vec k + \vec k_1- \vec k_2 - \vec k_3)d \vec k_1d \vec k_2 d \vec k_3.\\
\end{array}
\end{equation}
Подробная форма матричного элемента $T_{\vec k,\vec k_1,\vec k_2,\vec k_3}$
может быть найдена во множестве источников \cite{Zakharov-1966,
Zakharov-1967,Zakharov-1999}.
Функция $f_p (k)$ в (\ref{Kinetic_equation}) соответствует накачке волн, например,
ветром. Обычно накачка сосредоточена на больших масштабах.
Функция $f_d (k)$ представляет собой затухание волн благодаря вязкости
и обрушению волн. Эти функции известны лишь приблизительно. Однако
существует интервал между масштабами, в котором действие накачки и затухания
можно считать слабым и в котором и наблюдают Колмогоровские спектры в
эксперименте.

Рассмотрим стационарное решение уравнения (\ref{Kinetic_equation}), предполагая, что выполнены следующие утверждения:
\begin{itemize}
\item Среда изотропна по отношению к поворотам;
\item Дисперсионное соотношение является степенной функцией
$\omega=a k^\alpha$;
\item $T_{\vec k,\vec k_1,\vec k_2,\vec k_3}$ --- однородная функция:
$T_{\epsilon \vec k,\epsilon \vec k_1,\epsilon \vec k_2,\epsilon \vec k_3} = \epsilon^\beta T_{\vec k,\vec k_1,\vec k_2,\vec k_3}$.
\end{itemize}
В этих предположениях можно найти Колмогоровские решения \cite{ZLF92}:
\begin{equation}
\begin{array}{l}
\displaystyle
n_k^{(1)} = C_1 P^{1/3} k^{-\frac{2\beta}{3} - d},\\
\displaystyle
n_k^{(2)} = C_2 Q^{1/3} k^{-\frac{2\beta - \alpha}{3} - d}.\\
\end{array}
\end{equation}
Здесь $d$ --- это число измерений пространства ($d=2$ в случае гравитационных поверхностных волн).
Первое выражение представляет собой Колмогоровский спектр,
соответствующий постоянному потоку энергии $P$ в область малых масштабов
(прямой каскад энергии). Второе выражение --- это Колмогоровский спектр,
описывающий каскад волнового действия на большие масштабы, и $Q$ --- 
это поток действия. В обоих случаях $C_1$ и $C_2$ --- это безразмерные
<<Колмогоровские константы>>.

В случае гравитационных волн на глубокой воде $\omega=\sqrt{gk}$, а $\beta=3$.
Из работы \cite{Zakharov-1966} стало известно, что:
\begin{equation}
\label{direct_spectrum}
n_k^{(1)} = C_1 P^{1/3} k^{-4}.
\end{equation}
Действуя таким же образом \cite{Zakharov-1982}, получаем второй спектр:
\begin{equation}
n_k^{(2)} = C_2 Q^{1/3} k^{-23/6}.
\end{equation}

Эта глава посвящена исследованию спектра, отвечающего потоку энергии
(\ref{direct_spectrum}). Используя (\ref{I_k_expression}), можно вычислить:
\begin{equation}
\label{I_k_spectrum}
I_k = \frac{C_1 g^{1/2} P^{1/3}}{k^{7/2}}.
\end{equation}

\subsection{Численное моделирование}
\label{NumSim}
Для моделирования турбулентности поверхностных волн в уравнения (\ref{eta_psi_system})
необходимо ввести члены, отвечающие накачке и затуханию. Например, это можно сделать таким образом:
\begin{equation}
\begin{array}{l}
\displaystyle
\dot \eta_{\vec k} = ... + \gamma_k \eta_{\vec k},\\
\displaystyle
\dot \psi_{\vec k} = ... + \gamma_k \psi_{\vec k} + f_k.
\end{array}
\end{equation}
Такая модификация дает экспоненциальное затухание гармоник и линейную по времени накачку.
Модифицированные динамические уравнения (\ref{eta_psi_system}) решались
с использованием численной схемы, описанной в главе \ref{NumScheme}. Конкретный вид функций,
определяющих накачку и затухание, был таков:
\begin{equation}
\begin{array}{l}
\displaystyle
F_k = f_k e^{iR_{\vec k} (t)},\\
\displaystyle
f_k = 4 F_0 \frac{(k-k_{p1})(k_{p2}-k)}{(k_{p2} - k_{p1})^2};\\
\displaystyle
D_{\vec k} = \gamma_k \psi_{\vec k},\\
\displaystyle
\gamma_{k} = -\gamma_1, k \le k_{p1},\\
\displaystyle
\gamma_{k} = - \gamma_2 (k - k_d)^2, 
k > k_d.
\end{array}
\end{equation}
Здесь $R_{\vec k} (t)$ --- это равномерно распределенное случайное число в интервале $(0,2\pi)$.
Вычисления производились в области $L_x = L_y = 2\pi$ с периодическими граничными условиями. 
Параметры накачки были такими: $k_{p1} = 5,\; k_{p2} = 10,\; k_d = 64$.

Как было показано в параграфе \ref{GravRes}, при решении надо особо обратить внимание на
уровень гармоник, так как от этого зависит степень достоверности моделирования резонансных
взаимодействий волн. Экспериментально были найдены несколько подходящих уровней накачки.
Наилучшие результаты при приемлемом времени рассчетов получены для значений $F_0 = 2\times10^{-4}$
и $F_0 = 3\times10^{-4}$.

Значения параметров затухания адаптировались при рассчетах по следующим алгоритмам. Для
сохранения приемлемого качества решения, требовалось, чтобы в спектре наибольшая гармоника
отличалась от наименьшей в $10^{6} - 10^{7}$ раз. Обеспечивалось это посредством изменения
параметра $\gamma_2$. При таких условиях удалось избежать эффекта <<бутылочного горлышка>>,
описанного в \cite{Falkovich-1994}. 

Так как нас интересовал только спектр, соответствующий прямому
Колмогоровскому каскаду (перенос энергии в область с большими $k$), то единственное
требование, предъявлявшееся к спектру при $k \le k_{p1}$, сотояло в том, чтобы данные гармоники
были не больше, чем пиковая гармоника спектра. С этой целью проверялось выполнение данного условия
и при необходимости включалось затухание $\gamma_1$. В остальное время затухания в области малых
волновых чисел $k$ не было.

Получение стационарного состояния (то есть состояния, в котором спектр в целом не меняется)
производилось в несколько этапов. Очевидно, что наибольшее время занимает установление равновесия
в области относительно небольших волновых чисел, то есть в окрестности накачки. В таких условиях
считать на <<большой>> сетке было бы расточительно. Поэтому расчет начинался на сетке $128\times 128$,
при этом $k_d = 32$ до получения стационарного состояния. Затем происходил переход на сетку
$256\times 256$, при этом $k_d = 64$. Недостающие гармоники заполнялись нулями. После установления
равновесия был произведен переход на сетку $512\times 512$, при этом $k_d = 100$. После установления
равновесия и в этом режиме, было произведено усреднение проинтегрированного по углу квадрата спектра
(в силу симметрии системы) для получения коррелятора $<|\eta_{\vec k}|^2>$.

\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.3/inf.Ham.eps}
\label{Hamiltonian_fig}
\caption{Гамильтониан системы как функция времени.}
\end{figure}
На Рисунке~\ref{Hamiltonian_fig} приведен график Гамильтониана как функции времени для таких параметров модели:
$F_0 = 2\times10^{-4}, \gamma_1 = 10^{-3}, \gamma_2 \simeq 512, g=1$.

Корреляционная функция оказалась степенной в значительной части инерционного интервала, который составлял
область от $k_{p2}=10$ до $k_d=100$. Коррелятор показан на Рис. \ref{Correlator_fig}.
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.3/etacorr2.eps}
\label{Correlator_fig}
\caption{Корреляционная функция отклонения поверхности в двойном логарифмическом масштабе.}
\end{figure}
\begin{figure}
\includegraphics[angle=0,width=\PW]{Fig_4.3/comparespec3.eps}
\label{Correlator_mul}
\caption{Скомпенсированный коррелятор отклонения поверхности при различных степенях компенсации: $z=3.5$ сплошная линия
(теория слабой турбулентности), $z=4.0$ прерывистая линия (теория Филлипса).}
\end{figure}

Можно оценить показатель этой степенной функции. Скомпенсированные спектры
(для сравнения мы привели спектр Филлипса) показаны
на Рисунке~\ref{Correlator_mul}. Видно, что спектр, предсказанный теорией слабой турбулентности,
лучше подходит к результатам численного эксперимента.

На Рисунке~\ref{Broadening} видно уширение области, в которой наблюдается степенной спектр
в результате увеличения инерционного интервала при переходе на большую сетку.
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.3/spectra4.eps}
\label{Broadening}
\caption{Уширение области степенного спектра при увеличении количества точек на сетке.}
\end{figure}
Также представляет интерес проверить, существует ли зависимость показателя
в полученном степенном спектре от уровня накачки. Для этого был выполнен
расчет при уровне накачки $F_0 = 3\times 10^{-4}$ на сетке размером $256\times
256$. На Рис.~\ref{diffpump} хорошо
видно, что отличия показателей не наблюдается. Как и следовало ожидать,
наблюдается отличие потока $P$ в формуле (\ref{I_k_spectrum}), что
приводит лишь к параллельному смещению кривой.
\begin{figure}[!htbp]
\includegraphics[angle=0,width=\PW]{Fig_4.3/diffpump.eps}
\label{diffpump}
\caption{Инерционный интервал спектра при различных уровнях накачки.}
\end{figure}

 
 
\subsection{Заключение}
\label{GravTurbConcl}
Произведенное численное моделирование динамических уравнений волн на
поверхности жидкости в случае
гравитационных волн на глубокой воде показало, что в существенной области
инерционного интервала для коррелляционной функции отклонения поверхности наблюдается
Колмогоровский (степенной) спектр, подтверждающий предсказания теории
слабой турбулентности. Произведенные дополнительные численные
эксперименты
продемонстрировали, что степень убывания спектра не зависит от уровня
накачки и затухания в модели, что позволяет считать данный спектр
универсальным.
Таким образом, можно заключить, что предположения, в рамках которых получена теория слабой
турбулентности, являются верными и в случае гравитационных волн на поверхности глубокой жидкости.

