В 1952 году известный океанолог Лонге-Хиггинс ( Michael Longuet-Higgins) успешно доказал предположение, что волны в море отлично описываются суммой плоских волн бегущих в различных направлениях. Этот принцип очень хорошо демонстрируется картинкой расположенной ниже.

fig1

Как мы знаем волна в двумерном случае представляется в виде синуса

z(x) = sin(x).
Рассмотрим подробнее волну в трехмерном случае. В случае, когда наша волна зависит только от X, получаем следующую картинку:
sin_X
Как видите мы можем составить морское волнение из суммы волн, которые движутся либо в направлении X либо в направлении Y. Как же запустить волну так чтобы, она распространялась в произвольном направлении?
Для этого нам надо спроецировать каждую точку плоскости на нужное направление. А чтобы это сделать просто нам нужно вспомнить линейную алгебру и операцию скалярного произведения. Скалярное произведение — это произведение двух векторов, а результатом является не вектор, а лишь одно число.
Мы будем рассматривать два вектора:

  1. Направление распространения волны dir = (a,b),
  2. Позиция точки (X_i, Y_i).

Результатом скалярного произведения, назовем величину S.
S(x,y) = dir*(x,y);
Таким образом, мы спроецируем каждую точку плоскости на направление распространения волны и получим следующую формулу:
z(x,y) = sin(S(x,y));
Кроме того, наша волна будет иметь пространственную частоту, которую назовем k, эту переменную также называют «волновое число». Данная частота зависит от длины волны следующим образом:
k=2pi/wavelength;
Тогда формула преобразуется в следующий вид:
z(x,y) = sin(S(x,y)*k);
Попробуем это на деле, запрограммировав эту формулу в Matlab
dir = [0.5,0.5];
wavelength = 5;
x = [0:0.1:5*pi];
y = [0:0.1:5*pi];
z = zeros(length(x), length(y)); % выделяем память под переменную F
for ii=1:length(x)
for jj=1:length(y)
z(ii,jj)=sin(dir’*[x(ii), y(jj)]*2*pi/wavelength);% Matlab имеет очень мощные средства работы с векторами и матрицами поэтому, когда мы вектор (dir’) умножаем на вектор([x(ii), y(jj)] ), то умножение происходит по правилам скалярного произведения.
end
end
Surf(x,y,z); % отрисовка результата

