LBM with Computer Codes第五章学习笔记

从第五章才开始有代码，而且代码在书最后，看起来有点麻烦

一维扩散方程为如下形式

那个可以是任何需要求解的东西，比如温度，压力。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上的代码即可