二维巷道开挖模型文件

巷道开挖的数值模拟总让我想起小时候玩沙坑的感觉——挖掉一块区域，周围结构就会产生连锁反应。今天咱们用Python整一个二维巷道开挖的模型，看看土体应力怎么重新分布。

先上网格生成的核心代码：

import numpy as np
from scipy.sparse import lil_matrix

x = np.linspace(0, 50, 51)  # 50米巷道，1米网格
y = np.linspace(-20, 20, 41)
X, Y = np.meshgrid(x, y)

这相当于在巷道周围铺了张51x41的坐标网。注意Y轴范围是-20到20，这样上下各留20米缓冲带，避免边界效应影响计算结果。网格密度1米足够捕捉巷道周边的应力集中现象。

边界条件处理是重点：

def apply_boundary_conditions(K):
    # 固定底部边界
    bottom_nodes = np.where(Y == Y.min())
    for node in zip(bottom_nodes[0], bottom_nodes[1]):
        idx = node[0]*len(x) + node[1]
        K[idx*2:idx*2+2, :] = 0  # 固定x,y位移
        K[idx*2:idx*2+2, idx*2:idx*2+2] = 1e10
    return K

刚度矩阵处理这里用了刚度惩罚法，给边界节点赋极大刚度值。实际操作中发现直接消去自由度容易导致矩阵奇异，这种处理方式更稳定。注意节点编号是二维转一维的常见操作，这里用了网格坐标转索引的小技巧。

二维巷道开挖模型文件

开挖模拟的关键步骤：

excavated = (X >= 20) & (X <= 30) & (Y >= -5) & (Y <= 5)
E[excavated] = 1e-3  # 弹性模量骤降
sigma_yy[excavated] = 0  # 竖向应力释放

这里用了个取巧的方法：把开挖区域弹性模量降到接近零，同时释放竖向应力。比直接移除单元更稳定，特别适合小变形分析。实际项目中遇到大变形可能要改用单元生死法。

最后看个应力云图绘制实例：

plt.contourf(X, Y, sigma_yy, levels=20, cmap='jet')
plt.plot([20,30], [5,5], 'w--')  # 标出开挖边界
plt.plot([20,30], [-5,-5], 'w--')
plt.colorbar(label='Vertical Stress (kPa)')

这里用等高线填充图展示竖向应力分布，白色虚线标出开挖轮廓。注意应力单位换算，现场实测数据通常是MPa量级，但模型里用kPa更方便数值计算。

调试时发现个有趣现象：当巷道高宽比接近1:1时，顶板会出现对称的应力拱；但如果是狭长型巷道，应力会向两帮集中。这跟经典弹塑性理论预测的规律一致，说明咱们的简化模型抓住了主要矛盾。

完整代码大概200行左右，主要时间花在刚度矩阵组装和边界处理上。下次可以试试用GPU加速，或者引入节理面模拟岩层结构。数值模拟就像搭积木，每次改进都能发现新的应力传递路径，这大概就是岩土工程的魅力所在吧。