系综声子蒙特卡洛

电热输运

1D Ensemble Phonon Monte Carlo

写一写感觉对MC的认识又加深了一些. 实际上Ensemble MC相比于Kinetic MC，如果不考虑Phonon性质随温度的变化，那它在模拟过程中没有引入任何更多的信息，除了比Tracing慢以外给出的结果应该是一样的. 在纯phonon模拟中哪些性质需要我们知道局域格点温度呢？实际上就是散射时的频率抽样分布函数 \[ \frac{D(\omega, p) \hbar \omega}{\tau(\omega, p, T_{\text{loc}})}\left[f(\omega,T_\text{loc} ) - f(\omega, T_\text{eq})\right] \] 这一项里有两个地方出现了\(T_\text{loc}\)，一个是弛豫时间，一个是后面的玻色爱因斯坦分布，但是一旦我们采用差分近似 \[ \frac{D(\omega, p) \hbar \omega}{\tau(\omega, p, T_{\text{loc}})}\left[f(\omega,T_\text{loc} ) - f(\omega, T_\text{eq})\right] \approx \frac{D(\omega, p) \hbar \omega}{\tau(\omega, p, T_{\text{eq}})}\frac{\mathrm{d}f(\omega)}{\mathrm{d}T}\Big|_{T_\text{eq}} (T_\text{eq} - T_\text{loc}) \] 就不再需要统计局域温度了，因为此时同一格点内的\(T_\text{loc}\)都是一样的，在抽样分布函数归一化后最后的温度项就没了. 所以一旦采用了这种近似，Ensemble MC和Kinetic MC主要的区别就是对散射过程的处理.

对于局域温度统计过程，在MC中严格的统计是使用玻色爱因斯坦分布计算局域能量密度或者局域Pesudo能量密度， \[ \begin{aligned} (N^+ - N^-)\mathcal{E}_{\mathrm{eff}}^{d}/V &= \sum_p \int_{0}^{\omega_m} \hbar\omega D(\omega, p) \times \left[ f(\omega, T_\text{loc}) - f(\omega, T_{eq})\right] \mathrm{d}\omega \\ &= \sum_p \int_{0}^{\omega_m} \hbar\omega D(\omega, p) \times \left[\frac{1}{\exp(\frac{\hbar\omega}{k_bT_\text{loc}})- 1} - \frac{1}{\exp(\frac{\hbar\omega}{k_bT_{eq}}) - 1} \right] \mathrm{d}\omega \end{aligned} \] 当然我们也可以做同样的近似， \[ \begin{aligned} (N^+ - N^-) \mathcal{E}_{\mathrm{eff}}^{d}/V &\approx (T_\text{eq} - T_\text{loc})\sum_p \int_{0}^{\omega_m} \hbar\omega D(\omega, p)\frac{\mathrm{d}f(\omega)}{\mathrm{d}T}\Big|_{T_\text{eq}} \mathrm{d}\omega \end{aligned} \] 把温度从积分里提出来，从而简化温度的统计方式，但是这种近似与模拟过程无关，只与对模拟结果的统计有关.

Python写的1D薄膜导热的Ensemble phonon MC程序，还是使用Si的色散，和之前写过的Kinetic MC是一样的. 这里不考虑散射率随温度的变化并且对散射中的频率抽样分布函数进行差分近似，所以定义的抽样数组是基本一样的，但是在这里采用了严格的方式进行能量-温度转换，构建look-up table来对能量和温度进行映射.

导入一些库

1
2
3
4
5
6
7
8
9
10
11
12
13
from dataclasses import dataclass

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

# 如果使用Jupyter可以开启下两行魔法命令，第三行是科研绘图风格的SciencePlot库
# %matplotlib inline
# %config InlineBackend.figure_format = 'svg'
# plt.style.use('science')

读取Si的色散

