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

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

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

В этой формуле подстрочные индексы обозначают частные производные, а $$g\mathrm{=}\mathrm{}\mathrm{}{}^{\mathrm{}}$$9,81m/s2 - гравитационное ускорение.

Рассмотрим волну, пересекающую линейный наклон $$h(x$$) от области с постоянной глубиной $${}_{\mathrm{}}$$h2 до области с постоянной глубиной $${}_{\mathrm{}}{}_{\mathrm{}}$$h1≪h2. Преобразование Фурье относительно $$$$t превращает дифференциальное уравнение в частных производных водной волны в следующее обыкновенное дифференциальное уравнение для режима Фурье $$u(x,t)=U{}^{(x,}$$λ) eiobjectt.

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

Решение $${}_{\mathrm{u2}}(x,t{)}^{=\mathrm{ei\lambda }{(}_{\mathrm{t}}}+{\mathrm{x/c2}}^{)+R({\lambda }_{\mathrm{)}}}$$eiλ (t-x/c2) для глубоководной области является наложением двух волн:

Этот выбор $${}_{\mathrm{u2}}$$ удовлетворяет волновому уравнению в глубоководной области для любого $$R(\lambda$$).

Решение $${}_{\mathrm{u1}}(x,t)=T^{(\lambda ){\mathrm{ei\lambda }}_{\mathrm{}}(}$$t + x/c1) для мелководной области представляет собой передаваемую волну, проходящую влево с $${}_{\mathrm{постоянной}}\sqrt{{}_{\mathrm{}}}$$скоростью c1 = gh1. Этот $${}_{\mathrm{}}$$выбор u1 удовлетворяет волновому уравнению в области мелководья для любого $$$$коэффициента передачи T (λ).

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

Параметры и решения модели цунами в инструментарии символьной математики

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

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$$ представляет собой сложное выражение, включающее функции Бесселя. Он содержит две произвольные «константы», которые зависят от $$\lambda$$.

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

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

Для любого режима Фурье общее решение должно быть непрерывно дифференцируемой функцией $$$$x. Следовательно, значения функции и производные должны совпадать в точках шва $${}_{\mathrm{L1}}$$ и $${}_{\mathrm{L2}}$$. Это обеспечивает четыре линейных уравнения для $$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});

Невозможно непосредственно вычислить решение для $$\lambda \mathrm{=}$$ 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$$, перемещающийся влево с постоянной скоростью $${}_{\mathrm{c2}}$$ в глубоководной области.

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

Выберите $$\mathrm{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;

Выберите $$\mathrm{Nx}$$ точек образца в $$$$направлении x для каждой области. Создание точек выборки, $$\mathrm{X1}$$ для области мелководья, $$\mathrm{X2}$$ для области глубоководья и $$\mathrm{}\mathrm{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);

Вычислите преобразование Фурье входящего солитона на временной сетке $$\mathrm{из\; 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)));

Преобразование режимов Фурье отраженной волны в глубоководной области в числовые значения по сетке в $$\mathrm{пространстве\; (x},\lambda$$). Умножьте эти значения на коэффициенты Фурье в $$$$S и используйте функциюifft для вычисления отраженной волны в $$(x,t$$) пространстве. Заметим, что первая строка числовых данных $$$$R состоит из значений NaN, потому что правильная численная оценка символьных данных $$$$R для $$\lambda \mathrm{}$$= 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{}\mathrm{}\mathrm{}\mathrm{}\mathrm{gh\approx 700km/hour}$$.) Над глубоководьем амплитуда довольно небольшая, часто около 0,5 м и менее. При распространении на шельф, однако, цунами резко увеличивают свою высоту: сообщалось о амплитудах до 30 м и более.

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

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

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

Ссылки

[1] Дерек Г. Горинг и Ф. Райхлен, Tsunamis - Распространение длинных волн на шельф, Journal of Waterway, Port, Coasteral and Ocean Engineering 118 (1), 1992, pp. 41 - 63.

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