# Dam Break Simulation Code

Posted by Carlo
 Dam Break Simulation Code October 26, 2013 03:08PM Registered: 5 years ago Posts: 2
Hi everyone,
After trying to simulate a dam break problem with bi-phase LBM I decided to use a "simplest" model for shallow water cases.
I'm now dealing with some problem (I think related with the physical inputs conversion to LB units, but I'm not sure at all): the wave seems to move too slow in my opinion.

I've also tried to reproduce the results in "[ascelibrary.org]; but the model of mine collapse.

Can you search the MATLAB code for errors?

Any help will be greatly appreciated!

%% SHALLOW WATER SCRIPT
%Time definition
Nsteps = 1000;
N_out = 10;
% Simulation dimensions [m][second]
lx = 2;             %[m]
ly = 0.1;             %[m]

N = 50;            %[-]  NB: delta_x = round(lx/N)

delta_t = 10^-3    %[seconds]
delta_x = min(lx,ly) / N
c = delta_t/delta_x
viscosity = 10^-3;  %[m^2/s]
gravity = 9.81;     %[m/s^2]

% Evaluate lattice Units
nx = round(lx/delta_x);
ny = round(ly/delta_x);
g = gravity*delta_t^2/delta_x
c_s = 1/sqrt(3);
tau = (viscosity/c_s^2 * (delta_t/delta_x^2)) + 0.5

%return
h_max = 1;
h_min = 3;

omega = 1/tau;
[X,Y] = meshgrid(1:nx,1:ny);

% D2Q9 parameters
w =     [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36];
cx =    [0 1 0 -1 0 1 -1 -1 1];
cy =    [0 0 1 0 -1 1 1 -1 -1];
cxValue =    [0 c 0 -c 0 c -c -c c];
cyValue =    [0 0 c 0 -c c c -c -c];
opp =   [1 4 5 2 3 8 9 6 7];

%% INITIALIZZATION
% boundary matrix isWall
isWall = zeros(nx,ny);
isWall (1,2:ny-1) = 1;
isWall (nx,2:ny-1) = 1;
isWall (1:nx,ny) = 1;
isWall (1:nx,1) = 1;

%isWall (round(0.8*nx),2:round(0.3*ny)) = 1;
%isWall (round(0.8*nx),round(0.7*ny):ny) = 1;
%isWall (round(0.75*nx):round(0.8*nx),round(0.3*ny)) = 1;
%isWall (round(0.75*nx):round(0.8*nx),round(0.7*ny)) = 1;

%isWall (round(0.35*nx):round(0.45*nx),round(0.49*ny):round(0.51*ny)) = 1;

bbRegion = find(isWall);

%initialize height over bed for the domain
length_h_max = round(0.8*nx);
h = zeros(nx,ny);
h (1:length_h_max,:) = h_max;
h (length_h_max:nx,:) = h_min;

h (bbRegion) = 1;

%initialize bed geometry for the domain (zeros if there is no forcing term)
Z = zeros(nx,ny);

%initialize velocity of the sistem
ux = zeros(nx,ny);
uy = zeros(nx,ny);

%some other initializzation
fEq =  zeros(9,nx,ny);
fPop = zeros(9,nx,ny);
fX =   zeros(nx,ny);
fY =   zeros(nx,ny);
mEq =  zeros(9,nx,ny);
mIn = zeros(9,nx,ny);
figure1 = figure('Color',[1 1 1],'units','normalized','outerposition',[0 0 1 1]);
nz = max(h_min, h_max);
counterFrame = 1;

%% Construction of the initial equilibrium
fEq =   zeros(9,nx,ny);
pressure =  3 * g * h / (2 * c^2);
for l = 2:9
cu        = 3*(cxValue(l)*ux+cyValue(l)*uy) / (c^2);
fEq(l,:,:)   = h.* w(l) .* ...
( pressure + cu + 0.5*(cu.*cu) - 1.5*(ux.^2+uy.^2) / (c^2) );
end
% fEq expression for the central node
fEq(1,:,:) = h - reshape(sum(reshape(fEq,9,nx*ny)),nx,ny);
fIn = fEq;
fOut = fEq;

%% time cycle
for step = 1:Nsteps
step
%recover macroscopic variables
h = reshape(sum(reshape(fIn,9,nx*ny)),nx,ny);
h(bbRegion) = 0;
%Z(bbRegion) = 0;
ux  = reshape(cxValue*(reshape(fIn,9,nx*ny)),nx,ny)./h;
ux(bbRegion) = 0;
uy  = reshape(cyValue*(reshape(fIn,9,nx*ny)),nx,ny)./h;
uy(bbRegion) = 0;

