Welcome! Log In Create A New Profile

Advanced

Dam Break Simulation Code

Posted by CarloCarlo  
Dam Break Simulation Code
October 26, 2013 03:08PM
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
The function animazione3D is not defined in MATLAB (At line 153).
Re: Dam Break Simulation Code
November 01, 2013 08:26PM
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
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
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
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.