ciao a tutti
sto implementando un codice in Matlab in cui risolvo un problema di rottura della diga 1D con parete riflettente ma l'onda di Shock, arrivata alla fine del dominio, non rimbalza sulla parete anzi, già dall'inizio dell'esecuzione c'è un'altezza anomala alla fine del dominio spaziale; questo è il codice:
codice:
clear all
close all
clc
L = 1;
b = 0.25 * L;
global g
g = 9.81
hL = 2; %altezza iniziale a sinistra
hR = 0.1; %altezza iniziale a destra
uL = 0; %velocità iniziale a sinistra
uR = 0; %velocità iniziale a destra
%dominio spaziale
xL = -2 * L;
xR = L - b;
%dominio temporale
dt = 0.001;
t_in = dt;
t_out = 1;
IMAX = 100; % numero di punti per graficare il risultato
NMAX = 100000; % numero massimo di passi temporali
x = linspace(xL,xR,IMAX); %divido il dominio spaziale in IMAX celle
t=t_in; % pongo il tempo pari a t_in
for i=1:IMAX %condizioni iniziali
if(x(i)<0)
hi(i)=hL;
ui(i)=uL;
else
hi(i)=hR;
ui(i)=uR;
end
end
subplot(2,1,1) %grafico le condizioni iniziali
plot(x,hi,'r-')
subplot(2,1,2)
plot(x,ui,'r-')
drawnow
for j=1:NMAX %inizio ciclo temporale
if(t+dt>t_out)
dt=t_out-t;
end
if(t>=t_out)
break
end
for i=1:IMAX %inizio ciclo spaziale
xi = x(i)/t;
if(i==IMAX) %se ultima casella parete riflettente
[he(i),ue(i)]=ExactRiemannSWE(hL,uL,hR,-uR,-xi);%risolutore esatto del problema di Riemann
else %altrimenti
[he(i),ue(i)]=ExactRiemannSWE(hL,uL,hR,uR,xi); %risolutore esatto del problema di Riemann
end %termine if
end %termine ciclo spaziale
subplot(2,1,1) %grafico he ed ue
plot(x,he,'r-')
subplot(2,1,2)
plot(x,ue,'r-')
drawnow
t=t+dt; %incremento t
end %termine ciclo temporale
i risolutori di Riemann sono corretti, qualcuno può darmi una mano?
Grazie