数值计算 | 曲面(线)上法向量的散度是什么？

写在前面

在含有界面的流体流动模拟中，一般论文中常用的表面张力模型表达式为：

\Delta p=\gamma\nabla\cdot\mathbf{n}

其中，

\gamma

为表面张力系数，而在一些研究论文中通常将

\nabla\cdot\mathbf{n}

描述为曲率。

那么它究竟是什么呢？在本文中，将详细推导曲面（或曲线）上法向量的散度。

问题说明

在二维或三维空间中，定义曲线

h:\mathbb{R}\to\mathbb{R}

和曲面

f:\mathbb{R}^2\to\mathbb{R}

分别为

y=h(x)

和

z=f(x,y)

。假定这些函数具有连续二阶导数。

在二维空间中，曲线上法向量的散度等于曲线的曲率(Curvature)，也就是

\nabla\cdot\mathbf{n}=\kappa

；

在三维空间中，表面上法向量的散度等于主曲率(Principal Curvature)的和，主曲率的定义是包含法向量的平面与曲面交线的曲率的最大值和最小值，简单来说就是

\nabla\cdot\mathbf{n}=\kappa_1+\kappa_2

。

二维证明

根据第一节的介绍，我们可以定义曲线方程

F(x,y)=y-h(x)

，则对应（外）法向量为

\mathbf{n}=\frac{\nabla F(x,y)}{\left|\nabla F(x,y)\right|} = \frac{1}{\sqrt{1+h^2_x(x)}}\left[-h_x(x),1\right] \tag{2-1}

对法向量求散度可得

\begin{aligned} \nabla\cdot\mathbf{n} &= \frac{\partial\mathbf{n}_x}{\partial x} + \frac{\partial\mathbf{n}_y}{\partial y} \&= \frac{ -h_{xx}(x)\sqrt{1+h^2_x(x)} + h_x(x)\frac{1}{2}\frac{2h_x(x)h_{xx}(x)}{\sqrt{1+h^2_x(x)}} }{1+h^2_x(x)} + 0\&= \frac{-h_{xx}(x)}{\left[1+h^2_x(x)\right]^{\frac{3}{2}}} \end{aligned}\tag{2-2}

而多数微积分教科书中将曲线的曲率定义为正数

\kappa = \left|\frac{\mathrm{d}\phi}{\mathrm{d}s}\right|\tag{2-3}

其中，

\phi

为曲线切向量与水平线之间的夹角，

s

为曲线弧长。

但曲率值的正负表征的是曲线弯曲的方向，因此在曲率公式的推导过程中，将绝对值符号去掉，

\kappa = \frac{\mathrm{d}\phi}{\mathrm{d}s} \tag{2-4}

根据链式求导法则可得

\frac{\mathrm{d}\phi}{\mathrm{d}x} = \frac{\mathrm{d}\phi}{\mathrm{d}s}\frac{\mathrm{d}s}{\mathrm{d}x} \tag{2-5}

根据上图曲线的几何属性可知，夹角

\phi

为

\phi = \arctan\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right) = \arctan\left[h_x(x)\right] \tag{2-6}

将式（2-6）代入式（2-5）可得

\frac{\mathrm{d}\phi}{\mathrm{d}x} = \frac{h_{xx}(x)}{1+h^2_x(x)} \tag{2-7}

而弧长

s

关于坐标

x

的导数为

\frac{\mathrm{d}s}{\mathrm{d}x} = \frac{\mathrm{d}}{\mathrm{d}x} \int_0^x \sqrt{1+h^2_x(x)}\mathrm{d}\xi = \sqrt{1+h^2_x(x)} \tag{2-8}

综合式（2-7）和式（2-8），可得曲率表达式为

\kappa = \frac{{\mathrm{d}\phi}/{\mathrm{d}x}}{{\mathrm{d}s}/{\mathrm{d}x}} = \frac{h_{xx}(x)}{1+h^2_x(x)} \frac{1}{\sqrt{1+h^2_x(x)}} = \frac{h_{xx}(x)}{\left[1+h^2_x(x)\right]^{\frac{3}{2}}} \tag{2-9}

至此，式（2-2）与式（2-9）在形式上是一样的，只相差一个正负号。考虑将凸的曲线的曲率定义为正值，而这种曲线满足

h_{xx}(x) \leqslant 0

因此，完整定义为

\kappa = \frac{-h_{xx}(x)}{\left[1+h^2_x(x)\right]^{\frac{3}{2}}}\tag{2-10}

根据式(2-2)和式（2-10）可得二维情况下

\nabla\cdot\mathbf{n} = \kappa

三维证明

在三维空间中，曲面上任意一点

x_0

处存在多个曲率值。事实上在曲面上该点处，任意一个通过点

x_0

并包含该点处法向量的平面与

x_0

