二维巷道开挖模型文件巷道开挖的数值模拟总让我想起小时候玩沙坑的感觉——挖掉一块区域周围结构就会产生连锁反应。今天咱们用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*22, :] 0 # 固定x,y位移 K[idx*2:idx*22, idx*2:idx*22] 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, levels20, cmapjet) plt.plot([20,30], [5,5], w--) # 标出开挖边界 plt.plot([20,30], [-5,-5], w--) plt.colorbar(labelVertical Stress (kPa))这里用等高线填充图展示竖向应力分布白色虚线标出开挖轮廓。注意应力单位换算现场实测数据通常是MPa量级但模型里用kPa更方便数值计算。调试时发现个有趣现象当巷道高宽比接近1:1时顶板会出现对称的应力拱但如果是狭长型巷道应力会向两帮集中。这跟经典弹塑性理论预测的规律一致说明咱们的简化模型抓住了主要矛盾。完整代码大概200行左右主要时间花在刚度矩阵组装和边界处理上。下次可以试试用GPU加速或者引入节理面模拟岩层结构。数值模拟就像搭积木每次改进都能发现新的应力传递路径这大概就是岩土工程的魅力所在吧。
二维巷道开挖模型文件
发布时间：2026/6/27 22:32:55
二维巷道开挖模型文件巷道开挖的数值模拟总让我想起小时候玩沙坑的感觉——挖掉一块区域周围结构就会产生连锁反应。今天咱们用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*22, :] 0 # 固定x,y位移 K[idx*2:idx*22, idx*2:idx*22] 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, levels20, cmapjet) plt.plot([20,30], [5,5], w--) # 标出开挖边界 plt.plot([20,30], [-5,-5], w--) plt.colorbar(labelVertical Stress (kPa))这里用等高线填充图展示竖向应力分布白色虚线标出开挖轮廓。注意应力单位换算现场实测数据通常是MPa量级但模型里用kPa更方便数值计算。调试时发现个有趣现象当巷道高宽比接近1:1时顶板会出现对称的应力拱但如果是狭长型巷道应力会向两帮集中。这跟经典弹塑性理论预测的规律一致说明咱们的简化模型抓住了主要矛盾。完整代码大概200行左右主要时间花在刚度矩阵组装和边界处理上。下次可以试试用GPU加速或者引入节理面模拟岩层结构。数值模拟就像搭积木每次改进都能发现新的应力传递路径这大概就是岩土工程的魅力所在吧。
