Математика модели цунами

Уединенная волна (решение для солитона Уравнения Кортевега де Вриза) перемещается на постоянной скорости справа налево вдоль канала постоянной глубины. Это соответствует цунами, перемещающемуся по глубокому морю. В левом конце канала существует наклон, симулирующий континентальный шельф. После достижения наклона уединенная волна начинает увеличивать свою высоту. Когда вода становится очень мелкой, большая часть волны отражается назад в канал. Однако узкий, но высокий пик воды возникает в конце наклона и возобновляет сниженную скорость в первоначальном направлении инцидентной волны. Это - цунами, которое наконец обрушивается на берег, вызывая катастрофическое разрушение вдоль береговой линии. Скорость волны, приближающейся к берегу, сравнительно мала. Волна в конечном счете начинает повреждаться.

Используя линейную водную теорию без рассеяния, высоту $$u(x,t)$$ из волны свободной поверхности выше уровня воды без помех в одномерном канале различной глубины $$h(x)$$ решение следующего дифференциального уравнения с частными производными. (См. [2].)

В этой формуле индексы обозначают частные производные, и $$g=9.81m/s^2$$ гравитационное ускорение.

Рассмотрите волну, пересекающую линейный наклон $$h(x)$$ из области с постоянной глубиной $$h_2$$ в область с постоянной глубиной $$h_1\ll h_2$$. Преобразование Фурье относительно $$t$$ поворачивает дифференциальное уравнение с частными производными водной волны к следующему обыкновенному дифференциальному уравнению для режима Фурье $$u(x,t)=U(x,\omega )e^{i\omega t}$$.

Для областей с постоянной глубиной $$h$$, режимы Фурье перемещаются распространение волн в противоположных направлениях с постоянной скоростью $$c=\sqrt{gh}$$.

Решение $$u_2(x,t)=e^{i\omega (t+x/c_2)}+R(\omega )e^{i\omega (t-x/c_2)}$$ поскольку глубоководная область является суперпозицией двух волн:

Этот выбор $$u_2$$ удовлетворяет уравнению волны в глубоководной области для любого $$R(\omega )$$.

Решение $$u_1(x,t)=T(\omega )e^{i\omega (t+x/c_1)}$$ для мелководья область является переданной волной, перемещающейся налево с постоянной скоростью $$c_1=\sqrt{g{h}_1}$$. Этот выбор $$u_1$$ удовлетворяет уравнению волны в мелководной области для любого коэффициента передачи $$T(\omega )$$.

Для области перехода (наклон), использовать $$u(x,t)=U(x,w)e^{i\omega t}$$.

Параметры и решения модели цунами в Symbolic Math Toolbox

Задайте параметры модели цунами можно следующим образом. Игнорируйте зависимость от частоты $$\omega$$ в следующих обозначениях: $$R=R(\omega )$$, $$T=T(\omega )$$, $$U(x)=U(x,\omega )$$.

syms L H depthratio g positive
syms x t w T R U(x)

L1 = depthratio*L; 
L2 = L;

h1 = depthratio*H; 
h2 = H;
h(x) = x*H/L;

c1 = sqrt(g*h1);
c2 = sqrt(g*h2);

u(x,t)  = U(x)*exp(1i*w*t);
u1(x,t) = T*exp(1i*w*(t + x/c1));
u2(x,t) = exp(1i*w*(t + x/c2)) + R*exp(1i*w*(t - x/c2));

В области перехода по линейному наклону используйте dsolve решать ОДУ для преобразования Фурье $$U$$ из $$u$$.

wavePDE(x,t) = diff(u,t,t) - g*diff(h(x)*diff(u,x),x);
slopeODE(x) = wavePDE(x,0); 
U(x) = dsolve(slopeODE);

Решение $$U$$ сложное выражение, включающее Функции Бесселя. Это содержит две произвольных "константы", которые зависят от $$\omega$$.

Const = setdiff(symvar(U), sym([depthratio,g,H,L,x,w]))

Const = $$\left(\begin{array}{cc}{C}_1& C_2\end{array}\right)$$

