岩体动态响应分析程序的微机化及其应用

D0I:10.13374/j.issn1001-053x.1994.02.002 第16卷第2期 北京科技大学学报 Vol.16 No.2 1994 Journal of University of Science and Technology Beijing AT,1994 岩体动态响应分析程序的微机化及其应用+ 马明军 北京科技大学采矿系，北京100083 摘要开发微机化的数值分析程序是推广数值方法在岩体工程中应用的一条有效途径，本文首 先介绍一种微机数值分析程序的内存数据处理方法，然后介绍据此对岩体动态响应分析的二维非 线性有限元程序DYNPAK的微机化改造和应用实例. 关键词微机应用，数值方法，动态有限元，岩体工程 中图分类号TD325.2,TP311.52 The Improving of DYNPAK Program on Microcomputer and It's application Ma Mingjun Department of Mining and Mineral Engineering,USTB.Beijing,100083,PRC ABSTRACT Developing microcomputer programs have been regarded as a effective way to ex- pand the applications of numerical methods to rock engineering.This paper proposed a new core data-processing method on microcomputer programing of numerical methods so that core memory can be utilized as much as possible.Next,DYNPAK,a 2d FEM program on dy- namic analysis,was improved according to the proposed method on a microcomputer.Finally, a case study was given as a application of the improved DYNPAK. KEY WORDS mirocomputer application,numerical methods,dynamic FEM,rock engineering 应用数值方法解决岩体工程中的复杂问题已经引起了广泛重视，开发微机化的数值分析 程序则是推广数值方法在岩体工程中应用的一条有效途径山.数值方法应用于岩体工程以 来，许多的研究者在其理论、方法、程序和实践方面做出了积极的贡献，但是，这些程序大 都必须在大、中型计算机上运行，限制了其实际应用，应该看到，许多传统的程序都是花费 大量人力、财力而精心设计出来的，并且经过了值得信赖的实践考核，其核心部分和基本功 能是成熟的，因此，在微机上开发数值分析应用程序应该首先立足于传统程序的改造和发 展，以便收到事半功倍的效果·本文的工作正是基于这一指导思想而开展的， 199-02-11收稿第一作者男29岁博士讲师 +山东省黄金公司协作课题

第 卷 第 期 年 月 北 京 科 技 大 学 学 报 以山叮 吨 碑 岩体动态响应分析程序 的微机化及其应用 马 明 军 北京 科技大学采矿 系 ， 北京 以 摘要 开发微机化的数值分析程序是推广数值方法在岩 体工 程 中应 用 的 一 条有 效 途 径 本文 首 先介绍 一种微机数值分析程序 的 内存数据处理方法 ， 然 后介绍 据此对岩体动态响应分析的二维非 线性有 限元程序 卜 的微机化改造和应用 实例 关健词 微机应用 ， 数值方法 ， 动态有 限元 ， 岩体工程 中图分类号 ， ‘ 亩 众 泊 司 七 ， ， 氏石 ， 刃 ， 沈 卫 ℃ 心 浑 获名 俗币 桃 坛卫 在日 叮 涨过 一 姚 卫 ℃ 枉犷 刀 在剧 几旧 篮 饮 ， ， 面 】 ， 从 助川访 掀 止旧山 而 勿 ， 认吸 而 刊 万 姗 ， 吮 记 ， 五 ， 应用 数值方法解 决岩体工程 中的复杂 问题 已 经 引起 了广泛重 视 ， 开发微 机化 的数值分 析 程序则是 推广数值方法 在岩体工 程 中应用 的一条 有 效 途 径 汇’ 数 值 方 法 应 用 于 岩 体 工 程 以 来 ， 许多 的研究 者在其理论 、 方法 、 程序和 实践方 面做 出 了积极 的贡 献 但是 ， 这些 程 序大 都必须在 大 、 中型计算 机上 运行 ， 限制 了其实 际应用 应该看 到 ， 许多传 统 的程序 都是 花费 大量人力 、 财力 而 精 心设计 出来 的 ， 并且经过 了值得信赖 的实践考核 ， 其核 心部分 和 基本功 能是 成熟 的 因此 ， 在微机上 开发 数值分析 应 用 程 序 应 该 首 先 立 足 于 传 统程 序 的 改 造 和 发 展 ， 以便 收到 事半功倍 的效果 本文 的工 作正是基于 这一指 导思想 而 开展 的 一 一 收稿 第一作者 男 岁 博士 讲师 山 东省 黄金公 司协作课题 DOI ：10．13374／j ．issn1001－053x．1994．02．002

