声子BTE的离散坐标法

电热输运

声子BTE的离散坐标法

最近读了一篇动理学的文章，鲍华老师前两天开源的Gift-BTE也是采用确定性方法求解的声子BTE。确定性方法从程序实现上一个最大的好处就是，当把方程离散好了以后，整套求解流程是高度成熟的，可以利用开源的软件去建模、生成网格、做后处理，自己要写的内容主要就是一个求解器，而且其中的一些数值格式、迭代算法也可以借用很多成熟的方案，当然缺点就是占用的内存可能偏大，因为BTE的维度很高，最后得到的网格是很多个维度的乘积。对于蒙特卡洛模拟来说，整套求解逻辑是完全不同的，比如要模拟一个方块，确定性方法关注的是把这个方块离散成的一个个网格，而蒙特卡洛方法关注的是方块的六个边界对声子的反射过程，所以整套前后处理过程都需要自己亲自来写，当然它的好处就是可以处理复杂的界面，不需要考虑收敛性，因为不涉及网格生成的过程。同时它是用抽样来代替了实际对动量空间的离散，占用空间比较小。即使不懂BTE也可以用蒙特卡洛方法很快地模拟出一些结果，因为粒子的图像很清晰，就是很多的声子能量团在体系里面跑，哪里声子多，哪里就热一点。

这里以一维薄膜稳态导热过程为例，介绍一下离散坐标法（DOM）求解声子玻尔兹曼方程（BTE）的想法，大家还发展了很多先进的格式和算法，不过最基本的思路应该都差不多。

声子BTE方程组

弛豫时间近似（RTA）下稳态声子BTE可以写成 \[ \boldsymbol{v}_{k} \cdot \nabla_r f_k(\boldsymbol{r})= - \frac{f_k(\boldsymbol{r}) - f_k^{eq}(T_p(\boldsymbol{r}))}{\tau_k(T(\boldsymbol{r}))} \] 其中\(f_0\)是玻色-爱因斯坦分布， \[ f_k^{eq}(T_p) = \frac{1}{\exp(\hbar\omega_k/k_bT_p) - 1} \] 这个方程描述了波矢态为\(\vec{k}\)的声子的运动，\(f\)的含义就是该波矢态声子的数目分布，\(f\)再乘上\(\hbar\omega_k\)就是这个态的声子的能量。\(\vec{v}_k\)是这个波矢态的群速度，在\(\vec{k}\)确定后，群速度就固定了。波矢空间\(\vec{k}\)负责描述声子的性质，驱动声子在实空间\(\vec{r}\)中运动。其中\(T\)就是基于局域能量密度定义的温度，\(T_p\)是在散射项中为了确保能量守恒引入的假温度，它是通过能量守恒方程来定义的。

