考虑界面的法向导热

Y.C. Hua, B.Y. Cao. Slip boundary conditions in ballistic-diffusive heat transport in nanostructures. Nanoscale and Microscale Thermophysical Engineering, 2017, 21(3): 159-176

程序每填一个新功能，就得补一篇超哥文章的示例.. 智慧的源泉⛲

法向导热和面向导热的区别主要是分布函数的分布方向和温度梯度方向是否一致，我们在推导的时候采用了局域平衡近似，认为沿温度梯度方向\(\nabla f \approx \nabla f_0\)。弛豫时间近似下的稳态声子玻尔兹曼方程可以写成 \[ \vec{v}\cdot \nabla_r f = - \frac{f - f_0}{\tau} \] 对于面向导热过程，即平行于薄膜的输运过程，\(f\)是沿垂直于薄膜分布的（沿y方向），而\(f_0\)的梯度（温度梯度）是平行于薄膜的（沿x方向），此时玻尔兹曼方程可以简化成 \[ v_x \frac{\partial f_0}{\partial x} + v_y \frac{\partial f}{\partial y} = - \frac{f - f_0}{\tau} \] 即 \[ f + v_y\tau \frac{\partial f}{\partial y} = f_0 - v_x\tau\frac{\partial f_0}{\partial x} \equiv -S_0 \] 这里\(f_0\)是已知的，因此右边变成了常数，方程的左边既有\(f\)也有\(f\)的梯度，会解出来指数分布的解。

对于法向导热过程，\(f\)是垂直于薄膜分布的（沿x方向），\(f_0\)的梯度也是垂直于薄膜分布的，此时玻尔兹曼方程可以简化成 \[ v_x \frac{\partial f_0}{\partial x} = - \frac{f - f_0}{\tau} \] 即 \[ f = f_0 - v_x\tau \frac{\partial f_0}{\partial x} \] 这里\(f\)的梯度不见了，\(f\)可以直接用\(f_0\)来表示，从这个式子可以直接得到扩散导热定律， \[ \begin{aligned} q_{x} & =\frac{1}{4 \pi} \int_\omega \hbar \omega D(\omega) \mathrm{d} \omega \int_0^{2 \pi} \mathrm{d} \varphi \int_0^\pi f v_x \sin \theta \mathrm{d} \theta \\ & =\frac{1}{2} \int_\omega \hbar \omega D(\omega) \mathrm{d} \omega \int_0^\pi (-v_x\tau \frac{\partial f_0}{\partial x} + f_0) v_x \sin \theta \mathrm{d} \theta \\ &= -\frac{1}{2} \int_\omega \hbar \omega D(\omega) \mathrm{d} \omega \int_0^\pi v^2\cos\theta^2\tau \sin\theta \frac{\partial f_0}{\partial x} \mathrm{d}\theta \\ &= -\frac{1}{2} \int_\omega \hbar \omega D(\omega) \mathrm{d} \omega \int_{-1}^{1} v^2\mu^2\tau \frac{\partial f_0}{\partial x} \mathrm{d}\mu \\ &= - \frac{1}{3} \frac{\partial f_0}{\partial x}\int_\omega \hbar\omega D(\omega)v^2\tau\mathrm{d}\omega \\ &= - \frac{1}{3}\int C_\omega vl\mathrm{d}\omega \frac{\partial T}{\partial x} \approx - \frac{1}{3}Cvl \frac{\partial T}{\partial x} \end{aligned} \] 其中 \[ C_\omega = \hbar\omega D(\omega)\frac{\mathrm{d}f_0}{\mathrm{d}T} \] 比热表达式的物理含义很好理解，比热就是能量对温度的导数，\(\hbar\omega\)就是一个声子的能量，\(D(\omega)\)就是频率为\(\omega\)的声子态的数目，\(f_0\)就是平衡态下对应频率的声子数目，所以\(\hbar\omega D(\omega)f_0\)就是平衡态下频率为\(\omega\)的声子的总能量，这里只有\(f_0\)和温度有关，所以对温度的导数就长成了上面的样子。这和体材料热导率的推导是一样的，因为热导率的定义就是沿着温度梯度方向的热流大小和温度梯度的关系。但是在薄膜的法向导热问题理，\(f_0\)是未知的，我们需要用热流连续性方程得到\(f_0\)的表达式。