第2期 马明军等：岩体动态响应分析程序的微机化及其应用 .107. 1 内存空间的数据处理 与大、中型机相比，微机的主要缺点是内存容量较小，运算速度较低，因此，微机数值 分析应用程序的开发则必须适应微机的这些特点·在微机上实现运行大、中型算题的主要技 术途径可概括为以下几个方面：(1)优化内存空间的数据处理；(2)改变传统的程序结构 (变集中运行为分块运行)；(3)发展新的数值方法或耦合方法；(4)改进求解方法，下文 介绍一种内存空间的数据处理方法· 内存空间的数据处理要求解决以下几个问题：(1)在内存限制的条件下最大可能地利 用内存空间；(2)程序内模块间的数据通讯方便；(3)避免传统程序设计中多而杂的公用 语句，这些数组所占的存储空间应能随后被其他将要用到的数组取代， 1.1实施方法和主体子程序 PROGRAM MAIN COMMON MTOT,NP,IA(100000) MTOT=1000000. 其中，MAN为用户自定的程序名，IA为大型一维数组·在整个程序中（包括主程序 和子程序)，所有根据需要设定的分数组都存放于IA数组中，主要的数据处理操作将在IA 中进行.MTOT为LA的大小，NP为数组类型（整型、实型、字符型），它将用于各分数组 在IA中的定位，与主程序相对应，各子程序中的公用语句形式为： COMMON MTOT,NP,IA(1) 此外，在主程序中另外两个有名公用语句也是必要的（子程序可选择加人），形式为： COMMON/IOFILE/NF1,NF2 COMMON/CONTRL/NJNTS,NMBS 其中，IOFILE是输入输出设备通道号公用区，CONTRL是主控参数公用区. 各分数组在IA中的定位可通过下面的语句来实现，形式为： CALL DEFINI (NAME,NA,NR,NC)CALL DEFINE('NAME',NA,NR,NC) CALL DEFINH (NAMEE,NA,NR,NC) 其中，NAME为用户自定的数组名，NA为该数组在IA中的起始位置，NR为该数组的 行数，NC为该数组的列数.3个子程序（略）分别计算整型、实型和字符型分数组在A 中所占的长度和位置，并同时形成分数组目录，记录于IA的末端.分数组在IA中的定位 和存放形式如图1所示· 数组A 数组B 数组C 数组数组 数组 可以利用的剩余空间 目 存放区 存放区 存放区 目录 目录 图1 内存空间的数组动态定位 Fig.1 Dymamic lcation and dimension for arrays in core

