# Solve simple heat diffusion problem using LBM

Hi,
I tried to solve simple diffusion equaiton: alphaT’’(x) + x^2 = 0; T(0)=1;T(1) = 0 using D1Q3 but
I found I have to divide every term by a “Tune” factor to get the correct solution.
Although the two equations, "(alpha
T’’(x) + x^2 = 0)/Tune" and “alpha*T’’(x) + x^2 = 0”, should
have the same solution, that’s not the case in my code. Tune = 1; Tune = 4; Tune = 32 converge
to different solutions. When Tune >= 32, I can get the correct solution but I don’t know why it behaved like this.
Can anyone give me some suggestions ?

Thanks,
Min

## The code is attached below:


% solve heat diffusion equation : alpha*T''(x) + x^2 = 0; T(0)=1;T(1) = 0
% should have the same solution for (alpha*T''(x) + x^2 = 0)/Tune
clc;clear;
close all
NGrid = 101;                   % Grid number
n = NGrid;
x_exact = linspace(0,1,n);
Q = 3;                         % Number of Directions D1Q3
velocity = 0.0;                %  flow velocity
% --------------Initialization------------------------------------------%

pre_T(1:NGrid) = 0.0;
T(1:NGrid) = 0.0;

fT(1:n,1:Q) = 0.0;
T_L_wall = 1.0;         % T Left wall temperature
T_R_wall = 0.0;         % T Right wall temperature

wT_12 = 1/6;                  % for D1Q3
wT_0  = 4/6;                  % for D1Q3

% Initialize convergence criterion
sign = 0;                       % Exit Flag
e1 = 1E-10;                     % Convergence Criteria
mstep = 4E5;
% Tune = 1;   % won't work
Tune = 4;     % won't work
% Tune = 32;  % works !
alpha = 0.25;
alpha_lbm = alpha/Tune;
dx = 1/n;
dt = 1/n;
v = dx/dt;
cs_square  = 1.0/3*v^2; % for D1Q3
tau        = 1/2 + alpha_lbm/(cs_square*dt) ;

err1 = 1.0;
counter = 0;
while (sign==0)
counter = counter + 1 ;
pre_T = fT(:,1)+fT(:,2)+fT(:,3);
% -------  collision process D1Q3-------------------------------------------:
for i=1:n
fTeq(i,1)= wT_12 * T(i) * (1.0 + 1*velocity/cs_square);
fTeq(i,2)= wT_12 * T(i) * (1.0 - 1*velocity/cs_square);
fTeq(i,3)= wT_0  * T(i);

source =  x_exact(i)^2/Tune;

for k = 1:Q
if k ~=3
fT(i,k)=fT(i,k) + 1/tau*( fTeq(i,k)-fT(i,k) ) + dt*wT_12*source;
end

if k == 3
fT(i,k)=fT(i,k) + 1/tau*( fTeq(i,k)-fT(i,k) ) + dt*wT_0*source;
end

end
end
% -------  collision process D1Q3-------------------------------------------

% STREAMING
%===========================================================================%
for i=2:n
fT(n-i+2,1)=fT(n-i+1,1);
end

for i=1:n-1
fT(i,2)=fT(i+1,2);
end

% BOUNDARY CONDITIONS
%===============================================================================%
fT(1,1) = T_L_wall - fT(1,2) - fT(1,3); % for D1Q3
fT(n,2) = T_R_wall - fT(n,1) - fT(n,3); % for D1Q3

% CONVERGENCE
%===========================================================================%
% 4. calculate Macroscopic quantity (Non-dim_temperature)
T = fT(:,1)+fT(:,2)+fT(:,3);    % temperature for D1Q3
pre = sum(pre_T);
now = sum(T);
err1=abs((pre-now)/now);

if err1<e1
sign=1;
end

if counter>= mstep
break
end
end
T_exact = 1 - 0.666667*x_exact - 0.333333*x_exact.^4; % heat source = x^2
plot(x_exact,T,'--',x_exact,T_exact,'-')
legend('LBM','Exact',0);xlabel('x');ylabel('Temperature')
title(['heat source = 1/Tune*x^2' ', TL = ' num2str(T_L_wall) ', TR = ' num2str(T_R_wall) ', \tau = ' num2str(tau,'%3.2f')  10 ...
'\alpha'' = \alpha / Tune = ' num2str(alpha_lbm) ', ' ' \alpha  = ' num2str(alpha_lbm*Tune) ',' ' NGrid = ' num2str(n) ', Tune = ' num2str(Tune)
])