我们在方程左右两侧都乘上\(\hbar\omega_k\)，并引入参考温度，此时方程可以改写成， \[ \boldsymbol{v}_k \cdot \nabla_{\boldsymbol{x}} g_k(\boldsymbol{r})=\frac{g_k^{e q}(T_p(\boldsymbol{r}))-g_k(\boldsymbol{r})}{\tau_k(T(\boldsymbol{r}))}\quad\quad (1) \] 其中\(g_k = \hbar\omega_k \left[f_k - f_k^{eq}(T_\text{ref})\right], g_k^{eq} = \hbar\omega_k\left[f_k^{eq}(T(\vec{r})) - f_k^{eq}(T_\text{ref})\right]\)，这时候\(g_k\)就代表相对于参考温度的偏差能量了。我们进一步认为平衡分布函数之差可以解耦为\(g_k^{eq} = C_k (T(\vec{r}) - T_\text{ref})\)。采用分布函数描述的基本图像是，实空间中某一物理量的值，是所有波矢态的声子贡献叠加的总结果。所以局域能量密度的表达式就是 \[ u(\boldsymbol{r}) = \sum_\boldsymbol{k} g_k(\boldsymbol{r}) \] 局域温度\(T\)就是通过局域能量密度来定义的， \[ \sum_\boldsymbol{k}g_k(\boldsymbol{r}) = \sum_\boldsymbol{k}g_k^{eq}(\boldsymbol{r}) = \sum_{\boldsymbol{k}}C_k (T(\boldsymbol{r}) - T_\text{ref}) \] 于是 \[ T(\boldsymbol{r}) = \frac{\sum_\boldsymbol{k}g_k(\boldsymbol{r})}{\sum_{\boldsymbol{k}}C_k} + T_\text{ref} \quad\quad(2) \] 接下来我们获得\(T_p(\vec{r})\)的表达式。热流的表达式可以写成， \[ \boldsymbol{q}(\boldsymbol{r}) = \sum_{\boldsymbol{k}}\boldsymbol{v_k}g_k(\boldsymbol{r}) \] 如果没有内热源，根据能量守恒，热流的散度为0， \[ \nabla_r \cdot \boldsymbol{q}(\boldsymbol{r}) = \sum_k\left[g_k(\boldsymbol{r})\nabla_r \cdot \boldsymbol{v}_k +\boldsymbol{v}_k\cdot \nabla_r g_k(\boldsymbol{r}) \right] \] \(\vec{v}_k\)与空间坐标无关，因为BTE描述的是某一个\(\vec{k}\)的声子的运动方程，当\(\vec{k}\)固定后，\(\vec{v_k}\)就是一个常数了，它对空间坐标的偏导数为0，此时热流的散度变成 \[ \nabla_r \cdot \boldsymbol{q}(\boldsymbol{r}) = \sum_k\boldsymbol{v}_k\cdot \nabla_r g_k(\boldsymbol{r}) \] 我们再回去看看式 (1)，这实际上就是BTE左右两侧同时对波矢态求和， \[ \sum_{\boldsymbol{k}}\boldsymbol{v}_k \cdot \nabla_{\boldsymbol{x}} g_k(\boldsymbol{r})=\sum_{\boldsymbol{k}}\frac{g_k^{e q}(\omega_k, T_p(\boldsymbol{r}))-g_k(\boldsymbol{r})}{\tau_k(T(\boldsymbol{r}))} =\sum_{\boldsymbol{k}}\frac{C_k (T_p(\boldsymbol{r}) - T_\text{ref})-g_k(\boldsymbol{r})}{\tau_k(T(\boldsymbol{r}))} = 0 \] 于是 \[ T_p = \frac{\sum_{\boldsymbol{k}}g_k(\boldsymbol{r}) / \tau_k(T(\boldsymbol{r}))}{\sum_{\boldsymbol{k}}C_k / \tau_k(T(\boldsymbol{r}))} + T_\text{ref}\quad\quad (3) \] 所以我们就得到了采用确定性方法求解BTE的基本流程，首先，我们把声子的波矢空间离散成\(N\)个波矢态，这时我们就得到了\(N\)个方程 (1)，通过边界条件，我们可以确定边界处的一些声子态的分布函数的值。比如对于等温的黑体发射边界，垂直于边界的声子分布就是边界温度所处的玻色-爱因斯坦分布。我们假设体系的初始温度分布，用这个温度分布，可以计算每个空间格点的平衡分布\(g_k^{eq}\)。有了\(g_k^{eq}\)和边界处的\(g_k\)，我们就可以用一些数值方法去把方程 (1)去做离散，迭代出一轮新的\(g_k\)，有了新的\(g_k\)，就可以用式 (2)和式 (3)去算出新的\(T\)和\(T_p\)，并更新出一轮新的\(g_k^{eq}\)。全部更新结束后，就可以继续将离散化的方程 (1)去做迭代，直到收敛。

离散坐标法

在离散坐标法里，我们求解的并不是\(g_k\)，而是\(g_k\vec{v}\)，能量通量，或者称为辐射强度，它最早是用来解辐射输运方程的，通过光子和声子的类比，我们也可以用它来求解声子的输运方程。热辐射和声子输运的一个典型区别是，光子在空气中传播几乎是没有色散的，不同频率光的群速度几乎一样，所以本质上和灰体近似下的声子输运是一样的。前一节已经提到过了，在求解声子BTE的时候，我们需要把声子的波矢空间进行离散，我们需要求解每个态\((k_x, k_y, k_z)\)的声子BTE，然后再通过平衡分布把它们关联起来。

