前言

本学期选修了中国科学院大学开设的表面与界面的统计热力学,正好好久没写博客了,干脆把note放博客上面,笔记写的很乱,可能不严谨,一堆笔误之类的,凑合看

表界面热力学理论框架

该课程的理论主要包含四部分,分别是预备知识、 界面张力及涨落不稳定性、界面的浸润现象和界面相互作用。其中预备知识包含基本的热力学和统计物理知识,并会介绍部分微分几何知识(只包含经典微分几何,不含黎曼几何),第二部分主要关注瑞利-普拉托不稳定性;第三部分关注多相界面,第四部分主要讲述DLVO theory

第一章 预备知识

该部分主要介绍基本的热力学和统计物理,以及后续会用到的部分微分几何知识

热力学基础回顾

热力学定律是物理学中描述能量转换和热力学过程的基本原理,众所周知,当我们提到热力学时,通常会谈到热力学四大定律,它们分别是:

热力学第0定律:若两个热力学系统均与第三个系统处于热平衡状态,此两个系统也必互相处于热平衡;考虑三个热力学系统1、2、3;若系统1和3处于热平衡状态,即T_1=T_3;系统2和系统3也处于热平衡状态,即T_2=T_3,那么系统1和系统2也处于热平衡状态,即T_1=T_2;热力学第0定律揭示了热平衡的递移性,从而为温度的定义和测量提供了理论基础

热力学第1定律:该定律是能量守恒定律对非孤立系统的扩展,其可以写成:

\displaystyle dU=dW+dQ

这里我们均选取指向系统的方向为正方向,便于后续统一分析(即dW代表外界对系统做的功,dQ代表外界传入系统的热量)

热力学第2定律:该定律有克劳修斯表述和开尔文表述,为定量描述热力学第二定律,克劳修斯引入“熵”这一状态函数。其核心特征是:熵是广延量,仅由系统的平衡态决定,与过程无关。基于熵的定义,热力学第二定律可通过 “熵增原理” 转化为数学关系:

\displaystyle \Delta S \geq \frac{\Delta Q}{T}

移项得\Delta S-\frac{\Delta Q}{T}\geq0,将不等式左边定义为\Delta S^{tot},其中\Delta S^{tot}=\Delta S^{sys}+\Delta S^{res},sys代表系统,res代表热浴,其中\Delta S^{sys}=\Delta S\Delta S^{res}=-\frac{\Delta Q}{T}。其中不等式中的等于符号取在平衡态,可用于作为非平衡的判据

结合第1与第2定律,进行移项\Delta W=\Delta U-\Delta Q,之后利用\Delta W=\Delta U-T\Delta S+T\Delta S-\Delta Q,可以得到最大功原理,可以发现\Delta U-T\Delta S即为\Delta F,而后面部分其实就是T\Delta S^{tot},故而有\Delta W=\Delta F+T\Delta S^{tot},而第二定律又给出了这个表达式的边界,因此有\Delta W\geq \Delta F,改写成-\Delta W\leq -\Delta F,也就给出了最大功原理(亥姆霍兹自由能最小原理)

上述热力学第2定律还可以通过涨落定理导出,或者说涨落定理是热力学第2定律的统计证明(知乎上说的),由Crooks涨落定理:

\displaystyle <e^{-\Delta s^{tot}}>=1

此处\Delta s无量纲,玻尔兹曼常数取k_B=1,利用琴生不等式(对于凸函数,有<e^{x}>\geq e^{< x >}),可知<e^{-\Delta s^{tot}}>\geq e^{-< \Delta s^{tot} >},即1 \geq e^{-< \Delta s^{tot} >},由于<\Delta s^{tot}>=\Delta S^{tot},故得\Delta S^{tot}\geq0 (注意区分小s与大S,小s代表远离平衡位置的随机熵,大S则是总的熵)

同理,最大功原理可以通过Jarzynski equality导出,Jarzynski equality即:

\displaystyle <e^{-w}>=e^{-\Delta F}

根据上面类似的操作,容易得到e^{-\Delta F}\geq e^{-< w >},最后得到-\Delta W\leq -\Delta F

热力学第3定律:绝对零度不可达到

上面主要关注热力学第2定律及熵增定律,因为\Delta S^{tot}也即指明了熵产生,而\Delta S-\frac{\Delta Q}{T}\geq0则定义了平衡态,用于作为平衡态与非平衡态的判据。而对于自由能,直译便是“自由的能量”,中译中即热力学自由能可以用于对外做功,自由能是系统的内能减去不能用于做功的能量,因为前面提到了熵产生,并不是所有内能都可用于做功:

\displaystyle \begin{cases}
F=U-TS\\
T=\frac{\partial U}{\partial S}
\end{cases}

这也和前面提到的最大功原理-\Delta W\leq -\Delta F联系起来了

Legendre变换

经典热力学的目标是描述宏观性质的平均行为,而不是每个粒子或自由度的微观细节。需要注意,这里的Legendre变换是一个平衡态里才有的概念。考虑一个系统sys(S,V,N),其中S是熵,V是体积,N是粒子数;可以得到热力学第二定律:

\displaystyle
\begin{cases}
dU=TdS-PdV+\mu dN\\
(\frac{\partial U}{\partial S})_{V,N}\equiv T, (\frac{\partial U}{\partial V})_{S,N}\equiv -P,(\frac{\partial U}{\partial N})_{S,V}\equiv \mu
\end{cases}

当然这并非唯一的选择,由于熵难以测量,我们更希望用T作为变量而非S,因此我们可以做下面的操作:

\displaystyle \begin{cases}
d[U-TS] = -SdT - PdV + \mu dN = F(T,V,N)\\
T = \left( \frac{\partial U}{\partial S} \right)_{V,N}
\end{cases}

以此来得到F(T,V,N),同理类推你可以把V用类似的操作替换成P,把N替换成\mu等,但是这些热力学变量并不是相互独立的,譬如你利用d[U-TS+PV-\mu N]=-SdT+VdP-Nd\mu,你会发现你得到了Gibbs-Duhem equation,也就是说上面那一堆是等于0的,d[U-TS+PV-\mu N]=-SdT+VdP-Nd\mu=0,一种说法是宏观热力学是具有广延性的,而T,P,\mu都是强度量,S,V,N才是广延量。证明这个需要利用一个热力学势函数的性质,即热力学势函数都是欧拉一阶齐次函数