1
data = pd.read_csv('dataSi.txt', names=['frequency', 'DoS', 'velocity', 'frequency_width', 'relaxation_time', 'polarization'], sep=' ')
frequency DoS velocity frequency_width relaxation_time polarization
0 3.6875e+10 7.9627e+08 8253.5 7.375e+10 8.1497e-07 1
1 1.1063e+11 2.5122e+09 8250.7 7.375e+10 6.9282e-07 1
2 1.8438e+11 5.3781e+09 8247.8 7.375e+10 5.3303e-07 1
3 2.5813e+11 9.1975e+09 8245 7.375e+10 3.9602e-07 1
4 3.3188e+11 1.3909e+10 8242.1 7.375e+10 2.9494e-07 1
... ... ... ... ... ... ...

定义体系的基本参数

image-20220513135547791

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Teq = 300 # 参考平衡温度
Tl = 301 # 左侧高温热沉温度
Tr = 299 # 右侧低温热沉温度
L = 10e-7 # 体系的长度
Nx = 30 # 空间离散的cell数目
dt = 5e-12 # 时间步长
Ntt = 200 # 时间步数
Dx = L / Nx # cell的长度
xx = np.arange(Dx / 2, L + Dx / 2, Dx) # 各个cell的中心坐标
x_lines = np.arange(Dx, L+Dx, Dx) # cell的右侧边界线，这里主要是为了方便确定每个phonon所处的网格
Temperature = np.zeros((Ntt, Nx)) # 温度储存矩阵
PesudoTemperature = np.zeros((Ntt, Nx)) # 假温度储存矩阵
drift_projection = np.zeros(len(data)) # 统计一次随机行走各个phonon在x方向的投影长度，用于统计热流，并最终计算热导率，我们可以统计各个模态phonon的贡献

# 定义phonon对象，每个phonon由5个参数来描述，x为phonon的当前位置，next_x为phonon将要移动到的位置，这个量主要用于判断它是否撞到了边界. mu为phonon运动的方向，在1D问题里就是cos(theta). index为data数据的索引，用于确定当前phonon的性质. sign为phonon的符号，代表当前phonon相对平衡能量是正能量还是负能量 
@dataclass
class Phonon:
    x: float
    next_x: float
    mu: float
    index: int
    sign: int

构建能量和温度间的映射

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# Bose-enstein distribution
def f(data, T):
    return 1 / (np.exp(c.hbar * data.frequency / (c.k * T)) - 1)

# \int_0^{\omega_m} \hbar\omega D(\omega) f(T) - f(T_eq) \mathrm{d}\omega
def GetEnergyDeviation(data, T):
    return np.sum(c.hbar * data.frequency * data.DoS * (f(data, T) - f(data, Teq)) * data.frequency_width)

# \int_0^{\omega_m} \hbar\omega D(\omega)/\tau f(T) - f(T_eq) \mathrm{d}\omega
def GetPesudoEnergyDeviation(data, T):
    return np.sum(c.hbar * data.frequency * data.DoS * (f(data, T) - f(data, Teq)) / data.relaxation_time * data.frequency_width)

# 构建look-up table及插值函数的范围
T_range = np.linspace(298.8, 301.2, 30)

# 这两个产值函数分别把对应能量密度映射到温度上
EnergyToT = interp1d([GetEnergyDeviation(data, T) for T in T_range], T_range)
PesudoEnergyToT = interp1d([GetPesudoEnergyDeviation(data, T) for T in T_range], T_range)

构建抽样数组

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
de_dT = (c.hbar * data.frequency / Teq)**2 * np.exp(c.hbar * data.frequency / (c.k * Teq)) / c.k / (np.exp(c.hbar * data.frequency / (c.k * Teq)) - 1)**2 # hbar * omega * df/dT，注意这里多乘了一个hbar * omega，因为模拟的是声子能量团

base = data.DoS * de_dT * data.frequency_width
coll = base / data.relaxation_time
vel_hot = data.DoS * data.velocity * c.hbar * data.frequency * data.frequency_width * (f(data, Tl) - f(data, Teq))
vel_cold = data.DoS * data.velocity * c.hbar * data.frequency * data.frequency_width * (f(data, Teq) - f(data, Tr))