黑体边界

热流是某一空间坐标处的分布函数对整个速度空间（立体角空间+频率空间）做积分得到的值。对于薄膜内部某一节点的热流，就是上一节所做的积分值，是扩散定律给出的值 \[ q_x = - \frac{1}{3}Cvl \frac{\partial T}{\partial x} \] 对于边界节点的热流，一半立体角空间由热沉辐射得到，另一半立体角空间由对\(f\)做积分得到。高温热沉向半球空间辐射的总热流为， \[ \begin{aligned} q_h &= \frac{1}{4\pi} \int_{0}^{2\pi} \mathrm{d}\varphi\int \hbar\omega D(\omega) \mathrm{d}\omega \int_0^{\pi/2} f_0(T_h)v\cos\theta \sin\theta \mathrm{d}\theta \\ &= \frac{1}{2} \int \hbar\omega D(\omega) \mathrm{d} \omega \int_{0}^{1} f_0(T_h)v \mu \mathrm{d}\mu \\ &= \frac{1}{4} \int \hbar\omega D(\omega) f_0(T_h)v \mathrm{d}\omega \\ \end{aligned} \] 同理可以得到低温热沉向半球空间辐射的总热流为 \[ q_c = \frac{1}{4}\int \hbar\omega D(\omega) f_0(T_c)v\mathrm{d}\omega \] 边界节点向热沉方向的半球空间输运的总热流为 \[ \begin{aligned} q_\text{half} & =\frac{1}{4 \pi} \int_\omega \hbar \omega D(\omega) \mathrm{d} \omega \int_{0}^{2\pi} \mathrm{d} \varphi \int_{\pi/2}^\pi (-v_x\tau \frac{\partial f_0}{\partial x} + f_0) v_x \sin \theta \mathrm{d} \theta \\ &= \frac{1}{2} \int_\omega \hbar \omega D(\omega) \mathrm{d} \omega \int_{\pi/2}^\pi -v^2\cos\theta^2\tau \sin\theta \frac{\partial f_0}{\partial x} + f_0v\cos\theta\sin\theta \mathrm{d}\theta\\ &= \frac{1}{2}\int_\omega \hbar\omega D(\omega)\mathrm{d}\omega \int_{0}^{-1} v^2\mu^2\tau\frac{\partial f_0}{\partial x} - f_0 v\mu \mathrm{d}\mu \\ &= \frac{1}{2}\int_\omega \hbar\omega D(\omega)\mathrm{d}\omega \left[-\frac{1}{3}v^2\tau \frac{\partial f_0}{\partial x} - \frac{1}{2}f_0v\right]\\ &= -\frac{1}{6}Cvl \frac{\partial T}{\partial x} - \frac{1}{4}\int_\omega \hbar\omega D(\omega)f_0(T_\text{h, i}) v\mathrm{d}\omega \end{aligned} \] 于是可以得到 \[ q_x = q_h + q_\text{half} \] 即 \[ - \frac{1}{3}Cvl \frac{\partial T}{\partial x} = \frac{1}{4} \int_\omega \hbar\omega D(\omega) f_0(T_h)v \mathrm{d}\omega - \frac{1}{6}Cvl \frac{\partial T}{\partial x} - \frac{1}{4}\int_\omega \hbar\omega D(\omega)f_0(T_{h, i}) v\mathrm{d}\omega \] 于是 \[ \begin{aligned} -\frac{1}{6}Cvl\frac{\partial T}{\partial x} &= \frac{1}{4}\int_\omega \hbar\omega D(\omega )v (f_0(T_h) - f_0(T_{h, i})) \mathrm{d}\omega \\ &\approx \frac{1}{4}(T_h - T_{h, i})\int_\omega \hbar\omega D(\omega) v \frac{\mathrm{d}f_0}{\mathrm{d}T} \mathrm{d}\omega \\ &\approx \frac{1}{4}Cv(T_h - T_{h, i}) \end{aligned} \] 于是可以得到 \[ - \frac{2}{3}l \frac{\partial T}{\partial x} = T_h - T_{h, i} \] 同理可以得到 \[ -\frac{2}{3}l\frac{\partial T}{\partial x} = T_{c, i} - T_c \] 还有 \[ - \frac{\partial T}{\partial x} = \frac{T_{h, i } - T_{c, i}}{L} \] 于是 \[ \frac{4}{3}\frac{l}{L}(T_{h,i} - T_{c, i}) = (T_h - T_c) - (T_{h, i} - T_{c, i}) \]

