周报：2020.10.26 ~ 2020.11.1

2020.10.26 周一

已重构完SRT-LBM的顶盖驱动流C++代码（降低了一点常数复杂度），发现仍旧不能满足较大雷诺数的情况，查阅了相关文献后，发现使用MRT-LB模型（多松弛时间）似乎可以解决此问题，但仍旧较为复杂

研究了一晚上，发现高雷诺数的宏观大尺度流动模拟可太难了门槛和算力需求比较高，考虑临时更改模型为小尺度的自然对流流温耦合场模拟，这个应该会好做一些，而且可以应用在芯片散热这方面，目前的目标是先做一个简单的方腔自然对流流温耦合模型，大概就是何院士书中的如下案例：

image-20201026205644980

然而，关于边界条件和一些变量的初始化还是没弄明白，明天继续

2020.10.27

今天先熟悉了一下非平衡态外推格式，在顶盖驱动流内部加了两个小方块作为练习样例，如下图所示：

顶盖驱动 顶部方块

顶盖驱动 双方块

2020.10.28

完成了结题答辩PPT的理论基础板块制作，并且这两天成功实现了LBM的两个经典算例：Couette流和Poiseuille流，如下图所示：

  • Couette流：

    image-20201029235714276

Couette Flow 
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#include<bits/stdc++.h>
using namespace std;

const int Q = 9;
const int NX = 100;
const int NY = 100;
int e[Q][2] = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
double w[Q] = {4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36};
double rho[NX+1][NY+1], u[NX+1][NY+1][2], u0[NX+1][NY+1][2], f[NX+1][NY+1][Q], F[NX+1][NY+1][Q];
double U,Re,dx,dy,dt,Lx,Ly,rho0,tau,niu,c,cs;

inline double feq(int k,double rho,double u[2]){
    double res;
    double eu = e[k][0]*u[0]+e[k][1]*u[1];
    double uu = u[0]*u[0]+u[1]*u[1];
    res = w[k] * rho * (1 + eu/(cs*cs) + (eu*eu)/(2*cs*cs*cs*cs) - uu/(2*cs*cs));
    return res;
}

inline void init(){
    dx = 1.0;
    dy = 1.0;
    Lx = 1.0*dx*NX;
    Ly = 1.0*dy*NY;
    dt = 1.0;
    U = 0.1;
    rho0 = 1.0;
    Re = 1000;
    niu = U*Lx/Re;
    c = dx/dt;
    cs = c/sqrt(3.0);
    tau = niu/(cs*cs) + 0.5*dt;
    memset(u,0,sizeof(u));
    for(int i=0;i<=NX;i++){
        u[i][NY][0] = U; // lid speed = U
        for(int j=0;j<=NY;j++){
            rho[i][j] = rho0;
            for(int k=0;k<Q;k++){
                f[i][j][k] = feq(k,rho[i][j],u[i][j]);
            }
        }
    }
}

inline void evolution(){
    // collision & streaming
    for(int i=0;i<=NX;i++){
        for(int j=1;j<NY;j++){
            for(int k=0;k<Q;k++){
                int ip = i - e[k][0];
                int jp = j - e[k][1];
                // periodic boundary
                if(ip<0) ip=NX;
                if(ip>NX) ip=0;
                F[i][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp]) - f[ip][jp][k])/tau;
            }
        }
    }

    memcpy(u0,u,sizeof(u));
    memset(rho,0,sizeof(rho));
    memset(u,0,sizeof(u));

    // for(int j=1;j<NY;j++){
    //     for(int k=0;k<Q;k++){
    //         if(k==1||k==5||k==8){
    //             int ip = NX-e[k][0];
    //             int jp = j-e[k][1];
    //             F[NX][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp]) - f[ip][jp][k])/tau;
    //         }
    //         if(k==3||k==6||k==7){
    //             int ip = 0-e[k][0];
    //             int jp = j-e[k][1];
    //             F[0][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp]) - f[ip][jp][k])/tau;
    //         }
    //         if(k==1||k==5||k==8)
    //             F[0][j][k] = F[NX][j][k];
    //         if(k==3||k==6||k==7)
    //             F[NX][j][k] = F[0][j][k];
    //     }
    // }

    

    // get_macro_value
    for(int i=0;i<=NX;i++){
        for(int j=1;j<NY;j++){
            for(int k=0;k<Q;k++){
                f[i][j][k] = F[i][j][k];
                rho[i][j] += f[i][j][k];
                u[i][j][0] += e[k][0] * f[i][j][k];
                u[i][j][1] += e[k][1] * f[i][j][k];
            }
            u[i][j][0] /= rho[i][j];
            u[i][j][1] /= rho[i][j];
        }
    }

    // top & bottom boundary
    for(int i=0;i<=NX;i++){
        rho[i][0] = rho[i][1];
        rho[i][NY] = rho[i][NY-1];
        u[i][NY][0] = U;
        for(int k=0;k<Q;k++){
            f[i][0][k] = feq(k,rho[i][0],u[i][0]) + (1-1/tau) * (f[i][1][k] - feq(k,rho[i][1],u[i][1]));
            f[i][NY][k] = feq(k,rho[i][NY],u[i][NY]) + (1-1/tau) * (f[i][NY-1][k] - feq(k,rho[i][NY-1],u[i][NY-1]));
        }
    }
}

