HOW TO CREATE TSUNAMI MODEL

be discovered:

A = 5

X = 0 - 600

t = 0 - 60 minute

long = 0,5 - 30 kilometers

count the lamda, sigma and plot psi

[KSSV]

Do you have any reference?

[KSSV]

Share the reference documnet. How these fomulas are going to help?

[Sam Chak]

Hi @Fadilla Neisa

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]