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)
        ])