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, "(alphaT’’(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)
])