任意截面的导热过程

通解

温度梯度沿\(x\)方向，\(yz\)平面的形状任意但固定，\(x\)方向无限长。这个时候弛豫时间近似+小扰动近似下的偏差声子玻尔兹曼方程可以写成 \[ v_x\frac{\mathrm{d} f_0}{\mathrm{d} T} \frac{\mathrm{d}T}{\mathrm{d}x}+ v_y \frac{\partial g}{\partial y} + v_z \frac{\partial g}{\partial z} = - \frac{g}{\tau} \] 即 \[ v_y\tau\frac{\partial g}{\partial y} + v_z\tau \frac{\partial g}{\partial z} + g = - v_x\tau\frac{\mathrm{d}f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} = -S_0 \] 用特征线法来解这个一阶PDE，我们假设\(y\)和\(z\)是某个参数\(t\)的函数，并且存在 \[ \begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}t} &= v_y\tau, \frac{\mathrm{d}z}{\mathrm{d}t} = v_z\tau \end{aligned} \] 当\(t = 0\)时，\(y\)和\(z\)的值是确定的，我们假设\(y(0) = y_0, z(0) = z_0\)，于是我们可以得到 \[ y(t) = v_y\tau t + y_0, z(t) = v_z\tau t + z_0 \] 将特征线方程代回到原始微分方程中， \[ \frac{\mathrm{d}g}{\mathrm{d}t} + g = - S_0 \] 这和之前推导的面向导热过程一致了，我们可以得到这个一阶ODE的通解为 \[ g=-S_0 + C \exp \left(-t\right) \] 除了\(y\)和\(z\)，注意\(g\)还是方向的函数。假设边界的坐标可以表示为\(\vec{r_B} = (y_B, z_B)\)，若截面各边界是漫射边界条件，并且那么边界条件可以表示为边界处方向朝向结构内侧的\(g\)值为\(g(y_B, z_B)|_\text{inner}=0\)。"向内"这个方向关系用数学语言来描述就是\(y - y_B\)和\(v_y\)的符号是一样的。我们令\(t=0\)时\(y_0=y_B, z_0=z_B\)，此时特征线描述了分布函数从边界出发的传播过程，利用边界条件可以得到结构内侧的分布函数为 \[ g = S_0\left[\exp(- \left|\frac{y-y_B}{v_y\tau}\right|) - 1\right] = S_0\left[\exp(- \left|\frac{z-z_B}{v_z\tau}\right|) - 1\right] \] 我们还知道 \[ \begin{aligned} v_x &= v\cos\theta\\ v_y &= v\sin\theta\sin\varphi\\ v_z &= v\sin\theta\cos\varphi\\ \end{aligned} \] 于是可以把得到的分布函数整理成 \[ \begin{aligned} g &= S_0\left[\exp(- \frac{\sqrt{(y-y_B)^2 + (z-z_B)^2}}{(\sqrt{v_y^2 + v_z^2})\tau}) - 1\right]\\ &= S_0\left[\exp(- \frac{\left|\vec{r} - \vec{r_B}\right|}{l_0\sqrt{1 - \mu^2}}) - 1\right] \\ &= g_0 - S_0 \end{aligned} \] 这个表达式的物理含义是比较清晰的，\((\vec{r}, \vec{s})\)处的分布函数的值，由与其反向相交的边界处的分布函数的值确定，不同方向相交的边界是不一样的。分布函数从对应边界处传播到当前位置，经历了一个指数衰减过程。

