二维巷道开挖模型文件
二维巷道开挖模型文件
巷道开挖的数值模拟总让我想起小时候玩沙坑的感觉——挖掉一块区域，周围结构就会产生连锁反应。今天咱们用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加速，或者引入节理面模拟岩层结构。数值模拟就像搭积木，每次改进都能发现新的应力传递路径，这大概就是岩土工程的魅力所在吧。