Borrar filtros
Borrar filtros

HOW TO CREATE TSUNAMI MODEL

5 visualizaciones (últimos 30 días)
meoui
meoui el 17 de Mayo de 2024
Movida: Sam Chak el 17 de Mayo de 2024
be discovered:
A = 5
X = 0 - 600
t = 0 - 60 minute
long = 0,5 - 30 kilometers
count the lamda, sigma and plot psi
  4 comentarios
Sam Chak
Sam Chak el 17 de Mayo de 2024
Movida: Sam Chak el 17 de Mayo de 2024
Looking for this?
%% Tsunami model
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));
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]));
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});
%% Parameters
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;
%% incoming soliton wave
soliton = @(x,t) A.*sech(sqrt(3/4*g*A/H)*(x/c2+t)).^2;
%% creating time scale and sample points
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); % slope region
X1 = linspace(x_min, L1, Nx); % shallow water region
X2 = linspace(L2, x_max, Nx); % deep water region
%% Fourier transform of the incoming soliton
S = fft(soliton(- 0.8*TimeScale*c2, linspace(0, TimeScale, 2*(Nt/2) - 1)))';
S = repmat(S, 1, Nx);
%% Construct a traveling wave solution based on S
soliton = real(ifft(S.*exp(1i*W*X2/c2)));
%% Construct a reflected wave
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)));
%% Construct a transmitted wave
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)));
%% Construct the wave at the slope region
U12 = double(subs(subs(vpa(subs(U(x))), w, W), x, X12));
U12(1,:) = double(2/(1 + sqrt(depthratio)));
U12 = real(ifft(S.*U12));
%% Plot the Solution
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], 'linewidth', 2, 'Color', '#BEA256')
axis([x_min, x_max, -H-0.1, 0.6]), grid on
hold on
line1 = plot(X1, transmittedWave(1,:), 'linewidth', 4, 'Color', '#8CD0E1');
line12 = plot(X12, U12(1,:), 'linewidth', 3, 'Color', '#2AA7C0');
line2 = plot(X2, soliton(1,:) + reflectedWave(1,:), 'linewidth', 2, 'Color', '#0A7B88');
text(6, -0.6, 'Deep Water region')
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
meoui
meoui el 17 de Mayo de 2024
Movida: Sam Chak el 17 de Mayo de 2024
do you use this formula? for the coding

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Symbolic Computations in MATLAB en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by