多孔介质流动 多孔介质中的两相流动 多尺度模拟孔隙尺度建模Darcy-Brinkman-biot理论。 可以模拟粘性指进毛细管力驱动接触角研究。 模拟驱替和渗吸过程。多孔介质里的流动总带着点魔法色彩——你看不见的孔隙中两种流体像跳探戈似的相互推挤。实验室里做这类实验烧钱又费时数值模拟就成了勘探地下油藏、设计燃料电池的必备法器。今天咱们唠唠怎么用代码再现这些微观尺度上的博弈。先看个两相流驱替的典型场景。当注入流体把地层里的原油往外赶毛细管力会在孔隙咽喉处搞事情。用相场法模拟这个现象时得同时算压力场、速度场和相界面演化import numpy as np from scipy.sparse import linalg def solve_phase_field(): gamma 0.1 # 界面张力系数 epsilon 0.01 # 界面厚度 dt 1e-4 # 时间步长 # 构造压力泊松方程矩阵 A build_poisson_matrix(nx, ny) # 自定义的系数矩阵 p linalg.spsolve(A, div_u) # 解压力场 # 相场演化方程 c update_phase_field(c_old, u, gamma, epsilon, dt) return p, c这段伪代码里的压力求解用稀疏矩阵加速相场更新涉及对流-扩散方程。注意epsilon这个参数控制界面模糊程度——太小了数值震荡太大了物理失真一般在3-5个网格宽度较合适。说到多尺度耦合Darcy-Brinkman方程把宏观渗流和微观剪切粘性统一起来。下面这段FEniCS代码实现了该模型from fenics import * mesh UnitCubeMesh(24, 24, 24) V VectorFunctionSpace(mesh, P, 2) Q FunctionSpace(mesh, P, 1) u TrialFunction(V) p TrialFunction(Q) v TestFunction(V) q TestFunction(Q) # 达西项与纳维斯托克斯项耦合 a inner(grad(u), grad(v))*dx (mu/k)*inner(u, v)*dx L inner(p_div, v)*dx inner(f, v)*dx # 分裂式求解 solve(a L, u, bcs)这里把渗透率k作为场变量传入当k趋近无穷大时就退化为标准Navier-Stokes方程。注意动量方程里的mu/k项——它像刹车片一样孔隙越小k降低流体越难运动。多孔介质流动 多孔介质中的两相流动 多尺度模拟孔隙尺度建模Darcy-Brinkman-biot理论。 可以模拟粘性指进毛细管力驱动接触角研究。 模拟驱替和渗吸过程。粘性指进现象最能体现数值方法的健壮性。下图展示用格子玻尔兹曼方法模拟的非稳定驱替过程代码中边界条件处理很关键void LBM::streaming() { for (int i1; iNX-1; i) { // 避开边界 for (int j1; jNY-1; j) { for (int k0; kQ; k) { int ip i - cx[k]; int jp j - cy[k]; f_new[i][j][k] f[ip][jp][k]; // 迁移步骤 } } } apply_bounce_back(); // 固体边界反弹 }这里采用半步长反弹格式处理复杂孔隙边界。迁移时避开计算域边缘单独用applybounceback处理固体节点避免边界穿透问题。这种显式处理虽然计算量大但对GPU加速友好。接触角的影响更是个精细活。在相场法中需要将杨氏方程离散到固壁边界% 接触角边界条件处理 theta 60 * pi/180; # 目标接触角 grad_c compute_gradient(c); n_wall [0,1]; % 壁面法向量 grad_c_dot_n grad_c * n_wall; s grad_c_dot_n cos(theta)*norm(grad_c); c enforce_boundary(c, s); # 通过迭代修正c的值这里用表面法向梯度与接触角的余弦值构成约束条件。实际操作中常需要亚迭代才能收敛就像在刀尖上跳舞——参数敏感性极强。这些代码片段像乐高积木搭建成模拟多孔介质流动的虚拟实验室。当你在参数文件中把黏度比调到1000:1亲眼看到粘性指进如树枝分叉般生长时会真切感受到数值模拟的魔力——那些藏在岩石孔隙间的流体奥秘正被一行行代码逐步破译。