Hi everyone,
I’m new to the forum and to the LBM too. In this months I’ve being studying the theory behind the multiphase models (Swift binary liquid and Shan Chen). I’ve done some benchmark case in regard to the contact angle with a solid surface.
Now I’m trying to implement a bi-dimensional Dam Break model but I’m finding some difficulties. I was wondering if someone with more experience than me (everyone in here I presume 
) could tell me what model is the best in simulating such a problem.
I’ve already done some tests by modifying a script I found, but I’m unable to perform large-density-ratio simulation. I list below the matlab code I’ve created. Have anyone faced the same problem?
Any help is appreciated.
P.s. I’m a student in Civil Engineering by the way, be merciful
[i]%% THE MATLAB CODE
% dimension
lx = 1;
ly = 2;
N = 128;
N_iter = 3500;
N_out = 20;
% Set liquid propriety in LB system
radius = 25;
rhoLiquid=1.818
rhoGas=0.177
rhoWall=0.75
G=-4.9;
c_s = 1/sqrt(3);
%set LB Units
tau = 1;
delta_x = ly/N;
nx = round(lx/delta_x);
ny = round(ly/delta_x);
gravity = -0.0001
% lattice boltzmann initial velocities
ux = 0*ones(nx,ny);
uy = 0*ones(nx,ny);
% Some initialization
fX = zeros(nx,ny);
fY = zeros(nx,ny);
omega = 1/tau;
counterFrame = 1;
[X,Y] = meshgrid(1:nx,1:ny);
%% Liquid initial position
%---> walls
isWall = logical(zeros( nx, ny ));
isWall (1:nx,[1,ny]) = 1;
isWall ([1,nx],1:ny) = 1;
%---> Liquid
isLiquid = logical(zeros(nx,ny));
isLiquid (2:round(nx)/4,2:round(0.7*ny)) = 1;
%---> Gas
isGas = ~isLiquid - isWall;
% D3Q19 model grid initializzation
[ fIn, fOut, fEq, w, cx, cy, opp ] = initializeCells( nx, ny );
%% ------------------------> Simulation initializzation
rho = isLiquid.*rhoLiquid + isGas.*rhoGas + isWall.*rhoWall;
for l = 1:9
fEq(l,:,:)=w(l)*rho.*(1+3*(ux.*cx(l)+uy.*cy(l)) ...
+ 9/2*((cx(l)*cx(l)-1/3).*ux.*ux+2*cx(l)*cy(l)*ux.*uy+(cy(l)*cy(l)-1/3).*uy.*uy));
end
fIn = fEq;
fOut = fIn;
%% ------------------------> Time cicles
for step = 1:N_iter
%Compute macroscopic Variables
fIn = reshape(fIn,9, nx*ny);
rho = sum(fIn);
ux=(cx*fIn)./rho;
uy=(cy*fIn)./rho;
rho = reshape(rho,nx,ny);
ux = reshape(ux,nx,ny);
uy = reshape(uy,nx,ny);
fIn = reshape(fIn,9,nx,ny);
%Impose wall density
rho(isWall) = rhoWall;
ux (isWall) = 0;
uy (isWall) = 0;
%Force evaluation
force_sum_x=0.0;
force_sum_y=0.0;
for l=2:9
psi=1-exp(-circshift(rho, [cx(opp(l)),cy(opp(l))]));
force_sum_x=force_sum_x-G*w(l)*psi*cx(l);
force_sum_y=force_sum_y-G*w(l)*psi*cy(l);
end
fX=(1-exp(-rho)).*force_sum_x;
fY=(1-exp(-rho)).*force_sum_y + gravity;
%velocity after forcing
ux = ux + 0.5*fX./rho;
uy = uy + 0.5*fY./rho;
%collide & streaming
for l = 1:9
fEq(l,:,:)=w(l)*rho.*(1+3*(ux.*cx(l)+uy.*cy(l)) ...
+ 9/2*((cx(l)*cx(l)-1/3).*ux.*ux+2*cx(l)*cy(l)*ux.*uy+(cy(l)*cy(l)-1/3).*uy.*uy));
fPop(l,:,:)=w(l)*(1-0.5*omega)*((3*(cx(l)-ux)+9*cx(l)*(cx(l).*ux+cy(l).*uy)).*fX...
+(3*(cy(l)-uy)+9*cy(l).*(cx(l).*ux+cy(l).*uy)).*fY);
fOut(l,:,:) = fIn(l,:,:) - omega .* (fIn(l,:,:)-fEq(l,:,:)) + fPop(l,:,:); %collide
fIn(l,:,:) = circshift(fOut(l,:,:),[0,cx(l),cy(l)]); %stream
end
%Apply Bounce-Back no-slip
for x = 1:nx
for y=1:ny
for l = 1:9
newx=1+mod(x-1+cx(opp(l))+nx,nx);
newy=1+mod(y-1+cy(opp(l))+ny,ny);
if isWall(x,y) == 1
% boundary cell
fIn(opp(l),newx,newy) = fIn(l,x,y);
end
end
end
end
%Post processing utility
if mod(step, N_out) == 0
step
img = imagesc(flipud(rho'));
%img = quiver(X,Y,ux',uy');
axis equal
F(counterFrame) = getframe(gca);
counterFrame=counterFrame+1;
massaTot = sum(sum(rho))
end
end[/i]