Aztec在混凝土细观数值模拟中的应用研究

吴建平，赵军，宋君强，张卫民，马怀发

1.国防科技大学计算机学院，长沙410073

2.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京100038

1 介绍

混凝土性能研究是混凝土大坝地震破坏机理研究的关键技术之一，细观数值模拟作为一种重要研究手段，将混凝土试件看成由骨料、沙浆及其交界面组成的三相复合材料[1-3]。在模拟试件时，必须对各种材料的分布状况进行假设，假设骨料分布、骨料尺寸等参数均具有随机性，即随机骨料随机参数模型已经证明很有效[2-3]。该模型以有限元离散为基础，模拟过程中需要求解对应于全局刚度矩阵的稀疏线性方程组，以算出每一个加载步上的三个位移分量。该线性方程组的求解在总模拟时间中所占比重非常大，已经超过90%[4]。

对三维试件数值模拟中的稀疏线性方程组，采用直接解法时，不仅计算量很大，内存需求也很大。而且，为减少计算量和存储需求，常采用排序技术来减少矩阵分解过程中的填入元数量，但即使采用排序技术，对三维问题，计算量和内存需求依然相当庞大，难以接受。与直接解法相反，迭代法只要给定初始解向量和迭代计算规则，就可以依次不断计算出新的迭代值，在达到一定的停机准则时就可以停止迭代。迭代过程中，主要的计算量在矩阵向量乘上，这不会改变稀疏矩阵的非零元结构，所以，计算量与系数矩阵中的非零元个数成正比，且只需要存储系数矩阵与一些相应的向量，从而计算量和存储需求都相对较少[5]。因此，对三维问题，采用迭代求解是最佳选择。但迭代法面临收敛速度慢，甚至不收敛的问题。为此，需要引入预条件技术，以加速收敛过程。

Aztec是一个迭代求解线性方程组的软件包[6-7]，该软件包试图避开烦琐的并行编程细节。在创建分布矩阵、相应的右端向量和解向量，并进行相应的参数和选项设置后，就可以调用Aztec的求解器接口进行求解。Aztec包含多种K rylov子空间迭代法，包括CG、GMRES、BiCGSTAB等，这些方法可以与各种预条件结合使用，采用的预条件包括多项式预条件，在子区域内采用ILU的区域分解预条件，在子区域内采用GS迭代的预条件，以及基于最小二乘的预条件等。

虽然Aztec已经实现了并行计算，但是该软件包内部并未提供图分割接口。同时，由于混凝土试件采用的是无结构网格，而已有研究表明，采用M eTis软件包中基于嵌套二分的图分割算法最为有效[8]。因此，本文先利用图分割算法对混凝土试件进行图分割，之后在该分割下导出相应的分布矩阵和分布向量，相应的矩阵装配采用基于稀疏数据结构的高效并行算法[4]，之后将稀疏矩阵和右端向量导入到Aztec的求解器接口子程序中进行求解。

2 混凝土试件细观数值模拟

在混凝土细观数值模拟中，将混凝土试样看成是由骨料、灰浆以及其间界面组成的三相复合材料，这里以随机骨料随机参数模型为基础[1]。立方体试件用有限元离散，含骨料的区域用较小的六面体单元进行离散，只含灰浆的区域用较大的六面体单元进行离散，这两类离散区域用四面体单元连接。

所有有限元分为三类：骨料单元、灰浆单元和界面单元，分别对应于材料的三相。对每个有限元，要么有8个形函数，要么有4个形函数，这由所对应的单元是六面体单元还是四面体单元来确定。有限单元的每个节点对应于应力或应变的3个分量。对每个有限元应用损伤模型，并对每个节点应用简单的弹性关系，之后，就可以形成单元刚度矩阵。当每个节点上的荷载已知之后，节点位移可以通过求解对应于全局刚度矩阵的稀疏线性方程组来获得。

2.1 静载时的数值模拟

对静态加载，荷载一步一步地增加。在加载过程中，如果一个单元所受应力超过其所能承受的极限强度，单元刚度就会退化，即刚度是损伤参数的函数。因此，平衡方程是非线性的，采用增量方法进行求解。由于在荷载更大时产生了新的损伤，求解过程需要进行迭代。

用[Ki]、{Ui}、{Pi}分别表示单元刚度矩阵、位移和第i个加载步时的静荷载，则第i个加载步时的平衡方程可以写为：

且第i+1个加载步时的平衡方程可以写为：

因此

这里

且e为单元，e为单元的应变，d是对应于应变的损伤系数。这样，第i个加载步的求解过程可描述如下：

2.2 动态加载时的数值实验

在动态加载时，t时候与t+Dt时候的平衡方程可以分别写为：