第 期 马 明军等 岩体动态 响应分析程序 的微机化及其应 用 内存空 间的数据处理 与大 、 中型机相 比 ， 微机 的主要 缺 点是 内存容量较小 ， 运算 速度较低 ， 因此 ， 微机数值 分析应用程序 的开发则必须适应微机的这些 特点 在微机上 实现运行大 、 中型 算题 的主要技 术途径可概括 为 以下 几个方 面 优化 内存空 间的数据处理 改 变 传 统 的程 序 结构 变集 中运行为分块运行 发展新 的数值方法或藕合方 法 改进求解方法 下 文 介绍一种 内存空 间的数据处理方法 内存空 间的数据处理要求解 决 以下 几个 问题 在 内存 限 制 的 条件 下 最 大 可 能 地 利 用 内存空 间 程序 内模块 间的数据通讯方便 避 免传 统程 序 设 计 中多 而 杂 的公 用 语句 ， 这些数组所 占的存储空 间应能 随后被其他将要 用到 的数组取代 实施方法和主体子程序 洲 … ， 洲 其 中 ， 为用 户 自定 的程序名 ， 认 为大 型 一 维数 组 在 整 个程 序 中 包括 主 程 序 和子程序 ， 所有根据需要 设定 的分数组 都存放于 数组 中 ， 主要 的数据处理 操 作 将 在 中进行 为 的大小 ， 为数组类型 整 型 、 实型 、 字符型 ， 它将用 于各分数组 在 中的定位 与主程序相 对应 ， 各子程 序 中的公用语句 形式 为 ， ， 认 此外 ， 在 主程序 中另外 两个有 名公 用 语句也是 必要 的 子程 序 可 选 择加人 ， 形 式 为 ， … 丁 ， … 其 中 ， 是 输入输 出设备通道 号公用 区 ， 丁 是 主控参数公用 区 各分数组 在 中的定位 可 通过下 面 的语句来 实现 ， 形 式 为 卜 ’ ’， ， ， 刊 ’ ’ ， ， ， ， ， ， 其 中 ， 为用 户 自定 的 数组 名 ， 为 该 数组 在 中 的 起 始 位 置 ， 为 该 数 组 的 行数 ， 为该数组 的列 数 个 子 程 序 略 分 别 计算 整 型 、 实 型 和 字 符 型 分 数 组 在 中所 占的长度和位置 ， 并 同时形成分数 组 目录 ， 记 录 于 的末 端 分 数组 在 中 的定 位 和存放形式 如 图 所示 数组 数组 数组 可 以 利 用 的剩余空 间 存放 区 存放 区 存放 区 数组 目录 数组 目录 数组 目录 图 内存空 间的数组动态定位 ’ 功阿 位犯 日 成恤 叨 俪 日 ， 班℃

