Code of the Week #3: fvm::laplacian(nu, U)

本周,我们讨论一下这一个语句:

fvm::laplacian(nu, U);

这一项是动量方程

\frac{\partial U}{\partial t} + \nabla\cdot(UU) = -\nabla p + \nabla\cdot\tau

的右边第二项。

考虑不可压缩牛顿流体,则有

\tau = \nu ( \nabla U +(\nabla U)^T )

由于不可压缩牛顿流体有

\nabla U = (\nabla U)^T,\, \nabla\cdot U =0

因此, \tau 的散度为

\nabla\cdot\tau = \nu \nabla^2 U

这一步具体推导可见 wiki N-S 方程推导

现在,我们主要来看一下这一项是如何离散的。把 U 看成未知量 \varphi ,设 \nu=1

\iiint_V\nabla\cdot(\nabla \varphi) d V = \oint_{\partial V}\nabla \varphi d S

可以看到,这一项是对 \varphi 的梯度(是一向量)进行求散度,高斯积分后成为对 \varphi 的梯度求面积分,离散后为

\iiint_V\nabla\cdot(\nabla \varphi) d V = \sum \nabla \varphi d\vec{S}

设该项等于零,并以 x 方向为例进行说明,网格间距均为 \Delta x ,有

\sum \nabla \varphi d\vec{S} = \left(\frac{d \varphi}{d x}\right)_A \vec{S}_A + \left(\frac{d \varphi}{d x}\right)_B \vec{S}_B=0

考虑单个 cell,以 cell 4 进行解释,考虑 S_A = S_B = S

\sum \nabla \varphi d\vec{S} = S\left(\frac{\varphi_5-\varphi_4}{\Delta x} - \frac{\varphi_4 - \varphi_3}{\Delta x}\right) \\ = \frac{S}{\Delta x} ( -2\varphi_4 +\varphi_5 + \varphi_3)=0

设左右两面边界上的值分别为 \varphi^l_0\varphi^r_0 ,按此方法,对其余网格进行离散,形成矩阵乘法,有

\left[ \begin{array}% -3 & 1 & 0 & 0\\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1\\ 0 & 0 & 1 & -3\\ \end{array} \right] \times \left[ \begin{array}% \varphi_1\\ \varphi_2\\ \varphi_3\\ \varphi_4\\ \end{array} \right] = \left[ \begin{array}% -2\varphi^l_0\\ 0\\ 0\\ -2\varphi^r_0\\ \end{array} \right]

上式与动量方程的其他项相加,就形成动量方程的离散形式。

1 Like

提问:大神,最后的矩阵第一行和最后一行看不懂,即边界处,为啥是-3,-2啊,不应是-2,,1吗???

比如,网格 1

S\biggl(\underbrace{-\frac{\varphi_1 - \varphi_0^l}{\Delta x/2}}_{左面边界} + \underbrace{\frac{\varphi_2 - \varphi_1}{\Delta x}}_{网格 1 右面} \biggr)= 0

方向两边均乘以 2,有

\frac{S}{\Delta x} (-3\varphi_1+\varphi_2+2\varphi_0^l)=0
(-3\varphi_1+\varphi_2)=-2\varphi_0^l

懂了,我还以为边界值是假象一个单元中心的值,原来是面上值,谢谢博主

%E5%9B%BE%E7%89%87
没看明白,设该项等于零的依据和目的是什么?

应该是为了方便写成等式,只是一个假设,也可以是其他任意的数值。

1 Like