31、利用并行技术计算流体动力学过程

利用并行技术计算流体动力学过程

1. 数值算法

为求解相关方程组和附加条件，采用了按物理过程分裂的方法和有限体积网格法。分裂后的方程组形式如下：
- 方程（8）：$S_f^w / \rho_f^w [\varphi_f \rho_f^w]_t + (1 - S_f^w) / \rho_f^o [\varphi_f \rho_f^o]_t + DIG_f = 0$，其中$DIG_f = 1 / \rho_f^w div(\rho_f^w \vec{U}_f^w) + 1 / \rho_f^o div(\rho_f^o \vec{U}_f^o) + q_f^w / \rho_f^w + q_f^o / \rho_f^o$
- 方程（9）：$S_m^w / \rho_m^w [\varphi_m \rho_m^w]_t + (1 - S_m^w) / \rho_m^o [\varphi_m \rho_m^o]_t + DIG_m = 0$，其中$DIG_m = q_m^w / \rho_m^w + q_m^o / \rho_m^o$
- 方程（10）：$\varphi_f(S_f^w \rho_f^w \frac{\partial \varepsilon_f^w}{\partial t} + (1 - S_f^w) \rho_f^o \frac{\partial \varepsilon_f^o}{\partial t}) + \varphi_m(S_m^w \rho_m^w \frac{\partial \varepsilon_m^w}{\partial t} + (1 - S_m^w) \rho_m^o \frac{\partial \varepsilon_m^o}{\partial t}) + \frac{\partial}{\partial t}(1 - \varphi_f - \varphi_m) \rho_s \varepsilon_s + DIG_f^{\varepsilon} + DIG_m^{\varepsilon} + div \vec{W}_s = 0$

按物理过程分裂的算法，区分了方程（8）和（9），它们分别代表了裂缝和基质中，水和油质量组分组合在无饱和度转移情况下系统热力学参数的压导演化。而方程（10）则是整个系统（裂缝、孔隙基质和骨架传热）总能量的压导演化。

接着，使用有限体积法在合适的空间网格上对方程组进行近似。在时间上，采用具有恒定步长的隐式差分近似。对于二维情况，可在(r, z)几何下对计算域的径向层进行二维分层离散划分，从而对相关方程进行二维差分近似；对于一维几何，可使用已知的网格对应方法进行近似。最后，将离散方程线性化，得到如下形式的线性方程组：
[
\begin{cases}
C_0 y_0 - B_0 y_1 = \Phi_0, & k = 0 \
- A_k y_{k - 1} + C_k y_k - B_k y_{k + 1} = \Phi_k, & 1 \leq k \leq N - 1 \
- A_N y_{N - 1} + C_N y_N = \Phi_N, & k = 1
\end{cases}
]
其中，$A$、$B$、$C$是 2×2 矩阵（第一行负责压导过程，第二行负责热导过程），$y$是所需的压力$\delta P$和温度$\delta T$，$\Phi$是二维向量。

2. 矩阵扫描法

前面得到的线性代数方程组由块三对角矩阵定义：
[
A_k =
\begin{pmatrix}
A_{11}^{pk} & A_{12}^{pk} \
A_{21}^{\varepsilon k} & A_{22}^{\varepsilon k}
\end{pmatrix},
B_k =
\begin{pmatrix}
B_{11}^{pk} & B_{12}^{pk} \
B_{21}^{\varepsilon k} & B_{22}^{\varepsilon k}
\end{pmatrix},
C_k =
\begin{pmatrix}
C_{11}^{pk} & C_{12}^{pk} \
C_{21}^{\varepsilon k} & C_{22}^{\varepsilon k}
\end{pmatrix},
\Phi_k =
\begin{pmatrix}
\Phi_{pk} \
\Phi_{\varepsilon k}
\end{pmatrix}
]
这里$A_{12}^{pk} = A_{21}^{pk} = 0$，$B_{12}^{pk} = B_{21}^{pk} = 0$。

在初始测试模型和得到的矩阵系数时，当时间步长$\tau \to 0$和空间步长$h \to 0$，有$det ||C_k|| > 0$。若满足特定条件，模型会退化为之前描述的特殊情况。

为求解方程组（8） - （10），采用矩阵扫描法，该方法类似于等温情况下标量三点方程的扫描法。具体步骤如下：
- 前向步骤确定扫描系数：
- $\alpha_{k + 1} = (C_k - A_k \alpha_k)^{-1} B_k$，$k = 1, 2, \cdots, N - 1$，$\alpha_1 = C_0^{-1} B_0$
- $\beta_{k + 1} = (C_k - A_k \alpha_k)^{-1} (F_k + A_k \beta_k)$，$k = 1, 2, \cdots, N$，$\beta_1 = C_0^{-1} F_0$
- 后向步骤求解系统：
- $y_k = \alpha_{k + 1} y_{k + 1} + \beta_{k + 1}$，$k = N - 1, N - 2, \cdots, 0$，$y_N = \beta_{N + 1}$

