According to page 36 of the Book “LMB” by Mohamad, I have coded the example with Matlab. The example considers heat equation in a 1D bar with the length of 100. At t=0 (initial condition), temperature is set to zero. The left hand boundary has the temperature of 1 (boundary condition). The code simulates it successfully":
clear
clc
m=100;
f0=zeros(1,m+1);
f1=zeros(1,m+1);
f2=zeros(1,m+1);
tempr=zeros(1,m+1);
feq=zeros(1,m+1);
x=0:1:m;
dt=1;
dx=1;
omega=1.33333333;
mstep=200;
twall=1;
for kk=1:mstep
%collision process
for ii=1:(m+1)
tempr(ii)=f1(ii)+f2(ii);
feq(ii)=0.5*tempr(ii); %feq1=feq2
f1(ii)=(1-omega)f1(ii)+omegafeq(ii);
f2(ii)=(1-omega)f2(ii)+omegafeq(ii);
end
%streaming process
for jj=1:(m-1)
f1(m-jj+1)=f1(m-jj-1+1);
f2(jj-1+1)=f2(jj+1);
end
%Boundary Condition
f1(1)=twall-f2(1);
plot (x,tempr)
pause (0.03)
end
But lets change it a little. The initial (t=0) temprature follows the rule: T(x,0)=-0.024x^2+2.4x. Both left and right sides, have temprature of zero. But results at the right position (near x=100) are not correct. Please HELP. The code is:
clear
clc
m=100;
x=0:1:m;
tempr=-0.024x.^2+2.4x;
f0=zeros(1,m+1);
f1=0.5tempr;
f2=0.5tempr;
feq=0.5*tempr;
dt=1;
dx=1;
omega=1.333333333;
mstep=200;
twall=0;
for kk=1:mstep
%collision process
for ii=1:(m+1)
tempr(ii)=f1(ii)+f2(ii);
feq(ii)=0.5*tempr(ii); %feq1=feq2
f1(ii)=(1-omega)f1(ii)+omegafeq(ii);
f2(ii)=(1-omega)f2(ii)+omegafeq(ii);
end
%streaming process
for jj=1:(m-1)
f1(m-jj+1)=f1(m-jj-1+1);
f2(jj-1+1)=f2(jj+1);
end
%Boundary Condition
f1(1)=twall-f2(1);
f2(101)=twall-f1(101);
plot (x,tempr)
pause (0.01)
end