Heat Equation

Lattice Boltzmann Method LBM: Algorithms
#1

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)+omega*feq(ii);
    f2(ii)=(1-omega)*f2(ii)+omega*feq(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)+omega*feq(ii);
    f2(ii)=(1-omega)*f2(ii)+omega*feq(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