Для любого режима Фурье полное решение должно быть непрерывно дифференцируемой функцией $$x$$. Следовательно, значения функции и производные должны соответствовать в точках шва $$L_1$$ и $$L_2$$. Это обеспечивает четыре линейных уравнения для $$T$$, $$R$$, и эти две константы в $$U$$.

du1(x) = diff(u1(x,0),x);
du2(x) = diff(u2(x,0),x);
dU(x)  = diff(U(x),x);

eqs =  [ U(L1) == u1(L1,0), U(L2) == u2(L2,0),...
        dU(L1) == du1(L1), dU(L2) == du2(L2)];
unknowns = [Const(1),Const(2),R,T];

Решите эти уравнения.

[Cvalue1, Cvalue2, R, T] = solve(eqs, unknowns);

Замените результатами назад в $$R$$, $$T$$, и $$U$$.

U(x) = subs(U(x), {Const(1),Const(2)}, {Cvalue1,Cvalue2});

Вы не можете непосредственно оценить решение для $$\omega =0$$ потому что и числитель и знаменатель соответствующих выражений исчезают. Вместо этого найдите низкочастотные пределы этих выражений.

simplify(limit(U(x), w, 0)) 

ans = 
$$\frac{2}{\sqrt{\mathrm{depthratio}}+1}$$

simplify(limit(R, w, 0)) 

ans = 
$$-\frac{\sqrt{\mathrm{depthratio}}-1}{\sqrt{\mathrm{depthratio}}+1}$$

simplify(limit(T, w, 0))

ans = 
$$\frac{2}{\sqrt{\mathrm{depthratio}}+1}$$

Эти пределы удивительно просты. Они только зависят от отношения значений глубины, задающих наклон.

Замените символьными параметрами с числовыми значениями

Для следующих расчетов используйте эти численные значения в символьных параметрах.

g = 9.81;
L = 2;
H = 1;
depthratio = 0.04; 

h1 = depthratio*H;
h2 = H;

L1 = depthratio*L;
L2 = L; 

c1 = sqrt(g*h1);
c2 = sqrt(g*h2);

Задайте входящий солитон амплитуды $$A$$ перемещение налево с постоянной скоростью $$c_2$$ в глубоководной области.

A = 0.3;
soliton = @(x,t) A.*sech(sqrt(3/4*g*A/H)*(x/c2+t)).^2;

Выбрать $$Nt$$ точки выборки для $$t$$. Масштаб времени выбран в качестве кратного (временной) ширине входящего солитона. Сохраните соответствующие дискретизированные частоты преобразования Фурье в $$W$$.

Nt =  64;
TimeScale = 40*sqrt(4/3*H/A/g);
W = [0, 1:Nt/2 - 1, -(Nt/2 - 1):-1]'*2*pi/TimeScale;

Выбрать $$Nx$$ точки выборки в $$x$$ направление для каждой области. Создайте точки выборки $$X1$$ для мелководной области, $$X2$$ для глубоководной области, и $$X12$$ для наклонной области.

Nx = 100;
x_min = L1 - 4;
x_max = L2 + 12;
X12 = linspace(L1, L2, Nx);
X1  = linspace(x_min, L1, Nx);
X2  = linspace(L2, x_max, Nx);

Вычислите преобразование Фурье входящего солитона на сетке времени $$Nt$$ равноотстоящие точки выборки.

S = fft(soliton(-0.8*TimeScale*c2, linspace(0,TimeScale,2*(Nt/2)-1)))';
S = repmat(S,1,Nx);

Создайте перемещающееся решение волны в глубоководной области на основе данных Фурье в $$S$$.

soliton = real(ifft(S .* exp(1i*W*X2/c2)));

Преобразуйте режимы Фурье отраженной волны в глубоководной области к численным значениям по сетке в $$(x,\omega )$$ пробел. Умножьте эти значения с коэффициентами Фурье в $$S$$ и используйте функциональный ifft вычислить отраженную волну в $$(x,t)$$ пробел. Обратите внимание на то, что первая строка числовых данных $$R$$ состоит из значений NaN потому что соответствующая численная оценка символьных данных $$R$$ для $$\omega =0$$ не возможно. Задайте значения в первой строке $$R$$ как низкочастотные пределы.

