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方程为：

ParseError: KaTeX parse error: Undefined control sequence: \ at position 139: …bla p}{\rho}=0 \̲ ̲\operatorname{d…

$\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)式中，有

ParseError: KaTeX parse error: Undefined control sequence: \ at position 71: … p}{\partial t}\̲ ̲\operatorname{g…

那么

$\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$，

ParseError: KaTeX parse error: Undefined control sequence: \ at position 150: …div} \mathbf{v}\̲ ̲\varepsilon(\ma…

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

ParseError: KaTeX parse error: Undefined control sequence: \ at position 81: …f{v}+\nabla p=0\̲ ̲\epsilon\frac{\…

速度和压力的表达式

(14)式的三维形式为：

ParseError: KaTeX parse error: Undefined control sequence: \ at position 80: …{\partial x}=0 \̲ ̲\frac{\partial …

$\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}$

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

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

ParseError: KaTeX parse error: Undefined control sequence: \ at position 142: …t)}{\delta t}, \̲ ̲&\frac{\partial…

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

二阶精度：

$\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$)

ParseError: KaTeX parse error: Undefined control sequence: \ at position 219: …\partial y^4}, \̲ ̲& v_3=v_2+\frac…

由式(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)$

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

ParseError: KaTeX parse error: Undefined control sequence: \ at position 238: …x}\right) \tau \̲ ̲\tilde{v}_y=v_y…

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

求解渗透率

ParseError: KaTeX parse error: Undefined control sequence: \ at position 50: … LQ}{\Delta pS}\̲ ̲& K：渗透率\ & \mu：…

流量$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. 离散化时，如何确定体素网格的计算的定义域/如何根据图像确定定义域