Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

文献复现——基于有限差分的Stokes求解器

理论部分

多孔介质中单向流动模型

粘性不可压缩N-S方程

1. 质量守恒方程

已知可压缩流体的N-S质量守恒方程为：

\frac{\mathrm{D}\rho}{\mathrm{D}t}+\rho\nabla\cdot V=0,\tag{1}

对于不可压缩流，\rho_l=constant，因此不可压缩流的质量守恒方程为：\nabla\cdot V=0

2. 动量守恒方程

粘性可压缩流x方向动量方程为：

\rho\frac{\mathrm{D}u}{\mathrm{D}t}=-\frac{\partial p}{\partial x}+2\mu\frac{\partial^{2}u}{\partial x^{2}}+\mu\frac{\partial}{\partial y}\biggl(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\biggr)+\mu\frac{\partial}{\partial z}\biggl(\frac{\partial u}{\partial z}+\frac{\partial u}{\partial x}\biggr)+\rho f_{z}\tag{2}

结合不可压缩流的质量守恒方程，并忽略体积力，可得到粘性不可压缩流x方向动量方程：

\rho\frac{\mathrm{D}u}{\mathrm{D}t}=-\frac{\partial p}{\partial x}+\mu\nabla^{2}u\tag{3}

展开物质导数，

\frac{Du}{Dt}=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}

(3)式可变为：

\rho\frac{\partial \mathbf{u}}{\partial \mathbf{t}}+\rho u\frac{\partial \mathbf{u}}{\partial \mathbf{x}}={\mu} \Delta \mathbf{u}-{\nabla p}\tag{4}

整理即为

\rho\frac{\partial \mathbf{u}}{\partial \mathbf{t}}+\rho (\mathbf{u} \nabla) \mathbf{u}={\mu} \Delta \mathbf{u}-{\nabla p}\tag{5}

从而可得到粘性不可压缩流体的N-S方程为：

\frac{\partial \mathbf{v}}{\partial \mathbf{t}}+(\mathbf{v} \nabla) \mathbf{v}-\frac{\mu}{\rho} \Delta \mathbf{v}+\frac{\nabla p}{\rho}=0 \
\operatorname{div}\mathbf{v}=0\tag{6}

\mathbf{v}：速度场,\quad(\mathbf{v} \nabla)\mathbf{v}：对流项,\quad\mu：流体的动力粘度
,\quad\mu\Delta\mathbf{v}：粘性项,\quad p：压力场,\quad\nabla p：压力梯度

多孔介质中的流体流动，其属于小雷诺数流动，雷诺数Re=\frac{\rho \mathbf{U}L}{\mu}<0且<<1，为一个常数。

考虑到可压缩流体的连续性方程为：

\frac{\partial \rho}{\partial t}+\nabla \cdot (\rho u)=\frac{\partial \rho}{\partial t}+\rho \operatorname{div} \mathbf{v}+\mathbf{v}\operatorname{\nabla} \rho=0\tag{12}

假设\rho_0=1，将式(11)中密度的线性关系带入(12)式中，有

\frac{\partial \rho}{\partial t}=\epsilon\frac{\partial p}{\partial t}\
\operatorname{grad} \rho=\frac{\partial \rho}{\partial x} \mathbf{i}+\frac{\partial \rho}{\partial y} \mathbf{j}+\frac{\partial \rho}{\partial z} \mathbf{k}=\varepsilon\left(\frac{\partial p}{\partial x} \mathbf{i}+\frac{\partial p}{\partial y} \mathbf{j}+\frac{\partial p}{\partial z} \mathbf{k}\right)=\varepsilon \operatorname{\nabla} p

那么

\mathbf{v} \operatorname{\nabla} \rho=\varepsilon(\mathbf{v} \operatorname{\nabla} p)

将以上结果带入(12)中，得到

\epsilon\frac{\partial p}{\partial t}+\rho \operatorname{div} \mathbf{v}+\varepsilon(\mathbf{v} \operatorname{\nabla} p)=0\tag{13}

由于\epsilon<<1，\rho_0=1，

\rho \operatorname{div}\mathbf{v}=(\rho_0+\epsilon p)\operatorname{div} \mathbf{v}=\rho_0 \operatorname{div}\mathbf{v}=\operatorname{div} \mathbf{v}\
\varepsilon(\mathbf{v} \operatorname{\nabla} p)=0\
\rho \frac{\partial \mathbf{v}}{\partial \mathrm{t}}=(\rho_0+\epsilon p)\frac{\partial \mathbf{v}}{\partial \mathrm{t}}=\rho_0 \frac{\partial \mathbf{v}}{\partial \mathrm{t}}=\frac{\partial \mathbf{v}}{\partial \mathrm{t}}

(7)式可化简为(14)式，我们根据该方程组来求解速度和压力

