理解扩散项（二）：从梯度、散度到高维扩散
在上一篇关于扩散项的文章中，我们讨论了如何从二阶导数的定义来理解扩散项，并且介绍了如何将其离散化以进行数值模拟。在这篇文章中，我们会探讨一个更具有一般性的视角，从几何意义上理解扩散，以及解释一下为什么大家常说拉普拉斯算子是梯度的散度（divergence of gradient），我会尽量用直观的方式来解释这些概念。这篇文章我会使用二维有限体积的思路来证明，也许会损失一些数学上的严谨性，但我觉得这对于建立直观理解是非常有帮助的，这个思路也可以很自然的扩展到更高维的情况，希望能帮到大家更好的理解这些奇怪的“三角形”算子。
扩散与梯度
对于扩散现象，我们很直观的理解就是物质从高浓度区域向低浓度区域扩散，所以对于某个时刻的某个点，我们想要描述物质的流动方向如何，因此我们需要一个数学工具去描述这个现象。看到这里，聪明的你一定想到了，我们可以用梯度来描述嘛，因为梯度向量的方向指向函数值由低到高增加最快的方向（关于这里感兴趣的朋友可以看看我之前写的关于梯度的文章），因此对于浓度函数来说，我们只需要增加一个负号让方向指向浓度值由高到低的方向，并且添加上一个扩散系数$D$，就得到了菲克定律中的通量表达式，以2维为例：
其中 $\phi$ 是某种密度，只要一个物理量表现出扩散行为，我们都可以用这个结构来描述它的通量，这里通量是一个向量，其方向 $\frac{\vec{J}}{|\vec{J}|}$指向了物质流动的方向，而其大小$|\vec{J}|$ 则代表了单位时间内通过单位面积的物质量，而对于不同的现象，$\phi$ 代表的物理量也不一样：
- 比如我们讨论的浓度扩散，$\phi$ 就是浓度，一般记做$c$，因此通量就是 \(\vec{J} = -D \nabla c\)
- 而在热传导现象中，$\phi$ 就是温度，一般记做$T$，因此通量就是 \(\vec{q} = -k \nabla T\)（Fourier定律）
因此我们可以看到，虽然这些现象看起来很不一样，但它们的数学结构却是非常相似的，这也是为什么我们可以用同一个方程来描述它们的原因了。
我们再展开说说通量flux这个概念，通量表示单位时间内通过单位面积的物质量，它的量纲是$[物质的量]*[时间]^{-1}*[面积]^{-1}$，所以不太严谨地说当我们得知了某点的通量以后，我们用通量乘我们想要的时间和面积，就可以得到在这个时间内通过这个面积的物质量了。
变化的区域与物理量的守恒
有了通量$\vec{J}$之后，我们已经知道在某一点物质的流动方向与强度。但某一点的信息并非是我们关心的，我们真正关心的是：
一个小区域内的物质浓度如何随时间变化？
也就是由很多个点组成的区域才是我们感兴趣的对象，而根据守恒思想，区域内物质的变化只能来自边界的流入或流出，区域内部不会凭空产生或消失物质。因此问题转化为：
如何计算穿过区域边界的净流出通量？
如图所示，我们选定一个小矩形作为我们感兴趣的区域，黑色箭头向量则是我们在边界上每个点的通量，而边界则是围成矩形的黑色边了，一个向量总能分解为法向分量和切向分量，因此我们可以把通量$\vec{J}$分解为沿着边界的蓝色部分的切向分量和垂直于边界的红色部分的法向分量，而只有垂直于边界的法向分量才会真正穿过边界，改变区域内部的物质量，而沿着边界的切向分量则不会改变区域内部的物质量，因此我们需要的就只有$\vec{J}$在法向方向上的分量：
上式我们计算了某一点对于流入流出区域的贡献，接下来我们就算整个边界的通量，因此我们需要对边界进行积分：
\[\int_{\partial \Omega} \vec{J} \cdot \vec{n} \, dS\]千万不要害怕这个积分符号，我们可以把它理解成一个求和的过程，想象一下，我们把边界分成很多小块，每一小块上都有一个通量$\vec{J} \cdot \vec{n}$，我们把这些小块上的通量加起来，就得到了整个边界上的净流出通量了。就像上图我每条边都取了三个点，可以想象我们每条边取的点越多，我们的计算结果就越精确。这里的$\partial \Omega$表示的是这个小区域的所有边界，一般对于整个区域来说我们记作$\Omega$，而$\partial \Omega$(读作partial $\Omega$，和偏导数的partial derivative的出处完全一样)就是这个区域的边界了。
关于净流出通量（net outward flux）：一个值得点在于，为什么我们叫这个量净流出通量，是因为我们定义了法线是朝着区域的外部，因此当我们积分结果为正时，表示有净流出；当积分结果为负时，表示有净流入。因此如果你的法线定义是朝着区域内部的，那么当积分结果为正时，表示有净流入；当积分结果为负时，表示有净流出。这也是什么默认我们说散度大于0时，是发散的，而散度小于0时，是汇聚的原因了。一旦我们改变了法线的定义，我们计算结果的解释也需要改变了。
散度与高维扩散
上一步给出的，是针对某个具体区域的净流出通量。
然而我们真正希望得到的，并不是依赖于区域大小的整体信息，而是一个在空间中每一点都成立的描述。
因此，我们必须将这一净流出通量归一化到单位面积（或者体积，取决于问题维度）之上，使其不再依赖于所选区域的尺度，而成为一种局部的流出率。以2维为例，我们描述矩形的面积为$\Delta x \Delta y$，那么单位体积上的流出量就是：
\[\frac{1}{\Delta x \Delta y} \int_{\partial \Omega} \vec{J} \cdot \vec{n} \, dS\]当我们让$\Delta x$和$\Delta y$趋近于0时，这个表达式就变成了一个极限：
\[\lim_{\Delta x, \Delta y \to 0} \frac{1}{\Delta x \Delta y} \int_{\partial \Omega} \vec{J} \cdot \vec{n} \, dS\]这正是散度的定义！我们把通量$\vec{J}$替换为更为通用的向量场$\vec{F}$，以及不局限于二维空间，此处$\Delta V$表示我们所选区域的体积，在2维空间中就等于面积，也就是$\Delta x\Delta y$，我们就得到了散度的定义：
\[\text{div} \vec{F} := \lim_{\Delta V \to 0} \frac{1}{\Delta V} \int_{\partial \Omega} \vec{F} \cdot \vec{n} \, dS\]接下来我们会证明这个极限形式如何等价于我们最熟悉的那个散度表达式：
\[\text{div} \vec{F} := \lim_{\Delta V \to 0} \frac{1}{\Delta V} \int_{\partial \Omega} \vec{F} \cdot \vec{n} \, dS= \frac{\partial F_1}{\partial x} + \frac{\partial F_2}{\partial y}\]proof
不想看证明的朋友也可以跳过这一小节。
首先，我们将这里要处理的向量场$\vec{F}$替换为我们在扩散项中遇到的通量$\vec{J}$，因此我们要处理的积分变成了：
\[\int_{\partial \Omega} \vec{F} \cdot \vec{n} \, ds = \int_{\partial \Omega} \vec{J} \cdot \vec{n} \, ds\]
如图所示，我们就选取整一个向量场中的一小块矩形区域来做整个分析，我们分别计算每条边界的净流出通量，然后把它们加起来就得到了整个边界的净流出通量了：
与每条边的法向量点乘后，得到对净流出通量有贡献的部分：
\[\int_{\partial \Omega} \vec{J} \cdot \vec{n} \, ds =\int_{x_0}^{x_0+\Delta x}- J_2(x, y_0)\, dx +\int_{x_0}^{x_0+\Delta x}J_2(x, y_0+\Delta y)\, dx\] \[+ \int_{y_0}^{y_0+\Delta y} - J_1(x_0, y)\, dy + \int_{y_0}^{y_0+\Delta y} J_1(x_0+\Delta x, y)\, dy\]整理得：
\[\int_{\partial \Omega} \vec{J} \cdot \vec{n} \, ds=\int_{x_0}^{x_0+\Delta x} \big[ J_2(x, y_0+\Delta y) - J_2(x, y_0) \big] dx\] \[+\int_{y_0}^{y_0+\Delta y} \big[ J_1(x_0+\Delta x, y) - J_1(x_0, y) \big] dy\]取极限：
\[\lim_{\Delta x, \Delta y \to 0} \int_{\partial \Omega} \vec{J}\cdot \vec{n}\, ds =\lim_{\Delta x, \Delta y \to 0} \left( \int_{x_0}^{x_0+\Delta x} [ J_2(x, y_0+\Delta y) - J_2(x, y_0) ] dx + \int_{y_0}^{y_0+\Delta y} [ J_1(x_0+\Delta x, y) - J_1(x_0, y) ] dy \right)\]应用积分中值定理，以及当$\Delta x, \Delta y \to 0$时，$x$和$y$都趋近于$x_0$和$y_0$，我们就得到了：
\[=\lim_{\Delta x, \Delta y \to 0} \Big( [ J_2(x_0, y_0+\Delta y) - J_2(x_0, y_0) ] \Delta x +[ J_1(x_0+\Delta x, y_0) - J_1(x_0, y_0) ] \Delta y \Big)\]两边同时除以$\Delta x \Delta y$：
\[\lim_{\Delta x, \Delta y \to 0} \frac{1}{\Delta x \Delta y} \int_{\partial \Omega} \vec{J} \cdot \vec{n} \, ds = \underbrace{ \lim_{\Delta y \to 0} \frac{J_2(x_0, y_0+\Delta y) - J_2(x_0, y_0)}{\Delta y} }_{\text{definition of } \frac{\partial J_2}{\partial y}} + \underbrace{ \lim_{\Delta x \to 0} \frac{J_1(x_0+\Delta x, y_0) - J_1(x_0, y_0)}{\Delta x} }_{\text{definition of } \frac{\partial J_1}{\partial x}}\]最终，我们得到散度的定义：
\[\operatorname{div}\vec{J} := \lim_{\Delta x, \Delta y \to 0} \frac{1}{\Delta x \Delta y} \int_{\partial \Omega} \vec{J} \cdot \vec{n} \, ds = \frac{\partial J_1}{\partial x} + \frac{\partial J_2}{\partial y}\]最后我们再把通量$\vec{J}$替换为菲克定律中的表达式$-D \nabla \phi$，就得到了我们熟悉的扩散项的高维表达：
\[\text{div} \vec{J} = \text{div}(-D \nabla \phi) = -D \left( \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} \right) = -D \nabla^2 \phi\]扩散项与运输方程
假设我们需要计算的物理量为$\phi$，当我们只考虑扩散现象时，我们就得到了一个纯扩散方程：
\[\frac{\partial \phi}{\partial t} - D \nabla^2 \phi = 0\]在这个思路下，理解运输项就变得非常自然了，在扩散现象中，驱动物质流动的是一个与梯度相关的向量场，而在运输现象中，驱动物质流动的则是一个速度场$\vec{V}$，因此我们可以把运输项理解为一个速度场的散度：
\[\text{div} (\vec{V} \phi)\]因此我们就得到了同时考虑扩散和运输的方程：
\[\frac{\partial \phi}{\partial t} + \text{div}(\vec{V} \phi) - D \nabla^2 \phi = 0\]在推导运输方程时，我们通常将扩散项整理到等式右侧。这样一来，菲克定律中原本的负号在最终表达式中便不再显式出现。这也解释了一个常见疑问：为什么对流项以正的散度形式出现，而扩散项最终却变成了正的拉普拉斯算子。因为对流现象的通量为$\vec{V} \phi$，因此在守恒形式中写作$\text{div}(\vec{V}\phi)$。而扩散现象的通量根据菲克定律为$- D \nabla \phi.$因此守恒方程的本质形式其实是
\[\frac{\partial \phi}{\partial t} + \text{div}(\vec{V}\phi) = - \text{div}(-D\nabla\phi) + S.\]当将扩散通量代入并整理后，两个负号相互抵消，便得到常见的形式：
\[\frac{\partial \phi}{\partial t} + \text{div}(V\phi) = D \nabla^2 \phi + S.\]其中 \(S\) 表示源项。当运输对象由标量 \(\phi\) 替换为动量密度，并结合流体的应力构成关系时，我们就得到了大名鼎鼎的 Navier–Stokes 方程了（忽略连续性约束）。
Key takeaways
- 根据Fick定律，我们可以用梯度来描述扩散现象中的通量，扩散由梯度驱动，物质从高浓度区域向低浓度区域流动，则通量为$\vec{J} = -D \nabla \phi$。
- 根据守恒定律，区域内物质的变化只能来自边界的流入或流出，而只有法向分量的贡献才会改变区域内部的物质量，则边界的净流出通量为$ \int_{\partial \Omega} \vec{J} \cdot \vec{n} \, dS $。
- 通过边界积分计算净流出通量，并将其归一化到单位体积，我们得到了散度的定义，散度描述了一个向量场在空间中每一点的发散（流出流入多少）程度，$\text{div} \vec{F} := \lim_{\Delta V \to 0} \frac{1}{\Delta V} \int_{\partial \Omega} \vec{F} \cdot \vec{n} \, dS$。
- 最后，我们将通量表达式代入散度定义，得到了扩散项的高维表达式，证明了拉普拉斯算子与扩散现象之间的联系，$\text{div} \vec{J} = \text{div}(-D \nabla \phi) = -D \nabla^2 \phi$。这就是大家常说的，Laplacian是divergence of gradient的原因了.
☕ Happy learning journey~ 🛠️