\[ \therefore (1 + \frac{4}{3}Kn) (T_{h,i} - T_{c, i}) = T_h - T_c \]

透射反射边界

分析的对象还是这个一维薄膜法向导热过程，在经典的黑体边界条件下，从热沉发射的声子全部进入体系内，体系内部运动撞到热沉的声子全都被吸收。假设薄膜与热沉接触的边界是一个漫透射反射边界，此时从热沉发射到体系内部的声子数目等于黑体边界发射的声子数目乘上透射率；从体系内部撞到热沉的声子，以这一侧的透射率的概率被吸收，剩下的概率被反射回体系内部。

在这种情况下，体系内部的热流还是\(q_x = - \frac{1}{3}Cvl \frac{\partial T}{\partial x}\)，区别主要在边界上。假设高温热沉到系统内的透射率为\(t_{hf}\)，系统内部到高温热沉的反射率为\(r_{fh}\)，此时从边界射向系统的总热流包含了透射热流和反射热流两部分，另外透射部分和反射部分的材料也是不一样的， \[ \begin{aligned} q_h &= \frac{1}{4\pi} \int_{0}^{2\pi} \mathrm{d}\varphi\int \hbar\omega D(\omega) \mathrm{d}\omega \left[\int_0^{\pi/2} f_0(T_h)v t_{hf }\cos\theta \sin\theta \mathrm{d}\theta - \int_{\pi/2}^{\pi} fv r_{fh}\cos\theta\sin\theta \mathrm{d}\theta \right] \\ &= \frac{1}{2} \int \hbar\omega D(\omega) \mathrm{d}\omega \left[\int_0^{1} f_0(T_h)v t_{hf }\mu \mathrm{d}\mu - \int_{-1}^{0} (f_0 - v\mu\tau \frac{\partial f_0}{\partial x})v r_{fh}\mu \mathrm{d}\mu \right] \\ &= \frac{1}{2} \int \hbar\omega D(\omega) \mathrm{d}\omega \left[\int_0^{1} f_0(T_h)v t_{hf }\mu \mathrm{d}\mu - \int_{-1}^{0} f_0v r_{fh}\mu \mathrm{d}\mu + \int_{-1}^{0}\tau \frac{\partial f_0}{\partial x}v^2 r_{fh}\mu^2 \mathrm{d}\mu \right] \\ \end{aligned} \]

假设透射率和反射率和方向无关，此时正向热流可以写成

$$ \[\begin{aligned} q_h &= \frac{1}{4} t_{hf} \int \hbar\omega D(\omega) v f_0(T_h) \mathrm{d}\omega + \frac{1}{4}r_{fh} \int \hbar\omega D(\omega) v f_0\mathrm{d}\omega + \frac{1}{6}r_{fh}\int\hbar\omega D(\omega) v^2\tau \frac{\partial f_0}{\partial x} \mathrm{d}\omega \\ &= \frac{1}{4} t_{hf} \int \hbar\omega D(\omega) v f_0(T_h) \mathrm{d}\omega + \frac{1}{4}r_{fh} \int \hbar\omega D(\omega) v f_0\mathrm{d}\omega + \frac{1}{6}r_{fh}\frac{\partial T}{\partial x} Cvl \\ \end{aligned}\]

$$

同样根据热流连续性关系，可以得到

\[ q_x = q_h + q_\text{half} \]

也就是