# 内部初始化过程的phonon抽样数组，在本程序中不存在，因为假设体系初始就处于平衡温度
cumsum_base = np.cumsum(base) / np.sum(base)
# 散射phonon抽样数组
cumsum_coll = np.cumsum(coll) / np.sum(coll)
# 左侧高温热沉发射抽样数组
cumsum_hot = np.cumsum(vel_hot) / np.sum(vel_hot)
# 右侧低温热沉发射抽样数组
cumsum_cold = np.cumsum(vel_cold) / np.sum(vel_cold)

确定一颗声子能量团代表的能量

1
2
3
4
5
6
7
# 这个有什么明确规定吗？在这里的程序中，假设体系都处于高温热沉的温度时，计算此时体系的总能量偏差，用20000个phonon能量团来近似这个总能量偏差，从而得到一颗phonon能量团对应的能量，也可以确定每一个时间步高低温热沉发射的phonon能量团数目了. 按照这个能量大小，每一个dt时间步内高温热沉要发射268颗phonon
N_phonon = 20000
Eeff = GetEnergyDeviation(data, Tl) * L / N_phonon
q_hot = np.sum(data.DoS * data.velocity * c.hbar * data.frequency * (f(data, Tl) - f(data, Teq)) * data.frequency_width) / 4 * dt
q_cold = - np.sum(data.DoS * data.velocity * c.hbar * data.frequency * (f(data, Tr) - f(data, Teq)) * data.frequency_width) / 4* dt
q_hot_emitted_per_dt = int(q_hot / Eeff)
q_cold_emitted_per_dt = int(q_cold / Eeff)

开始模拟

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
# 这里采取了累积概率函数抽样的形式，np.random.choice可以完成相同的需求
@nb.jit(nopython=True)
def rand_choice_nb(arr, cumsum_prob):
    """Sample an element in arr with the probability given in prob.
    """
    return arr[np.searchsorted(cumsum_prob, np.random.random(), side="right")]

# 在这里用phonons列表储存系统内模拟的所有phonon
phonons = []