获得了分布函数以后，可以获得\((y, z)\)处沿\(x\)方向的热流： \[ \begin{aligned} q_x(y, z) &= \frac{1}{4\pi} \int \hbar \omega D(\omega) d \omega \int_0^{2 \pi} \int_{-1}^1 g v_x \mathrm{d} \mu \mathrm{d} \varphi \\ &= \frac{1}{4\pi} \int \hbar \omega D(\omega) d \omega \int_0^{2 \pi} \int_{-1}^1 g v \mu\mathrm{d}\mu \mathrm{d} \varphi \\ &= \frac{1}{4\pi} \int \hbar \omega D(\omega) d \omega \int_0^{2 \pi} \int_{-1}^1 (g_0 - S_0) v \mu\mathrm{d}\mu \mathrm{d} \varphi \end{aligned} \] 积分得到整个截面处沿\(x\)方向的热流，\(A\)是截面的面积， \[ Q_x = \int_A q_x(y, z) \mathrm{d}s \] 其中\(S_0\)这一项与位置无关，我们先把它积掉 \[ \begin{aligned} Q_{S_0} &= \frac{A}{4\pi} \int \hbar \omega D(\omega) d \omega \int_0^{2 \pi} \int_{-1}^1 S_0 v \mu\mathrm{d}\mu \mathrm{d} \varphi \\ &= \frac{\mathrm{d}f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} A \int \hbar \omega D(\omega) v^2\tau d \omega \int_{0}^1 \mu^2\mathrm{d}\mu \\ &= A\frac{\mathrm{d}T}{\mathrm{d}x} \frac{1}{3} \int \frac{\mathrm{d}f_0}{\mathrm{d}T}\hbar\omega D(\omega)v^2\tau \mathrm{d}\omega \\ & = A\frac{\mathrm{d}T}{\mathrm{d}x}\frac{1}{3} \int Cvl \mathrm{d}\omega \end{aligned} \] 这就是体材料热导率导热这一项，所以前面的 \[ g_0 = S_0 \exp(- \frac{\left|\vec{r} - \vec{r_B}\right|}{l_0\sqrt{1 - \mu^2}}) \] 这一项体现了尺寸效应的影响，不同的结构\(g_0\)的表达式不同，这里要清楚\(\vec{r}_B\)是和分布函数的方向有关的。

矩形截面

假设我们研究的是一个纳米线中的导热过程，那么截面是一个长方形，我们可以分析一下这种情况。假设截面\(z\)方向的长度为\(L_z\)，\(y\)方向的长度为\(L_y\)，某一点的坐标是\((y, z)\)。这是一个对称的结构，因此我们只要分析方向的反向延长线和两个边界相交的情况就可以了，这把这个结果乘2就是最终结果了，我们这里分析左侧和上边界相交的情况。

当坐标处于\((y, z)\)的时候，从左侧边界出发的角度，就是从\(\varphi_1\)到\(\varphi_2\)这一段，这一段的\(g_0\)可以写成 \[ g_0 = S_0 \exp(- \frac{y}{v\sin\theta\sin\varphi}) \] 从上方边界出发的角度，就是\(\varphi_2\)到\(\varphi_3\)这一段， \[ g_0 = S_0 \exp(- \frac{L_z - z}{v\sin\theta\cos\varphi}) \] 其中 \[ \varphi_1 = \arctan\left[y / z\right], \varphi_2 = \pi/2 + \arctan\left[(L_z - z) / y\right], \varphi_3 = \pi + \arctan\left[(L_y - y)/(L_z - z)\right] \] 于是这一项对总热流的贡献可以写成一个五重积分的形式： \[ \begin{aligned} Q_{g_0} &= - \frac{\mathrm{d}f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x}\frac{1}{\pi} \int_0^{L_y}\mathrm{d}y \int_0^{L_z} \mathrm{d}z\int \hbar \omega D(\omega)v^2\tau d \omega \int_{0}^1 \mu^2\mathrm{d}\mu \left[ \int_{\varphi_1}^{\varphi_2}\exp(- \frac{y}{l_0\sqrt{1-\mu^2}\sin\varphi}) \mathrm{d}\varphi + \int_{\varphi_2}^{\varphi_3} \exp(- \frac{L_z - z}{l_0\sqrt{1-\mu^2}\cos\varphi}) \mathrm{d} \varphi \right]\\ \end{aligned} \]