\frac{\partial \mathbf{v}}{\partial \mathrm{t}}-\mu \Delta \mathbf{v}+\nabla p=0\
\epsilon\frac{\partial p}{\partial t}+\operatorname{div} \mathbf{v}=0
\tag{14}

速度和压力的表达式

(14)式的三维形式为：

\frac{\partial v_x}{\partial t}-\mu \Delta v_x+\frac{\partial p}{\partial x}=0 \
\frac{\partial v_y}{\partial t}-\mu \Delta v_y+\frac{\partial p}{\partial y}=0 \
\frac{\partial v_z}{\partial t}-\mu \Delta v_z+\frac{\partial p}{\partial z}=0
\tag{15}

\Delta v_x=\frac{\partial^2 v_x}{\partial x^2}+\frac{\partial^2 v_x}{\partial y^2}+\frac{\partial^2 v_x}{\partial z^2}

需要对速度的偏导数和压力的偏导数作有限差分，从而解出数值解。

分别对速度和压力进行差分(一阶精度)，

\begin{gathered}
&\frac{\partial v_i}{\partial t}\left(t_0\right) \approx \frac{v_i\left(t_0+\delta t\right)-v_i\left(t_0\right)}{\delta t}, \
&\frac{\partial p}{\partial t}\left(t_0\right) \approx \frac{p\left(t_0+\delta t\right)-p\left(t_0\right)}{\delta t}, \
&\frac{\partial p}{\partial x}\left(x_0\right) \approx \frac{p\left(x_0+\delta x\right)-p\left(x_0\right)}{\delta x}, \frac{\partial p}{\partial y} \text { 和 } \frac{\partial p}{\partial z}\text { 类似之 }\
&\operatorname{div} \mathbf{v}\left(\mathbf{x}_0\right)=\frac{\partial v_x}{\partial x}+\frac{\partial v_y}{\partial y}+\frac{\partial \nu_z}{\partial z} \approx \
&\approx \frac{v_x\left(x_0+\delta x\right)-v_x\left(x_0\right)}{\delta x}+\frac{v_y\left(y_0+\delta y\right)-v_y\left(y_0\right)}{\delta y}+\frac{v_z\left(z_0+\delta z\right)-v_z\left(z_0\right)}{\delta z}
\end{gathered}

速度二阶导的差分(几阶精度差分方程截断误差就是几阶，推导在于凑数消去截断误差)：

二阶精度：

\frac{\partial^2 v_x}{\partial x^2}=\frac{v_x(x_0-\delta x)-2v_x(x_0)+v_x(x_0+\delta x)}{(\delta x)^2}\tag{16}

四阶精度：

\frac{\partial^2 v_x}{\partial x^2}=\frac{-v_x(x_0-2\delta x)+16v_x(x_0-\delta x)-30v_x(x_0)+16v_x(x_0+\delta x)-v_x(x_0+2\delta x)}{12(\delta x)^2}\tag{17}

由上面的二阶精度和四阶精度方案我们可以知道，要求解速度的二阶导，对于二阶精度方案需要临近的两个体素，对于四阶精度方案需要临近的四个体素。

上图为体素化后利用有限差分法计算速度二阶导的例子，图a为四阶精度，图b为二阶精度，图c为从A孔隙体素位置计算速度二阶导的例子，对于二阶精度有6种组合，四阶精度有15种组合

黑色的代表固相体素，绿色的代表正在被计算的孔隙体素，蓝色的代表当前没有被计算的孔隙体素

对于图c中A位置的体素，其速度为v_2，其左侧为代表固相的体素，其和固相体素交界面的上的速度为v_1，其右侧有四个孔隙体素，分别对应速度v_3,v_4,v_5,v_6，在四阶精度下，从A处分别作差分求出其他位置的速度(\delta x=\delta y=\delta z=1)

\begin{aligned}
& v_1=v_2-\frac{1}{2} \frac{\partial v_2}{\partial y}+\frac{1}{8} \frac{\partial^2 v_2}{\partial y^2}-\frac{1}{48} \frac{\partial^3 v_2}{\partial y^3}+\frac{1}{384} \frac{\partial^4 v_2}{\partial y^4}, \
& v_3=v_2+\frac{\partial v_2}{\partial y}+\frac{1}{2} \frac{\partial^2 v_2}{\partial y^2}+\frac{1}{6} \frac{\partial^3 v_2}{\partial y^3}+\frac{1}{24} \frac{\partial^4 v_2}{\partial y^4}, \
& v_4=v_2+2 \frac{\partial v_2}{\partial y}+2 \frac{\partial^2 v_2}{\partial y^2}+\frac{4}{3} \frac{\partial^3 v_2}{\partial y^3}+\frac{2}{3} \frac{\partial^4 v_2}{\partial y^4},\
& v_5=v_2+3 \frac{\partial v_2}{\partial y}+\frac{9}{2} \frac{\partial^2 v_2}{\partial y^2}+\frac{9}{6} \frac{\partial^3 v_2}{\partial y^3}+\frac{27}{8} \frac{\partial^4 v_2}{\partial y^4}, \
& v_6=v_2+4 \frac{\partial v_2}{\partial y}+8 \frac{\partial^2 v_2}{\partial y^2}+\frac{32}{3} \frac{\partial^3 v_2}{\partial y^3}+\frac{32}{3} \frac{\partial^4 v_2}{\partial y^4} .
\end{aligned}

