外，为满足矿岩自然崩落模拟的要求，自行开发增设了崩落判据程序。整个程序源代码近1万行，总框图如图5—4所示。
图5-5网格自动划分程序框图 图5-4程序流程总框图
整个粘弹塑性有限元程序求解过程主要分为如下七个不同的阶段：
(1) 有限元网格的自动形成
有限元网格的自动划分程序主要针对岩土工程问题编写的。BOPACE 3D源程序中的有限元前处理极不完善，是将所有节点和单元按一定的顺序进行一个一个的数据输入。考虑到工程实际计算模拟模型一般都较大，节点和单元少则几千，多则几万甚至几十万的数据输入量。如此繁琐而艰巨的数据录入工作，给前期建模带来了极大的难度，且随着数据输入量的增大极易出错。有时尽管进行了反复的校核，但仍可能留下人为的错误而无法觉察到，这种错误的存在势必导致不正确的计算结果。如果在计算过程中发现计算结果不尽合理而要改变网格的划分，这时全部数据必须重新填写和输入，这样将导致工作量成倍的增加。
鉴于BOPACE 3D前处理的数据录入极不方便和难以确保数据的正确性，本文针对三维有限元划分网格的特点，自行开发了基于空间六面体8结点等参单元三维有限元网格自动划分的岩土工程通用程序，解决了BOPACE程序前处理难题，大大地提高了前期数据录入工作效率和数据的可靠度。网格自动划分程序的框图如图5-5所示。
(2) 有限元的图形显示
有限元网格及边界条件的图形显示很重要，它可以直观显示有限元网格自动形成的结果，直观了解网格形成的正确性及合理性，指导使用人员修改有限元网格自动形成程序的输入参数。
本文使用FORTRAN90加油站的图形功能及OpenGL三维图形函数库来实现图形显示，在FORTRAN90加油站编译系统上编译, 有限元网格的图形显示程序框图如图5-6所示。
尽管AutoCAD具有很强的绘图功能，能实现图形的显示和绘制，但得通过数据文件转换图形数据，然后使用Lisp语言或C语言绘图，如此处理存在许多问题。在FORTRAN语言内部通过
30
OpenGL三维图形函数库绘图具有很大的方便性，不必传递图形数据，动态修改图形也很方便。其次工程设计人员比较熟悉FORTRAN语言，使用简单、方便。
OpenGL是硬件和图形的软件接口，实际上就是一个三维图形和模型库。由于它在三维真实感图形制作中具有优秀的性能，所以，诸如Microsoft、SGI、LBM、DEC、SUN、HP等在计算机市场占主导地位的大公司，都将它作为自己的图形标准，从而使之成为新一代的三维图形工业标准。
OpenGL不仅仅是一个图形库，更是一个API(Application Programming Interface)。它本身是一个与硬件无关的编程接口，可以在不同的硬件平台上得到实现，也正是因为如此，OpenGL中没有包含处理窗口和用户输入的命令。在OpenGL中不提供三维造型的高级命令，虽然0penGL也是通过基本的几何图元——点、线、多边形来建立物体模型，但更确切地说它应该被称为新一代的三维图形开发标准。现在有很多优秀的三维图形软件，如3Dmax等可以方便地建立模型，但难以对其进行控制。把这些模型转化为用Fortran语言编写的OpenGL程序，就可以随心所欲地控制这些模型，制作CAD，制作三维动画，实现虚拟仿真，这些都使我们制作三维真实感图形更方便、更真实。
OpenGL可以制作各种各样的三维图形，方便地实现三维图形的交互操作。OpenGL提供的图形操作有：绘制物体、、变换、着色、光照、反走样、混合、雾、位图和图像、纹理映射和交互操作和动画。
C语言及FORTRAN语言可以很方便地调用OpenGL图形函数，实现图形显示，特别是0penGL提供了许多坐标变换函数，绘制三维图形比较方便。因此，用FORTRAN语言绘制三维图形，选用OpenGL图形库是惟一较好的途径。而仅仅使用FORTRAN语言内部图形函数绘制平面图形是可以的，但绘制三维图形就非常困难。
图5-6网格自动划分程序框图
(3) 形成节点及单元数据文件
网格自动生成后，对网格节点进行编号，形成有限元分析所需要的节点及单元数据文件。 (4) 常数结构矩阵的装配
在进行KU?R求解之前，计算单元刚度矩阵，并组装成总刚度矩阵K，并贮存在硬盘上。此外，计算并贮存有效线性结构刚度矩阵。
(5) 载荷矢量计算
31
对每个时间步，计算外部施加的荷载矩阵R并贮存到硬盘上。
(6) 逐步求解
采用迭代的方法计算位移、应变和应力程序，该程序是在BOPACE3D程序的基础上开发而成的，是有限元程序的重要组成部分，其功能强大、适用范围广。经过对源程序和理论文本的分
图5-7 迭代求解位移、应变和应力程序框图 图5-8 崩落判据程序框图
析和理解，得出位移、应力与应变有限元程序框图如图5-7所示。
(7) 崩落判据程序框图
崩落判据程序是根据对三矿区岩体特性进行研究后，经理论推导，得出相应公式，并在此基础上将其转化为程序的，崩落判据程序框图如图5-8所示。
5.3 Flac3D原理简介
5.3.1 空间导数的有限差分近似
三维快速拉格朗日法采用了混合离散方法，区域被划分为常应变六面体单元的集合体，而在计算过程中，程序内部又将每个六面体分为以六面体角点为角点的常应变四面体的集合体，变量均在四面体上进行计算，六面体单元的应力、应变取值为其内四面体的体积加权平均。如图5-9所示一四面体，节点编号为1到4，
0第ｎ面表示与节点ｎ相对的面，设其内任一点的速率分量为ui，则可由高斯公式得：
图5-9 四面体
32
00?Vui,jdV??SuinjdS （5-4）
其中V为四面体的体积，S为四面体的外表面，ni为外表面的单位法向向量分量。对于常
0应变单元，ui为线性分布，ni在每个面上为常量，由式(5-4)可得：
0ui,j??13V4?l?10uinjl(l)S(l) （5-5）
其中，上标l表示节点l的变量，(l)表示面l的变量。
5.3.2 运动方程
三维快速拉格朗日法以节点为计算对象，将力和质量均集中在节点上，然后通过运动方程在时域内进行求解。节点运动方程可表示为如下形式：
0?ui?tll?Fi(t)mll (5-6)
其中Fi(t)为在ｔ时刻l节点的在ｉ方向的不平衡力分量，可由虚功原理导出。ml为l节点的集中质量，在分析动态问题时采用实际的集中质量，而在分析静态问题时则采用虚拟质量以保证数值稳定，对于每个四面体，其节点的虚拟质量为：
m?la19Vmax{?n(i)iS(l)?,i?1,3} （5-7）
2其中a1?K?43G，Ｋ为体积模量，Ｇ为剪切模量。式(5-7)的前提是计算时步?t=1。
任一节点的虚拟质量为包含该节点的所有四面体对该节点的贡献之和。将上式左端用中心差分来近似，则可得到：
0ui(t?l?t20)?ui(t?l?t2)?Fi(t)mll??t （5-8）
5.3.3 应变、应力及节点不平衡力
三维快速拉格朗日法由速度来求某一时步的单元应变增量，如下式：
??ij?1200j，i(ui，j?u)??t （5-9）
有了应变增量，即可由本构方程求出应力增量：
??ij?Hij(?ij，??ij)???ij （5-10）
其中Ｈ为已知的本构方程，??ij为大变形情况下对应力所作的旋转修正：
??ij?(?ik?kj??ik?kj)??t （5-11）
其中：?ij?1200j，i(ui，j?u)
33