灰体近似，同时用体材料热导率归一化，可以得到 \[ \begin{aligned} (Q_{g_0} + Q_{S_0}) / Q_{S_0} &= k_\text{eff} / k_\text{bulk}\\ &= 1 - \frac{3}{\pi L_y L_z} \int_0^{L_y}\mathrm{d}y \int_0^{L_z} \mathrm{d}z\int_{0}^1 \mu^2\mathrm{d}\mu \left[ \int_{\varphi_1}^{\varphi_2}\exp(- \frac{y}{l_0\sqrt{1-\mu^2}\sin\varphi}) \mathrm{d}\varphi + \int_{\varphi_2}^{\varphi_3} \exp(- \frac{L_z - z}{l_0\sqrt{1-\mu^2}\cos\varphi}) \mathrm{d} \varphi \right]\\ \end{aligned} \] ### 有限长度纳米线

Hua Y C, Cao B Y. Ballistic-diffusive heat conduction in multiply-constrained nanostructures[J]. International Journal of Thermal Sciences, 2016, 101: 126-132.

假设长度有限，这里还会涉及到法向的边界影响，这时候可以近似认为法向边界和侧向边界的影响是独立的，采用马西森定则来降低推导的困难程度。假如载流子在运动过程中受到两种散射机制A和B的作用，他们的散射率分别是\(\Gamma_A=1/\tau_A\)和\(\Gamma_B = 1/\tau_B\)。如果A、B两种机制的影响是相互独立的，那么总的散射率\(\Gamma\)就可以表示成\(\Gamma = \Gamma_A + \Gamma_B = 1/\tau_A + 1/\tau_B\)，也就是说我们可以分别研究在单独某一种散射机制影响下载流子的运动。邵成博士和Shiomi开发的P-TRANS中给出的所谓"马西森定则模拟方法"就是这种想法的应用，近似认为边界散射的影响可以独立于声子本征散射，此时总的散射率可以写成

\[ \frac{1}{\Lambda_\text{eff}(\omega)} = \frac{1}{\Lambda_\text{pp}(\omega)} + \frac{1}{\Lambda_\text{bdy}(\omega)} \] 这种分离从模拟的角度理解，就代表在仿真的时候可以只考虑边界散射而不再考虑本征散射，粒子只有撞到界面的时候才会发生透射/反射，其他时刻始终沿直线运动，这其实就和图形学里的光线追踪差不多了。通过这种模拟获得了\(\Lambda_\text{bdy}(\omega)\)后，可以直接通过马西森定则计算得到总散射率，最后再通过积分的方式获得结构的等效热导率。由于不需要在仿真的时候考虑本征散射过程，这大幅降低了模拟的计算需求。从解析的角度，我们也可以分别推导法向导热过程和侧向导热过程的抑制函数，最后通过马西森定则得到多约束的总计影响， \[ l_t = \frac{1}{l_0^{-1} + l_x^{-1} + l_c^{-1}} = \frac{l_0}{1 + l_0/l_x + l_0/l_c} \] 对于法向导热过程，我们已经知道了等效热导率的表达式： \[ k_\text{eff} / k_\text{bulk} = \frac{1}{\frac{4}{3}\frac{l_0}{L_x} + 1} \] 于是 \[ k_\text{eff} = \frac{1}{3} Cv^2 \frac{1}{\tau_0^{-1} + \tau_x^{-1}} = \frac{1}{3}Cv \frac{l_0}{1 + l_0/l_x} = \frac{k_\text{bulk}}{1 + l_0/l_x} \] 可以得到 \[ l_0 / l_x = \frac{4}{3}\frac{l_0}{L_x} \]

