C语言改变封闭图形颜色,LBM的封闭方腔自然对流新手代码求修改

有谁用LBM做流体的，本人新手，写了一个双分布耦合的程序，结果有问题，请教高手改正,感谢！本人用matlab写的

function LBGK_D2Q9_closedcavityflow

% Lattice Boltzmann LBE, geometry: D2Q9, model: CLBGK

% Density distribution function, f, feq[I][J][N].

% Temperature distribution function, g, geq[I][J][N].

% Velocity field, u[I][J], v[I][J].

% Temperrature field, T[I][J].

% Space node: i,j (I,J),      Discrete velocity : n (1:9).

% Space step: dx (1/I)   dy (1/J)

% Alterable parameter: Pr,Ra,Ma

% Relaxation time, tau_f,tau_T

I=10;J=10;N=9;K=4;

dx=1;dy=1;dt=1;%c=1;

H=dx*I;H=dy*J;

T0=19;T1=1;Tm=(T0+T1)/2;P0=0;

Pr=0.71;Ma=0.1;

Ra=10000;

niu=H*Ma*(Pr/(3*Ra))^(1/2);

tau_f=3.0*niu+0.5;

chi=H*Ma/(3*Pr*Ra)^(1/2);

tau_T=2.0*chi+0.5;

e=[1,0;0,1;-1,0;0,-1;1,1;-1,1;-1,-1;1,-1;0,0];

%Initialization

P=P0*ones(I+1,J+1);

for j=1:J+1

for i=1:I+1

u(i,j)=0;v(i,j)=0;uu(i,j)=0;

T(i,j)=Tm;P(i,j)=P0;

for n=1:N

U=[u(i,j),v(i,j)];

p1=P(i,j);

f(i,j,n)=feq(U,p1,n);

end

for k=1:K

U=[u(i,j),v(i,j)];

Tt=T(i,j);

g(i,j,k)=geq(U,Tt,k);

end

end

T(1,j)=T0;

T(I+1,j)=T1;

end

% evolution

for iter=1:1000

for j=2:J

for i=2:I

for k=1:K

ip=i-e(k,1);

jp=j-e(k,2);

G(i,j,k)=g(ip,jp,k)+(geq([u(ip,jp),v(ip,jp)],T(ip,jp),k)-g(ip,jp,k))/tau_T;

end

for n=1:N

ip=i-e(n,1);

jp=j-e(n,2);

F(i,j,n)=f(ip,jp,n)+(feq([u(ip,jp),v(ip,jp)],P(ip,jp),n)-f(ip,jp,n))/tau_f+Force(T(ip,jp),n);

end

end

end

%计算宏观量

for i=2:I

for j=2:J

P(i,j)=0;

T(i,j)=0;

u(i,j)=0;v(i,j)=0;

for k=1:K

g(i,j,k)=G(i,j,k);

T(i,j)=T(i,j)+g(i,j,k);

end

for n=1:N

f(i,j,n)=F(i,j,n);

u(i,j)=e(n,1)*f(i,j,n)+u(i,j);

v(i,j)=e(n,2)*f(i,j,n)+v(i,j);

end

for n=1:N-1

P(i,j)=f(i,j,n)+P(i,j);

end

uu(i,j)=sqrt(u(i,j).^2+v(i,j).^2);

P(i,j)=0.6*P(i,j)-0.4*uu(i,j).*uu(i,j);

end

end

%边界处理

% 1. 左右外推式边界

for j=1:J+1

for n=1:N

P(I+1,j)=P(I,j);

f(I+1,j,n)=feq([u(I+1,j),v(I+1,j)],P(I,j),n)+f(I,j,n)-feq([u(I,j),v(I,j)],P(I,j),n);

P(1,j)=P(2,j);

f(1,j,n)=feq([u(1,j),v(1,j)],P(2,j),n)+f(2,j,n)-feq([u(2,j),v(2,j)],P(2,j),n);

end

for k=1:K

g(I+1,j,k)=geq([u(I+1,j),v(I+1,j)],T(I,j),k)+g(I,j,k)-geq([u(I,j),v(I,j)],T(I,j),k);

g(1,j,k)=geq([u(1,j),v(1,j)],T(2,j),k)+g(2,j,k)-geq([u(2,j),v(2,j)],T(2,j),k);

end

T(1,j)=T0;

T(I+1,j)=T1;

end

%上下外推式边界

for i=2:I

for n=1:N

P(i,1)=P(i,2);

f(i,1,n)=feq([u(i,1),v(i,1)],P(i,2),n)+f(i,2,n)-feq([u(i,2),v(i,2)],P(i,2),n);

P(i,J+1)=P(i,J);

f(i,J+1,n)=feq([u(i,J+1),v(i,J+1)],P(i,J),n)+f(i,J,n)-feq([u(i,J),v(i,J)],P(i,J),n);

end

for k=1:K

g(i,1,k)=geq([u(i,1),v(i,1)],T(i,2),k)+f(i,2,k)-geq([u(i,2),v(i,2)],T(i,2),k);

g(i,J+1,k)=geq([u(i,J+1),v(i,J+1)],T(i,J),k)+f(i,J,k)-geq([u(i,J),v(i,J)],T(i,J),k);

end

end

end

%输出图形

x=0:dx:H;

y=0:dy:H;

%   mesh(x,y,uu)

contour(y,x,uu)

%     quiver(x,y,u,v)

function result1 = feq(U,p1,n)

w=[1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36,4/9];

e=[1,0;0,1;-1,0;0,-1;1,1;-1,1;-1,-1;1,-1;0,0];

eu=e(n,1).*U(1)+e(n,2).*U(2);

uv=U(1)^2+U(2)^2;

if n==9

result1=w(n)*(3*eu+4.5*eu.*eu-1.5*uv)-5*p1/3;

elseif  n<5

result1=w(n)*(3*eu+4.5*eu.*eu-1.5*uv)+p1/3;

else

result1=w(n)*(3*eu+4.5*eu.*eu-1.5*uv)+p1/12;

end

function result2 = geq(U,Tt,n)

e=[1,0;0,1;-1,0;0,-1;1,1;-1,1;-1,-1;1,-1;0,0];

et=e(n,1).*U(1)+e(n,2).*U(2);

result2=Tt*(1+2*et)/4;

function result3 = Force(T,n)    %只受重力

e_F=[0,0;0,1;0,0;0,-1;0,0;0,0;0,0;0,0;0,0];

a=[0;-1];Tm=10;Ma=0.1;H=10;

ef=e_F(n,1).*a(1)+e_F(n,2).*a(2);

result3=-ef*Ma^2*(T-Tm)/(6*H*2*Tm);