因此

其中

[M]为单元质量矩阵，[C]为单元阻尼矩阵，[Kd]为单元动态刚度矩阵，[Ks]为单元静态刚度矩阵，{Pd}为动态荷载，[Ps]为静态荷载。选取δ≥0.5，γ=0.25(0.5+δ)2，并记

则方程（3）可以近似为：

方程（4）是非线性的，可以利用下述迭代进行近似

其中

除了需要求解形如式（5）的线性方程组外，对[M]和[C]的计算还需要计算[K]的某些特征值，这也需要求解相应的线性方程组。

3 Aztec软件包

Aztec是一个迭代软件包，用于并行求解稀疏线性方程组

其中，A是n×n稀疏矩阵，b是n维向量，x是要求解的n维向量。Aztec试图避开烦琐的并行编程细节，因此最复杂的环节是为特定应用问题进行矩阵分布状态的指定。一旦创建了分布矩阵，就可以在任意运行Aztec的并行计算机上进行计算。Aztec提供了两种稀疏矩阵格式，一种是逐点修正型稀疏行DMSR格式，另一种是逐块行DVBR格式，分别是MSR格式和VBR格式在分布存储环境下的版本。由于现有混凝土细观数值模拟程序中采用CSR格式，且CSR格式转换到MSR格式代价很低，因此，文中采用DMSR格式进行分布矩阵指定。

任意长为n的向量通过某种分割方法，将元素分配到特定的处理器上。当计算向量y=Ax中的分量时，一个处理器仅计算分配到其上的y分量，这些分量显式地存储在该处理器上，分量对应的全局标号存储在update索引集中。update又分为两个子集，即internal和border。集合internal中每个索引对应的分量只需要利用本处理器上的信息进行更新。集合border中含有的索引对应于要在矩阵向量乘中进行更新，但需要从其他处理器获取x分量值的y分量。在矩阵向量乘中，不在本处理器上，但又需要用来更新border中y分量的x分量指标存储在索引集external中。这些x分量显式存储在其他处理器上，且在进行矩阵向量乘时，需要从这些处理器经由通信来获得。Aztec在对向量进行存储时，先存储internal对应的分量，之后存储border对应的分量，最后存储external对应的分量。同时，所有从同一个处理器接收到的external分量进行连续存储。

Aztec包含许多K rylov子空间迭代法，包括CG、GMRES、BiCGSTAB、CGS和TFQMR。这些K rylov方法可以与各种预条件结合使用，可用的预条件包括Jacobi迭代，Neumann展开多项式，最小二乘多项式，非重叠区域分解k步对称化GS和区域分解预条件，在区域分解预条件中，每个子区域上的局部预条件可以选用LU分解、ILUT、ILU(k)、RILU(k，ω)、BILU(k)和ICC(k)，而对子区域解到全局解的组合可以采用由拥有者确定最终分量值的标准方式和对分量分别各自进行相加的对称化方式两种方式，区域分解中的重叠度可以采用参数来进行指定。

4 数值实验

针对非结构网格上混凝土细观数值模拟所得稀疏线性方程组，就Aztec软件包中提供的预条件和迭代法组合进行实验。由于这里的稀疏线性方程组全都对称正定，因此迭代法全部选用CG法。同时，为利用预条件迭代进行求解，预条件本身必须具有对称性，因此仅对Aztec中具有对称性且相对比较有效的预条件进行测试，包括区域分解预条件，非重叠区域分解k步对称化GS即SGS（k）和最小二乘k次多项式即LS（k）。区域分解预条件中采用对分量各自分别进行相加的对称化方式，对重叠度0即无重叠的DDM（0）与重叠度1的DDM（1）进行测试，且对局部是否采用RCM排序进行测试，在子区域局部对具有对称性的ICC（0）和ICC（1）进行测试。在所有实验中，初始猜测都选为0向量，且迭代的停机准则选为残量2范数下降10-4倍。

所有测试结果均在一机群系统上得到，该系统采用Intel®Xeon®CPU X 7350@2.93GHz（Cache 4 096 KB），且用Infiniband进行互连，操作系统为Red Hat 4.1.2-44版的Linux，消息传递接口采用MPICH2 1.4.0rc2，所用编译器为Intel Fortran 10.1.018。

实验针对两个混凝土试件进行，离散网格形状参见文献[3]，试件1是一个三级配混凝土，其大小为1 100mm×300 mm×300mm。离散节点有44 117个，有限单元有53 200个。两条支撑轴分别位于x=-0.36m且z=0m，x=0.54m且z=0m处，即处于底部且分别离左右边界0.1m处。两条加载轴分别位于x=-0.06m且z=0.3m，x=0.24m且z=0.3m处，即处于试件顶部且分别距左右边界0.4m。试件2是一湿筛混凝土，其大小为550mm×150 mm×150mm。离散节点数为71 013，有限单元数为78 800。两支撑轴分别位于x=-0.15m且z=0m，x=0.3m且z=0m处，两加载轴分别位于x=0m且z=0.15m，x=0.15m且z=0.15m处。