欧拉n阶齐次函数是指:

\displaystyle f(\lambda x_1,\lambda x_2,\cdots,\lambda x_N)=\lambda^n f(x_1,x_2,\cdots,x_N)

其具有这样的性质:

\displaystyle nf(x_1,x_2,\cdots,x_N)=\sum_{i=1}^{N} x_i\frac{\partial f}{\partial x_i}

证明过程只需要对其定义式两边对\lambda做偏导即可:

\displaystyle \sum_{i=1}^{N}\frac{\partial f}{\partial (\lambda x_i)}\frac{\partial (\lambda x_i)}{\partial \lambda}=n\lambda^{n-1}f(x_1,x_2,\cdots,x_N)

此时令\lambda=1即得到nf(x_1,x_2,\cdots,x_N)=\sum_{i=1}^{N} x_i\frac{\partial f}{\partial x_i}

回到我们先前讨论的Gibbs-Duhem equation上来,U(S,V,N)=\frac{\partial U}{\partial S}S+\frac{\partial U}{\partial V}V+\frac{\partial U}{\partial N}N=TS-PV+\mu N,自然而然的有d[U-TS+PV-\mu N]=0,也即证明了Gibbs-Duhem equation

需要注意的是这是一个对于宏观热力学系统才成立的方程,对于纳米系统如分子层次的微小结构,此时不可使用该定律,方程右边可能不为0,因为此时的热力学势函数可能不再是欧拉一阶齐次函数

统计力学基础

统计力学是理论物理的一个大分支,它利用概率论和统计学的方法来研究由大量粒子组成的宏观系统;由于该部分的讲述针对完全未学过统计力学的学生,没从Boltzmann那一套讲法出发,于是这部分直接copy上课板书

首先,在统计力学中,我们引入状态变量(state variable),状态变量是指在动态系统中,可以描述系统数学状态的一组变量,它应能确定系统未来的演化行为,这里我们记状态变量为\omega(t)

有了状态变量后,我们就有了概率,写为p(\omega)d\omega,其中p(\omega)是概率密度函数,在全空间中,应当有\int p(\omega)d\omega=1,之后我们有可观察量(某一个具体的物理量):A(\omega),则其平均值< A >=\int A(\omega)p(\omega)d\omega,其实这就是统计力学的核心,定义状态变量,获得概率密度函数,然后拿去算你需要算的物理量

我们落实到具体的physical quantities来讲,在一个孤立的系统中,哈密顿量通常代表系统的内能,即<\hat{H}(\omega)>=U,注意,这里的\hat{H}并非哈密顿算子,只是写成这样方便区分,根据玻尔兹曼熵公式有S=<-k_B\log{p(\omega)}>,这里的\log其实是\ln,为了和板书一致我就没有更改。这里的S便是Gibbs-Shannon entropy,而\hat{s}=-k_B\log{p(\omega)}则是Stochastic entropy

接着有F=U-TS=<\hat{H}(\omega)>+k_BT<\log{p(\omega)}>=<\hat{H}(\omega)+k_BT\log{p(\omega)}>

这一堆东西还可以写成<\hat{f}(\omega)>,其中\hat{f}(\omega)=\hat{H}(\omega)+k_BT\log{p(\omega)}=\hat{H}(\omega)-T\hat{s}(\omega),又称Stochastic free energy

接着讲了Free energy(non-equilibrium),其可以根据下面的公式进行计算:

