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 “http://ascelibrary.org/doi/pdf/10.1061/41036(342)338” 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)