R = double(subs(subs(vpa(subs(R)), w, W), x ,X2));
R(1,:) = double((1-sqrt(depthratio)) / (1+sqrt(depthratio)));
reflectedWave = real(ifft(S .* R .* exp(-1i*W*X2/c2)));

Используйте тот же подход в переданной волне в мелководной области.

T = double(subs(subs(vpa(subs(T)),w,W),x,X1));
T(1,:) = double(2/(1+sqrt(depthratio)));
transmittedWave = real(ifft(S .* T .* exp(1i*W*X1/c1)));

Кроме того, используйте этот подход в наклонной области.

U12 = double(subs(subs(vpa(subs(U(x))),w,W),x,X12));
U12(1,:) = double(2/(1+sqrt(depthratio)));
U12 = real(ifft(S .* U12));

Постройте решение

Для более сглаженной анимации сгенерируйте дополнительные точки выборки с помощью тригонометрической интерполяции вдоль столбцов данных о графике.

soliton = interpft(soliton, 10*Nt);
reflectedWave = interpft(reflectedWave, 10*Nt);
U12 = interpft(U12, 10*Nt);
transmittedWave = interpft(transmittedWave, 10*Nt);

Создайте анимированный график решения что показы в отдельном окне рисунка.

figure('Visible', 'on');
plot([x_min, L1, L2, x_max], [-h1, -h1, -h2, -h2], 'Color', 'black')
axis([x_min, x_max, -H-0.1, 0.6])
hold on

line1 = plot(X1,transmittedWave(1,:), 'Color', 'blue');
line12 = plot(X12,U12(1,:), 'Color', 'blue');
line2 = plot(X2,soliton(1,:) + reflectedWave(1,:), 'Color', 'blue');

for t = 2 : size(soliton, 1)*0.35
  line1.YData = transmittedWave(t,:);
  line12.YData = U12(t,:);
  line2.YData = soliton(t,:) + reflectedWave(t,:);
  pause(0.1)
end

Больше о цунами

В реальной жизни цунами имеют длину волны сотен километров, часто перемещающихся на скоростях больше чем 500 км/час. (Обратите внимание на то, что средняя глубина океана составляет приблизительно 4 км, соответствуя скорости $$\sqrt{gh}\approx 700km/hour$$.) По глубокому морю амплитуда скорее мала, часто приблизительно 0,5 м или меньше. При распространении на полку, однако, цунами увеличивают свою высоту существенно: амплитуды до 30 м и о большем сообщили.

Одно интересное явление - то, что несмотря на то, что цунами обычно приближаются к береговой линии как к расширению фронта волны для сотен перпендикуляра километров к направлению, в котором они перемещаются, они не наносят универсальный ущерб вдоль побережья. В некоторых точках они вызывают бедствия, тогда как только умеренные явления волны наблюдаются в других местах. Это вызывается различными наклонами от морского дна до континентального шельфа. На самом деле очень крутые склоны заставляют большую часть цунами быть отраженной назад в область глубоководных, тогда как маленькие наклоны отражают меньше волны, передавая узкую, но высокую волну, несущую много энергии.

Запустите симуляцию для различных значений $$L$$, которые соответствуют различным наклонам. Чем более крутой наклон, тем ниже и менее мощный волна, которая передается.

Обратите внимание на то, что эта модель игнорирует эффекты трения и дисперсия. На полке симуляция теряет свой физический смысл. Здесь, эффекты трения важны, вызывая повреждение волн.

Ссылки

[1] Дерек Г. Горинг и Ф. Рэйчлен, Цунами - Распространение Длинных волн на Полку, Журнал Водного пути, Порта, Прибрежной и Океанской Разработки 118 (1), 1992, стр 41 - 63.

[2] Х. Лэмб, гидродинамика, Дувр, 1932.