对于面向导热过程，上面我们也已经推出了等效热导率的表达式，可以整理出来 \[ \begin{aligned} k_\text{eff} / k_\text{bulk} &= \frac{l_0^{-1}}{l_c^{-1} + l_0^{-1}}\\ &= \frac{1}{l_0/l_c + 1}\\ &= 1 - \frac{3}{\pi L_y L_z} \int_0^{L_y}\mathrm{d}y \int_0^{L_z} \mathrm{d}z\int_{0}^1 \mu^2\mathrm{d}\mu \left[ \int_{\varphi_1}^{\varphi_2}\exp(- \frac{y}{l_0\sqrt{1-\mu^2}\sin\varphi}) \mathrm{d}\varphi + \int_{\varphi_2}^{\varphi_3} \exp(- \frac{L_z - z}{l_0\sqrt{1-\mu^2}\cos\varphi}) \mathrm{d} \varphi \right] \end{aligned} \] 所以 \[ \frac{1}{1 - \frac{3}{\pi L_y L_z} \int ...} - 1 = l_0 / l_c \] 所以 \[ k_\text{wire} / k_\text{bulk} = 1 / (1 + \frac{4}{3}\frac{l_0}{L_x} + \frac{1}{1 - \frac{3}{\pi L_y L_z} \int ...} - 1) = 1 / (\frac{4}{3}\frac{l_0}{L_x} + \frac{1}{1 - \frac{3}{\pi L_y L_z} \int ...}) \] 这个结果没法解析求解，只能够数值求解。

import numpy as np
from numba import njit


@njit
def compute_integral(L_y, L_z, l_0, n_y, n_z, n_mu, n_phi):
    # 创建离散点，避免端点
    y_vals = np.linspace(0, L_y, n_y)[1:-1]
    z_vals = np.linspace(0, L_z, n_z)[1:-1]
    mu_vals = np.linspace(0, 1, n_mu)[1:-1]
    
    dy = y_vals[1] - y_vals[0]
    dz = z_vals[1] - z_vals[0]
    dmu = mu_vals[1] - mu_vals[0]

    # 初始化积分结果为0
    integral = 0.0
    
    # 执行四重循环来计算离散积分
    for i in range(len(y_vals)):
        for j in range(len(z_vals)):
            for k in range(len(mu_vals)):
                
                # 计算φ的限制
                phi_1 = np.arctan((y_vals[i]) / (z_vals[j]))
                phi_2 = np.pi / 2 + np.arctan((L_z - z_vals[j]) / (y_vals[i]))
                phi_3 = np.pi + np.arctan((L_y - y_vals[i]) / (L_z - z_vals[j]))
                
                phi_vals1 = np.linspace(phi_1, phi_2, n_phi)[1:-1]
                phi_vals2 = np.linspace(phi_2, phi_3, n_phi)[1:-1]

                dphi1 = phi_vals1[1] - phi_vals1[0]
                dphi2 = phi_vals2[1] - phi_vals2[0]
                
                for l in range(len(phi_vals1)):
                    # 计算两个不同φ区间的积分部分
                    term1 = np.exp(- y_vals[i] / (np.sin(phi_vals1[l]) * np.sqrt(1 - mu_vals[k]**2) * l_0)) * mu_vals[k]**2
                    integral += term1 * dphi1
    
                for l in range(len(phi_vals2)):
                    # 计算两个不同φ区间的积分部分
                    term2 = np.exp((L_z - z_vals[j]) / (np.cos(phi_vals2[l]) * np.sqrt(1 - mu_vals[k]**2) * l_0)) * mu_vals[k]**2
                    integral += term2 * dphi2

    return integral

k = []

for L in [0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:
    # 定义常量
    L_y = L
    L_z = L
    L_x = 10.0
    l_0 = 1.0

    # 定义离散点的数量
    n_y = 100
    n_z = 100
    n_mu = 100
    n_phi = 100

    # 计算离散积分的最终结果
    integral = compute_integral(L_y, L_z, l_0, n_y, n_z, n_mu, n_phi)
    integral *= (3 / (np.pi * L_y * L_z)) * (L_y / n_y) * (L_z / n_z) * (1 / n_mu) 

    # 计算G_rec^-1的最终结果
    G_rec_inv = 1 - integral
    G = 1 / G_rec_inv
    F = 1 + 4 * l_0 / (3 * L_x)

    k.append(1 / (F + G - 1))