所在曲面的交线都是一条二维曲线，也都有相应的曲率。我们期望计算特定方向上的曲率，通过改变方向计算得到最大和最小的曲率值，也就是主曲率。直接计算在数学上会非常麻烦，为了简化代数运算，我们首先记住一点，下面等式在笛卡尔直角坐标系中与坐标轴位置无关。因此，我们可以通过旋转和平移坐标系以方便分析。

我们先在想要计算曲率和法向量散度的地方固定一个点，记为

x_0=(x,y,f(x,y))

，并在该处建立局部坐标系

x_0-uvw

。接下来我们通过变换坐标轴使得

ux_0v

平面与

x_0

处切平面重合，并且

x_0w

轴法向量同向（此时也重合了），如下图所示。

因而，在局部坐标系

x_0-uvw

下有如下等式成立：

\begin{aligned} f(0,0) &= 0 \f_u(0,0) &=0 \f_v(0,0) &= 0\\end{aligned} \tag{3-1}

其中，坐标系

x_0-uvw

中点

(0,0,0)

对应坐标系

O-xyz

中

x_0

点。

因此，通过法向量（也就是

x_0w

轴）和

ux_0v

平面上的一条直线

v=ru

可以定义一个平面，该平面与曲面

f

的交线在局部坐标系

x_0-uvw

中的表达式为

w = f(u,ru) \tag{3-2}

其中，

r

代表了我们所选的方向。

要计算式（3-2）所标示曲线的曲率，我们必须计算曲面

f

在参数

r

所确定方向上的微分，我们可以用线段长度参数

t

来表示曲面

f

，令

t^2 = u^2 + v^2 = u^2\left(1+r^2\right) \tag{3-3}

则可以得到

\begin{aligned} u^2 &= \frac{t^2}{1+r^2}\\frac{\mathrm{d}u}{\mathrm{d}t} &= \frac{1}{\sqrt{1+r^2}}\w = f\left(u,ru\right) &= f\left(t/\sqrt{1+r^2},\ rt/\sqrt{1+r^2}\right)\\end{aligned} \tag{3-4}

与上一节原理类似，我们同样根据式（2-4）计算参数

r

方向对应的曲率，也就是

\frac{\mathrm{d}\phi}{\mathrm{d}t} = \frac{\mathrm{d}\phi}{\mathrm{d}s}\frac{\mathrm{d}s}{\mathrm{d}t} =\kappa_r \frac{\mathrm{d}s}{\mathrm{d}t}

通过几何关系可以知道夹角

\phi

的表达式为

\phi = \arctan\left(\frac{\mathrm{d}w}{\mathrm{d}t}\right) = \arctan\left[\frac{f_u(u,ru)+rf_v(u,ru)}{\sqrt{1+r^2}}\right] \tag{3-5}

故，夹角

\phi

关于线段长度

t

的导数为

\begin{aligned} \frac{\mathrm{d}\phi}{\mathrm{d}t} &= \frac{1}{1+\frac{\left[f_u(u,ru)+rf_v(u,ru)\right]^2}{1+r^2}} \frac{f_{uu}(u,ru)+2rf_{uv}(u,ru)+r^2f_{vv}(u,ru)}{1+r^2}\&= \frac{f_{uu}(u,ru)+2rf_{uv}(u,ru)+r^2f_{vv}(u,ru)}{1+r^2+\left[f_u(u,ru)+rf_v(u,ru)\right]^2}\\end{aligned}

而交线式（3-2）上弧长

s

关于

t

的导数为

\begin{aligned} \frac{\mathrm{d}s}{\mathrm{d}t} &= \frac{\mathrm{d}}{\mathrm{d}t} \int_0^u \sqrt{ 1+\frac{\left[f_u(u,ru)+rf_v(u,ru)\right]^2}{1+r^2} }\sqrt{1+r^2}\mathrm{d}\xi\&= \sqrt{ 1+\frac{\left[f_u(u,ru)+rf_v(u,ru)\right]^2}{1+r^2} }\\end{aligned}

综合上述两式可得在参数

r

方向对应的曲率的表达式为

\begin{aligned} \kappa_r =& \frac{\mathrm{d}\phi/\mathrm{d}t}{\mathrm{d}s/\mathrm{d}t} \\= & \frac{f_{uu}(u,ru)+2rf_{uv}(u,ru)+r^2f_{vv}(u,ru)}{1+r^2+\left[f_u(u,ru)+rf_v(u,ru)\right]^2}\times\\& \sqrt{ \frac{1+r^2}{1+r^2+\left[f_u(u,ru)+rf_v(u,ru)\right]^2} } \end{aligned}

而我们想要求

x_0

处处的曲率，因此将

t=0

代入，并根据式（3-1）可得

\kappa_r = \frac{f_{uu}(0,0)+2rf_{uv}(0,0)+r^2f_{vv}(0,0)}{1+r^2} \tag{3-6}