\[ - \frac{1}{3}Cvl \frac{\partial T}{\partial x} = \frac{1}{4} t_{hf} \int \hbar\omega D(\omega) v f_0(T_h) \mathrm{d}\omega + \frac{1}{4}r_{fh} \int \hbar\omega D(\omega) v f_0\mathrm{d}\omega + \frac{1}{6}r_{fh}\frac{\partial T}{\partial x} Cvl - \frac{1}{6}Cvl \frac{\partial T}{\partial x} - \frac{1}{4}\int_\omega \hbar\omega D(\omega)f_0 v\mathrm{d}\omega \]

可以得到

\[ - \frac{1}{6}(1 + r_{fh})Cvl\frac{\partial T}{\partial x} = \frac{1}{4} t_{hf} \int \hbar\omega D(\omega) v f_0(T_h) \mathrm{d}\omega + \frac{1}{4}(r_{fh}-1) \int \hbar\omega D(\omega) v f_0\mathrm{d}\omega \]

选取\(0\,\text{K}\)作为参考温度

\[ - \frac{1}{6}(1 + r_{fh})Cvl\frac{\partial T}{\partial x} = \frac{1}{4} t_{hf} \int \hbar\omega D(\omega) v (f_0(T_h) - 0) \mathrm{d}\omega + \frac{1}{4}(r_{fh}-1) \int \hbar\omega D(\omega) v (f_0-0)\mathrm{d}\omega \]

在灰体近似下，只有一个态，此时\(\hbar\omega D(\omega)f_0(T_h)\)就代表从0K到\(T_h\)所吸收的热量，于是上式可以改写成

\[ - \frac{1}{6}(1 + r_{fh})Cvl\frac{\partial T}{\partial x} = \frac{1}{4} t_{hf} C_hv_hT_h + \frac{1}{4}(r_{fh}-1) CvT_{h, i} \]

于是

\[ - \frac{2}{3}\frac{1 + r_{fh}}{1-r_{fh}}l\frac{\partial T}{\partial x} = \frac{t_{hf} C_hv_h}{Cv(1-r_{fh})}T_h -T_{h, i} \]

选取DMM模型来描述界面的透射率，

\[ t_{hf} = r_{fh} = \frac{Cv}{C_hv_h + Cv} \] 于是

\[ \frac{t_{hf} C_hv_h}{Cv(1-r_{fh})} = \frac{Cv/(C_hv_h + Cv) \times C_hv_h}{Cv (C_hv_h / (C_hv_h+Cv))} = \frac{CvC_hv_h}{CvC_hv_h} = 1 \] 于是

\[ T_h - T_{h,i} = - \frac{2}{3}\frac{1 + r_{fh}}{1-r_{fh}}l\frac{\partial T}{\partial x} \]

同理可以得到 \[ T_{c,i} - T_c = - \frac{2}{3}\frac{1 + r_{fc}}{1-r_{fc}}l\frac{\partial T}{\partial x} \] 于是

$$ T_h - T_c - (T_{h,i} - T_{c,i}) = Kn(T_{h,i} - T_{c,i}) ( + )

$$

\[ \therefore \left[1 + \frac{2}{3}Kn\left(\frac{1+r_{fh}}{1-r_{fh}} + \frac{1+r_{fc}}{1-r_{fc}}\right)\right] (T_{h,i} - T_{c, i}) = T_h - T_c \]

---

假如一个硅薄膜两边贴着锗的热沉， \[ C_\text{si} = 0.93\times 10^6 \,\mathrm{J/m^2K}, v_\text{si} = 1804\,\mathrm{m/s} \]

\[ C_\text{Ge} = 0.87\times 10^6 \, \mathrm{J/m^3K}, v_\text{Ge} = 1042\, \mathrm{m/s} \]

所以用DMM模型得到从硅到锗的穿透率就是0.35，从锗到硅的穿透率就是0.65。这个在蒙特卡洛仿真里也比较好实现，定义一个吸收-漫反射边界就可以了，如果随机数小于穿透率，就被界面吸收，否则被弹回来。发射按照原有的黑体发射就可以了，只是在最后的热流统计的时候用锗的性质来计算后在乘上锗到硅的穿透率。