Результат выполнения этого кода:
sin_F_1
Как видим действительно наша волна располагается в направлении вектора (0.5,0.5).
Далее нам необходимо заставить волну двигаться, для этого основное уравнение необходимо дополнить следующими параметрами
z(x,y,t) = sin(S(x,y)*k+t*w);
Как видим уравнение дополнилось членом t*w, где w – временная частота, которая связана с периодом волны как w = 2*pi/T, здесь T – период в секундах. Вообще говоря, период и длина волны связаны как связаны частота (w) и волновое число (k). Отношение w/k называют «фазовой скоростью» волны. В зависимости от различных условий фазовая скорость волны будет разная, однако установлено, что существует зависимость между w и k, такая зависимость называется дисперсионным соотношением. Например, для волн на конечной глубине зависимость будет следующая:
w(k) = sqrt(g*k*tanh(k*h));
Выглядит немного страшно, но зато это соотношение точное и не имеет никаких приближений.
Выразим это в виде кода на Matlab. Поскольку теперь наша волна может распространяться во времени, необходимо представить результат моделирования не одним графиком, а видеофайлом. Наш код существенно изменился, поэтому обратите внимание на комментарии в коде.
dir = [0.5, 0.5]; %направление распространения волны
wavelength = 5; %длина волны в метрах
Gg = 9.81; %гравитационная постоянная
H = 15; %глубина в метрах
x = [0:0.1:5*pi]; %определяем границы изменения пространственной переменной X
y = [0:0.1:5*pi]; %определяем границы изменения пространственной переменной Y
z = zeros(length(x), length(y)); % выделяем память под переменную F
fig=figure;% переменная в которой будет храниться окно графика
aviobj = VideoWriter('wave.avi'); %создаем объект связанный с видеофайлом
open(aviobj); %открываем объект для записи видеофреймов
k = 2*pi/wavelength; % расчет волнового числа
w = sqrt(k*Gg*tanh(k*H)); % расчет частоты по дисперсионному соотношению
for t = 0:0.1:3 %цикл по времени
for ii=1:length(x) %цикл по пространственной переменной Х
for jj=1:length(y) %цикл по пространственной переменной Y
z(ii,jj)=sin(dir*[x(ii), y(jj)]'*k+t*w); %расчет высоты в каждой точке домена
end
end
view(180, 80) %определение позиции с которой будет показано движение волны
surf(x,y,z,'LineStyle','none'); % отрисовка результата, параметр LineStyle указан для того чтобы сетка на поверхности не отрисовывалась
frame = getframe(fig); % получение фрейма из графика
writeVideo(aviobj,frame); % запись фрейма в видеофайл
end %завершение цикла по времени
close(fig); %закрытие графика
close(aviobj); %завершение потока, ведущего запись в файл

Видеофайл с результатом моделирования можно скачать здесь. Либо можно посмотреть по этой ссылке на youtube.
Сейчас волновое поле представляет из себя одиночную волну распространяющуюся в любом заданном направлении. Для того чтобы представить реалистичную волновую поверхность необходимо представить сумму плоских волн, как и говорилось ранее, чем больше таких волн будет тем реалистичнее будет поверхность. Модифицируем наш код добавив расчет одновременно 10 волн а также некоторую случайность в направлении их распространения и высоты.


dir = rand(10,2); %направление распространения волны
wavelength = 2*rand(10,1)+3; %длина волны в метрах
A = rand(10,1)+2;
Gg = 9.81; %гравитационная постоянная
H = 15; %глубина в метрах
x = [0:0.1:10*pi]; %определяем границы изменения пространственной переменной X в метрах
y = [0:0.1:10*pi]; %определяем границы изменения пространственной переменной Y в метрах
z = zeros(length(x), length(y)); % выделяем память под переменную F
fig=figure;% переменная в которой будет храниться окно графика
aviobj = VideoWriter('wave2.avi'); %создаем объект связанный с видеофайлом
open(aviobj); %открываем объект для записи видеофреймов
k = 2*pi./wavelength; % расчет волнового числа
w = sqrt(k*Gg.*tanh(k*H)); % расчет частоты по дисперсионному соотношению
for t = 0:0.1:10 %цикл по времени
for ii=1:length(x) %цикл по пространственной переменной Х
for jj=1:length(y) %цикл по пространственной переменной Y
for m=1:10
z(ii,jj)=z(ii,jj)+A(m)*sin(dir(m,:)*[x(ii), y(jj)]'*k(m)+t*w(m)); %расчет высоты в каждой точке домена
end
z(ii,jj) = z(ii,jj)/10;
end
end
view(180, 80) %определение позиции с которой будет показано движение волны
surf(x,y,z,'LineStyle','none'); % отрисовка результата, параметр LineStyle указан для того чтобы сетка на поверхности не отрисовывалась
set(gca,'ZLim',[-1*H max(A)]);
frame = getframe(fig); % получение фрейма из графика
writeVideo(aviobj,frame); % запись фрейма в видеофайл
end %завершение цикла по времени
close(fig); %закрытие графика
close(aviobj); %завершение потока ведущего запись в файл

Результат этого моделирования можно скачать тут. Также можно посмотреть он-лайн на youtube
Безусловно такое моделирование, несмотря на страшную формулу дисперсионного соотношения, не совсем физично, однако в следующих статьях мы постараемся сделать существенно более физичным и реалистичным.

Спасибо за внимание!