Dear every one, I’m now simulating Multiphase (water and vapor) problem. I’m using P-R(Peng-Robinson) EOS, all parameter are defined with reference in (Equation of state in a lattice Boltzmann model, Peng Yuan 2006). The problem is in the fifth or sixth loop, the density will be negative and the simulation will be fail and density plot cannot be finished. Anyone would help me to check the problem ? It’s quite urgent, thanks and really appreciate in advance! The code is in Matlab
clear
clc
global yDim; global xDim;
xDim = 50;
yDim = 50;
G = -1;
omega = 1;
omega_psi = 0.344; % it is used to calculate psi in P-R EOS
tau = 1/omega;
T = 0.018; % T and Tc are used to compute psi,in SC model, T is 1/3
Tc = 0.03;
maxT = 300;
tPlot = 1; % how many iteration it draw a plot
% allocate space
rho(yDim,xDim) = zeros;
uux(9,yDim,xDim) = zeros;
uuy(9,yDim,xDim) = zeros;
f(yDim,xDim,9) =zeros;
F(yDim,xDim,2) = zeros;
u_eq(yDim,xDim,2) = zeros;
f_eq(yDim,xDim,9) = zeros;
% D2Q9 lattice constants
t = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
v(1:9,1) = [ 0, 1, 0, -1, 0, 1, -1, -1, 1];
v(1:9,2) = [ 0, 0, 1, 0, -1, 1, 1, -1, -1];
% initialization of density (a bubble in the center)
radius = yDim/5;
for i = 1:xDim
for j = 1:yDim
if ((i-xDim/2)^2 + (j-yDim/2)^2) ^ 0.5 <= radius
rho(j,i) = 8; % density inside of bubble
else
rho(j,i) = 0.01; % denisty outside around the bubble
end
end
end
% initialization of distribution function
for i=1:9
f(:,:,i) = t(i).*rho;
end
% main loop
for cycle = 1:maxT
cycle
% Macroscopic properties
for x=1:xDim
for y=1:yDim
rho(y,x) = f(y,x,1) + f(y,x,2) +f(y,x,3) + f(y,x,4)+ f(y,x,5)...
+ f(y,x,6) + f(y,x,7) +f(y,x,8) + f(y,x,9);
end
end
[alpha,p,psi] = func_PR_compute_psi(T,Tc,omega_psi,rho,G);
% uDot calculation
for x = 1: xDim
for y = 1: yDim
for i = 1:9
uux(i,y,x) = f(y,x,i) .* v(i,1);
uuy(i,y,x) = f(y,x,i) .* v(i,2);
end
end
end
uux_sum = sum(uux);
uuy_sum = sum(uuy);
uux_sum2(:,:) = uux_sum(1,:,:); % convert 3D uux_sum to 2D uux_sum2 matrix
uuy_sum2(:,:) = uuy_sum(1,:,:);
uDot(:,:,1) = uux_sum2 ./ rho;
uDot(:,:,2) = uuy_sum2 ./ rho;
% stream of one phase (stream of density)
[Ftemp] = func_SCMP_stream_density(psi,t,v);
% force
for x = 1:xDim
for y = 1:yDim
F(y,x,1) = -G .* psi(y,x) .* Ftemp(y,x,1);
F(y,x,2) = -G .* psi(y,x) .* Ftemp(y,x,2);
end
end
% equilibrium velocity
for x= 1:xDim
for y=1:yDim
u_eq(y,x,1) = uDot(y,x,1) - F(y,x,1).*tau./rho(y,x);
u_eq(y,x,2) = uDot(y,x,2) - F(y,x,2).*tau./rho(y,x);
end
end
u = u_eq;
% uSqr and feq
for x = 1:xDim
for y = 1: yDim
u_Sqr(y,x) = u(y,x,1).*u(y,x,1) + u(y,x,2).*u(y,x,2);
end
end
for i=1:9
for x =1:xDim
for y=1:yDim
u_xy(y,x,i) = u(y,x,1)*v(i,1) + u(y,x,2)*v(i,2);
f_eq(y,x,i) = t(i)*rho(y,x)*(1 + 3*u_xy(y,x,i) + 4.5*u_xy(y,x,i)*u_xy(y,x,i) - 1.5*u_Sqr(y,x));
end
end
end
% collide
for i=1:9 % collision only happend without obstacle
for x=1:xDim
for y=1:yDim
f(y,x,i) = (1-omega) * f(y,x,i) + omega * f_eq(y,x,i);
end
end
end
% streaming
[f] = func_SCMP_stream_periodic(f);
% visulization
if(mod(cycle,tPlot)==0)
imagesc(rho); colorbar
title('Heavy phase (liquid) density');
axis equal off; drawnow
end
end
function [alpha,p,psi] = func_PR_compute_psi(T,Tc,omega,rho,G)
Cs = (1/3)^0.5;
C0 = 6.0;
R = 1;
a = 2/49;
b = 2/21;
alpha = (1 + (0.37464 + 1.54226omega - 0.26992omega^2) * (1-(T/Tc)^(1/2)) )^2;
p = rho.* R .* T ./ (1-b.*rho) - a.*alpha.*rho.^2./ (1+2.*b.*rho - b.^2.*rho.^2);
psi = (2.*(p - Cs^2.*rho) ./ C0 ./ G ).^0.5;
end
function [Ftemp] = func_SCMP_stream_density(psi,t,v)
global xDim; global yDim;
Ftemp(yDim,xDim,2) = zeros;
for x=1:xDim
if x > 1
xn = x-1; else xn = xDim; end
if x < xDim
xp = x+1; else xp = 1; end
for y=1:yDim
if y > 1
yn = y-1; else yn = yDim; end
if y < yDim
yp = y+1; else yp = 1; end
Fx_temp = t(2).* v(2,1).* psi(y,xp) + t(3).* v(3,1).* psi(yp,x)+...
t(4).* v(4,1).* psi(y,xn) + t(5).* v(5,1).* psi(yn,x)+...
t(6).* v(6,1).* psi(yp,xp) + t(7).* v(7,1).* psi(yp,xn)+...
t(8).* v(8,1).* psi(yn,xn) + t(9).* v(9,1).* psi(yn,xp);
Fy_temp = t(2).* v(2,2).* psi(y,xp) + t(3).* v(3,2).* psi(yp,x)+...
t(4).* v(4,2).* psi(y,xn) + t(5).* v(5,2).* psi(yn,x)+...
t(6).* v(6,2).* psi(yp,xp) + t(7).* v(7,2).* psi(yp,xn)+...
t(8).* v(8,2).* psi(yn,xn) + t(9).* v(9,2).* psi(yn,xp);
Ftemp(y,x,1) = Fx_temp;
Ftemp(y,x,2) = Fy_temp;
end
end
function [fout] = func_SCMP_stream_periodic(fin)
global xDim; global yDim;
fout(yDim,xDim,9) = zeros;
for x=1:xDim
if x > 1
xn = x-1;
else
xn = xDim; end
if x < xDim
xp = x+1;
else
xp = 1; end
for y=1:yDim
if y > 1
yn = y-1;
else
yn = yDim; end
if y < yDim
yp = y+1;
else
yp = 1; end
fout(y,x,1) = fin(y,x,1);
fout(y,xp,2) = fin(y,x,2);
fout(yp,x,3) = fin(y,x,3);
fout(y,xn,4) = fin(y,x,4);
fout(yn,x,5) = fin(y,x,5);
fout(yp,xp,6) = fin(y,x,6);
fout(yp,xn,7) = fin(y,x,7);
fout(yn,xn,8) = fin(y,x,8);
fout(yn,xp,9) = fin(y,x,9);
end
end