Azzera filtri
Azzera filtri

HOW TO CREATE TSUNAMI MODEL

6 visualizzazioni (ultimi 30 giorni)
meoui
meoui il 17 Mag 2024
Spostato: Sam Chak il 17 Mag 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 Commenti
Sam Chak
Sam Chak il 17 Mag 2024
Spostato: Sam Chak il 17 Mag 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 il 17 Mag 2024
Spostato: Sam Chak il 17 Mag 2024
do you use this formula? for the coding

Accedi per commentare.

Risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by