\displaystyle F=\int d\omega p(\omega)[\beta \hat{H}(\omega)+\log{p(\omega)}]k_BT\\
=k_BT\int d\omega p(\omega)\log\frac{p(\omega)}{\frac{1}{Z}e^{-\beta\hat{H}(\omega)}}-k_BT\log{Z}\\
=k_BT\int d\omega p(\omega)\log{\frac{p(\omega)}{p_s(\omega)}}-k_BT\log{Z}\\
=k_BTD_{KL}(p(\omega)//p_s(\omega))-k_BT\log{Z}

其中\beta=\frac{1}{k_BT},另外有Z=\int d\omega e^{-\beta\hat{H}(\omega)},并且\frac{1}{Z}e^{-\beta\hat{H}(\omega)}=p_s(\omega),且\int d\omega p_s(\omega)=1

D_{KL}是Kullback-Leibler Divergence,也可以叫Relative entropy,其具体定义为,给定两个概率密度函数PDF,p_1(\omega)p_2(\omega),其相对熵可以写成:

\displaystyle D_{KL}(p_1(\omega)//p_2(\omega))=<\log{\frac{p_1(\omega)}{p_2(\omega)}}>_{p_1(\omega)}\\
=\int d\omega p_1(\omega)\log{\frac{p_1(\omega)}{p_2(\omega)}}

下标处的p_1(\omega)代表以p_1作为概率测度。最后自由能可以写成:

\displaystyle F-(-k_BT\log{Z})=D_{KL}(p(\omega)//p_s(\omega))

其中F_{eq}=(-k_BT\log{Z})为平衡态自由能,而D_{KL}(p(\omega)//p_s(\omega))为凸函数,且D_{KL}(p(\omega)//p_s(\omega))\geq0


接下来先review一下之前的东西,上次我们说到对于Nano system,宏观的系统广延性可能会丢失,然后我们接着开启了统计力学的部分

对于统计力学,首先第一步是确定状态变量,这里举两个例子,经典:相空间中的任意一点;量子:波函数。需要注意的是,状态空间有连续和离散两种区别,需要分别处理

在给定状态变量\omega后,我们就有了Probability:

\displaystyle p(\omega)=\begin{cases}
probability,discrete \ space\\
probability \ density \ function,continuous \ space
\end{cases}

此处可以看出,若是连续空间,那此时某一点的p(\omega)=0,是无意义的

之后给定一个可观测量:<\hat{A}(\omega)>,那么显然对于随机变量\hat{A}(\omega)有:

\displaystyle <\hat{A}(\omega)>=\begin{cases}
\sum_{\omega}p(\omega)\hat{A}(\omega),discrete\\
\int d\omega p(\omega)\hat{A}(\omega),continuous
\end{cases}

再次回到具体的physical quantities上来,在一个孤立的系统中,哈密顿量通常代表系统的内能,即<\hat{H}(\omega)>=U,熵则可以写成S=<\hat{s}(\omega)>,其中\hat{s}是stochastic entropy,有\hat{s}=-k_B\log{p(\omega)}

这里对\hat{s}做一点额外的解释,熵是信息的度量,一个事件发生的概率越小则其所蕴含的信息量越大,具体写为:

\displaystyle \hat{s}=-k_B\log p(\omega)=k_B\log{\frac{1}{p(\omega)}}

显然有p(\omega)越小,\hat{s}越大,而对于玻尔兹曼熵变,则有:

\displaystyle S_{Boltzmann}=k_B\log{W}=k_B\log{\frac{1}{\frac{1}{W}}}

其中\frac{1}{W}=p,如果一个微观状态的概率很小,意味着该状态越不常见,系统越可能处于其他状态,这增加了系统的混乱程度,因此熵就越大。反之,如果一个微观状态的概率很大,意味着该状态越常见,系统越可能稳定地处于该状态,熵就越小。比如离散的情况下:

\displaystyle <k_B\log{W}>_p=\sum_{\omega}p(\omega)k_B\log{W}

此外,关于为什么熵公式中含有\log函数,我们进行了如下的讨论:用语言解释即是熵的广延性(可加性);在物理上,我们先令S=k_Bf(\omega),其中f(\omega)是一个未知的函数,此时考虑一个系统sys,其具有微观状态数W,此时将该系统沿着中间切成两半,分成系统1和系统2,并分别具有微观状态数W_1W_2,二者只在界面处有交换,令该热力学系统为无穷大,此时交换可以忽略,那么二者的微观状态数应该满足:

\displaystyle W_{sys}=W_1\cdot W_2

而熵应该满足:

\displaystyle S_{sys}=S_1+S_2

而只有\log函数具有这样的性质,因为:

\displaystyle S_{sys}=S_{1+2}=k_B\log{W}=k_B\log{W_1W_2}\\
=k_B\log{W_1}+k_B\log{W_2}=S_1+S_2

接着是自由能F=U-TS,stochastic free energy可以写成\hat{f}(\omega)=\hat{H}(\omega)-T\hat{s}(\omega),然后有:

\displaystyle F=<\hat{f}(\omega)>=<\hat{H}(\omega)-T\hat{s}(\omega)>

还可以写成:

\displaystyle F=k_BTD_{KL}(p(\omega)//p_s(\omega))+F_{eq}

其中F_{eq}=-k_BT\log{Z},为平衡态自由能,Z为配分函数,Z=\int d\omega e^{-\beta \hat{H}(\omega)},并且有:

\displaystyle D_{KL}(p(\omega)//p_s(\omega))=\int d\omega p(\omega)\log{\frac{p(\omega)}{p_{s}(\omega)}}

p_s(\omega)=\frac{1}{Z}e^{-\beta\hat{H}(\omega)},也即Gibbs-Boltzmann probability

KL散度

前面我们提到过:

\displaystyle k_BTD_{KL}(p(\omega)//p_s(\omega))=F-F_{eq}

接着详细讨论一下KL散度,KL散度还具有一些别的性质,比如非负性,即:

\displaystyle D_{KL}(p(\omega)//p_s(\omega))\geq 0

若我们把F写成p(\omega)的泛函,则有F[p(\omega)]-F_{eq}=D_{KL}(p(\omega)//p_{eq}),认真观察一下可以发现我们得到了热力学中的极小原理:最小自由能原理。我们进行简单的证明,首先,D_{KL}是凸函数,我们将其按照数学定义写为:

\displaystyle D_{KL}(p(\omega)//p_{eq}(\omega))=\mathbb{E}^{p(\omega)}[\log{\frac{p(\omega)}{p_{eq}(\omega)}}]\\
= \mathbb{E}^{p(\omega)}[-\log{\frac{p_{eq}(\omega)}{p(\omega)}}]

其中\mathbb{E}代表数学期望,上标代表以p(\omega)为概率测度。由于-\log{x}是凸函数,因此我们利用琴生不等式,有:

\displaystyle \mathbb{E}[-\log{x}]\geq -\log{\mathbb{E}(x)}

进而可以得到:

\displaystyle \mathbb{E}^{p(\omega)}[-\log{\frac{p_{eq}(\omega)}{p(\omega)}}]\geq -\log{\mathbb{E}^{p(\omega)}[\frac{p_{eq}(\omega)}{p(\omega)}]}

其中等于号取在:

\displaystyle -\log{\int d\omega p(\omega)\frac{p_{eq}(\omega)}{p(\omega)}}=-\log{\int d\omega p_{eq}(\omega)}=0

这里利用了\int d\omega p_{eq}(\omega)=1,这也就说明了,无论系统从哪一个态出发,最后会演化到唯一的平衡态。这里其实还引出了Lyapunov function,不过此处不再进行深入讨论了

动力学

接下来我们研究系统的动力学,即系统从非平衡态到稳态(Steady state)的过程。考虑一个具有遍历性(Ergodicity)的系统,我们需要概率p(\omega,t)的演化方程,其中该系统全概率守恒,即:

\displaystyle \sum_{\omega}p(\omega,t)=1

而所谓遍历性,在数学中指动力系统的一种性质,它表明系统在长时间内可以访问其状态空间的任意一部分,表现为系统性质的不变性。接下来我们引入主方程:

\displaystyle \frac{d}{dt}p_{\omega}(t)=\sum_{\omega'}W_{\omega\omega'}p_{\omega'}(t)

这里假设系统的未来发展仅取决于其当前状态,而不是其过去的历史,即做了马尔可夫近似,另外为了简单我们仅考虑了离散的情形,方程中W_{\omega\omega'}是转移速率矩阵,给定时间t,系统由\omega'态演化到\omega态,这其中,由于\sum_{\omega}p(\omega,t)=1,因此对于任意p_{\omega'}(t),可以得到:

\displaystyle\frac{d}{dt}\sum_{\omega}p_{\omega}(t)=\sum_{\omega'}(\sum_{\omega}W_{\omega\omega'})p_{\omega'}(t)=0

显然有\sum_{\omega}W_{\omega\omega'}=0,这说明了该转移速率矩阵在任一列的加和为0,继续写即:

\displaystyle W_{\omega'\omega'}+\sum_{\omega\neq\omega'}W_{\omega\omega'}=0

最后有W_{\omega'\omega'}=-\sum_{\omega\neq\omega'}W_{\omega\omega'},也即该矩阵的任意一个对角元等于负的该列其他元素的和。物理上,这表示对于任意状态\omega',其概率流失速率(由对角元W_{\omega'\omega'}描述)在数值上等于从该状态转移到所有其他状态的总速率(由非对角元之和\sum_{\omega\neq\omega'}W_{\omega\omega'}给出)

主方程的解的形式为:

\displaystyle p_{\omega}(t)=\sum_{\lambda}C_{\lambda}e^{-\lambda t}\phi_{\lambda}(\omega)

其中\lambda为特征值且\lambda >0\phi_{\lambda}(\omega)对应特征值\lambda的特征函数,而C_{\lambda}为展开系数,继续推导还可以得到:

\displaystyle \sum_{\omega'}W_{\omega\omega'}\phi_{\lambda}(\omega')=-\lambda\phi_{\lambda}(\omega)

这里如果我们考虑具体的例子,比如考虑\omega代表位置,即\omega\rightarrow\vec{x},则对于连续情况,我们有:

\displaystyle \frac{\partial}{\partial t}p(\vec{x},t)=-\nabla\cdot\vec{J}(\vec{x},t)

其中流矢量写成:

\displaystyle\vec{J}(\vec{x},t)=-D\nabla p(\vec{x},t)+\vec{u}(\vec{x},t)p(\vec{x},t)

如果只看流矢量的第一部分,就得到了扩散方程,或者说Fick定律,这部分也称为扩散流,D为扩散系数。而流矢量的第二部分为漂移流,两部分一起考虑,就得到了Fokker–Planck equation,为了研究Local velocity,我们把流矢量写成:

\displaystyle\vec{J}(\vec{x},t)=p(\vec{x},t)\vec{\gamma}(\vec{x},t)\\
=p(\vec{x},t)[\vec{u}(\vec{x},t)-D\frac{1}{p(\vec{x},t)}\nabla p(\vec{x},t)]\\
=p(\vec{x},t)[\vec{u}(\vec{x},t)-D\nabla\log{p(\vec{x},t)}]

因此有:

\displaystyle\vec{\gamma}(\vec{x},t)=\vec{u}(\vec{x},t)-D\nabla\log{p(\vec{x},t)}

\vec{u}(\vec{x},t)不含时,即\vec{u}(\vec{x},t)=\vec{u}(\vec{x}),我们可以把它写成\vec{u}(\vec{x})=\kappa\vec{f}(\vec{x})=-\kappa\nabla\hat{H}(\vec{x})的形式,其中\vec{f}=-\nabla\hat{H}(\vec{x})\kappa为迁移率

继续推导我们有:

\displaystyle \vec{\gamma}(\vec{x},t)=-\kappa\nabla\hat{H}(\vec{x})-D\nabla\log{p(\vec{x},t)}\\
=-\kappa\nabla[\hat{H}(\vec{x})+\frac{D}{\kappa}\log{p(\vec{x},t)}]

此时若令D=k_BT\kappa=\frac{k_BT}{\xi},其中\xi=\kappa^{-1},则我们可以得到爱因斯坦关系(Einstein relation)如下:

\displaystyle \vec{\gamma}(\vec{x},t)=-\kappa\nabla[\hat{H}(\vec{x})+k_BT\log{p(\vec{x},t)}]\\
=-\kappa\nabla\hat{f}(\vec{x},t)

\hat{f}即前面提到的Stochastic free energy,其中D为扩散系数,而\xi为摩擦系数,得到爱因斯坦关系的前提要求是t\rightarrow\infty并且\vec{\gamma}(\vec{x},t)=0,也即-\hat{H}(\vec{x})=\frac{D}{\kappa}\log{p(\vec{x},t)},或者写成-\frac{\kappa}{D}\hat{H}(\vec{x})=\log{p(\vec{x},t)},此时p(\vec{x},t)=e^{-\frac{\kappa}{D}\hat{H}(\vec{x})}

其实仔细回顾这部分推导还可以看出,没有外界驱动的情况下,系统会演化到玻尔兹曼分布,其实就是解Fokker–Planck equation的定态解,很显然,平衡态时,密度分布不显含时间

这里我自己再啰嗦两句,在物理学的诸多领域中,我们常会观察到一种极具共性的数学结构:许多核心的物理矢量场,都可以通过某个标量函数的负梯度来定义。比如F=-\nabla V,亦或是E=-\nabla \phi,以及这里的\vec{f}=-\nabla\hat{H}(\vec{x}),第一个公式揭示了保守力场中势能与保守力的关系,第二个公式揭示了电场中电场强度与静电势的关系,而值得注意的是,磁场中并不存在这样的数学结构,这一本质差异,恰恰从数学层面印证了磁场的一个关键物理特性:它不存在平衡态,或者说它无法达到真正的平衡态

Rewrite一下主方程:

\displaystyle \frac{d}{dt}p_{\omega}(t)=\sum_{\omega'}W_{\omega\omega'}p_{\omega'}(t)=\sum_{\omega\neq\omega'}W_{\omega\omega'}p_{\omega'}(t)-\sum_{\omega\neq\omega'}W_{\omega'\omega}p_{\omega}(t)

上面式子可以合并为

\displaystyle \frac{d}{dt}p_{\omega}(t)=\sum_{\omega\neq\omega'}[W_{\omega\omega'}p_{\omega'}(t)-W_{\omega'\omega}p_{\omega}(t)]

两部分分别表示获得和流失的概率,我们可以定义从状态\omega\omega'的概率流:

\displaystyle J_{\omega\omega'}=W_{\omega\omega'}p_{\omega'}(t)-W_{\omega'\omega}p_{\omega}(t)=-J_{\omega'\omega}

此时主方程写为:

\displaystyle \frac{d}{dt}p_{\omega}(t)=\sum_{\omega\neq\omega'}J_{\omega\omega'}=-\sum_{\omega\neq\omega'}J_{\omega'\omega}

当系统达到稳态(steady state)时,若所有状态的概率不再随时间变化,即\frac{d}{dt}p_{\omega}(t)=0,也即对于任意的状态\omega\omega',有J_{\omega\omega'}=0,此时称系统满足细致平衡(detailed balance)条件。此时,每条转移路径上的正向流与反向流相互抵消,系统处于平衡态(equilibrium state);反之,若稳态下某些概率流非零,则系统处于非平衡稳态。平衡态是稳态的一种,但是要求更严格

由于时间关系,跳过涨落定理这部分,后面有时间再补充......

经典微分几何

界面的维度和我们所考虑的维度是n-1的关系,比如,考虑二维的体系,那么只会有一维的界面存在,要清楚的表述界面就需要一些微分几何,首先需要参数化,对于m维界面,我们引入m个实参数,将界面表示为从参数域到物理空间的映射,例如一维界面可用单参数{r}(u)表示,二维界面可用{r}(u,v)表示;以一维界面为例,当u\rightarrow u+du时,\vec{r}(u)\rightarrow \vec{r}(u+du)d\vec{r}=\vec{r}(u+du)-\vec{r}(u)=\vec{r_u}du,其中\vec{r_u}=\frac{\partial \vec{r}}{\partial u},是切向量,它指向曲线在u处的切线方向,而弧长微元ds=|\vec{r_u}|du;定义单位切向量\hat{t}(u)=\frac{\vec{r_u}}{|\vec{r_u}|}=\frac{d\vec{r}}{ds},之后定义\frac{d\hat{t}(s)}{ds}=\kappa_c\hat{n}(s),其中\kappa_c是曲率,用于衡量曲线在该点偏离直线的程度(弯曲的剧烈程度),\hat{n}(s)是单位主法向量,它指向曲线的凹侧(即弯曲的圆心方向),并且与切向量垂直,对于二维体系中的一维界面,定义了这些量便足以描述

对于二维界面,则可用{r}(u,v)表示,考虑一个二维涨落界面,其参数化表示\vec{r}(u,v),当参数发生微小变化u\rightarrow u+du,v\rightarrow v+dv时,位置矢量的微分为:

\displaystyle d\vec{r}(u,v)=\frac{\partial \vec{r}}{\partial u}du+\frac{\partial \vec{r}}{\partial v}dv=\vec{r_u}du+\vec{r_v}dv

其中\vec{r_u}\vec{r_v}是曲面的切向量,它们张成该点处的切平面,平面的方向是法向,法向量定义为:

\displaystyle \hat{n}(u,v)=\frac{\vec{r_u}\times\vec{r_v}}{|\vec{r_u}\times\vec{r_v}|}

它垂直于切平面,表征曲面的局部取向;切向量\vec{r_u}\vec{r_v}张成的平行四边形面积即为面积微元:

\displaystyle dA=|\vec{r_u}\times\vec{r_v}|dudv

弧长微元ds由第一基本形式给出:

\displaystyle (ds)^2=d\vec{r}\cdot d\vec{r}=E(du)^2+2Fdudv+G(dv)^2

其中E=\vec{r_u}\cdot\vec{r_u},F=\vec{r_u}\cdot\vec{r_v},G=\vec{r_v}\cdot\vec{r_v},弧长微元写成矩阵形式为:

\begin{pmatrix} du&dv \end{pmatrix}
\begin{pmatrix} E&F\\F&G \end{pmatrix}
\begin{pmatrix} du\\dv \end{pmatrix}

中间的矩阵即I矩阵,也称为微分几何的第一基本形式,其行列式通常记为g或者直接写为det(I),其与面积微元的关系为:

\displaystyle dA=\sqrt{g}dudv

曲面的弯曲由第二基本形式描述,即:

\displaystyle L(du)^2+2Mdudv+N(dv)^2

其中L=\hat{n}\cdot\vec{r_{uu}},M=\hat{n}\cdot\vec{r_{uv}},N=\hat{n}\cdot\vec{r_{vv}},写成矩阵形式为:

\begin{pmatrix} du&dv \end{pmatrix}
\begin{pmatrix} L&M\\M&N \end{pmatrix}
\begin{pmatrix} du\\dv \end{pmatrix}

中间的矩阵即II矩阵,也称为微分几何的第二基本形式,给定切方向\vec{r_u}du+\vec{r_v}dv后,法曲率为:

\displaystyle \frac{II}{I}=\frac{L(du)^2+2Mdudv+N(dv)^2}{E(du)^2+2Fdudv+G(dv)^2}

法曲率\kappa(u,v)依赖于切方向,若\kappa(u,v)>0,则代表曲面沿切方向朝法向量的正侧弯曲(凹),若\kappa(u,v)<0,则代表曲面沿切方向朝法向量的负侧弯曲(凸);所有切方向中,法曲率的极值称为主曲率\kappa_1\kappa_2,其方向互相垂直,;而高斯曲率K=\kappa_1\kappa_2和平均曲率H=\frac{1}{2}(\kappa_1+\kappa_2)则是点上的内蕴或不变量

在曲面某点,形状算子是一个线性变换,它将切向量映射到切向量,描述了曲面在该点沿不同方向的弯曲程度;形状算子的特征值就是主曲率\kappa_1\kappa_2,对应的特征方向就是主方向;方程det(II-\kappa I)=0是形状算子的特征方程在参数表示下的形式。它的两个解给出主曲率,通过求解久期行列式得到两个主曲率后就可以利用第一和第二基本形式表示高斯曲率和平均曲率:

\displaystyle K=\kappa_1 \kappa_2=\frac{LN-M^2}{EG-F^2}

\displaystyle H=\frac{1}{2}(\kappa_1+\kappa_2)=\frac{EN+GL-2FM}{2(EG-F^2)}

第二章 界面问题

其实前面还有好多没写完,有空来补充

Monge参数化

Monge gauge,可以叫Monge参数化或Monge规范形式或者蒙日表象,具体是指,一个曲面被表示为某个坐标函数(通常是高度)关于另外两个坐标的显式函数。最常见的形式是:

\displaystyle \vec{r}=\vec{r}(x,y,h(x,y))

或者写成分量形式:

\displaystyle \begin{cases}
x=u\\
y=v\\
z=h(u,v)
\end{cases}

曲面的位置向量由两个水平坐标(x,y)和它们所决定的高度z=h(x,y)唯一确定,Monge形式极大地简化了曲面的第一和第二基本形式的计算,让几何量(如曲率)的表达式变得非常简洁,切向量非常简单,可以写为:

\displaystyle \begin{cases}
r_x=(1,0,h_x)\\
r_y=(0,1,h_y)
\end{cases}

第一基本形式系数为:

\displaystyle \begin{cases}
E=r_x\cdot r_x=1+h_x^2\\
F=r_x\cdot r_y=h_x h_y\\
G=r_y\cdot r_y=1+h_y^2
\end{cases}

第一基本形式的行列式值为:

\displaystyle g=det(I)=EG-F^2=1+h_x^2+h_y^2=1+(\vec{\nabla}h)^2

单位法向量为:

\displaystyle \hat{n}=\frac{r_x\times r_y}{\sqrt{g}}=\frac{(-h_x,-h_y,1)}{\sqrt{1+h_x^2+h_y^2}}

第二基本形式系数为:

\displaystyle \begin{cases}
L=\hat{n}\cdot r_{xx}=\frac{h_{xx}}{\sqrt{g}}\\
M=\hat{n}\cdot r_{xy}=\frac{h_{xy}}{\sqrt{g}}\\
N=\hat{n}\cdot r_{yy}=\frac{h_{yy}}{\sqrt{g}}
\end{cases}

将上述简化的系数代入一般的高斯曲率和平均曲率公式可以得到:

\displaystyle K=\kappa_1 \kappa_2=\frac{LN-M^2}{EG-F^2}=\frac{\frac{h_{xx}h_{yy}}{g}-\frac{h_{xy}^2}{g}}{g}=\frac{h_{xx}h_{yy}-h_{xy}^2}{g^2}

\displaystyle H=\frac{EN+GL-2FM}{2(EG-F^2)}=\frac{(1+h_x^2)h_{yy}+(1+h_y^2)h_{xx}-2h_x h_y h_{xy}}{2g\sqrt{g}}

若此时考虑h_x\ll 1,h_y\ll 1,即曲面非常平坦,斜率很小;我们可以对曲率公式进行线性近似(一阶近似),二阶小量可以忽略,此时有:

\displaystyle K\approx h_{xx}h_{yy}-h_{xy}^2,H\approx \frac{h_{xx}+h_{yy}}{2}

毛细不稳定性

毛细不稳定性(Rayleigh-Plateau不稳定性)是指一个圆柱状的液体射流(或细丝)在表面张力作用下,会自发地断裂成一串液滴的现象,其背后的驱动力源于Laplace Pressure,我们通过几种不同方式来对其进行推导和解释,首先是基于能量最小化原理的简化静力学论证,考虑初始圆柱半径R_0,长度为L,侧表面积(界面面积)为A_1=2\pi R_0 L,假设破碎后液滴断裂成n个相同的球状液滴,每个半径为r_0,总表面积为A_2=n4\pi r_{0}^2,引入体积守恒,即:

\displaystyle \pi R_0^2 L=n\frac{4}{3}\pi r_{0}^3

得到n=\frac{3R_0^2L}{4r_0^3},代入总表面积:

\displaystyle A_2=n4\pi r_{0}^2=\frac{3\pi R_0^2 L}{r_0}

比较面积,若A_2<A_1则破碎后表面积减小,能量降低,导出r_0>\frac{3}{2}R_0,其物理含义为表面张力驱使系统趋向表面积最小的状态。当r_0>\frac{3}{2}R_0时,液滴总表面积小于圆柱侧面积,破碎在热力学上有利

接下来是基于统计热力学的线性稳定性分析,首先考虑一个柱状界面,其参数化为\vec{r}=(x,y,z)=(h\cos{\theta},h\sin{\theta},z)h(z)是柱状界面在轴向位置z处的半径,第一基本形式系数:E=r_z\cdot r_z=1+h_z^2G=r_{\theta}\cdot r_{\theta}=h^2F=0,面积元dA=\sqrt{g}dzd\theta=h\sqrt{1+h_z^2}dzd\theta,总表面积能:

\displaystyle H[h]=\gamma\int dA=\gamma\int_0^L dz\int_0^{2\pi} d\theta h\sqrt{1+h_z^2}

做小斜率展开,当h_z \ll1时,有:

\displaystyle H[h]\approx\gamma\int_0^L dz\int_0^{2\pi} d\theta h[1+\frac{1}{2}h_z^2]

继续可以写成:

\displaystyle \gamma[\int_0^L dz\int_0^{2\pi} d\theta h+\frac{1}{2}\int_0^L dz\int_0^{2\pi} d\theta hh_z^2]

引入平均半径,由平均半径定义\bar{h}=\frac{1}{L}\int_0^L h(z)dz,因此:

\displaystyle H[h]=2\pi \gamma L\bar{h}+\frac{1}{2}\int_0^L dz\int_0^{2\pi} d\theta hh_z^2

由于完美柱状界面的能量为H_0=2\pi\gamma Lh_0,因此能量变化为:

\displaystyle \Delta H[h]=2\pi \gamma L(\bar{h}-h_0)+\frac{1}{2}\int_0^L dz\int_0^{2\pi} d\theta hh_z^2

引入体积守恒V=\int_0^L \pi h^2dz=\pi h_0^2 L,有\int_0^L h^2dz= h_0^2 L,由于现在的界面是已经开始有涨落的界面,我们将h(z)展开,代表已经有涨落后的界面,设h(z)=h_0+\delta h,其中\delta h是微小扰动,代入体积守恒:

\displaystyle \int_0^L[h_0^2+2h_0\delta h+(\delta h)^2]dz=h_0^2L

\displaystyle h_0^2L+2h_0\int_0^L\delta hdz+\int_0^L(\delta h)^2dz=h_0^2L

\displaystyle 2h_{0}\int_{0}^{L} \delta hdz + \int_{0}^{L} (\delta h)^{2} dz = 0

可以导出:

\displaystyle \int_{0}^{L}\delta hdz=-\frac{1}{2h_{0}}\int_{0}^{L}(\delta h)^{2}dz

由于平均半径可写为:

\displaystyle \bar{h}=\frac{1}{L}\int_{0}^{L}h(z)dz=h_{0}+\frac{1}{L}\int_{0}^{L}\delta hdz

所以:

\displaystyle \bar{h}=h_0-\frac{1}{2h_0L}\int_0^L(\delta h)^2dz

最后发现:

\displaystyle \bar{h}-h_{0}=-\frac{1}{2h_{0}L}\int_{0}^{L}(\delta h)^{2}dz<0

再次回到这个式子中来:

\displaystyle \Delta H[h]=2\pi \gamma L(\bar{h}-h_0)+\frac{1}{2}\int_0^L dz\int_0^{2\pi} d\theta hh_z^2

由于\bar{h}-h_{0}<0,所以第一项为负,而第二项被积函数均为正,所以第二项为正;因此能量变化\Delta H是正是负取决于哪个项占主导,换言之,第一项主导,则会导致breaking,若第二项主导,则是稳定的

当然,这里是在轴对称假设下,h只是z的函数,可以对\theta积分,在一般情况下,hz\theta的函数,不能直接积分掉\theta,但线性稳定性分析显示,只有轴对称模式(m=0)可能导致失稳,所以研究轴对称扰动是合理的简化,如果考虑更一般的情况,即hz\theta的函数,那么此时的表面积能会变化成:

\displaystyle H[h]=\gamma\int dA=\gamma\int_0^L dz\int_0^{2\pi} d\theta h\sqrt{1+h_z^2+\frac{1}{h^2}h_{\theta}^2}

做小斜率展开,当h_z \ll1,h_{\theta} \ll1时,有:

\displaystyle H[h]=\gamma\int dA=\gamma\int_0^L dz\int_0^{2\pi}d\theta h[1+\frac{1}{2}h_z^2+\frac{1}{2}\frac{1}{h^2}h_{\theta}^2]

\displaystyle H[h]=\gamma 2\pi L\bar{h}+\frac{\gamma}{2}\int_0^L dz\int_0^{2\pi}d\theta h[h_z^2+\frac{1}{h^2}h_{\theta}^2]

引入体积守恒条件\bar{h^2}=h_0^2,前面我们是通过微扰变分法引入的,其实都是一样的,体积守恒要求平均平方半径\bar{h^2}等于未扰动的圆柱半径的平方h_0^2,其中平均被定义为:

\displaystyle \bar{(\cdot)}=\frac{1}{2\pi L}\int_0^L \int_0^{2\pi} dzd\theta (\cdot)

将界面高度展开为傅里叶级数:

\displaystyle h(z,\theta)=\sum_{m=-\infty}^{+\infty}\sum_qe^{iqz}e^{im\theta}\tilde{h}_m(q)

其中q是轴向波数,是离散值,m是整数方位角波数,系数为:

\displaystyle \tilde{h}_m(q)=\frac{1}{2\pi L}\int_0^L\int_0^{2\pi}e^{-iqz-im\theta}h(z,\theta)d\theta dz

计算\overline{h^{2}}

\displaystyle \overline{h^{2}} = \frac{1}{2\pi L} \int_{0}^{L} \int_{0}^{2\pi} h^{2}d\theta dz

代入傅里叶展开:

\displaystyle h^{2} = \sum_{q_{1},q_{2}} \sum_{m_{1},m_{2}} e^{i(q_{1}+q_{2})z} e^{i(m_{1}+m_{2})\theta} \tilde{h}_{m_{1}}(q_{1}) \tilde{h}_{m_{2}}(q_{2})

利用正交性:

\displaystyle \int_{0}^{L} e^{i(q_{1}+q_{2})z}dz = L \delta_{q_{1},-q_{2}}, \int_{0}^{2\pi} e^{i(m_{1}+m_{2})\theta}d\theta = 2\pi \delta_{m_{1},-m_{2}}

得到:

\displaystyle \bar{h^2} = \sum_{q,m} \tilde{h}_{m}(q) \tilde{h}_{-m}(-q) = \sum_{q,m} |\tilde{h}_{m}(q)|^{2}

这个求和包括所有波数qm,特别包含零模q=0,m=0,其系数\tilde{h_0}(0)=\bar{h}是平均半径;因此,我们可以将零模项分离出来:

\displaystyle \bar{h^2} = \bar{h}^{2} + \sum_{(q,m) \neq (0,0)} |\tilde{h}_{m}(q)|^{2}

由体积守恒\bar{h^2}=h_0^2,有:

\displaystyle h_{0}^{2} = \bar{h}^{2} + \sum_{q,m} |\tilde{h}_{m}(q)|^{2}

解得:

\displaystyle h_{0}=\sqrt{\overline{h}^{2}+\sum_{q,m}|\tilde{h}_{m}(q)|^{2}}=\overline{h}\sqrt{1+\frac{1}{\overline{h}^{2}}\sum_{q,m}|\tilde{h}_{m}(q)|^{2}}

假设扰动小,近似为:

\displaystyle h_{0}\approx\overline{h}\left(1+\frac{1}{2\overline{h}^{2}}\sum_{q,m}|\tilde{h}_{m}(q)|^{2}\right)=\overline{h}+\frac{1}{2\overline{h}}\sum_{q,m}|\tilde{h}_{m}(q)|^{2}

因此:

\displaystyle \overline{h}-h_{0}=-\frac{1}{2\overline{h}}\sum_{q,m}|\tilde{h}_{m}(q)|^{2}

线性近似(高斯近似)下取\bar{h} \approx h_{0},有:

\displaystyle \bar{h} - h_{0} \approx -\frac{1}{2h_{0}} \sum_{q,m} |\tilde{h}_{m}(q)|^{2}<0

面积能H=\gamma\int dA,很容易写出柱状界面的面积元并得到总面积(前面写过这里不再重复),总面积:

\displaystyle A=\int dA\approx\int_{0}^{L}\int_{0}^{2\pi}\left(h+\frac{h}{2}h_{z}^{2}+\frac{1}{2h}h_{\theta}^{2}\right)d\theta dz

改写成:

\displaystyle A=2\pi L\bar{h}+\int_{0}^{L}\int_{0}^{2\pi}\left(\frac{h}{2}h_{z}^{2}+\frac{1}{2h}h_{\theta}^{2}\right)d\theta dz

性近似(高斯近似)下:

\displaystyle A=2\pi L\bar{h}+\int_{0}^{L}\int_{0}^{2\pi}\left(\frac{h_0}{2}h_{z}^{2}+\frac{1}{2h_0}h_{\theta}^{2}\right)d\theta dz

进行傅里叶变换:

\displaystyle \int h_z^2d\theta dz=2\pi L\sum_{q,m}q^2|\tilde{h}_m(q)|^2,\int h_\theta^2d\theta dz=2\pi L\sum_{q,m}m^2|\tilde{h}_m(q)|^2

总面积:

\displaystyle A=2\pi L\bar{h}+\pi Lh_0\sum_{q,m}\left(q^2+\frac{m^2}{h_0^2}\right)|\tilde{h}_m(q)|^2

初始未扰动面积为A_0=2\pi h_0 L,故面积变化:

\displaystyle \Delta A=A-A_0=2\pi L(\bar{h}-h_0)+\pi Lh_0\sum_{q,m}\left(q^2+\frac{m^2}{h_0^2}\right)|\tilde{h}_m(q)|^2

代入\bar{h} - h_{0} = -\frac{1}{2h_{0}} \sum_{q,m} |\tilde{h}_{m}(q)|^{2},可以得到:

\displaystyle \Delta A=2\pi L\left(-\frac{1}{2h_0}\sum_{q,m}|\tilde{h}_m(q)|^2\right)+\pi Lh_0\sum_{q,m}\left(q^2+\frac{m^2}{h_0^2}\right)|\tilde{h}_m(q)|^2

合并求和项:

\displaystyle \Delta A=\frac{\pi L}{h_0}\sum_{q,m}\left(h_0^2q^2+m^2-1\right)|\tilde{h}_m(q)|^2

能量变化\Delta H=\gamma \Delta A,故:

\displaystyle \frac{\Delta H}{\gamma}=\frac{\pi L}{h_0}\sum_{q,m}\left(h_0^2q^2+m^2-1\right)|\tilde{h}_m(q)|^2

\displaystyle \Delta H=\gamma\frac{\pi L}{h_0}\sum_{q,m}\left(h_0^2q^2+m^2-1\right)|\tilde{h}_m(q)|^2

将式子写为:

\displaystyle \Delta H=\gamma\pi Lh_0\sum_{q,m}\left(q^2+\frac{1}{h_0^2}m^2-\frac{1}{h_0^2}\right)|\tilde{h}_m(q)|^2

观察稳定性系数,每个模式(m,q)对能量变化的贡献系数为:

\displaystyle L(m,q)=(q^2+\frac{1}{h_0^2}m^2-\frac{1}{h_0^2})

稳定性判据:若L(m,q)>0,模式稳定(扰动增加能量),反之若L(m,q)<0,模式不稳定(扰动降低能量),观察发现若m\neq0,由于m^2\geq1,有L(m,q)\geq0,界面总是稳定或中性,当m=0(轴对称扰动)时,L(0,q)=(q^2-\frac{1}{h_0^2}),不稳定条件为L(0,q)<0,即:

\displaystyle q^2<\frac{1}{h_0^2} \Rightarrow |q|<\frac{1}{h_0}

对应波长\lambda=2\pi/|q|>2\pi h_0,长波扰动降低系统表面能,导致柱状射流断裂成液滴;更深入的讨论是,波数q决定了轴向弯曲的剧烈程度,对应于\kappa_1的变化尺度,h_0^{-1}(初始圆周曲率) 则是系统的固有弯曲尺度,对应于\kappa_2的基准值,稳定性条件q^2<\frac{1}{h_0^2}意味着:只有当轴向的弯曲变化比圆周的固有弯曲更平缓(长波)时,系统才会失稳。如果轴向变化太剧烈(短波),增加的曲率能会抵消表面积减少的好处

润湿现象

第三章 输运与界面动力学方程

Langevin方程

Fokker–Planck方程

考虑概率密度p(x,t)的演化方程:

\displaystyle \frac{\partial}{\partial t}p(x,t)=-\nabla\cdot J(x,t)

其中概率流密度J(x,t)为:

\displaystyle J(x,t)=-D\nabla p+u(x)p=p\left[u(x)-D\nabla\log p\right]

定义局域速度v(x,t)=u(x)-D\nabla\log p,则J=pv,假设漂移速度u(x)可写成梯度形式:

\displaystyle u(x)=-\kappa\nabla H(x)

其中\kappa为迁移率,H(x)为势能函数(或者你可以理解成哈密顿量),将u(x)=-\kappa\nabla H(x)代入J得到:

\displaystyle J=-D\nabla p-\kappa\nabla H(x)p

利用爱因斯坦关系D=k_BT\kappa代入得:

\displaystyle J=-\kappa\left(k_BT\nabla p+p\nabla H\right)

提取公因子p

\displaystyle J=-\kappa p\left(\frac{k_BT}{p}\nabla p+\nabla H\right)=-\kappa p\nabla\left(H+k_BT\log p\right)

定义随机自由能(其实是有效自由能密度):

\displaystyle \hat{f}(x,t)=H(x)+k_BT\log p(x,t)

则概率流简化为:

\displaystyle J(x,t)=-p\kappa\nabla\hat{f}

于是Fokker-Planck方程可重写为:

\displaystyle \frac{\partial p}{\partial t}=-\nabla\cdot J=\nabla\cdot\left(p\kappa\nabla\hat{f}\right)

引入稳态条件,t\rightarrow\infty时,流的散度为0,即\nabla\cdot J=0,由于\kappap均大于0,故J=0等价于\nabla \hat{f}=0,因此\hat{f}=H+k_BT\log p为常数:

\displaystyle H+k_BT\log p=\text{const}\Rightarrow p(x)\propto e^{-H(x)/(k_BT)}

\frac{1}{k_BT}通常可以写为\beta,因此p(x)\propto e^{-\beta H},这就是玻尔兹曼分布

Edwards-Wilkinson方程

第四章 补充知识

随机热力学

此作者没有提供个人介绍
最后更新于 2026-01-04