二维钻孔封孔效果模拟案例
钻孔封孔效果模拟这事挺有意思的。咱今天拿个简化版的二维模型练手，用Python搞个渗流场可视化。先说场景：地下50米有个直径0.2米的钻孔，现在要往里面注水泥浆，得看看封堵效果咋样。
先整网格。用numpy搞个非均匀网格，靠近钻孔的位置网格加密：
import numpy as np radius = 0.2 # 钻孔半径 domain_size = 5.0 # 模拟区域边长 r = np.concatenate([np.linspace(radius, 0.5, 30), np.linspace(0.5, domain_size, 70)[1:]]) theta = np.linspace(0, 2*np.pi, 120) R, Theta = np.meshgrid(r, theta)这里用极坐标系更贴合钻孔形状。靠近钻孔的0.2米到0.5米区间安排了30层网格，外围则用较稀疏的70层，总网格数控制在合理范围。注意[1:]的用法是为了避免端点重复。
接着上物理模型，考虑浆液扩散的非牛顿流体特性：
def calc_viscosity(shear_rate): # 赫歇尔-巴尔克利模型 tau_0 = 120 # 屈服应力 k = 0.5 # 稠度系数 n = 0.3 # 流动指数 return tau_0/(shear_rate + 1e-6) + k*(shear_rate)**(n-1)这个粘度模型是关键，避免了牛顿流体假设带来的误差。1e-6的小量防止除零，实际工程中常见处理方式。参数取值根据实际水泥浆配比调整，这里用典型值示意。
边界条件设置要讲究，注浆压力随时间变化：
injection_pressure = np.interp(time, [0, 10, 30], [5e6, 8e6, 6e6]) # 分段线性加压模拟真实注浆过程——初始阶段逐步升压，中期维持高压，后期缓慢降压。用np.interp做线性插值比写循环更高效，适合向量化计算。
求解器部分用显式差分，注意稳定性条件：
dt = 0.1 * (dx**2) * rho_max / (viscosity_max) # CFL条件 for _ in range(steps): new_p = p.copy() # 极坐标下的扩散方程离散 new_p[1:-1,1:-1] = p[1:-1,1:-1] + dt*( ... ) p = new_p这里时间步长根据最大粘度和网格尺寸动态调整，避免计算发散。扩散项的具体离散需要考虑极坐标的1/r项，代码中简写为(...)，实际需要展开拉普拉斯算子。
最后来个动态可视化：
plt.figure(figsize=(10,6)) contour = plt.contourf(X, Y, pressure_field, levels=20, cmap='jet') plt.plot(wellbore_x, wellbore_y, 'w-', linewidth=2) # 绘制钻孔 plt.colorbar(label='Pressure (Pa)') plt.title(f'Grout propagation t={time:.1f}s')用极坐标转笛卡尔坐标后的X,Y做显示更直观。颜色映射选jet突出压力梯度，白线标出钻孔位置。动态更新时只需要更新contour的数据，比反复创建新图快得多。
跑完模拟发现个有趣现象：在注浆中期，浆液在钻孔北侧形成指进现象。检查代码发现是网格各向异性导致的，后来改用自适应网格后改善。做工程仿真时，这种数值假象要和物理真实现象区分开，得反复验证。