在对两试件的静载实验中，每步增加的荷载为0.25 kN。当在某步出现损伤单元时，在该步上就需要连续求解多个线性方程组，以校正由刚度矩阵退化引起的位移。对试件1，在第439步出现损伤单元，并在第567步完全失效，总共需要求解696个线性方程组。对试件2，在第59步出现损伤，且在第94步完全失效，总共需要求解178个线性方程组。在对试件2的动载实验中，时间步长取为0.001 s，且荷载线性地以增量0.8 kN的方式增加，在第38个加载步完全失效，由于实际计算之前，需要计算某些特征值以确定相关参数，这也需要求解线性方程组，所以模拟过程总共需要求解221个线性方程组。

对试件1、试件2的静载实验，以及对试件2的动载实验结果分别列于表1、表2、与表3，其中时间单位为s，“—”表示迭代未能正常结束。从表中实验结果可见，采用RCM排序时，从迭代次数上看，对LS和SGS预条件，与不采用RCM排序时完全相同。对基于区域分解的并行ICC预条件，迭代次数一般稍有减少。从计算时间上看，仅在局部采用ICC（1）且采用重叠区域分解时才具有改进效果，在其他情况下求解时间一般反而较长，这可能是因为并行预条件的有效性需要全局协调一致才能充分体现所致，即使局部预条件很有效，如果区域分解时策略适当，导致子区域间对应的舍弃元素相对较大，将对全局计算效果的改进形成很大限制。

此外，从实验结果中还可以看到，当处理器个数增加1倍时，有时出现求解时间减少1倍多，即超线性加速比现象，这应是Cache的影响所致。同时，可以注意到，在很多情况下，在采用相同个数处理器时，ICC（0）结合无重叠区域分解时的求解时间在所有预条件方法中最短，这为以后进行针对混凝土材料细观数值模拟时，进行预条件技术的研究提供了比较基准，即当与Aztec软件包进行求解性能比较时，只要与其中不采用RCM排序的无重叠区域分解型并行ICC（0）预条件CG迭代进行比较。

表1 对试件1进行静载模拟时平均每个线性方程组的迭代次数和求解时间

表2 对试件2进行静载模拟时平均每个线性方程组的迭代次数和求解时间

表3 对试件2进行动载模拟时平均每个线性方程组的迭代次数和求解时间

5 结论

本文利用Aztec软件包中的各种预条件技术，对目前混凝土细观数值模拟中大型对称正定稀疏线性方程组的求解，进行了实验比较研究。结果表明，在基于区域分解的并行ICC、无重叠对称化GS迭代、最小二乘等预条件技术中，采用无重叠区域分解型并行ICC时求解时间最短，且采用RCM排序一般不仅不能缩短求解时间，反而常导致求解时间更长，因此，在本应用中不采用RCM排序更有利。

[1]唐春安，朱万成.混凝土损伤与断裂－数值模拟[M].北京：科学出版社，2003.

[2]马怀发，陈厚群.全级配大坝混凝土动态损伤破坏机理研究及其细观力学分析方法[M].北京：中国水利水电出版社，2008.

[3]马怀发.全级配大坝混凝土动态性能细观力学分析研究[D].北京：中国水利水电科学研究院，2005.

[4]吴建平，王正华，朱星明，等.混凝土细观力学分析程序中的快速算法与并行算法设计[J].计算力学学报，2008，25（3）：352-358.

[5]吴建平，王正华，李晓梅.稀疏线性方程组的高效求解与并行计算[M].长沙：湖南科技出版社，2004.

[6]Wu Jianping，Zhao Jun，Song Junqiang，et al.Impact of two factors on several domain decomposition based parallel incomplete factorizations for the meso-scale simulation of concrete[C]//Proceedings of the 3rd International Conference on Information and Computing Science（ICIC2010），Wuxi，China，2010.

[7]Shadid J N，Tum inaro R S.Aztec-a parallel preconditioned K rylov solver library：algorithm description version 1.0，Technical Report Sand95[R].Sandia National Laboratories，Albuquerque NM，87185，1995.

[8]Tum inaro R S，Heroux M，Hutchinson S A，et al.Official Aztec user’s guide，Version 2.1，Technical Report Sand99-8801J[R].Sandia National Laboratories，Albuquerque NM，87185，1999.