对流扩散方程

了解了一维扩散方程的求解原理后，就可以逐渐将其推广至二维、三维以及流动问题。首先加入最简单的流动，本文所指的对流是流体粒子整体运动（平流），即所有粒子的宏观运动速度都相同。

在笛卡尔坐标系中，一维对流扩散方程可以表示为

这个方程描述的是对流和扩散同时发生的过程，比如在恒定速度的水流中滴一滴墨水，墨水浓度随时间变化的情况。

格子Boltzmann法求解

对流扩散方程的格子玻尔兹曼方程和扩散方程的相同，即

唯一的不同是在平衡分布函数上增加了表示对流效应的外力项

其中，$ \vec{u} $为对流速度矢量，对于二维问题，$ \vec{u}=ui+vj $，i和j分别为x方向和y方向上的单位矢量，$ c_k $为沿着流动方向的单位矢量，$ c_k = \frac{\Delta x}{\Delta t}i+\frac{\Delta y}{\Delta t}j $。对于D2Q9模型，当$c_k=1$时，声速$c_s$为$ \frac{1}{\sqrt{3}} $。

举例

100*100平板上，仅有底部温度保持为1，其他三边均为0。速度分别为u=0.1和v=0.4，介质的热扩散系数为1，模拟时间步长为400。

matlab代码实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
% D2Q9tem
clear all
n = 100;                % 计算区域长度
m = 100;                % 计算区域宽度
f = zeros(9, n, m);     % 分布函数初始化
feq = zeros(1, 9);      % 平衡分布函数
rho = zeros(n, m);      % 密度

u = 0.1;                % 向上速度分量
v = 0.4;                % 向右速度分量
dt = 1;                 % 时间步长
dx = 1;                 % 空间步长
dy = dx;

alpha = 1;              % 扩散系数
ck = dx/dt;             % 格子运动速度
csq = ck*ck;

omega = 1.0/(3.*alpha/(csq*dt)+0.5);        % 松弛频率
mstep = 400;                                % 迭代次数

w = [4/9.0, 1/9.0, 1/9.0, 1/9.0, 1/9.0, 1/36.0, 1/36.0, 1/36.0, 1/36.0];
density = 0;                                % 初始值
tw = 1;                 % 上边界温度

cu = [0, 1, 0, -1, 0, 1, -1, -1, 1];        % 速度方向x
cv = [0, 0, 1, 0, -1, 1, 1, -1, -1];        % 速度方向y

% 初始化分布函数
for j = 1:m
    for i = 1:n
        for k = 1:9
            f(k,i,j) = w(k) * density;
            if i==0 % 上边界保持恒温
                f(k, i, j) = w(k) * tw;
            end
        end
    end
end

% 主循环
for kk = 1:mstep
    
    % 计算rho
    for j = 1:m
        for i = 1:n
            sum = 0;
            for k = 1:9
                sum = sum + f(k, i, j);
            end
            rho(i,j) = sum;
        end
    end
    
    % collision
    for j = 1:m
        for i = 1:n
            for k = 1:9
                feq(k) = w(k)*rho(i,j)*(1+3*(cu(k)*u + cv(k)*v) / ck);
                f(k, i, j) = omega * feq(k) + (1 - omega)*f(k,i,j);
            end
        end
    end
    
    % streaming
    for j = m:-1:2
        for i = 1:n-1
            f(3, i, j) = f(3, i, j-1);
            f(7, i, j) = f(7, i+1, j-1);
        end
    end
    for j = m:-1:2
        for i = n:-1:2
            f(2, i, j) = f(2, i-1, j);
            f(6, i, j) = f(6, i-1, j-1);
        end
    end  
    for j = 1:m-1
        for i = n:-1:2
            f(5, i, j) = f(5, i, j+1);
            f(9, i, j) = f(9, i-1, j+1);
        end
    end 
    for j = 1:m-1
        for i = 1:n-1
            f(4, i, j) = f(4, i+1, j);
            f(8, i, j) = f(8, i+1, j+1);
        end
    end 
    
    % boundary
    % bottom
    for j = 1:m
        f(2, 1, j) = w(2)*tw + w(4)*tw - f(4, 1, j);
        f(6, 1, j) = w(6)*tw + w(8)*tw - f(8, 1, j);
        f(9, 1, j) = w(9)*tw + w(7)*tw - f(7, 1, j);
    end
    % top
    for j = 1:m
        f(7, n, j) = - f(9, n, j);
        f(4, n, j) = - f(2, n, j);
        f(8, n, j) = - f(6, n, j);
        f(3, n, j) = - f(5, n, j);
        f(1, n, j) = 0;
    end
    % right
    for i = 1:n
        f(9, i, m) = - f(7, i, m);
        f(8, i, m) = - f(6, i, m);
        f(5, i, m) = - f(3, i, m);
        f(2, i, m) = - f(4, i, m);
        f(1, i, m) = 0;
    end
    % left
    for i = 1:n
        f(3, i, 1) = - f(5, i, 1);
        f(7, i, 1) = - f(9, i, 1);
        f(6, i, 1) = - f(8, i, 1);
        f(2, i, 1) = - f(4, i, 1);
        f(1, i, 1) = 0;
    end
end

for j = 1:m
    for i = 1:n
        sum = 0;
        for k = 1:9
            sum = sum+f(k,i,j);
        end
        rho(i,j) = sum;
    end
end

结果图