# 在Kinetic MC中，外循环是对phonon循环，内循环才是对时间（自由程）循环，因为追踪不同phonon是独立的，在Ensemble MC中因为每一步都需要统计各个网格内的phonon分布，所以时间循环放在外面
for t in range(Ntt):
    ############# 内部声子运动 #############
    
    remove_indexes = []
    drift_projection = np.zeros(len(data))

    # Drift-process
    for i, phonon in enumerate(phonons):
        phonon.next_x = phonon.x + data.loc[phonon.index, "velocity"] * phonon.mu * dt
        # 判断该phonon是否撞到了热沉边界，在1D问题中用位置来确定就好了，在不规则几何内部可以用射线法进行碰撞检测
        if phonon.next_x < 0:
            phonon.next_x = 0
            remove_indexes.append(i)
        elif phonon.next_x > L:
            phonon.next_x = L
            remove_indexes.append(i)
        # 统计一次随机行走各个phonon在x方向的投影长度
        drift_projection[phonon.index] += phonon.sign * (phonon.next_x - phonon.x)

    # 删除掉撞到撞到热沉边界上的phonon，这里必须从后向前删除phonon，保证phonons数组中的index不变
    for index in remove_indexes[::-1]:
        phonons.pop(index)

    ############# 边界发射 #############

    # 从热边界发射
    for i in range(q_hot_emitted_per_dt):
        mu = np.sqrt(np.random.rand())
        phonon_index = rand_choice_nb(data.index.values, cumsum_hot.values)
        phonon = Phonon(0, 0, mu, phonon_index, 1)
        time = dt * np.random.rand()
        phonon.next_x = phonon.x + data.loc[phonon.index, "velocity"] * phonon.mu * time
        if phonon.next_x >= L:
            phonon.next_x = L
            drift_projection[phonon.index] += phonon.sign * (phonon.next_x - phonon.x)
        else:
            drift_projection[phonon.index] += phonon.sign * (phonon.next_x - phonon.x)
            phonons.append(phonon)

    # 从冷边界发射
    for i in range(q_cold_emitted_per_dt):
        mu = -np.sqrt(np.random.rand())
        phonon_index = rand_choice_nb(data.index.values, cumsum_cold.values)
        phonon = Phonon(L, 0, mu, phonon_index, -1)
        time = dt * np.random.rand()
        phonon.next_x = phonon.x + data.loc[phonon.index, "velocity"] * phonon.mu * time
        if phonon.next_x <= 0:
            phonon.next_x = 0
            drift_projection[phonon.index] += phonon.sign * (phonon.next_x - phonon.x)
        else:
            drift_projection[phonon.index] += phonon.sign * (phonon.next_x - phonon.x)
            phonons.append(phonon)

    # 更新phonon位置
    for phonon in phonons:
        phonon.x = phonon.next_x

    # 统计热流
    heat_flux = np.sum(drift_projection) / dt * Eeff / L

    ############# 统计这一步的结果 #############
    scatter_index_arrays = [[] for _ in range(len(xx))]  # 这个数组用于确定每个网格中是在phonons数组中index为哪些的phonon发生了散射
    Energy = np.zeros(Nx) # 能量储存矩阵
    PesudoEnergy = np.zeros(Nx) # 假能量储存矩阵
    for index, phonon in enumerate(phonons):
        # 判断当前phonon所处于的网格
        position_index = np.searchsorted(x_lines, phonon.x)
        Energy[position_index] += phonon.sign * Eeff / Dx # 因为是能量密度，所以除以网格尺寸Dx
        PesudoEnergy[position_index] += (
            phonon.sign * Eeff / data.loc[phonon.index, "relaxation_time"] / Dx
        )
        Temperature[t] = EnergyToT(Energy)
        PesudoTemperature[t] = PesudoEnergyToT(PesudoEnergy)
        # 顺便判断一下这个phonon是否发生了散射
        if np.random.rand() < 1 - np.exp(
            -dt / data.loc[phonon.index, "relaxation_time"]
        ):
            scatter_index_arrays[position_index].append(index)

    ############# 开始散射 #############
    # 这部分的逻辑可参考Peraud和Hadjiconstantinou的论文，核心就是假设当前网格内有N+个正phonon和N-个负phonon发生了散射（上面一步中进行了判断，判断结果储存在scatter_index_arrays中），最终程序中实际散射的phonon只有|N+ - N-|个，这些散射后的声子的符号由当前格点的温度决定，最终再删除其他选中的phonon
    remove_indexes = []
    for position_index, scatter_index_array in enumerate(scatter_index_arrays):
        net_scatter_phonon_number = abs(
            int(np.sum([phonons[i].sign for i in scatter_index_arrays[position_index]]))
        )
        phonon_sign = np.sign(Temperature[t, position_index] - Teq)
        # scatter_phonons里记录的是程序中实际发生散射的|N+ - N-|个声子在phonons数组中的index
        scattered_phonons = np.random.choice(
            scatter_index_array, net_scatter_phonon_number, replace=False
        )
        # 抽样发生散射后phonon的方向及性质
        for phonon_index in scattered_phonons:
            phonons[phonon_index] = Phonon(
                phonons[phonon_index].x,
                phonons[phonon_index].next_x,
                2 * np.random.rand() - 1,
                rand_choice_nb(data.index.values, cumsum_coll.values),
                phonon_sign,
            )
        # 统计要删除掉的声子，这里只能最后一起删，否则列表的index会出现问题
        remove_indexes.append(set(scatter_index_array) - set(scattered_phonons))
    remove_set = set()
    for subset in remove_indexes:
        remove_set = remove_set | subset
    for index in list(sorted(remove_set))[::-1]:
        phonons.pop(index)
    print('Effective k is {} W/mK'.format(heat_flux / (2/L)))

看一下结果

jojo

1
print('Effective k is {} W/mK'.format(heat_flux / (2 / L)))

Effective k is 68.82739201872596 W/mK