img

而此时如果我们假设声子的性质是各向同性的，那么第一布里渊区里的波矢态就不需要用笛卡尔坐标\((k_x, k_y, k_z)\)来描述和离散，而是可以用球坐标\((k, \theta, \phi)\)。这时波矢大小\(k\)用于描述声子的性质，比如频率、比热容、群速度大小，\(k\)一样的声子性质完全一样，除了群速度方向不同，群速度方向则用\((\theta, \phi)\)来描述。这样在离散波矢态的时候，我们可以把波矢大小和立体角分开离散，这样会节约对波矢空间离散的成本，当然这里会出现一些坐标转换时的系数。进一步地，如果采用了灰体近似，那么声子的性质沿球的半径方向也是没有分布的，这时只对立体角进行离散就可以了。在各向同性的近似下，群速度方向和\((\theta, \phi)\)所指向的方向是一样的。但一定要清楚当使用DFT计算的全带声子性质的时候并不是这样的，我们引入球坐标、引入立体角，本质上还是为了描述确定某一个声子态。也就是说虽然\(k\)相同的声子态的性质相同，但是它们本质上并不是同一个声子态，它们的分布函数是不同的。

slide-16

对于薄膜的一维稳态导热问题，对于各向同性声子性质，方程 (1)可以写成， \[ v_k \mu\frac{\mathrm{d}g_{k, \mu}(x)}{dx} + \frac{g_{k, \mu}(x)v_k }{\tau_kv_k}=\frac{v_k g_k^{e q}(T_p(x))}{\tau_kv_k} \] 其中\(\mu = \cos\theta\)，此时某个声子态就可用\((k, \mu)\)来描述。\(g_{k, \mu}v_k\)就是能流密度，或者说是辐射强度，我们把它定义成\(I_{k, \mu}\)，此时方程可以写成 \[ \frac{\mathrm{d}I_{k, \mu} (x)}{\mathrm{d}x} + \frac{I_{k, \mu}(x)}{\tau_kv_k\mu} = \frac{I^{eq}_k(T_p(x))}{\tau_kv_k\mu} \] 采用迎风格式对对流项进行离散，让\(I_{k, \mu}(x)\)接收到来自上游方向边界辐射的信息， \[ \frac{I_{k, \mu} (x_i) - I_{k, \mu}(x_{i-1})}{\Delta x} + \frac{I_{k, \mu}(x_i)}{\tau_kv_k\mu} = \frac{I^{eq}_k(T_p(x_i))}{\tau_kv_k\mu} \] 于是可以得到\(I_{k, \mu}(x_i)\)的更新方式， \[ I_{k, \mu}(x_i) = \left[\frac{I^{eq}_k(T_p(x_i))\Delta x}{\tau_kv_k\mu} + I_{k, \mu}(x_{i-1})\right]/ (1 + \frac{\Delta x}{\tau_k v_k \mu}) \] 假设高温热沉在下面，低温热沉在上面，那么对于\(\mu>0\)的都按照上面这个式子去更新，对于\(\mu<0\)的都按照低温热沉的方向去更新。那么在获得了\(I_{k, \mu}(x_i)\)以后，怎么计算当地的平衡分布和温度呢？还是用式 (2)和式 (3)，只不过这时候由于我们的变量是球坐标下的辐射强度，所以表达式略有不同，到这里才到了离散坐标法。离散坐标法的基本思路是：假设辐射强度在某个立体角单元内不随方向变化，则整个立体角空间可被划分为多个控制角度，分别求解每个控制角度上的能流强度，最后利用加权求和得到当地总能流密度，\(\int I \mathrm{d}\Omega = \sum_{i=1}^N w_i I_i\)， \[ \sum_k\sum_{i=1}^{N}w_i I_{k, \mu}(x_i) / 4\pi \mathrm{d}k= \sum_k I_{k, \mu}^{eq}(x_i)\mathrm{d}k \] 其中\(I_{k, \mu} = g_{k, \mu} v_k\)其实两侧还是对所有的声子态求和，只不过这时候我们是把总的波矢态拆分成了波矢大小\(k\)和立体角\(\Omega\)，左侧里层求和是对立体角求和，外层是对波矢大小求和，加到一起就是对所有的波矢态求和。但是由于我们采用了离散坐标法，实际上我们相当于完成了一次对立体角的积分，所以除以\(4\pi\)之后才是实际得到的能流密度。如果我们只对\(\mu\)进行离散，那么除以2就可以了，因为我们没有对\(\phi\)方向的\(2\pi\)进行积分。通过上面的等式，我们就可以分别计算假温度和温度了，在积分累加的时候我们需要乘上这一段的微元大小，

