深度解析三维细观力学数值分析的程序设计

在第8章中已经给出了混凝土试件损伤破坏数值模拟的基本方程，这些方程包括应变率强化关系、损伤本构关系及其演化法则以及混凝土试件损伤破坏静、动力学方程，并介绍了数值计算方法和相应的计算程序。与二维混凝土细观力学分析程序的步骤类似，在进行三维混凝土细观力学计算时，首先按照混凝土实际级配的骨料含量和比例生成球形（多面体）骨料模型（马怀发，2005；李运成等，2006）。将所生成的随机骨料模型投影到混凝土试件的计算区域，按照第7章给出的混凝土细观有限元剖分方法对整个试件进行有限元剖分，并判断单元属性，生成单元编码和结点编码数据块。三维问题求解自由度巨大，目前普通PC机的内存量远不能满足如此大规模方程的求解要求，为了使程序能够在普通PC机上运行，在三维程序设计中采用了整体刚度矩阵分块求解算法，求解过程的中间变量采用外存储读写方式，其结构特点和各个程序块间的相互依赖关系如图12.1所示。

图12.1　程序中各予程序的相互依赖关系

由随机骨料模型及剖分程序生成的结点编码、单元编码数据和其他材料参数、荷载信息、约束情况以及程序控制参数的输入通过input子程序实现。在数据输入完成后，利用stochparam子程序，分别对混凝土细观各相单元随机地生成其中每个单元的初始弹性模量与静拉强度，使得这两个参数都遵循对数正态分布，具体实现算法参见第11章。(www.xing528.com)

在进行加载过程数值模拟之前，先利用sdtk子程序为每个单元形成各自的初始单元刚度矩阵，对四面体单元调用stiff4来形成，对六面体单元调用stiff8来形成。由于在该数值计算模型中，有限元网格始终不变，所以整体刚度矩阵的稀疏结构也始终不变，而原程序中采用变带宽算法求解对应于整体刚度矩阵的稀疏线性方程组，从而在主程序中调用hband模块来计算整体刚度矩阵每行的半带宽以及全局半带宽。同时，由于在利用变带宽算法求解稀疏线性方程组时，存储需求相当大，所以当计算规模比较大时，很容易出现内存不足的问题。为解决这个问题，在整体刚度矩阵按变带宽形式形成与进行Cholesky分解时，采用分块技术来实现，将整个矩阵分成多个块，一次形成该块中的元素并对与之对应的Cholesky分解因子进行计算，每算完一个块，就将其存储到磁盘上。在Cholesky分解完成之后，利用三角线性方程组来求解整个稀疏线性方程组时，再逐次调用相应的数据文件来进行计算。对整个刚度矩阵形成分块信息通过调用divbl来实现。

计算的主体由solution模块来实现，在该模块中，首先调用load子程序来输入承受集中荷载的结点号与这种结点上的集中荷载值及加载时比例系数对应的向量，调用gravity施加重力荷载，调用subparam来输入损伤模型中需要用到的参数。进行静、动力学分析的过程包括两个阶段，首先是静载加载过程，之后是动载加载过程，从而可以进行静载与动载的混合加载模拟，在静载加载到一定程度后再加动载。

在进行静载加载过程的模拟时，按一定增量加载，每次的加载增量用一个比例系数来控制。每增加一次荷载，调用asesk利用分块技术形成并分解整体刚度矩阵，调用snodeforce来计算由于位移引起的结点力并累计到结点荷载上，调用foba来利用上下三角线性方程组计算结点位移增量，调用substrain来计算各单元的应变增量，调用dstress来计算各单元的应力增量，调用stressmain计算主应力。调用judge，根据最大拉应力准则判断在新的应力水平下是否有新单元发生拉伸损伤，并记录单元在加载历史上主应变的最大值与应力达到抗拉强度时的主拉应变，在单元有新的损伤时更新刚度矩阵，再求解结点位移增量与应力增量，这个迭代过程一直进行到没有新的单元损伤时再增加荷载。在每一次荷载增加时，判断是否由于损伤单元过多而导致试件大变形失稳，如果失稳则计算终止。

在对动载加载过程进行模拟时，结点总位移等于静载位移与动载位移之和，静载位移是静载加载完毕时引起的位移。动载计算是一个时间相关问题，在时域上采用Newmark方法进行离散。程序首先调用dyncoeff形成Newmark系数，调用formdm求出结点阻尼系数，调用eigv求得结构频率，之后，在每个时间步上进行计算，过程与静载计算过程相似。与静载不同的是，需要在结点荷载上叠加由于损伤静位移所引起的转移结点荷载，以及对应于结点加速度与结点速度的广义结点荷载，这两个过程分别通过调用dsforce与gpmc来实现。dsforce计算由损伤产生的荷载增量，gpmc形成广义结点荷载。asesk和dasesk分别调用demove子程序生成刚度和广义刚度的分块矩阵。