从第五章才开始有代码，而且代码在书最后，看起来有点麻烦
一维扩散方程为如下形式
那个可以是任何需要求解的东西，比如温度，压力。a对温度来说就是热传递参数，对动量传输来说就是动力学粘度。可以看到对于时间t只有一阶导，因为待求解的东西只取决于之前的时间。比如墨水在清水中扩散，某个点是否会有墨水，在空间上来说可能和这个点的上方和下方是否有墨水有关，但在时间上只与这个时间之前是否有墨水有关，而不可能取决于未来是否有墨水。
比如正宗的热传递方程
这时候a就等于 k / pC，于是简化为
写成机器能看懂的形式
然后开搞1阶LBM，我们看看D1Q2，也就是只计算这点左右两边的点，不包括中心点。首先是一些参数
% LBM 1D diffusion equation D1Q2 clear m = 101; dx = 1.0; rho = zeros(m);%这玩意就是温度T f1 = zeros(m); flux = zeros(m); x = zeros(m); x(1) = 0.0; for i = 1:m-1 x(i+1)=x(i)+dx; end alpha = 0.25; omega = 1/(alpha + 0.5); twall = 1.0;%左边界温度为1.0 nstep = 200; for i = 1:m f1(i) = 0.5*rho(i); f2(i) = 0.5*rho(i); end首先设置omega，由
得
LBM包含两个步骤，首先是Collision，然后是Streaming
%Collision for k1 = 1:nstep for i = 1:m feq = 0.5 * rho(i); f1(i) = (1-omega)*f1(i) + omega*feq; f2(i) = (1-omega)*f2(i) + omega*feq; end然后Streaming
%Streaming for i = 1:m-1 f1(m-i+1) = f1(m-i);% f1(m) = f1(m-1) , f1(m-1) = f1(m-2) ... f1(2) = f1(1) f2(i) = f2(i+1); % f2(1) = f2(3) , f2(2) = f2(3) ... f2(m-1) = f2(m) end边界条件
%Boundary condition f1(1) = twall - f2(1); f1(m) = f1(m-1); f2(m) = f2(m-1);然后由于Chapman–Enskog Expansion
for j = 1:m rho(j) = f1(j) + f2(j); end end画图
figure(1) plot(x,rho) title('Tempture') xlabel('X') ylabel('T')然后热流
%Flux for k = 1:m flux(k) = omega*(f1(k) - f2(k)); end figure(2) plot(x,flux,'o') title('Flux, time step = 200') xlabel('X') ylabel('Flux')然后把nstep设置成2，看看只迭代两次有多么不精确
然后开搞D1Q3
首先添加一些参数
w0 = 4./6.; w1 = 1./6.; w2 = w1; c2 = 1./3.; omega = 1/(3 * alpha + 0.5); f0 = zeros(m);然后Collision
feq0 = w0 * rho(i); feq = w1 * rho(i); f0(i) = (1-omega)*f1(0) + omega*feq0; f1(i) = (1-omega)*f1(i) + omega*feq; f2(i) = (1-omega)*f2(i) + omega*feq;边界条件
%Boundary condition f1(1) = twall - f2(1) - f0(1); f1(m) = f1(m-1); f2(m) = f2(m-1); f0(m) = f0(m-1);最终计算温度T
rho(j) = f1(j) + f2(j) + f0(j);flux也需要改一下，看看我的这个库吧，看看提交记录就知道了
https://github.com/clatterrr/cfdcode
D1Q3最终效果
然后看看D2Q5，也就是上下左右加中心点，没什么重要的变换，看看github上的代码即可
