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

发布于 2023-06-09  472 次阅读


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

  1. 动量守恒方程

粘性可压缩流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}

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

1693916668886.png

上图为体素化后利用有限差分法计算速度二阶导的例子,图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孔隙几何形状的预处理,第二部分是对流动建模。预处理部分会显示哪些体素中的压力与速度需要重新计算。其中对压力的重新计算只在流体和孔隙体素中进行。

三维结构的预处理

扫描图导入

1686317595984.png

问题

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

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

Everything not saved will be lost.
最后更新于 2023-09-05