温度： \[ \begin{aligned} &\sum_{i=1}^{N}\frac{w_i I_{k, \mu}(x_i) / 4\pi}{v_k} \mathrm{d}k = \frac{I_{k, \mu}^{eq}(x_i)}{v_k} \\ &4\pi\sum_k \frac{I_{k, \mu}^{eq}(x_i)}{v_k} \mathrm{d}k= (T(x_i) - T_\text{ref}) \int_k C_k \mathrm{d}k \end{aligned} \] 假温度： \[ \begin{aligned} &\sum_{i=1}^{N}\frac{w_i I_{k, \mu}(x_i) / 4\pi}{v_k\tau_k} \mathrm{d}k = \frac{I_{k, \mu}^{pesudo, eq}(x_i)}{v_k\tau_k} \\ &4\pi\sum_k \frac{I_{k, \mu}^{pesudo, eq}(x_i)}{v_k \tau_k} \mathrm{d}k=(T_p(x_i) - T_\text{ref}) \int_k C_k / \tau_k \mathrm{d}k \end{aligned} \] 当然，在灰体近似下就更加简单了，没有了对\(k\)求和的部分，整个BTE求解的过程就少掉了一个维度，同时假温度和温度是没有区别的。

一维薄膜稳态导热

这里分别给出灰体近似和各向同性色散条件下采用DOM求解一维稳态薄膜导热过程的程序，考虑色散时，声子的性质来自于https://github.com/jeanphilippeperaud/1D_KMC_matlab/tree/master，灰体近似下，声子的性质采用

image-20230815200221356

jojo

色散版本：

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
import os
import sys

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import constants as c
import numba as nb

# 薄膜两侧温度
T_h = 301
T_c = 299
T_ref = 300
T_init = 300

# 读取硅的声子性质
data = pd.read_csv('dataSi.txt', names=['frequency', 'DoS', 'velocity', 'frequency_width', 'relaxation_time', 'polarization'], sep=' ')
data['de_dT'] = (c.hbar * data.frequency / T_ref)**2 * np.exp(c.hbar * data.frequency / (c.k * T_ref)) / c.k / (np.exp(c.hbar * data.frequency / (c.k * T_ref)) - 1)**2
data['C'] = data.de_dT * data.DoS
C = np.sum(data.C * data.frequency_width)
C_div_tau = np.sum(data.C / data.relaxation_time * data.frequency_width)
k = np.sum(data.C * data.relaxation_time * data.velocity**2 * data.frequency_width / 3)

# 立体角离散
mu_size = 48
mu, w = np.polynomial.legendre.leggauss(mu_size)

# 空间离散
L = 100e-9
N = 101
L_grids = np.linspace(0, L, N)
dx = L / (N-1)

# 初始化
radiation_density = np.zeros((len(data), mu_size, N))
equilibrium_radiation_density = np.ones((len(data), N))
for x in range(N):
    equilibrium_radiation_density[:, x] = data.C * data.velocity * (T_init - T_ref) / (4 * np.pi) * data.frequency_width

temperature = np.ones(N) * T_ref
pesudo_temperature = np.ones(N) * T_ref

tau = data.relaxation_time.to_numpy()
v = data.velocity.to_numpy()
d_omega = data.frequency_width.to_numpy()
Ck = data.C.to_numpy()
number_of_k = len(data)