为方便计算，我们引入常量符号

a=f_{uu}(0,0)

b=2f_{uv}(0,0)

，

c=f_{vv}(0,0)

，则此时计算

\kappa_r

的导数为

\begin{aligned} \frac{\mathrm{d}\kappa_r}{\mathrm{d}r} &= \frac{\mathrm{d}}{\mathrm{d}r} \left(\frac{a+br+cr^2}{1+r^2}\right) \& = \frac{b+2cr}{1+r^2} - \frac{2r\left(a+br+cr^2\right)}{\left(1+r^2\right)^2}\\ &= \frac{ b+2cr+br^2+2cr^3 -2ar-2br^2-2cr^3}{\left(1+r^2\right)^2} \&= \frac{-br^2+\left(2c-2a\right)r+b}{\left(1+r^2\right)^2}\\end{aligned}

根据上式可知该式的临界点（Critical Point）即为曲率的最大值点和最小值点，于是我们将得到下面关于

r

的一元二次方程

-br^2+\left(2c-2a\right)r+b = 0 \tag{3-7}

假设

r_0

满足该方程，那么我们将

-1/r_0

代入可得

\begin{aligned} &\frac{-b}{r_0^2} - \frac{2c-2a}{r_0} + b \\=& \frac{-1}{r_0^2}\left[b+\left(2c-2a\right)r_0-br_0^2\right] \\=& 0 \end{aligned}

根据上式推导可知，

-1/r_0

也满足该方程。这也表明主曲率所在的方向相互垂直。我们直接将

r

和

-1/r

代入可得

\begin{aligned} \kappa_1 + \kappa_2 &= \frac{a+br+cr^2}{1+r^2} + \frac{a-b/r+c/r^2}{1/r^2+1}\&= \frac{a+br+cr^2+ar^2-br+c}{1+r^2} \&= \frac{\left(a+c\right)\left(1+r^2\right)}{1+r^2} \&= a+c = f_{uu}(0,0) + f_{vv}(0,0)\\end{aligned}

与二维情况相同，我们也需要选择一个符号来表示凸的曲面的正曲率，与二维同理可得

\kappa_1 + \kappa_2 = -f_{uu}(0,0) - f_{vv}(0,0) \tag{3-8}

现在，我们在局部坐标系

x_0-uvw

下计算面法向量的散度为

\begin{aligned} \nabla \cdot \mathbf{n} &= \nabla \cdot \left\{\frac{\left(-f_u,-f_v,1\right)}{ \sqrt{f^2_u+f^2_v+1} }\right\}\&= \frac{ -f_{uu}-f_{vv}}{ \sqrt{f^2_u+f^2_v+1} } +\left(-f_u,-f_v,1\right)\cdot \frac{-\left(f_uf_{uu}+f_vf_{vu},f_uf_{uv}+f_vf_{vv},0\right)}{\left(f^2_u+f^2_v+1\right)^{\frac{3}{2}}}\&=\frac{ -f_{uu}-f_{vv}}{ \sqrt{f^2_u+f^2_v+1} } +\frac{f_u\left(f_uf_{uu}+f_vf_{vu}\right)+f_v\left(f_uf_{uv}+f_vf_{vv}\right)}{\left(f^2_u+f^2_v+1\right)^{\frac{3}{2}}} \\end{aligned} %\tag{3-13}

因此，根据式（3-1）所述性质，在

(0,0)

处有

\nabla\cdot \mathbf{n} = -f_{uu}(0,0)-f_{vv}(0,0) \tag{3-9}

根据式（3-8）式（3-9）可得三维情况下

\nabla \cdot \mathbf{n} = \kappa_1 + \kappa_2

证毕。

参考资料

A Proof that the Divergence of a Surface Normal Is Equal to the sum of the Principal Curvatures
https://willson.cm.utexas.edu/Research/Sub_Files/Surface_Phenomena/Spring%202000/surface_normal_proof.pdf#:~:text=Then%20in%20t%20w%20o%20dimensions%2C%20the%20div,to%20surface%20is%20equal%20sum%20principal%20curv%20atures.

Nicholas M. Patrikalakis, Takashi Maekawa, Wonjoon Cho. Shape Interrogation for Computer Aided Design and Manufacturing:Principal normal and curvature
https://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node23.html

Nicholas M. Patrikalakis, Takashi Maekawa, Wonjoon Cho. Shape Interrogation for Computer Aided Design and Manufacturing:Curvature and curvature vector
https://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node119.html#sec:dgint_trans_curv

请问如何根据二阶导判断函数凹凸性?
https://www.zhihu.com/question/349019716

Relation between divergence of unit normal and radius of curvature
https://physics.stackexchange.com/questions/698035/relation-between-divergence-of-unit-normal-and-radius-of-curvature

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。