# 基于有限差分法的用于求解三维孔隙模型的N-S求解器

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

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

## 理论部分

### 多孔介质中单向流动模型

#### 粘性不可压缩N-S方程

1. 质量守恒方程

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

1. 动量守恒方程

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

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

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

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

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

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

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

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

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

# 问题

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

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

Everything not saved will be lost.