# large density ratio in Shan Chen by using P-R EOS

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)

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

Hello
I suggest changing this line:
psi = (2.*(p - Cs^2.rho) ./ C0 ./ G ).^0.5;
with this:
psi = (6.
(p - 1.*rho) ./ G ).^0.5;
and then decrease your initial values for densities. I mean 8 and .01 to about 3 and .5. Generally, it depends on your temperature.

Best regards

Dear hello_yjl_1,

could you solve your problem ?

Best regards
Markus