@nb.jit(nopython=True)
def dispersion_func(radiation_density, equilibrium_radiation_density, pesudo_temperature, temperature):

    for it in range(100):
        for k in range(number_of_k):
            for solid_angle in range(int(mu_size/2), mu_size):
                radiation_density[k, solid_angle, 0] = Ck[k] * (T_h - T_ref) * v[k] / (4 * np.pi)
                for x in range(1, N):
                    radiation_density[k, solid_angle, x] = (radiation_density[k, solid_angle, x-1] + equilibrium_radiation_density[k, x] * dx / (tau[k] * v[k] * mu[solid_angle])) / (1  + dx / (tau[k] * v[k] * mu[solid_angle])) 

            for solid_angle in range(int(mu_size/2)):
                radiation_density[k, solid_angle, -1] = Ck[k] * (T_c - T_ref) * v[k] / (4 * np.pi)
                for x in range(N-2, -1, -1):
                    radiation_density[k, solid_angle, x] = (- radiation_density[k, solid_angle, x+1] + equilibrium_radiation_density[k, x] * dx / (tau[k] * v[k] * mu[solid_angle])) / (-1 + dx / (tau[k] * v[k] * mu[solid_angle]))
        
        for x in range(N):
            pesudo_energy = 0
            energy = 0
            for k in range(number_of_k):
                pesudo_energy += np.sum(w * radiation_density[k, :, x]) / 2 * d_omega[k] / v[k] / tau[k]
                energy += np.sum(w * radiation_density[k, :, x]) / 2 * d_omega[k] / v[k]
            pesudo_temperature[x] = T_ref + pesudo_energy * np.pi * 4 / C_div_tau
            temperature[x] = T_ref + energy * np.pi * 4 / C
            equilibrium_radiation_density[:, x] = Ck * v * (pesudo_temperature[x] - T_ref) / (4 * np.pi)
        
dispersion_func(radiation_density, equilibrium_radiation_density, pesudo_temperature, temperature)

灰体版本：

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
import os
import sys

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import constants as c
import numba as nb

# 灰体性质
C = 690 * 2330
v = 6400
l = 43.7e-9
tau = l / v
k = 1/3 * C * v * l

# 薄膜两侧温度
T_h = 301
T_c = 299
T_ref = 300
T_init = 300

# 立体角离散
mu_size = 48
mu, w = np.polynomial.legendre.leggauss(mu_size)

# 空间离散
L = 50e-9
N = 11
L_grids = np.linspace(0, L, N)
dx = L / (N-1)

# 初始化
radiation_density = np.zeros((mu_size, N))
equilibrium_radiation_density = np.ones((N)) * C * v * (T_init - T_ref) / (4 * np.pi)

temperature = np.ones(N) * T_ref
pesudo_temperature = np.ones(N) * T_ref

@nb.jit(nopython=True)
def gray_func(radiation_density, equilibrium_radiation_density, pesudo_temperature):
    for i in range(1000):

        for solid_angle in range(int(mu_size/2), mu_size):
            radiation_density[solid_angle, 0] = C * (T_h - T_ref) * v / (4 * np.pi)
            for x in range(1, N):
                radiation_density[solid_angle, x] = (radiation_density[solid_angle, x-1] + equilibrium_radiation_density[x] * dx / (tau * v * mu[solid_angle])) / (1  + dx / (tau * v * mu[solid_angle])) 

        for solid_angle in range(int(mu_size/2)):
            radiation_density[solid_angle, -1] = C * (T_c - T_ref) * v / (4 * np.pi)
            for x in range(N-2, -1, -1):
                radiation_density[solid_angle, x] = (- radiation_density[solid_angle, x+1] + equilibrium_radiation_density[x] * dx / (tau * v * mu[solid_angle])) / (-1 + dx / (tau * v * mu[solid_angle]))

        for x in range(N):
            equilibrium_radiation_density[x] = np.sum(w * radiation_density[:, x]) / 2

        pesudo_temperature = T_ref + equilibrium_radiation_density * np.pi * 4 / C / v

p_t = gray_func(radiation_density, equilibrium_radiation_density, pesudo_temperature)