Этот пример моделирует явление волны цунами с помощью Symbolic Math Toolbox™ для решения дифференциальных уравнений.
Эта симуляция является упрощенной визуализацией явления и основана на статье Горинга и Райхлена [1].
Одиночная волна (солитонное решение уравнения Кортевега-де Фриса) движется с постоянной скоростью справа налево по каналу постоянной глубины. Это соответствует цунами, перемещающемуся по глубокому морю. В левом конце канала расположен склон, имитирующий континентальный шельф. После достижения склона одиночная волна начинает увеличивать его высоту. Когда вода становится очень мелкой, большая часть волны отражается назад в канал. Однако узкий, но высокий пик воды возникает в конце склона и протекает с пониженной скоростью в исходном направлении падающей волны. Это цунами, которое наконец-то обрушивается на берег, вызывая катастрофические разрушения вдоль береговой линии. Скорость волны, приближающейся к берегу, сравнительно невелика. Волна в конце концов начинает ломаться.
Используя линейную теорию рассеянной воды, высота волны свободной поверхности выше уровня ненарушенной воды в одномерном канале различной глубины является решением следующего дифференциального уравнения с частными производными. (См. [2].)
В этой формуле нижние индексы обозначают частные производные, и - ускорение свободного падения.
Рассмотрим волну, пересекающую линейный уклон из области с постоянной глубиной в область с постоянной глубиной . Преобразование Фурье относительно превращает дифференциальное уравнение с частными производными водной волны в следующее обыкновенное дифференциальное уравнение для режима Фурье .
Для областей с постоянной глубиной , режимами Фурье являются бегущие волны, распространяющиеся в противоположных направлениях с постоянной скоростью .
Решение для глубоководной области является суперпозиция двух волн:
волна, идущая налево с постоянной скоростью
волна, идущая вправо с амплитудой, заданной частотно-зависимым коэффициентом отражения
Этот выбор удовлетворяет волновому уравнению в глубоководной области для любой .
Решение для мелководной области является переданная волна, перемещающаяся налево с постоянной скоростью . Этот выбор удовлетворяет волновому уравнению в мелководной области для любого коэффициента передачи .
Для области перехода (уклон) используйте .
Определите параметры модели цунами следующим образом. Игнорируйте зависимость от частоты в следующих обозначениях: , , .
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
чтобы решить ОДУ для преобразования Фурье из .
wavePDE(x,t) = diff(u,t,t) - g*diff(h(x)*diff(u,x),x); slopeODE(x) = wavePDE(x,0); U(x) = dsolve(slopeODE);
Решение является сложным выражением с участием функций Бесселя. Он содержит две произвольные «константы», которые зависят от .
Const = setdiff(symvar(U), sym([depthratio,g,H,L,x,w]))
Const =
Для любого режима Фурье полное решение должно быть непрерывно дифференцируемой функцией . Следовательно, значения функции и производные должны совпадать в точках шва и . Это обеспечивает четыре линейных уравнения для , , и две константы в .
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);
Замените результаты обратно на , , и .
U(x) = subs(U(x), {Const(1),Const(2)}, {Cvalue1,Cvalue2});
Вы не можете непосредственно оценить решение для потому что и числитель, и знаменатель соответствующих выражений исчезают. Вместо этого найдите пределы низких частот этих выражений.
simplify(limit(U(x), w, 0))
ans =
simplify(limit(R, w, 0))
ans =
simplify(limit(T, w, 0))
ans =
Эти пределы удивительно просты. Они зависят только от отношения значений глубин, определяющих уклон.
Для следующих расчетов используйте эти числовые значения для символьных параметров.
Ускорение свободного падения:
Глубина канала:
Коэффициент глубины между мелкими и глубокими областями:
Длина области уклона:
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 = 0.3; soliton = @(x,t) A.*sech(sqrt(3/4*g*A/H)*(x/c2+t)).^2;
Выбирать точки выборки для . Шкала времени выбирается как кратная (временной) ширине входящего солитона. Сохраните соответствующие дискретизированные частоты преобразования Фурье в .
Nt = 64; TimeScale = 40*sqrt(4/3*H/A/g); W = [0, 1:Nt/2 - 1, -(Nt/2 - 1):-1]'*2*pi/TimeScale;
Выбирать точки выборки в направление для каждой области. Создайте точки выборки для мелководной области, для глубоководной области, и для области уклона.
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);
Вычислите преобразование Фурье входящего солитона на временной сетке равноудаленные точки выборки.
S = fft(soliton(-0.8*TimeScale*c2, linspace(0,TimeScale,2*(Nt/2)-1)))'; S = repmat(S,1,Nx);
Создайте решение бегущей волны в глубоководной области на основе данных Фурье в .
soliton = real(ifft(S .* exp(1i*W*X2/c2)));
Преобразуйте режимы Фурье отраженной волны в глубоководной области в числовые значения по сетке в пространство. Умножите эти значения коэффициентами Фурье в и используйте функцию ifft
для вычисления отраженной волны в пространство. Обратите внимание, что первая строка числовых данных состоит из значений NaN, поскольку надлежащая числовая оценка символьных данных для невозможно. Задайте значения в первой строке как пределы низкой частоты.
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 км, что соответствует скорости .) Над глубоким морем амплитуда довольно небольшая, часто около 0,5 м и менее. При распространении на шельф, однако, цунами резко увеличивают их высоту: сообщалось об амплитудах до 30 м и более.
Одно из интересных явлений заключается в том, что, хотя цунами обычно приближаются к береговой линии как к фронту волны, простирающемуся на сотни километров перпендикулярно направлению их движения, они не наносят равномерного ущерба вдоль побережья. В некоторых точках они вызывают бедствия, в то время как в других местах наблюдаются только умеренные волновые явления. Это вызвано различными склонами от морского дна до континентального шельфа. На самом деле очень крутые склоны заставляют большую часть цунами отражаться назад в область глубокой воды, в то время как небольшие склоны отражают меньше волны, передавая узкую, но высокую волну, несущую много энергии.
Запустите симуляцию для различных значений , которые соответствуют разным склонам. Чем круче склон, тем ниже и менее мощна волна, которая передается.
Обратите внимание, что эта модель игнорирует эффекты дисперсии и трения. На полке симуляция теряет свой физический смысл. Здесь важны эффекты трения, вызывающие разрушение волн.
[1] Дерек Г. Горинг и Ф. Райхлен, Цунами - Распространение длинных волн на шельф, Journal of Waterway, Port, Coastal and Ocean Engineering 118 (1), 1992, pp. 41 - 63.
[2] H. Lamb, Hydrodynamics, Dover, 1932.