实现该方法需要矩阵的加法、乘法、转置等操作，这些操作被实现为单独的函数，矩阵扫描法本身也被封装为一个单独的函数，以便在计算的每一层（迭代）中调用。

矩阵扫描并行化算法基于一个已知的算法，主要计算在 K100 超级计算机上进行，使用 C 语言和 MPI 标准实现。使用的 MPI 函数包括：MPI Init、MPI Comm size、MPI Comm rank、MPI Get processor name、MPI Bcast（模型参数分布）、MyRange（计算域分解）、MPI Isend/MPI Irecv（相邻进程间的数据交换）、MPI Allreduce。

并行化的一个重要步骤是解决一个短问题，所有 MPI 进程会进行系数的集体交换。与并行标量扫描相比，现在的系数由矩阵和向量表示，因此需要为进程间集体交互的变量分配更多内存。

3. 计算结果

程序模块在以下任务中进行了测试：长时间停机后，油井以井底恒定压力投入生产。假设油藏边界保持恒定压力，且未观察到相邻油井的影响。对于油井周围的径向油藏（由基质空间和裂缝网络组成），给定了一系列参数：
|参数|数值|
| ---- | ---- |
|$P_f$|$250 \times 10^5$ (Pa)|
|$P_m$|$250 \times 10^5$ (Pa)|
|$\varphi_f$|0.01|
|$\varphi_m$|0.1|
|$\rho_w$|1118 (kg/m³)|
|$\rho_o$|730 (kg/m³)|
|$S_{\alpha}^w$|0.36|
|$\sigma$|0.12 (1/m²)|
|$k_f$|$1 \times 10^{-12}$ (m²)|
|$k_m$|$1 \times 10^{-16}$ (m²)|
|$k_{rw}$|$7.7(S_{\alpha}^w)^4 - 12.07(S_{\alpha}^w)^3 + 6.9(S_{\alpha}^w)^2 - 1.8(S_{\alpha}^w) + 0.2$|
|$k_{ro}$|$(S_{\alpha}^w)^2 \times 0.03 + 0.002(S_{\alpha}^w) + 0.0002$|
|$\mu_o$|$0.67 \times 10^{-3}$ (Pa·s)|
|$\mu_w$|$0.3 \times 10^{-3}$ (Pa·s)|

计算得到了以下结果：
- 压力曲线 ：不同裂缝绝对渗透率（$1 \times 10^{-10}$ m²、$1 \times 10^{-11}$ m²、$1 \times 10^{-12}$ m²）下，压力随距油井距离的变化曲线。结果表明，裂缝渗透率越高，油井周围的压降漏斗越大且出现得越快。
- 温度曲线 ：不同裂缝绝对渗透率下，温度随时间的变化曲线。裂缝渗透率越高，温度下降越快，这与实验观测数据相符。
- 水饱和度曲线 ：在不同时间点（300 s、600 s），裂缝网络和基质中，水饱和度随空间的变化曲线。裂缝中水饱和度的增加比孔隙部分更快，这与流体仅沿裂缝网络流动有关。随着油井生产时间的延长，水饱和度增加，导致油井产量下降。

对于并行算法的分析，使用以下标准公式计算加速比和效率：
- 加速比$S_m = T_1 / T_m$
- 效率$E_m = S_m / m \times 100\%$

当将研究半径划分为 1000 和 5000 个点时，绘制了加速比和效率的图表。结果显示，对于 1000 和 5000 个计算点，m = 14 和 16 分别是最优的并行进程数，此后加速比下降。随着计算点数的增加，需要增加并行进程数。同时，在特定计算中，存在有限最大加速比的效应，即对于特定的计算网格和固定数量的网格节点，存在一个最优的节点数。

在油藏网络的情况下，并行化效率取决于网络图的总分支数和边的最小长度。总分支数决定了所需的并行进程数，边的最小长度影响最终效率。单个油藏越长或越宽，其数值分析所需的计算网格点数就越多，并行效果也就越明显。

4. 结论

通过建立基于双孔隙度模型的非等温数学模型，采用按物理过程分裂的算法构建了带时间权重的差分格式，确保了天然裂缝系统和油藏孔隙部分之间通量的正确性和一致性。使用并行化的矩阵扫描法求解方程组，并通过计算实验验证了该数值算法及其并行实现的有效性。得到了并行算法的加速比和效率随进程数变化的图表，确定了特定问题下的最优进程数（14 和 16）。绘制了压力、温度和饱和度曲线，展示了生产井周围这些参数的变化情况，为调整油井的生产模式和设置提供了依据。

