J-P. Péraud and N. G. Hadjiconstantinou, "An alternative approach to efficient simulation of micro/nanoscale phonon transport", Applied Physics Letters, 101, 153114, 2012.
# 这里采取了累积概率函数抽样的形式，np.random.choice可以完成相同的需求，根据cumsum_prob这个累积概率分布函数去抽样arr数组. @nb.jit(nopython=True) defrand_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")]
for particle inrange(N): Ri = np.random.rand() # 声子可能从三个源发射，左侧高温热沉，右侧低温热沉，以及体系内部。由于在这个模拟中体系的初始能量为0，所以实际上并不会在内部发射。在(0,1)间随机抽样一个数，确定声子从哪里发射。 if Ri < enrgInit / enrgTot: # 在这种情况下从内部发射 x0 = L * np.random.rand() # 确定发射的初始位置 phonon_index = rand_choice_nb(data.index.values, cumsum_base.values) # 在KMC里温度偏差用de_dT近似了，所以这里在cumsum_base里抽样就可以了 R = 2 * np.random.rand() - 1# 确定发射的仰角，R = cos(theta) vx = data.loc[phonon_index, 'velocity'] * R # x方向群速度的投影 psign = np.sign(Tinit - Teq) # 这颗能量团相对于参考温度的能量是正还是负 elif Ri > enrgInit / enrgTot and Ri < 1 - enrgRight / enrgTot: # 在这种情况下从左墙发射 x0 = 0 phonon_index = rand_choice_nb(data.index.values, cumsum_vel.values) R = np.sqrt(np.random.rand()) # R=cos(theta)，这时是面发射，所以累积分布函数和内部发射存在差别，包括仰角的积分上限以及可见面积 vx = data.loc[phonon_index, 'velocity'] * R psign = np.sign(Tl - Teq) else: # 此时从右墙发射 x0 = L phonon_index = rand_choice_nb(data.index.values, cumsum_vel.values) R = -np.sqrt(np.random.rand()) vx = data.loc[phonon_index, 'velocity'] * R psign = np.sign(Tr - Teq) t0 = np.random.rand() * tmax # 这个particle的初始时间.每个particle的初始时间是随机抽样的，最大为tmax，这个是为了描述时间相关的非稳态演化过程. finished = False# False代表着这个particle的生命周期还没结束，我们还要继续追踪它的运动. im = np.where(tt >= t0)[0][0] # 找到t0时间对应的时间数组的index，在fortran和matlab里数组的index从0开始，这和其他的语言不太一样》 whilenot finished: Delta_t = - data.loc[phonon_index, 'relaxation_time'] * np.log(np.random.rand()) # Delta_t就是到下一次散射前该particle运动的时间 t1 = t0 + Delta_t x1 = x0 + Delta_t * vx # 统计当前particle的贡献，这个逻辑是在两次散射间粒子可能运动了很长时间，跨过了很多个时间数组节点，因此我们要统计每个节点处particle的位置. while (im < Ntt and t0 <= tt[im] and t1 > tt[im]): x_ = x0 + (tt[im] - t0) * vx indx = int(np.floor(x_ / L * Nx)) if (indx < Nx) and (indx >= 0): T[im, indx] = T[im, indx] + psign * Eeff / C / Dx # 认为在小温差下比热是不变的，此时能量偏差和温度的对应关系是很简单的 Qx[im, indx] = Qx[im, indx] + psign * Eeff * vx / Dx im += 1 # 抽样散射后的声子团频率 phonon_index = rand_choice_nb(data.index.values, cumsum_col.values) # 更新声子团频率、速度、位置、时间、运动方向 R = 2 * np.random.rand() - 1 vx = data.loc[phonon_index, 'velocity'] * R x0 = x1 t0 = t1 if (t0 > tmax) or (x0 < 0) or (x0 > L): finished = True
