沿薄膜的导热过程
沿薄膜的导热过程
陈刚. 纳米尺度能量输运和转换: 对电子, 分子, 声子和光子的统一处理[M]. 清华大学出版社, 2014.
华钰超. 微观与宏观结合的声子弹道扩散导热研究[D]. 清华大学, 2018.
Y.C. Hua, B.Y. Cao. Ballistic-diffusive heat conduction in multiply-constrained nanostructures. International Journal of Thermal Sciences, 2016, 101: 126-132
面向导热过程
沿薄膜的导热过程，又被称为面向导热过程，和泊肃叶流的情况差不多。温度梯度沿\(x\)方向，\(y\)方向是一个有限长度的薄膜，受到边界散射的影响，薄膜内的热流沿厚度方向的分布不再是均匀的，出现了热流滑移现象，这将导致薄膜的等效热导率低于其体材料的值。等效热导率降低的原因是在声子经过边界散射后，沿\(x\)方向输运的能量部分丢失了。如果上下两个边界都是镜面反射边界，那么载流子在和边界发生碰撞时沿\(x\)方向的分量将不会改变，此时等效热导率不会发生改变。而当边界部分发生漫反射时载流子碰撞到边界后的会重新按照黑体发射抽样方向，失去碰撞前向右的速度分量，导致热流减小，在相同的温度梯度下等效热导率降低。这里推导完全漫反射下的薄膜等效热导率的表达式。
在弛豫时间近似下，稳态声子玻尔兹曼方程可以写为 \[ \vec{v}\cdot \nabla_r f = - \frac{f - f_0}{\tau} \] 其中，\(f_0\)代表平衡态分布，\(f_0\)通过温度\(T\)随空间坐标而变化。对于薄膜内的二维导热问题，玻尔兹曼方程可写成 \[ v_x \frac{\partial f}{\partial x} + v_y\frac{\partial f}{\partial y} = - \frac{f - f_0}{\tau} \] 其中\(v_x = v\sin\theta\cos\varphi, v_y = v\cos\theta\)，这里极角\(\theta\)是从\(v_y\)方向算起。引入一个偏差函数\(g\)，定义为 \[ g = f - f_0 \] 此时玻尔兹曼方程就可以写成 \[ v_x \frac{\partial f_0}{\partial x} + v_x \frac{\partial g}{\partial x} + v_y\frac{\partial f_0}{\partial y} + v_y\frac{\partial g}{\partial y} = - \frac{g}{\tau} \] 式中，\(\partial f_0 / \partial y = 0\)，这是因为对沿薄膜方向的输运而言，\(f_0\)仅仅是\(x\)的函数。在\(x\)方向上，我们认为\(f\)的变化相比于\(f_0\)的变化很小，只保留\(\partial f_0/\partial x\)这一项，忽略\(\partial g / \partial x\)这一项。以上近似导出 \[ v_x\frac{\partial f_0}{\partial x} + v_y \frac{\partial g}{\partial y} = - \frac{g}{\tau} \] 或 \[ \tau v \mu \frac{\partial g}{\partial y} + g = - \tau v_x \frac{\partial f_0}{\partial x} \equiv - S_0 \] 这里\(\mu = \cos\theta\)。我们称\(S_0\)为源函数。在玻尔兹曼方程中，分布函数是三维空间坐标和三维动量空间坐标的函数。在这个二维导热问题中，由于对称性，\(g(y, \mu)\)是\(y\)和\(\mu\)的函数。与\(y\)相关是因为\(y\)方向有边界存在，所以分布函数\(f\)或\(g\)与\(y\)有关。
化简后的表达式是一个一阶微分方程，只需要一个边界条件。如果是这样，那薄膜的两个表面将怎样影响输运过程？答案在于我们需要对所有速度分量，即对表达式中所有的\(\mu\)方向取一个边界条件。在\(y=0\)取一个边界点，我们同时需要载流子流向边界和离开边界的边界条件。因为流向边界\(y=0\)的载流子可能来自另一个边界\(y=d\)，这些载流子的分布是未知的。因此，为流向边界的载流子给定已知的边界条件通常是很困难的。相反，我们通常在\(y=0\)和\(y=d\)处仅为离开边界的载流子给定边界条件，每一个都覆盖一半的角度空间，这与在一个边界给定整个空间的边界条件是相当的。这一点在随后的论述中将更加清晰。
一阶微分方程
\[ \tau v \mu \frac{\partial g}{\partial y} + g = - S_0 \]
定义\(\eta = y / (\tau v)\)，于是可以得到 \[ \mu \frac{\partial g}{\partial \eta} + g = - S_0 \] 可以用常数变易法来获得一阶微分方程的通解，先解 \[ \mu \frac{\partial g}{\partial \eta} + g = 0 \] 得到\(g = c \exp(-\frac{\eta}{\mu})\)，我们假设\(g\)的通解的形式为\(g = c(\eta)\exp(- \frac{\eta}{\mu})\)，代回到原微分方程中，可以得到 \[ c^\prime \exp(-\frac{\eta}{\mu}) = -S_0 \] 于是得到 \[ c(\eta) = - S_0 \exp(\frac{\eta}{\mu}) + c_1 \] 代回到通解中可以得到 \[ g = -S_0 + c_1\exp(-\frac{\eta}{\mu}) \] 上面提到了，一阶微分方程只需要一个边界条件，但是我们需要对所有速度分量，即对表达式中所有的\(\mu\)方向取一个边界条件。为了方便利用两个边界的信息，我们把通解按照方向拆成两个， \[ \begin{aligned} & g^{+}(\eta, \mu)=g(\eta, \mu)=g_1 \exp \left(-\frac{\eta}{\mu}\right)-S_0 \quad(0<\mu<1) \\ & g^{-}(\eta, \mu)=g(\eta, \mu)=g_2 \exp \left(\frac{\xi-\eta}{\mu}\right)-S_0 \quad(-1<\mu<0) \end{aligned} \] 其中\(\xi = d/v\tau\)是无量纲厚度，引入它只是为了后面代入边界条件的时候好看一点。当边界漫反射式地散射声子，并且没有声子透射出边界时，离开边界的声子呈各向同性分布，并遵守局域玻色-爱因斯坦分布，即\(f = f_0\)。我们得到 \[ g^+(y = 0) = 0\quad \text{和}\quad g^-(y = d) = 0 \] 代入通解中，可以得到 \[ g_1 = g_2 = S_0 \] 于是我们得到了分布函数的解 \[ \begin{aligned} & g^{+}(\eta, \mu)=g(\eta, \mu)=S_0\left[\exp \left(-\frac{\eta}{\mu}\right)-1\right] \quad(0 \leqslant \mu \leqslant 1) \\ & g^{-}(\eta, \mu)=g(\eta, \mu)=S_0\left[\exp \left(\frac{\xi-\eta}{\mu}\right)-1\right] \quad(-1 \leqslant \mu \leqslant 0) \end{aligned} \] 其中 \[ S_0 = - \tau v_x \frac{\partial f_0}{\partial x} \]
等效热导率
利用得到的解，我们可以计算沿\(x\)方向的总热流 \[ q = \int_0^d\int_0^{\omega_m}\hbar\omega\frac{D(\omega)}{4\pi}\int_{0}^{2\pi}\int_0^{\pi}\sin\theta f v_x\mathrm{d}\theta\mathrm{d}\varphi\mathrm{d}\omega\mathrm{d}y \] 把\(f = g + f_0\)代进去，发现\(f_0\)各向同性积分后就消掉了，于是 \[ q = \frac{1}{4\pi} \int_0^d \int_0^{\omega_m} \hbar\omega D(\omega) \int_0^{2\pi}\mathrm{d}\varphi \int_{0}^1 \left[g^+(y, \mu) + g^-(y, -\mu)\right] v_x\mathrm{d}\mu\mathrm{d}\omega\mathrm{d}y \] 代入\(g\)的表达式以及\(v_x = v\sin\theta\cos\varphi\)， \[ \begin{aligned} q &=- \frac{1}{4\pi} \int_0^d \int_0^{\omega_m} \hbar\omega D(\omega) \int_0^{2\pi}\mathrm{d}\varphi \int_0^1 \tau v_x \frac{\partial f_0}{\partial x}\left[\exp \left(-\frac{\eta}{\mu}\right)-1 + \exp \left(\frac{\xi-\eta}{-\mu}\right)-1\right] v_x\mathrm{d}\mu \mathrm{d}\omega\mathrm{d}y\\ &= - \frac{1}{4\pi} \int_0^d \int_0^{\omega_m} \hbar\omega D(\omega) \int_0^{2\pi}\mathrm{d}\varphi \int_0^1 \tau v^2 \sin^2\theta\cos^2\varphi \frac{\partial f_0}{\partial x}\left[\exp \left(-\frac{\eta}{\mu}\right) + \exp \left(\frac{\xi-\eta}{-\mu}\right)-2\right] \mathrm{d}\mu \mathrm{d}\omega\mathrm{d}y\\ &= - \frac{1}{4}\tau v^2 \frac{\mathrm{d} f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} \int_0^d \int_0^{\omega_m} \hbar\omega D(\omega) \int_0^1 (1 - \mu^2) \left[\exp \left(-\frac{\eta}{\mu}\right) + \exp \left(\frac{\xi-\eta}{-\mu}\right)-2\right] \mathrm{d}\mu \mathrm{d}\omega\mathrm{d}y\\ &= - \frac{1}{4}\tau^2 v^3 \frac{\mathrm{d} f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} \int_0^\xi \int_0^{\omega_m} \hbar\omega D(\omega) \int_0^1 (1 - \mu^2) \left[\exp \left(-\frac{\eta}{\mu}\right) + \exp \left(\frac{\xi-\eta}{-\mu}\right)-2\right] \mathrm{d}\mu \mathrm{d}\omega\mathrm{d}\eta\\ \end{aligned} \] 各个积分的积分限都是无关的，所以可以交换积分次序，先把厚度方向进行积分，跟\(\eta\)有关的部分是 \[ \begin{aligned} &\int_0^\xi\left[\exp \left(-\frac{\eta}{\mu}\right) + \exp \left(\frac{\xi-\eta}{-\mu}\right)-2\right]\mathrm{d}\eta \\ &= \int_0^\xi \mathrm{d}\left[-\mu \exp(-\frac{\eta}{\mu}) + \mu \exp(\frac{\xi - \eta}{-\mu}) - 2\eta\right] \\ &= -\mu \exp(-\frac{\xi}{\mu}) + \mu + \mu - \mu\exp(- \frac{\xi}{\mu}) - 2\xi\\ &= -2 \mu \left(\exp(-\frac{\xi}{\mu}) - 1\right) - 2 \xi \end{aligned} \] 代回到上面的积分中 \[ q = \frac{1}{2}\tau^2 v^3 \frac{\mathrm{d} f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} \int_0^{\omega_m} \hbar\omega D(\omega) \mathrm{d}\omega \int_0^1 (1 - \mu^2) \left[ \mu \left(\exp(-\frac{\xi}{\mu}) - 1\right) + \xi\right] \mathrm{d}\mu \\ \] 其中， \[ \int_0^{1} (1-\mu^2) \mathrm{d}\mu = \frac{2}{3} \] 于是 \[ q = \frac{1}{3} \tau^2v^3 \frac{\mathrm{d} f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} \xi \int_0^{\omega_m} \hbar\omega D(\omega) \mathrm{d}\omega+ \frac{1}{2}\tau^2 v^3 \frac{\mathrm{d} f_0}{\mathrm{d}T}\frac{\mathrm{d}T}{\mathrm{d}x} \int_0^{\omega_m} \hbar\omega D(\omega) \mathrm{d}\omega \int_0^1 (1 - \mu^2) \mu \left(\exp(-\frac{\xi}{\mu}) - 1\right) \mathrm{d}\mu \\ \] \(v\)和\(\tau\)是依赖于频率的，但是在灰体近似下，我们可以得到下面的表达式 \[ \frac{q/d}{\mathrm{d}T/\mathrm{d}x}= k_\text{eff} = \frac{1}{3} Cvl - \frac{1}{2}C v\frac{l^2}{d} \int_0^1 (1 - \mu^2) \mu \left(1 - \exp(-\frac{\xi}{\mu})\right) \mathrm{d}\mu \\ \] 体材料热导率 \[ k_\text{bulk} = \frac{1}{3}Cvl \] 于是 \[ k_\text{eff} / k_\text{bulk} = 1 - \frac{3l}{2d} \int_0^{1}(1 - \mu^2) \mu \left(1 - \exp(-\frac{\xi}{\mu})\right) \mathrm{d}\mu \] 文献里可能会有其他的表示形势，比如在这里给出的等效热导率的表达式是
Y.C. Hua, B.Y. Cao. Ballistic-diffusive heat conduction in multiply-constrained nanostructures. International Journal of Thermal Sciences, 2016, 101: 126-132
\[ G_{i n p}^{-1}=1-\frac{3 l_0}{2 L_y} \int_0^1\left(1-\exp \left(-\frac{L_y}{\sqrt{1-\mu^2} l_0}\right)\right) \mu^3 d \mu \]
这里是\(\theta\)极角选取的坐标系不同导致的，可以定义\(1-\mu^2 =x^2\)，然后代回我们得到的表达式里 \[ \begin{aligned}k_\text{eff} / k_\text{bulk} &= 1 - \frac{3l}{2d} \int_0^{1}(1 - \mu^2) \mu \left(1 - \exp(-\frac{\xi}{\mu})\right) \mathrm{d}\mu \\ &= 1 - \frac{3l}{2d} \int_0^{1}x^2 \sqrt{1-x^2} \left(1 - \exp(-\frac{\xi}{\sqrt{1-x^2}})\right) \mathrm{d}\sqrt{1-x^2} \\ &= 1 - \frac{3l}{2d} \int_0^{1}x^3 \left(1 - \exp(-\frac{\xi}{\sqrt{1-x^2}})\right) \mathrm{d}x \\ \end{aligned} \] 另一个更常见的表达形式是
陈刚. 纳米尺度能量输运和转换: 对电子, 分子, 声子和光子的统一处理[M]. 清华大学出版社, 2014.
\[ \frac{k_\text{eff}}{k_\text{bulk}}=1-\frac{3}{8 \xi}\left[1-4\left(E_3(\xi)-E_5(\xi)\right)\right] \]
其中\(E_n(x)\)为\(n\)阶指数积分 \[ E_n(x)=\int_0^1 \mu^{n-2} \exp \left(-\frac{x}{\mu}\right) \mathrm{d} \mu \] 同样可以从我们得到的表达式整理出这个结果，其实就是把各个项都分别算出来 \[ \begin{aligned}k_\text{eff} / k_\text{bulk} &= 1 - \frac{3l}{2d} \int_0^{1}(1 - \mu^2) \mu \left(1 - \exp(-\frac{\xi}{\mu})\right) \mathrm{d}\mu \\ &= 1 - \frac{3l}{2d} \left[\frac{1}{4}- \int_0^1 \mu\exp(-\frac{\xi}{\mu}) \mathrm{d}\mu+ \int_0^1 \mu^3 \exp(-\frac{\xi}{\mu}) \mathrm{d}\mu\right] \\ &= 1-\frac{3}{8 \xi}\left[1-4\left(E_3(\xi)-E_5(\xi)\right)\right] \end{aligned} \]
一些讨论
从上面的推导可以很清楚地看出来，在非傅里叶导热过程中，结构的等效热导率并只是几何结构的函数，而是和声子发射方式强相关的，不同的加热方式会导致不同的结果，因为本质上热导率这个概念是用傅里叶定律定义出来的，用在微纳尺度下只是为了直观地表现一些结果。所以很显然用法向热导率+面向等效热导率是没法完全描述一个薄膜内部的声子输运过程的。比如对于下面这个GaN HEMTs的仿真，即使对于各个薄膜都用上两个等效热导率，最后得到的结温也远远低于直接求解声子玻尔兹曼方程得到的结果，因为这时候热源集中在很小的一个区域，这是一个热扩展过程，而不是单纯地沿着某一个方向进行输运。
Hao Q, Zhao H, Xiao Y, et al. Electrothermal studies of GaN-based high electron mobility transistors with improved thermal designs[J]. International Journal of Heat and Mass Transfer, 2018, 116: 496-506.