利用并行技术计算流体动力学过程

5. 技术要点总结

为了更清晰地理解整个计算过程，下面对关键技术要点进行总结：
- 数值算法流程
1. 采用按物理过程分裂的方法和有限体积网格法对方程组进行处理。
2. 在合适的空间网格上使用有限体积法进行近似，时间上采用隐式差分近似。
3. 对离散方程进行线性化，得到线性方程组。
- 矩阵扫描法步骤
1. 前向步骤确定扫描系数：$\alpha_{k + 1} = (C_k - A_k \alpha_k)^{-1} B_k$，$\beta_{k + 1} = (C_k - A_k \alpha_k)^{-1} (F_k + A_k \beta_k)$。
2. 后向步骤求解系统：$y_k = \alpha_{k + 1} y_{k + 1} + \beta_{k + 1}$。
- 并行化实现要点
1. 使用 C 语言和 MPI 标准在 K100 超级计算机上实现。
2. 利用 MPI 函数进行模型参数分布、计算域分解、数据交换等操作。
3. 解决短问题时，所有 MPI 进程进行系数的集体交换。

下面用 mermaid 格式的流程图展示整个计算流程：

graph LR
    A[开始] --> B[数值算法]
    B --> B1[按物理过程分裂方程组]
    B --> B2[有限体积法空间近似]
    B --> B3[隐式差分时间近似]
    B --> B4[离散方程线性化]
    B4 --> C[矩阵扫描法]
    C --> C1[前向确定扫描系数]
    C --> C2[后向求解系统]
    C --> D[并行化实现]
    D --> D1[使用 C 语言和 MPI 标准]
    D --> D2[MPI 函数操作]
    D --> D3[解决短问题]
    D --> E[计算结果分析]
    E --> E1[压力曲线分析]
    E --> E2[温度曲线分析]
    E --> E3[水饱和度曲线分析]
    E --> E4[并行算法加速比和效率分析]
    E4 --> F[结论]
    F --> F1[验证算法有效性]
    F --> F2[确定最优进程数]
    F --> F3[调整油井生产模式]
    F --> G[结束]

6. 实际应用与展望

在实际的油藏开发中，这些计算结果和算法具有重要的应用价值。例如，通过分析压力、温度和饱和度曲线，可以更好地了解油藏的动态变化，优化油井的生产策略，提高油井的产量和效率。同时，并行算法的应用可以显著缩短计算时间，提高计算效率，使得对复杂油藏模型的分析成为可能。

未来的研究方向可以包括：
- 模型的进一步优化 ：考虑更多的物理因素，如化学反应、多相流等，以提高模型的准确性和适用性。
- 并行算法的改进 ：探索更高效的并行算法，进一步提高计算效率和加速比。
- 与实际数据的结合 ：将计算结果与实际油藏数据进行对比和验证，不断完善模型和算法。

总之，利用并行技术计算流体动力学过程为油藏开发提供了一种有效的工具和方法，具有广阔的应用前景和研究价值。通过不断的研究和改进，有望在油藏开发领域取得更好的成果。

7. 常见问题解答

为了帮助读者更好地理解和应用上述内容，下面列出一些常见问题及解答：
|问题|解答|
| ---- | ---- |
|为什么要采用按物理过程分裂的方法？|按物理过程分裂的方法可以将复杂的方程组分解为相对简单的部分，便于分别处理和分析，同时可以更清晰地描述系统的物理过程。|
|矩阵扫描法与标量三点方程扫描法有什么区别？|矩阵扫描法的系数由矩阵和向量表示，而标量三点方程扫描法的系数为标量。因此，矩阵扫描法在计算时需要进行更多的矩阵运算，并且在并行化时需要为进程间集体交互的变量分配更多内存。|
|并行算法的加速比和效率受哪些因素影响？|并行算法的加速比和效率受计算点数、并行进程数、网络图的总分支数和边的最小长度等因素影响。随着计算点数的增加，需要增加并行进程数；网络图的总分支数决定了所需的并行进程数，边的最小长度影响最终效率。|
|如何确定最优的并行进程数？|可以通过计算不同进程数下的加速比和效率，绘制相应的图表，找到加速比最大时对应的进程数，即为最优的并行进程数。|

通过以上解答，希望能帮助读者更好地理解和应用相关技术。在实际应用中，如有其他问题，可以进一步深入研究和探索。