LBM气体吸附流动程序求教

zouguokafeiwu

LBM初学者，基于李庆老师气泡沸腾那篇文章的程序进行修改的，我研究的是超临界气体在狭缝孔隙的吸附流动，所以只是单一孔隙单相气体流动，应该是很简单的，但是我想要考虑壁面吸附，在加入Shan-Chen模型时能够产生正常的抛物线形的速度剖面，但是加上考虑壁面相互作用的的Fads后就很快出现错误，但经过最近的尝试，我发现是不是shanchen只能体现壁面密度有增加的现象，并不能真正模拟壁面的吸附作用呢？
另外，还有一个问题，就是下面这个子程序中的S( 2 ,  和 S( 3 ,   中在李庆老师的子程序中只考虑了流体间相互作用力，但是没有考虑浮力，所以在我现在这个程序中是否需要将它替换为总力呢？

global lx ly lxy rho cc c_squ g tau_e tau_t ux uy ex ey sigm obst;

G = -1;
w = single([4/3 1/3 1/3 1/3 1/3 1/12 1/12 1/12 1/12]);
S = single(zeros(9, lxy));
psx_re = zeros(9, lx, ly);

psx = sqrt(abs(2.*(p - rho.*c_squ)./G./cc./cc));
psx = reshape(psx, lx, ly);
for i = 1:9
       psx_re(i, : , = circshift(psx, [-ex(i), -ey(i)]);
end
psx = reshape(psx, 1, lxy);
psx_re = reshape(psx_re, 9, lxy);
Fmx = -G.*psx.*(w*(psx_re.*(ex'*ones(1, lxy))));
Fmy = -G.*psx.*(w*(psx_re.*(ey'*ones(1, lxy))));
Fmx(obst) = 0;
Fmy(obst) = 0;

Fbx = zeros(1, lxy);
Fby = (rho - mean(rho)).*(-g);

Fx = Fmx + Fbx;
Fy = Fmy + Fby;

ux     = (ex*reshape(ff, 9, lxy) + Fx/2)./rho;   
uy     = (ey*reshape(ff, 9, lxy) + Fy/2)./rho;
uF      = ux.*Fx + uy.*Fy;
Fm2     = Fmx.^2 + Fmy.^2;
psx2    = psx.^2;
S(1, :) = 0;
S(2, :) = 6*uF + sigm*Fm2./psx2./(tau_e - 0.5);
S(3, :) = -6*uF - sigm*Fm2./psx2./(tau_t - 0.5);
S(4, :) = Fx;
S(5, :) = -Fx;
S(6, :) = Fy;
S(7, :) = -Fy;
S(8, :) = 2*(ux.*Fx - uy.*Fy);
S(9, :) = ux.*Fy + uy.*Fx;

end