Решение дифференциальных уравнений с частными производными

Этот пример симулирует явление волны цунами при помощи Symbolic Math Toolbox™, чтобы решить дифференциальные уравнения.

Эта симуляция является упрощенной визуализацией явления и основана на статье путем Бодания и Raichlen [1].

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

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

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

utt=g(hux)x

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

Рассмотрите волну, пересекающую линейный наклон h(x) из области с постоянной глубиной h2 в область с постоянной глубиной h1h2. Преобразование Фурье относительно t поворачивает дифференциальное уравнение с частными производными водной волны к следующему обыкновенному дифференциальному уравнению для режима Фурье u(x,t)=U(x,ω)eiωt.

-ω2U=g(hUx)x

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

u(x,t)=C1(ω)eiω(t+x/c)+C2(ω)eiω(t-x/c)

Решение u2(x,t)=eiω(t+x/c2)+R(ω)eiω(t-x/c2) поскольку глубоководная область является суперпозицией двух волн:

  • волна, перемещающаяся налево с постоянной скоростью c2=gh2

  • волна, перемещающаяся направо с амплитудой, данной зависимым частоты отражательный коэффициент R(ω)

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

Решение u1(x,t)=T(ω)eiω(t+x/c1) для мелководья область является переданной волной, перемещающейся налево с постоянной скоростью c1=gh1. Этот выбор u1 удовлетворяет уравнению волны в мелководной области для любого коэффициента передачи T(ω).

Для области перехода (наклон), использовать u(x,t)=U(x,w)eiωt.

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

Задайте параметры модели цунами можно следующим образом. Игнорируйте зависимость от частоты ω в следующих обозначениях: R=R(ω), T=T(ω), 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 сложное выражение, включающее Функции Бесселя. Это содержит две произвольных "константы", которые зависят от ω.

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

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

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

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

2depthratio+12/(sqrt (depthratio) + 1)

simplify(limit(R, w, 0)) 
ans = 

-depthratio-1depthratio+1- (sqrt (depthratio) - 1) / (sqrt (depthratio) + 1)

simplify(limit(T, w, 0))
ans = 

2depthratio+12/(sqrt (depthratio) + 1)

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

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

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

  • Гравитационное ускорение: g=9.81m/sec2

  • Глубина канала: H=1m

  • Отношение глубины между отмелью и глубокими областями: depthratio=0.04

  • Длина наклонной области: L=2m

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 перемещение налево с постоянной скоростью c2 в глубоководной области.

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,ω) пробел. Умножьте эти значения с коэффициентами Фурье в S и используйте функциональный ifft вычислить отраженную волну в (x,t) пробел. Обратите внимание на то, что первая строка числовых данных R состоит из значений NaN потому что соответствующая численная оценка символьных данных R для ω=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 км, соответствуя скорости gh700km/hour.) По глубокому морю амплитуда скорее мала, часто приблизительно 0,5 м или меньше. При распространении на полку, однако, цунами увеличивают свою высоту существенно: амплитуды до 30 м и о большем сообщили.

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

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

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

Ссылки

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

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

Для просмотра документации необходимо авторизоваться на сайте