108 北京科技大学学报 第16卷 对于某些只起过渡作用的数组可用下面的语句来释放，形式为：CALL DELETE(NAME). 该子程序（略）释放名为NAME的数组及其目录在IA中所占的空间，以便被其他数组所 取代，从而有效利用内存空间， 1.2示例 现在，用下面的例子来简要说明内存空间数据处理的过程， PROGRAM DATAGN COMMON MTOT,NP,IA(100000) MTOT=100000 +。， CALL CORDFN (KXCOR,KYCOR,KZCOR,NJNTS) CALL XYZDAT(KXCOR,KYCOR,KZCOR) STOP END SUBROUTINE CORDFN(KXCOR,KYCOR,KZCOR,NJNTS) COMMON MTOT,NP IA(1) CALL DEFINE (XCOR',KXCOR,1,NJNTS CALL DEFINE (YCOR',KYCOR,1,NJNTS) CALL DEFINE ('ZCOR',KZCOR,1,NINTS RETURN END SUBROUTINE XYZDAT (KXCOR KYCOR,KZCOR) COMMON MTOT,NP,IA(1) … RETURN END 在该例中，NNTS表示节点数，由主程序部分赋值.调用CORDFN后返回KXCOR, KYCOR、KZCOR设定了3个实型数组在IA中的位置，并已经在IA中设置了3个数组的长 度和相应的目录信息，调用XYZDAT后，节点的x,y,z坐标依次存放在IA数组中，供后 续的计算调用或处理， 传统的处理方法通常是在各程序单位中硬性规定数组的大小并由有名公用语句传递·例 如：COMMON/COORD/X(1000),Y(1000),Z(1000). 由于大型应用程序的数组很多，使得COMMON语句很多，既不简明，也不能充分利 用有限的内存空间，而且修改一个数组，就会牵一动百，株连其他.使用DIMENSIOND定

北 京 科 技 大 学 学 报 第 卷 对于某 些只起过 渡作 用 的数组 可 用 下 面 的语句来释放 ， 形式为二 叹 ， 该子程序 略 释放名 为 的数组及 其 目录 在 认 中所 占的空 间 ， 以 便 被 其他 数 组 所 取代 ， 从而 有效利 用 内存 空 间 示例 现在 ， 用 下 面 的例子来 简要说 明 内存 空 间数据处理 的过程 五 ， ， ， ， ， ， ， 」 刀 旧 ， ， ， 丑们 ， ’ ， ， ， 丁 ‘ ， ， ， ， ‘ ‘ ， ， ， 丁 旧 欠 ， ， ， 在 该 例 中 ， 表 示 节 点 数 ， 由主 程 序 部 分 赋 值 调 用 后 返 回 、 、 设定 了 个 实型 数组 在 中的位置 ， 并 已 经在 中设置 了 个数组 的长 度 和相 应 的 目录信 息 ， 调 用 后 ， 节点 的 ， ， 坐标 依 次存 放 在 数组 中 ， 供后 续 的计算调 用 或处理 传 统 的处理 方法 通 常是在各程 序 单位 中硬性规定数组 的大 小并 由有 名公用语句 传递 例 如 ， ， 刃 … 由于 大 型应 用 程序 的数组很多 ， 使 得 语 句 很 多 ， 既 不 简 明 ， 也 不 能 充 分利 用有 限 的 内存 空 间 ， 而且修 改一 个数组 ， 就 会牵一 动百 ， 株 连其 他 使 用 定

第2期 马明军等：岩体动态响应分析程序的微机化及其应用 .109. 数组也会遇到同样的问题，而且大量数组的虚实结合还会降低运行速度， 2 DYNPAK程序的微机化改造 二维非线性动态有限元程序DYNPAK是用FORTRAN语言写成的，原始程序由Owen 和Hinton于1980年给出).该程序可模拟材料的物理非线性和几何非线性，以及对地震载 荷或脉冲载荷的响应·微机化改造的主要内容是引进前面讨论的内存空间数据处理方 法，以便：(1)该程序能在普通微机上运行；(2)扩大程序的解题能力；(3)增强程序的 可读性，便于修改和扩充，此外，还对该程序的输入系统做了改进，将传统的固定格式输 入改为可标识的自由格式输人，大大地方便了输入数据的准备，下面是改造后的部分主程序， PROGRAM DYNPAK IMPLICIT REAL (A-H,O-Z) COMMON MTOT,NP,IA (100 000) COMMON/CONTRL/NPOIN,NELEM,NDOFN,NMAT,NVFIX,... COMMON/IOFILE/IPR,NF1,NF2,NF3,. … CALL ARRDFN (KCORD,KIFPR,KMATN,. CALL INPUTD (KCORD,KIFPR,KMATN,.. …… CALL LUMASS (KCORD,KMATN,.. … STOP END SUBROUTINE ARRDFN(KCORD,KIFPR,KMATN,. IMPLICIT REAL (A-H,O-Z) COMMON MTOT,NP,IA (1) COMMON/CONTRL/NPOIN,NELEM,NDOFN,NMAT,NVFIX,.. CALL DEFINE (CORD',KCORD,NDOFN,NPOIN) CALL DEFINE (IFPR,KIFPR,NDOFN,NPOIN CALL DEFINE (MATN',KMATN,1,NELEM ”+ RETURN END 改造后的DYNPAK程序可在普通IBM-PC/XT/AT/286/386微机上运行，由于保留 了原始程序的核心部分（语法按MS-FORTRAN77修改）和基本功能，示例计算的结果与 Owen和Hinton给出的结果吻合很好，限于篇幅，在此不一一给出

第 期 马 明军等 岩体动态响应分析程序 的微机化及其应用 数组也会遇到 同样 的 问题 ， 而 且大量 数组 的虚 实结合还 会 降低 运行 速 度 程序 的微机化 改造 二 维非 线性 动态有 限元程序 是 用 人 语言 写 成 的 ， 原 始 程 序 由 和 于 年给 出 「】 该程 序可模拟 材料 的物理 非 线性 和几何 非 线性 ， 以 及 对 地 震 载 荷 或脉 冲 载 荷 的 响应 微 机 化 改 造 的 主 要 内容 是 引 进 前 面 讨 论 的 内 存 空 间 数 据 处 理 方 法 ， 以便 该程序能在普通微 机上 运行 扩大程序 的解题 能 力 增 强 程 序 的 可读性 ， 便于修改和 扩充 此外 ， 还 对该 程 序 的输 入 系 统 做 了 改 进 ， 将 传 统 的 固 定 格 式输 人改 为可 标识 的 自由格式输人 ， 大大地方便 了输人数据 的准备 下面是 改造后 的部分主程序 仃 一 ， 一 ， ， 丁 刊田 呵 ， ， ， ， 』 ， … 刀〕 正 ， ， 于 ， ， ， … ， ， ， … ， ， ， … ， ， … 叮 ， ， ， … 吐 】 ， ， 一 均 ， ， 陌 ， ， ， ， ， … ’ ， ， ， ， ， ， ， 卜丁 卜 ’ ‘ ， ， ， 改造后 的 程 序 可 在 普 通 一 叉丁 戌 ， 微 机 上 运 行 由于 保 留 了原始程序 的核心部分 语法按 一 修 改 和基本功 能 ， 示 例计 算 的结 果 与 叭七 和 给 出的结果 吻合很 好 ， 限于 篇 幅 ， 在此不 一 一 给 出

·110 北京科技大学学报 第16卷 3应用实例 下面介绍运用微机化的DYNPAK程序分析某选矿厂地表对露天矿爆破地震的动态响 应，计算模型取表土和一部分基岩，高100m,宽300m,分为315个等参单元.表土和基 岩都设为弹塑性材料，采用Mohr-Coulomb屈服准则.岩、土体动力学参数如表1所示. 有限元网格如图2所示.图中的A、B两点之间为选矿厂所处的部位， 在动态分析之前，先进行静态有限元分析，由静态分析的结果作为动态分析的初应力输 入·动载荷主要考虑来自水平方向的爆破激振，激振输入的加速度一时间曲线由现场的爆 破地震测试提供， 表1岩、土体动力学参数 Table 1 Dynamic parameters of rock and soil 岩性弹模泊松粘结内靡容重 比 力擦角 (MPa) (MPa)(）kg/m) 表士40 0.3 215 1.8 基岩1650000.2524.548 2.85 图2边坡断面的有限元网格 Fig.2 FEM mesh for slope section of some opencast mine 该分析所考察目标是选矿厂部位的地表振速、位移及变形情况（下面的结果中未计入地 表的静态位移与变形).图3示出了A、B两点水平位移的时间历程；图4示出了厂区地表倾 斜和水平变形的时间历程.A,B两点的最大振速分别为3.24cm/s和7.56cm/s.按照国家有 关标准，地表建筑物的安全性属需要警惕的范围， 20 10 0 0.002 0.40.60.81.0 t/s 图3A、B两点的水平位移历程 Fig.3 Displacement histories for points A.B

北 京 科 技 大 学 学 报 第 卷 应用 实例 下 面介 绍运 用 微 机 化 的 卜 程 序 分 析 某 选 矿 厂 地 表 对露 天 矿 爆破 地 震 的 动 态 响 应 计算模型取表 土和 一部分基岩 ， 高 ， 宽 ， 分 为 个 等 参单元 表 土 和 基 岩都设为 弹 塑性材料 ， 采 用 一 屈 服 准 则 岩 、 土 体动 力 学 参数 如 表 所 示 有 限元 网格 如 图 所示 图 中的 、 两点之 间为选矿厂 所处的部位 在动态分析之前 ， 先进行静态有 限元分析 ， 由静态分析 的结果作 为 动态分 析 的初 应力 输 人 动载荷 主要 考虑来 自水平 方 向的爆破激振 ， 激振输入 的加 速 度 一 时 间 曲线 由现 场 的爆 破地震测 试提供 表 岩 、 土体动 力学参数 址山 】妙 让 脚皿搜七招 健 比 川 ，司 岩性 弹 模 泊 松 比 粘 结 力 巧 内 摩 擦角 容 重 洲 傀 耐 巨于井 二 三王三王三， ，’ ， 月 月 布】〕二 厂 一 丁 二」二二工二 田 － 一 一 万口月『以月〕 二 厂 一 厂二一 二川 刀 刀 盯 〔 了叮「二二二 二 厂 口门口 川口 厂门 曰口 川口 巨二二 口口口 川口 二刃 口习三 川 二三二 川 二三」 川 表 土 基岩 《刀 抖 图 边坡断面的有限元网格 万触 日她 由 俪 幽脾 “ ℃山扣 成 ，，世 卿曰眨以 口血此 该分析所考察 目标是 选矿厂部位 的地表振速 、 位移及 变形情况 下 面 的结果 中未计人 地 表 的静态位 移 与变形 图 示 出 了 、 两点水平位移 的 时 间历程 图 示 出了厂 区地表倾 斜和水平 变形 的时 间历程 、 两点 的最大振速分别为 和 按照 国家有 关标准 ， 地表建筑物 的安全性 属需要警惕 的范 围 一 馋攀︸口目 一 ︸ 几户 才矛 图 、 两点 的水平位移历程 瑰 晚自“ 盯犯成 肠演面七 阮 ，血妇

第2期 马明军等：岩体动态响应分析程序的微机化及其应用 .111· 0.100ra b 0.050 -0.000 -0.050 -0.100 图4厂区地表的倾斜和水平变形历程 Fig.4 Surface inclination and deformation histories for workshop 4结论 从上述讨论可知：(1)本文介绍的内存空间数据处理方法是先进的，它不仅可以充分 利用内存空间，而且语句简明易读，易于修改和扩充，尤其适用于在内存空间有限的微机上开 发数值分析程序；(2)用该法改造传统的数值分析程序，只需对有关语句作少量修改，而 且可以扩大原程序的解题能力；(3)DYNPAK程序的微机化改造是成功的，实例计算的成 果已提供给有关生产单位参考；(4)DYNPAK程序只给出了自由边界及固定边界的处理， 这将给计算带来一定的误差，应该在该程序中进一步引人人工边界处理方法， 参考文献 1刘怀恒.岩土力学，1989,10(2)：15~18 2蔡跃军，微机有限元弹塑性分析程序及其在矿业工程中的应用.[博士论文]，长沙：中南工业大学，1990 3 Owen DR J,Hinton E.Finite Element in Plasticity.In:Swansea UK eds.Pineridge Press,1980.377~429 4廖振鹏.地震工程和工程振动，1982,2(1)：11~14

第 期 马 明军等 岩体动 态 响应分析程序 的微机化及其应用 一日 · 习 欲俐降关、昌 一 旧 一 一畜草︸日昌 一 。 、 图 厂 区地表的倾斜和水平变形历程 瑰 反『白 沁由 自犯 回 山而门圈位扣 地血的巴 助 佩刘 脚 结 论 从上 述讨论可 知 本文介绍 的 内存空 间数据 处理 方 法 是 先 进 的 ， 它 不 仅 可 以 充 分 利用 内存空 间 ， 而且语句简明易读 ， 易于修改和 扩充 ， 尤其适用于 在 内存空 间有 限 的微机上开 发数值分析程序 用该法 改造传 统 的数值分析 程 序 ， 只 需 对有 关 语句 作 少 量 修改 ， 而 且可 以扩大原程序 的解题 能力 程序 的微机化改造是成功 的 ， 实例 计 算 的成 果 已提供给有 关生产单位参考 卜 程序 只 给 出 了 自由 边 界 及 固定 边 界 的处 理 ， 这将给计算带来一定 的误 差 ， 应 该 在 该 程 序 中进 一 步 引 人 人 工 边 界 处理 方法 参 考 文 献 刘 怀恒 岩 土力学 ， ， 巧 一 蔡跃军 微机有 限元弹塑性分析程序及 其在 矿业工程 中的应用 博士论文 长沙 中南工 业大学 ， 男 阴 ， 而 此 眺石 钾 玩 泥 ‘ 已 辱 亡洛 ， 一 廖振鹏 地震工程和 工 程振动 ， ， 一