%Construction of the equilibrium
fEq =   zeros(9,nx,ny);
pressure = 3 * g * h /  (2 * c^2);
for l = 2:9
cu        = 3*(cxValue(l)*ux+cyValue(l)*uy) / (c^2);
fEq(l,:,:)   = h.* w(l) .* ...
( pressure + cu + 0.5*(cu.*cu) - 1.5*(ux.^2+uy.^2) / (c^2) );
end
fEq(1,:,:) = h - reshape(sum(reshape(fEq,9,nx*ny)),nx,ny);

% Evaluation of symmetric and antisimmetric part of fIn & fEq (for TRT)
%     for l = 1:9
%         fIn_sym(l,:,:) = ( fIn(l,:,:) + fIn(opp(l),:,:) )/2;
%         fIn_skw(l,:,:) = ( fIn(l,:,:) - fIn(opp(l),:,:) )/2;
%         fEq_sym(l,:,:) = ( fEq(l,:,:) + fEq(opp(l),:,:) )/2;
%         fEq_skw(l,:,:) = ( fEq(l,:,:) - fEq(opp(l),:,:) )/2;
%     end
%     tau_skw = tau;
%     tau_sym = 0.5 + 1/(12*(tau_skw - 0.5));
%     omega_skw = 1/tau_skw;
%     omega_sym = 1/tau_sym;

%effettuo le collisioni
fX = 0;
fY = 0;
for l = 1:9
fPop(l,:,:)=w(l)*(1-0.5*omega).*((3*(cxValue(l)-ux)+9*cxValue(l)*(cxValue(l).*ux+cyValue(l).*uy)).*fX...
+(3*(cyValue(l)-uy)+9*cyValue(l)*(cxValue(l).*ux+cyValue(l).*uy)).*fY);

fOut (l,:,:) = fIn(l,:,:)*(1 - omega) + omega*fEq(l,:,:) + (delta_t/6*c^2)*fPop(l,:,:);
%fOut (l,:,:) = fIn(l,:,:) - omega_sym*(fIn_sym(l,:,:) - fEq_sym(l,:,:)) - omega_skw*(fIn_skw(l,:,:) - fEq_skw(l,:,:)) + (delta_t/6*c^2)*fPop(l,:,:);
end

%Apply Bounce-Back no-slip
for l=1:9
fOut(l,bbRegion) = fIn(opp(l),bbRegion);
end
% stream DF
for l = 1:9
fIn(l,:,:) = circshift(fOut(l,:,:),[0,cx(l),cy(l)]);
end

% Post processing utilitys
if mod(step, N_out) == 0
% 3D animation
axes_animazione3D = subplot(2,2,1);
animazione3D (X*delta_x, Y*delta_x, flipud(Z'+ h'), flipud(Z'), nx*delta_x, ny*delta_x, nz, axes_animazione3D);
axes_topView = subplot(2,2,2);
%top view
topView( flipud(Z + h), nx, ny, axes_topView )
F(counterFrame) = getframe(gcf);
time = delta_t*step;
time_str = strcat('elapsed time: ', num2str(time), ' [sec]');
suptitle(time_str);
% middle sextion
subplot(2,2,3);
plot(1:nx, (Z(:,round(ny/2))' + h(:,round(ny/2))'));
drawnow;

hTot(counterFrame) = sum(sum(h));
counterFrame=counterFrame+1;
end
end

figure()
plot(1:counterFrame-1,hTot)
 Re: Dam Break Simulation Code October 29, 2013 03:57AM Registered: 6 years ago Posts: 10
The function animazione3D is not defined in MATLAB (At line 153).
 Re: Dam Break Simulation Code November 01, 2013 08:26PM Registered: 6 years ago Posts: 7
Dear Carlo
How do you select delta_t=0.001??
It seems that it should be selected based on delta_x and your boundary condition not spurious
good luck
 Re: Dam Break Simulation Code November 09, 2013 08:45AM Registered: 5 years ago Posts: 1
Hello
Please, could you give me a (serial) fortran code with LBM??
I need it for any natural convection problem
 Re: Dam Break Simulation Code September 10, 2014 10:09AM Registered: 5 years ago Posts: 14
Dear all,
did someone make any progress with shallow water equations?
I was wondering if anyone thought of implementing wave inlets as BCs (such as Stokes I, Stokes II, irregular waves...) to be used in hydrodynamics simulations.
I think it could be an interesting area of development.
Christian
 Re: Dam Break Simulation Code May 13, 2017 12:33AM Registered: 8 years ago Posts: 2
Hello Carlo. I was running your code and i found two inconsistencies, one of them is the velocity, i think it must be the inverse the one you declared in the code, it will make the fluid move faster. And the other one is the coefficient of force is not correct. Can you tell me the bibliographic source that inspired you to use that term? Thanks
Sorry, you do not have permission to post/reply in this forum.