dambreak

%A simple example to simulate a dambreak in 1D
%basic instruction to write a VolcFlow file
nrow = 3;   %equations are in 2D. At least 3 lignes (or columns) should be defined, even for 1D simulation
ncol = 500; %number of column of the simulation
dx_horiz = 3/ncol;   %defines the space step (in m) as a function
                     %of the domain length and the number of cells chosen
dy_horiz = dx_horiz; %space step in y
dt=0.001;         %time step (in s)
dtplot = dt*50;   %time step for plot
tmax = 0.3;         %duration of the simulation (s)
tstep_adjust = 0; %0 = the time step is fixed by the user, 1 = its varies in time using a CFL condition
g=9.81;                  %gravity (m/s2)
delta_int = 0/180*pi;    %internal friction angle (for Coulomb behaviour). Set to 0 here. (rad)
delta_bed = 0/180*pi;  %external friction angle (for Coulomb behaviour). Set to 0 here. (rad)
cohesion =    0;         %cohesion or yield strength (depending of the rheology used), Pa
rho = 1000;              %density (kg/m3)
viscosity = 0;           %dynamic viscosity in Pa s
representation = ‘cla; patch( »xdata », [x(1) x x(end)], »ydata », [0 h(2,:) 0],  »facecolor », [0 0.8 1]); patch( »xdata », [x(1) x x(end)], »ydata », [-0.5 h(2,:)*0 -0.5], »facecolor », [0.5 0.5 0.5]); axis([x(1) x(ncol) -0.5 2]);‘; %function used for the visualisation of results (at each dtplot)
f_avi =  »;  %avi file done with the plots.  » = no movie is done
f_data =  »; %save various variables at each dtplot.  » = no data is saved
%boundary conditions should be defined for 1D simulation
bound_cond = ‘h2(1,:)=h2(2,:); h2(nrow,:)=h2(2,:); u_xx(:,1)=0; u_xx(:,ncol-1)=0; u_yy(1,:)=0; u_yy(nrow-1,:)=0;’;
%definition of the initial thickness
h(1:nrow, 1:ncol)=0.5;
h(1:nrow, 1:ncol/2)=1.5;
%initial topography: flat, z=0, same size than h
z = h*0;
%vectors x and y – not needed by VolcFlow (for the visualisation of results only)
x = [dx_horiz/2:dx_horiz:ncol*dx_horiz-dx_horiz/2];
y = [dy_horiz/2:dy_horiz:nrow*dy_horiz-dy_horiz/2];
%initial velocity – 0 here
u_xx = z(1:nrow, 1:ncol-1)*0;
u_xy = z(1:nrow, 1:ncol-1)*0;
u_yy = z(1:nrow-1, 1:ncol)*0;
u_yx = z(1:nrow-1, 1:ncol)*0;