Hi everyone,
I translated the fortran code from Haibo Huang’s book (Multiphase Lattice Boltzmann Method) into Matlab, but somehow I didn’t get the parabolic velocity profile as it supposed to be. Please do help me if I have mistakes in writing my code or the input parameter for my simulatiion is making the run unstable. I will also appreciate if someone can share multiphase color gradient code with me (billal.m.aslam@gmail.com)
Warm Regards,
Billal
Here are my Matlab code :
% Code Purpose :
% Simulate 2D multiphase poiseuille flow using LB Method (D2q9q2)
% Model : Rothman-Keller Color Gradient Model (improved, Huang(2014))
% Fluids: wetting (red), nonwetting (blue)
% Lattice Numbering
% c6 c2 c5
% \ | /
% c3 -c9 - c1
% / | \
% c7 c4 c8
%***********************************************************************
clear; clc;
% Computational Domain & Lattice Variables*****************
global Nx Ny cx cy cs w
Nx = 2;
Ny = 101;
u_x=zeros(Nx,Ny);
u_y=zeros(Nx,Ny);
rho=zeros(Nx,Ny);
fequ1=zeros(9,1);
fequ2=zeros(9,1);
feq=zeros(9,1);
ff=zeros(9,Nx,Ny); %pDF with 9 vel, rest velocity is last
cx = [1 0 -1 0 1 -1 -1 1 0]; % velocities, x components
cy = [0 1 0 -1 1 1 -1 -1 0]; % velocities, y components
w = [1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36 4/9]; % weights
Bi = [2/27 2/27 2/27 2/27 5/108 5/108 5/108 5/108 -4/27]; %2nd collision wt
r0=0; r1=1; r2=sqrt(2);
rsq= [r1 r1 r1 r1 r2 r2 r2 r2 r0]; %length of velocity vector
cc = 1; %lattice speed dx*=1, dt*=1
cs = cc*cc/3; %speed of sound squared
%RK model parameters*******************************
beta = 0.5; %tunable parameter, affect IFT 0<beta<1
A = 0.0001; %tunable parameter, affect IFT
aa1 = 25; %width of inner fluid in latice nodes (red)
alfa_r = 0.8;
alfa_b = 0.4;
rho_ri = 1-alfa_b;
rho_bi = 1-alfa_r;
P0_r = 3*(1-alfa_r)/5;
P0_b = 3*(1-alfa_b)/5;
%%% SS convergence parameter and timestep
t_max = 10000; %max time step
teval=100;
tol=1e-10;
u_old=u_x;
tau_r = 0.8;
tau_b = 1;
tau=0; %modified at interface
%%% relaxation parameter at interface
delta=0.98;
s1=2tau_rtau_b/(tau_r+tau_b);
t1=s1;
s2=2*(tau_r-s1)/delta;
s3=-s2/(2delta);
t2=2(t1-tau_b)/delta;
t3=t2/(2*delta);
%unwanted terms variables
pr=zeros(Ny,1); %used in density fd (red)
par=zeros(Ny,1);
pb=zeros(Ny,1); %used in density fd (blue)
pab=zeros(Ny,1);
cr=zeros(Ny,1);
cb=zeros(Ny,1);
% Forcing terms
df=0.00000015;
G0=zeros(Ny,1);
%**************************************************
%Flow Obstacles/Conduit Geometry
bound=zeros(Nx,Ny);
for j=1:Ny
for i=1:Nx
if j==1||j==Ny
bound(i,j)=1; %2D pipe
end
end
end
%Initialize macroscopic variables
f_r1=zeros(9,1);
f_b1=zeros(9,1);
f_r=zeros(9,Nx,Ny);
f_b=zeros(9,Nx,Ny);
Fx=zeros(Nx,Ny);
Fy=zeros(Nx,Ny);
coslai=zeros(9,1);
rho_r=zeros(Nx,Ny);
rho_b=zeros(Nx,Ny);
for j=1:Ny
for i=1:Nx
rho_b(i,j)=rho_bi;
end
end
%Initialize fluid distribution and feq
for j=1:Ny
for i=1:Nx
if abs(j-(1+Ny)/2)<=aa1
rho_r(i,j)=rho_ri; %red fluid at center
rho_b(i,j)=0;
end
rho(i,j)=rho_r(i,j)+rho_b(i,j);
f_r1=get_feq(rho_r(i,j),u_x(i,j),u_y(i,j),alfa_r,f_r1);
f_b1=get_feq(rho_b(i,j),u_x(i,j),u_y(i,j),alfa_b,f_b1);
for k=1:9
f_r(k,i,j)=f_r1(k);
f_b(k,i,j)=f_b1(k);
ff(k,i,j)=f_r1(k)+f_b1(k);
end
if bound(i,j)==1
rho_r(i,j)=0;
rho_b(i,j)=1;
end
end
end
tic;
%% Main Loop**********************************************************
for t=1:t_max
disp([‘LBM-RK is running = ’ num2str(t) ’ timestep’]);
%calculate error after t timestep and write output (to be updated)
% if mod(t,teval)==1
% conv = abs(mean(u_x(:))/mean(u_old(:))-1);
% disp(['relative change = ’ num2str(conv)]);
% if conv<tol
% break
% else
% u_old = u_x;
% end
% end
%streaming step
f_r=stream(f_r);
f_b=stream(f_b);
%calculate macro variables
for j=1:Ny
for i=1:Nx
if bound(i,j)==0
%density
rho_r(i,j)=f_r(1,i,j)+f_r(2,i,j)+f_r(3,i,j)...
+f_r(4,i,j)+f_r(5,i,j)+f_r(6,i,j)...
+f_r(7,i,j)+f_r(8,i,j)+f_r(9,i,j);
rho_b(i,j)=f_b(1,i,j)+f_b(2,i,j)+f_b(3,i,j)...
+f_b(4,i,j)+f_b(5,i,j)+f_b(6,i,j)...
+f_b(7,i,j)+f_b(8,i,j)+f_b(9,i,j);
rho(i,j)=rho_r(i,j)+rho_b(i,j);
%x-velocity
u_x(i,j)=(f_b(1,i,j)+f_b(5,i,j)+f_b(8,i,j))...
-(f_b(3,i,j)+f_b(6,i,j)+f_b(7,i,j))...
+(f_r(1,i,j)+f_r(5,i,j)+f_r(8,i,j))-...
(f_r(3,i,j)+f_r(6,i,j)+f_r(7,i,j));
%y-velocity
u_y(i,j)=(f_b(2,i,j)+f_b(5,i,j)+f_b(6,i,j))...
-(f_b(4,i,j)+f_b(7,i,j)+f_b(8,i,j))...
+(f_r(2,i,j)+f_r(5,i,j)+f_r(6,i,j))-...
(f_r(4,i,j)+f_r(7,i,j)+f_r(8,i,j));
else
rho_r(i,j)=0;
rho_b(i,j)=rho_bi; %blue as wetting fluid
end
end
end
%collision step
%%% get unwanted term in RK model (read p102-103 Hb.Huang's book)
for j=1:Ny
if j>1 && j<Ny %do we need to calculate at each x?
pr(j)=0.5*(rho_r(1,j+1)-rho_r(1,j-1));
pb(j)=0.5*(rho_b(1,j+1)-rho_b(1,j-1));
end
end
for j=1:Ny
if j>1 && j<Ny %do we need to calculate at each x?
par(j)=0.5*(pr(j+1)*u_x(1,j+1)-pr(j-1)*u_x(1,j-1));
pab(j)=0.5*(pb(j+1)*u_x(1,j+1)-pb(j-1)*u_x(1,j-1));
end
end
for j=1:Ny
for i=1:Nx
if bound(i,j)==0
fequ1=get_feq(rho_r(i,j), u_x(i,j), u_y(i,j), alfa_r,fequ1);
fequ2=get_feq(rho_b(i,j), u_x(i,j), u_y(i,j), alfa_b,fequ2);
psi=(rho_r(i,j)-rho_b(i,j))/(rho_r(i,j)+rho_b(i,j));
%calculate tau at interface
if psi>delta
tau=tau_r;
elseif psi>0 && psi<=delta
tau=s1+s2*psi+s3*psi*psi;
elseif psi<=0 && psi>=-delta
tau=t1+t2*psi+t3*psi*psi;
elseif psi<-delta
tau=tau_b;
end
%adding body force
if abs(j-(1+Ny)/2)<=aa1
G0(j)=df;
end
%eliminating unwanted extra terms (forcing in x dir)
cr(j)=(tau_r-0.5)*((1/3)-P0_r)*par(j);
cb(j)=(tau_b-0.5)*((1/3)-P0_b)*pab(j);
for k=1:9
% f_r(k,i,j)=fequ1(k)+(1-1/tau)*(f_r(k,i,j)-fequ1(k))…
% -cr(j)*w(k)cx(k)/cs;
% f_b(k,i,j)=fequ2(k)+(1-1/tau)(f_b(k,i,j)-fequ2(k))…
% -cb(j)*w(k)*cx(k)/cs;
%body force added as source term
ff(k,i,j)=f_r(k,i,j)+f_b(k,i,j)+w(k)*cx(k)*G0(j)/cs;
end
end
end
end
%recoloring step
for j=1:Ny
for i=1:Nx
if bound(i,j)==0
y_n=mod(j,Ny)+1;
x_e=mod(i,Nx)+1;
y_s=Ny-mod(Ny+1-j,Ny);
x_w=Nx-mod(Nx+1-i,Nx);
%calculate color gradient vector
Fx(i,j) = cx(1) *(rho_r(x_e,j) - rho_b(x_e,j))...
+ cx(5) *(rho_r(x_e,y_n) - rho_b(x_e,y_n))...
+ cx(8) *(rho_r(x_e,y_s) - rho_b(x_e,y_s))...
+ cx(3) *(rho_r(x_w,j) - rho_b(x_w,j))...
+ cx(6) *(rho_r(x_w,y_n) - rho_b(x_w,y_n))...
+ cx(7) *(rho_r(x_w,y_s) - rho_b(x_w,y_s));
Fy(i,j) = cy(2) *(rho_r(i,y_n) - rho_b(i,y_n))...
+ cy(5) *(rho_r(x_e,y_n) - rho_b(x_e,y_n))...
+ cy(6) *(rho_r(x_w,y_n) - rho_b(x_w,y_n))...
+ cy(4) *(rho_r(i,y_s) - rho_b(i,y_s))...
+ cy(7) *(rho_r(x_w,y_s) - rho_b(x_w,y_s))...
+ cy(8) *(rho_r(x_e,y_s) - rho_b(x_e,y_s));
fm=sqrt(Fx(i,j)*Fx(i,j)+Fy(i,j)*Fy(i,j));
if fm<0.00000001 %avoid zero denominator
for k=1:9
f_r(k,i,j)=rho_r(i,j)/rho(i,j)*ff(k,i,j);
f_b(k,i,j)=rho_b(i,j)/rho(i,j)*ff(k,i,j);
end
else
%add second collision operator
k=9;
ff(k,i,j)=ff(k,i,j)+A*fm*(-Bi(k));
for k=1:8
coslai(k)=(cx(k)*Fx(i,j)+cy(k)*Fy(i,j))/rsq(k)/fm;
ff(k,i,j)=ff(k,i,j)+...
A*fm*(w(k)*coslai(k)*coslai(k)*rsq(k)*rsq(k)-Bi(k));
end
temp=rho_b(i,j)*rho_r(i,j)/rho(i,j)/rho(i,j);
%separation of two fluids distribution
f_r(9,i,j)=(rho_r(i,j)/rho(i,j))*ff(9,i,j);
f_b(9,i,j)=(rho_b(i,j)/rho(i,j))*ff(9,i,j);
for k=1:8
feq(k)=w(k)*rho(i,j);
f_r(k,i,j)=(rho_r(i,j)/rho(i,j))*ff(k,i,j)...
+beta*temp*feq(k)*coslai(k);
f_b(k,i,j)=(rho_b(i,j)/rho(i,j))*ff(k,i,j)...
-beta*temp*feq(k)*coslai(k);
end
end
end
end
end
%bounceback (no-slip boundary ops) step
for i=1:Nx
for j=1:Ny
if bound(i,j)==1
temp=ff(1,i,j);
ff(1,i,j)=ff(3,i,j);
ff(3,i,j)=temp;
temp=ff(2,i,j);
ff(2,i,j)=ff(4,i,j);
ff(4,i,j)=temp;
temp=ff(5,i,j);
ff(5,i,j)=ff(7,i,j);
ff(7,i,j)=temp;
temp=ff(6,i,j);
ff(6,i,j)=ff(8,i,j);
ff(8,i,j)=temp;
end
end
end
end
T=toc;
disp([‘Elapsed Running time = ’ num2str(T) ’ seconds’]);
X=1:Ny;
plot(X,u_x(1,:));
figure
plot(X,u_y(1,:));
%functions
function fequ=get_feq(rho, ux, uy, alfa, fequ)
global cx cy cs w
u_n=zeros(9,1);
u_squ=ux.^2+uy.^2;
tmp1=0.5*u_squ/cs;
for k=1:8
u_n(k)=cx(k)*ux+cy(k)*uy;
end
fequ(9)=w(9)rho((u_n(9)/cs)+(0.5u_n(9)u_n(9)/cs^2)-tmp1)+rhoalfa;
for k=1:4
fequ(k)=w(k)rho((u_n(k)/cs)+(0.5u_n(k)u_n(k)/cs^2)-tmp1)…
+rho(1-alfa)/5;
end
for k=5:8
fequ(k)=w(k)rho((u_n(k)/cs)+(0.5*u_n(k)u_n(k)/cs^2)-tmp1)…
+rho(1-alfa)/20;
end
end
function ff=stream(ff)
global Nx Ny
fprop=zeros(9,Nx,Ny);
for y=1:Ny
for x=1:Nx
% Streaming step (Periodic streaming of whole domain)
y_n=mod(y,Ny)+1;
x_e=mod(x,Nx)+1;
y_s=Ny-mod(Ny+1-y, Ny);
x_w=Nx-mod(Nx+1-x, Nx);
fprop(1,x_e,y)=ff(1,x,y);
fprop(2,x,y_n)=ff(2,x,y);
fprop(3,x_w,y)=ff(3,x,y);
fprop(4,x,y_s)=ff(4,x,y);
fprop(5,x_e,y_n)=ff(5,x,y);
fprop(6,x_w,y_n)=ff(6,x,y);
fprop(7,x_w,y_s)=ff(7,x,y);
fprop(8,x_e,y_s)=ff(8,x,y);
end
end
for y=1:Ny
for x=1:Nx
for k=1:8
ff(k,x,y)=fprop(k,x,y);
end
end
end
end