由式(10)可知边界速度为0，结合四阶精度方案的二阶导数，导出 \frac{\partial^2 v_x}{\partial y^2}

\frac{\partial^2 v_2}{\partial y^2}=-\frac{65}{12} v_2+\frac{22}{9} v_3-\frac{1}{2} v_4+\frac{2}{21} v_5-\frac{1}{108} v_6
写成有限差分形式：

\frac{\partial^2 v_x\left(x_0\right)}{\partial y^2} \approx-\frac{65}{12} v_x\left(x_0\right)+\frac{22}{9} v_x\left(x_0+\delta y\right)-\frac{1}{2} v_x\left(x_0+2 \delta y\right)+\frac{2}{21} v_x\left(x_0+3 \delta y\right)-\frac{1}{108} v_x\left(x_0+4 \delta y\right)

现在，把上面求出来的速度的偏导数、压力的偏导数的差分形式带入到连续性方程中：

\begin{gathered}
\tilde{v}_x=v_x+\left(\mu\left(\frac{\partial^2 v_x}{\partial x^2}+\frac{\partial^2 v_x}{\partial y^2}+\frac{\partial^2 v_x}{\partial z^2}\right)-\frac{p\left(x_0+\delta x\right)-p\left(x_0\right)}{\delta x}\right) \tau \
\tilde{v}_y=v_y+\left(\mu\left(\frac{\partial^2 v_y}{\partial x^2}+\frac{\partial^2 v_y}{\partial y^2}+\frac{\partial^2 v_y}{\partial z^2}\right)-\frac{p\left(y_0+\delta y\right)-p\left(y_0\right)}{\delta y}\right) \tau \
\tilde{v}_z=v_z+\left(\mu\left(\frac{\partial^2 v_z}{\partial x^2}+\frac{\partial^2 v_z}{\partial y^2}+\frac{\partial^2 v_z}{\partial z^2}\right)-\frac{p\left(z_0+\delta z\right)-p\left(z_0\right)}{\delta z}\right) \tau \
\tilde{p}=p-\frac{1}{\varepsilon}\left(\frac{v_x\left(x_0+\delta x\right)-v_x\left(x_0\right)}{\delta x}+\frac{v_y\left(y_0+\delta y\right)-v_y\left(y_0\right)}{\delta y}+\frac{v_z\left(z_0+\delta z\right)-v_z\left(z_0\right)}{\delta z}\right) \tau
\end{gathered}

\tilde{v}_x ,\tilde{v}_y,\tilde{v}_z,\tilde{p}是下一个时间步数求出的v_x,v_y,v_z,p

求解渗透率

\begin{aligned}
& 达西定律：K=\frac{\mu LQ}{\Delta pS}\
& K：渗透率\
& \mu：动力粘度\
& L：压力差作用范围\
& Q：流过截面的流量\
\end{aligned}

流量Q可表示为：Q=\langle \mathbf{v} \rangle S_{eff}，\langle \mathbf{v} \rangle为每单位面积的流体的平均速度，S_{eff}为截面对应的有效面积

则达西定律可改写为K=\frac{\mu L\langle \mathbf{v} \rangle S_{eff}}{\Delta pS}，这样，在求解出速度场和已知压力边界条件的情况下，利用达西定律就可以求出渗透率

对于不可压缩流体，其流速处处相等。如果给定压力边界条件\frac{\Delta p}{L}=grad \ p=1，达西定律可进一步简化为：K=\frac{\mu \langle \mathbf{v} \rangle S_{eff}}{S}。

对于实际求解，我们将面积用体积替换，K=\frac{\mu \langle \mathbf{v} \rangle S_{eff}}{S}，这样可以平滑计算时的噪声。

算法结构和误差准则

算法分为两个部分，第一部分是3D孔隙几何形状的预处理，第二部分是对流动建模。预处理部分会显示哪些体素中的压力与速度需要重新计算。其中对压力的重新计算只在流体和孔隙体素中进行。

三维结构的预处理

扫描图导入

问题

1. 如何实现扫描图的导入并复现结构

2. 离散化时，如何确定体素网格的计算的定义域/如何根据图像确定定义域