inline double ERR(){
    double tp1=0,tp2=0;
    for(int i=1;i<NX;i++){
        for(int j=1;j<NY;j++){
            tp1 += ((u[i][j][0]-u0[i][j][0])*(u[i][j][0]-u0[i][j][0])+(u[i][j][1]-u0[i][j][1])*(u[i][j][1]-u0[i][j][1]));
            tp2 += (u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]);
        }
    }
    tp1 = sqrt(tp1);
    tp2 = sqrt(tp2);
    return tp1/(tp2+1e-30);
}

inline void output(int num){
    ostringstream name;
    name<<"couette_flow"<<num<<".dat";
    ofstream out(name.str().c_str());
    out<<"Title= \"LBM Couette Flow\"\n"<<"VARIABLES=\"X\",\"Y\",\"U\",\"V\"\n"<<"ZONE T=\"BOX\",I="<<NX+1<<",J="<<NY+1<<",F=POINT"<<endl;
    for(int i=0;i<=NX;i++){
        for(int j=0;j<=NY;j++){
            out<<double(i)/Lx<<" "<<double(j)/Ly<<" "<<u[i][j][0]<<" "<<u[i][j][1]<<endl;
        }
    }
}

int main(){
    init();
    for(int i=1;;i++){
        evolution();
        if(i%100==0){
            double error = ERR();
            printf("%dth error = %e\n",i,error);
            if(i>=1000){
                if(i%1000==0) output(i);
                if(error<1.0e-5) break;
            }
        }
    }
}

  • Poiseuille流：

    image-20201029235758244

Poiseuille Flow 
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
#include<bits/stdc++.h>
using namespace std;

const int Q = 9;
const int NX = 100;
const int NY = 100;
int e[Q][2] = {{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
double w[Q] = {4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36};
double gamma[Q] = {0,1.0/3,1.0/3,1.0/3,1.0/3,1.0/12,1.0/12,1.0/12,1.0/12};
double rho[NX+1][NY+1], u[NX+1][NY+1][2], u0[NX+1][NY+1][2], f[NX+1][NY+1][Q], F[NX+1][NY+1][Q];
double U,Re,dx,dy,dt,Lx,Ly,rho0,tau,niu,c,cs,f_x,f_y;

struct rectangle
{
    int x1,y1,x2,y2;
}rec[1];

inline double feq(int k,double rho,double u[2]){
    double res;
    double eu = e[k][0]*u[0]+e[k][1]*u[1];
    double uu = u[0]*u[0]+u[1]*u[1];
    res = w[k] * rho * (1 + eu/(cs*cs) + (eu*eu)/(2*cs*cs*cs*cs) - uu/(2*cs*cs));
    return res;
}

inline bool inrec(int x,int y,int id){
    return x>=rec[id].x1&&x<=rec[id].x2&&y>=rec[id].y1&&y<=rec[id].y2;
}

inline void set_single_rec_geo(int id,int x1,int y1,int x2,int y2){
    rec[id] = rectangle{x1,y1,x2,y2};
}

inline void init(){
    dx = 1.0;
    dy = 1.0;
    Lx = 1.0*dx*NX;
    Ly = 1.0*dy*NY;
    dt = 1.0;
    rho0 = 1.0;
    f_x = 1e-5;
    f_y = 0;
    Re = 3000;
    niu = 0.01;
    c = dx/dt;
    cs = c/sqrt(3.0);
    tau = niu/(cs*cs) + 0.5*dt;

    set_single_rec_geo(0,30,30,60,60);

    memset(u,0,sizeof(u));
    for(int i=0;i<=NX;i++){
        for(int j=0;j<=NY;j++){
            rho[i][j] = rho0;
            for(int k=0;k<Q;k++){
                f[i][j][k] = feq(k,rho[i][j],u[i][j]);
            }
        }
    }
}

inline void evolution(){
    // collision & streaming
    for(int i=0;i<=NX;i++){
        for(int j=1;j<NY;j++){
            if(inrec(i,j,0)) continue;
            for(int k=0;k<Q;k++){
                int ip = i - e[k][0];
                int jp = j - e[k][1];
                // periodic boundary
                if(ip<0) ip=NX;
                if(ip>NX) ip=0;
                F[i][j][k] = f[ip][jp][k] + (feq(k,rho[ip][jp],u[ip][jp]) - f[ip][jp][k])/tau;
            }
        }
    }

    memcpy(u0,u,sizeof(u));
    memset(rho,0,sizeof(rho));
    memset(u,0,sizeof(u));

    // get_macro_value
    for(int i=0;i<=NX;i++){
        for(int j=1;j<NY;j++){
            if(inrec(i,j,0)) rho[i][j] = 1.0;
            else{
                for(int k=0;k<Q;k++){
                    f[i][j][k] = F[i][j][k] + gamma[k] * (f_x * e[k][0] + f_y * e[k][1]);
                    rho[i][j] += f[i][j][k];
                    u[i][j][0] += e[k][0] * f[i][j][k];
                    u[i][j][1] += e[k][1] * f[i][j][k];
                }
            }
            u[i][j][0] /= rho[i][j];
            u[i][j][1] /= rho[i][j];
        }
    }

    // top & bottom boundary
    for(int i=0;i<=NX;i++){
        rho[i][0] = rho[i][1];
        rho[i][NY] = rho[i][NY-1];
        for(int k=0;k<Q;k++){
            f[i][0][k] = feq(k,rho[i][0],u[i][0]) + (1-1/tau) * (f[i][1][k] - feq(k,rho[i][1],u[i][1]));
            f[i][NY][k] = feq(k,rho[i][NY],u[i][NY]) + (1-1/tau) * (f[i][NY-1][k] - feq(k,rho[i][NY-1],u[i][NY-1]));
        }
    }
}

inline double ERR(){
    double tp1=0,tp2=0;
    for(int i=1;i<NX;i++){
        for(int j=1;j<NY;j++){
            tp1 += ((u[i][j][0]-u0[i][j][0])*(u[i][j][0]-u0[i][j][0])+(u[i][j][1]-u0[i][j][1])*(u[i][j][1]-u0[i][j][1]));
            tp2 += (u[i][j][0]*u[i][j][0]+u[i][j][1]*u[i][j][1]);
        }
    }
    tp1 = sqrt(tp1);
    tp2 = sqrt(tp2);
    return tp1/(tp2+1e-30);
}

inline void output(int num){
    ostringstream name;
    name<<"poiseuille_flow_with_rec"<<num<<".dat";
    ofstream out(name.str().c_str());
    out<<"Title= \"LBM Couette Flow\"\n"<<"VARIABLES=\"X\",\"Y\",\"U\",\"V\"\n"<<"ZONE T=\"BOX\",I="<<NX+1<<",J="<<NY+1<<",F=POINT"<<endl;
    for(int i=0;i<=NX;i++){
        for(int j=0;j<=NY;j++){
            out<<double(i)/Lx<<" "<<double(j)/Ly<<" "<<u[i][j][0]<<" "<<u[i][j][1]<<endl;
        }
    }
}

int main(){
    init();
    for(int i=1;;i++){
        evolution();
        if(i%100==0){
            double error = ERR();
            printf("%dth error = %e\n",i,error);
            if(i>=1000){
                if(i%1000==0) output(i);
                if(error<1.0e-5) break;
            }
        }
    }
}

两个基本算例主要是熟悉两种边界条件处理格式：

  1. 周期性边界格式
  2. 压力边界格式

其中周期性边界处理只需要巧妙地利用模运算实现左壁和右壁边界分布函数的碰撞和迁移，即左壁面的右方向离散速度由右壁面确定，反之亦然。

2020.10.29

主要实现了对cpu芯片背后半导体排布的简单建模，大概就是模拟下面这个小方框里面的一些半导体块受横向流体作用时的流场速度分布

cpu

模拟结果如下图所示：

cpu_fx1e-5_25w lattice_2wdt_streamline

模拟使用的参数为：500*500的网格，格子运动粘度0.01，x正方向压力1e-5

可以看到在空当比较大的位置流体能很顺利地流过，而在中间空隙比较小的地方流体也可以通过，不过速度较小，类似于被”汽封“了；并且，在一些半导体块尾侧会出现卡门涡街，这属于可以预料到的情况。

遇到的问题

2020.10.30 ~ 2020.10.31

基本上LBM的东西算是告一段落了，我不太想再搞下去了，已经约定好周日去交大见蔡老师了，目前LBM阶段性任务算是肝完了，先放一放了，这两天临时抱佛脚再仔细看看NLOS吧。

2020.11.1

啊~总算是给蔡老师一个交代了，我以为我一个月没啥进度，结果老师说我进度还不错，嘿嘿，关键让我高兴的是，我不用再做这个NLOS的东西啦，老师给我讲述了一下俺们课题组近年来的研究方向，并给我发了33篇关于计算光谱学的文献，从今天起，我要过上每天看文献的日子啦，我还主动跟老师说每周什么时候汇报一下，初步定了是每周一晚上，通过腾讯会议的方式，nice！终于有适合我的活干了，干巴得！

项目记录
周报

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！