岩土工程有限元快速解法的研究与应用

岩土工程有限元快速解法的研究与应用吴梦喜1余进1（1．中国科学院力学研究所，北京100190）摘要：岩土工程有限元计算往往需要求解很大规模的问题，这对有限元软件的求解规模和求解速度提出了更高的要求。线性方程组的求解方法，在很大程度上决定了有限元的求解规模和速度。本文对线性方程组快速求解器进行了研究，重点研究了稀疏直接求解器和迭代求解器在岩土工程有限元（包括渗流有限元和应力变形有限元）计算中的应用，对它们的优缺点进行了探讨。研究表明，稀疏直接求解器具有稳定可靠的优点，无论对于渗流问题，还是应力变形问题，都可以高效快速的求解，可以求解的规模也相当大。对于超大规模的问题，宜根据具体的问题采用不同的迭代解法。对于渗流问题，形成的线性方程组的系数矩阵是对称正定的，采用无填充不完全分解预处理共轭梯度法（ICCG（0））或者对称超松弛预处理共轭梯度法（SSOR-CG），可以高效的进行求解。对于应力变形问题，形成的线性方程组的系数矩阵可能是对称不定的，使用ICCG（0）和SSOR-CG不能保证收敛；采用双门槛不完全分解预处理Krylov子空间法可以求解一些问题，但是依然存在参数难以选取的问题。在以上基础上，将快速求解器引入到课题组开发的有限元计算软件中，并应用于实际工程的分析，达到了显著提高求解速度和求解规模的目的。关键词：有限元；稀疏直接求解器；Krylov子空间法；预处理技术0前言随着有限元技术的发展，需要求解实际问题的规模越来越大。有限元法中形成的大型线性方程组的快速求解是有限元计算的关键技术，一直是研究的热点。研究人员常常使用自主开发的有限元软件分析实际工程问题。在方程组求解方面，这些软件很多都还在采用变带宽一维数组存储方程组并使用传统的三角分解法求解。这限制了软件所能求解问题的规模。对于规模稍大的问题，往往求解速度过慢，动辄需要计算几天才能得出结果。对于非线性或动力问题，这种求解方法的劣势更加明显。因此，研究求解规模更大、速度更快的的线性方程组的求解器，并将其引入到自主开发的有限元软件中，是十分迫切和必要的。近年来，线性方程组快速求解技术得到了很大的发展，很多的商业有限元软件都采用快速求解技术，求解规模和求解速度得以大大提高。消化吸收这些先进的技术，并将其引入到自主开发的有限元软件中来，是十分必要的，也是完全可行的。有限元计算形成的线性方程组一般来说是大型稀疏对称的，针对不同的问题，方程组的系数矩阵可能是正定的，也可能不正定。针对线性方程组稀疏的特点，采用全稀疏存储策略[1]，只存储系数矩阵非零元素的位置和数值，使大型稀疏线性方程组的存储量得以大大减小。使用稀疏直接解法或者迭代解法，可以实现大型线性方程组的快速求解。稀疏直接解法是20世纪70年代发展起来的一种线性方程组的直接求解方法。无论系数矩阵病态与否、正定与否，都可以得到稳定可靠的结果。稀疏直接解法也采用三角分解，这样在分解过程中在一部分零元位置将产生非零元，这些非零元叫做填充元，要增加存储空间。在进行三角分解之前有一个填充元优化的过程。利用图论的原理，采用最小度算法[2]或者嵌套剖分法[3]-[5]，使分解过程中产生的填充元尽可能的少。然后通过符号分解（SymbolicFactorization），确定分解因子的非零结构，并按照这个结构进行数值分解（NumericFactorization），然后回代求出变量。相比于传统的三角分解法和波前法，稀疏直接解法可以求解的问题规模更大，求解的速度更快。稀疏直接解法的缺点也是明显的，虽然填充元优化可以减少很多的填充元，但是三角分解会破环系数矩阵的稀疏性，分解因子的非零元远比系数矩阵要多，一般是系数矩阵非零元的10倍以上，依然需要占用大量的存储空间。因此当求解问题的规模很大时，迭代解法会更有优势，有时候由于内存的限制，甚至别无选择。Krylov子空间法是一类基于Krylov子空间的迭代方法，适合求解大型稀疏线性方程组。包括适用于求解对称正定方程组的共轭梯度法(ConjugateGradient，简称CG)，适合于求解对称不定问题的极小残量法(MinimalResidua1，简称MINRES)，适合求解一般问题的广义极小残量法(GeneralizedMinimalResidua1，简称GMRES)、双共轭梯度法(Bi-ConjugateGradient，简称Bi-CG）、拟极小残量法(Quasi-minimalResidualMethod，简称QMR)等。Krylov子空间法收敛速度依赖于系数矩阵的特征值分布，分布越集中，收敛越快。有限元计算形成的线性方程组的特征值分布往往很分散，因而常常需要采用预处理技术，对原方程组进行同解变换以改善特征值分布，从而加快迭代收敛速度。结合预处理技术的Krylov子空间法具有良好的收敛特性和较高的数值稳定性，已成为目前求解高维稀疏线性方程组的主要方法。在岩土工程有限元分析中，存在渗流、应力变形及其耦合计算。渗流计算中形成的线性方程组往往是对称正定的，应力变形有限元计算中，由于土与混凝土结构物接触界面的存在，形成的线性方程组有可能是对称不定的。针对这两类问题的具体特点，将线性方程组快速求解器引入到岩土工程有限元计算中，提高有限元软件的计算规模和计算速度，是本文要介绍的主要内容。1稀疏直接求解器稀疏直接求解器（SparseDirectSolver）是求解规模不太大的方程组的一种有效的求解方法。它充分利用刚度矩阵的稀疏性，运用以图论为基础的稀疏矩阵重排序技术。改变稀疏矩阵的行和列，达到减少填充元的目的。实现一个高效的稀疏直接求解器往往是很复杂的，涉及很多方面的知识，往往不是一个工程师能够在短期内实现的。有很多的研究小组，开发出很多高效的稀疏直接求解器。这些求解器采用了复杂的填充元优化算法，调用了高性能数值计算库，具有很高的求解效率。这些求解器有商业的、免费不开源的和免费开源的。它们在不同的软件开发集成环境中提供了简单的接口，很容易让软件开发者调用。直接调用这些求解器，一方面可以避免重复的工作，一方面这些求解器有很多的实际使用经验，具有很高的可靠性。表1列出了常用的一些稀疏直接求解器的情况。表1几种稀疏直接求解器数值类型接口语言版本矩阵类型性质求解器名称实数复数FortranC串行并行对称一般免费开源DSS√√√√√√PARDISO√√√√√√√√√MUMPS√√√√√√√√√√HSL√√√√√√√√SuperLU√√√√√√√√√Umfpack√√√√√√√√TAUCS√√√√√√√√GSS√√√√√√√为了将合适的稀疏求解器集成到岩土工程有限元计算软件之中，本文选取了三个有代表性的求解器进行测试：CXML的DSS、Intel的PARDISO、MUMPS。CXML的DSS可以在CompaqVisualFortran编译器里直接调用；IntelPARDISO可以在IntelVisualFortran编译器里面直接调用；MUMPS主软件使用F90编写，部分代码使用C编写，可以在任何可以混合编程的平台下（如VisualStudio）编译成链接库然后在该平台下调用。数值测试平台为微机，其CPU是Pentium(R)D3.40GHz，1G内存，WindowsXP操作系统。选取两组测试矩阵，第一组矩阵来自岩土工程应力变形有限元计算，第二组来自TimDavis的UFSparseMatrixCollection，可以在其首页上[6]下载得到，这些矩阵都是实际的各个领域有限元计算中形成的。两组矩阵的描述见表2、表3。表2第一组测试矩阵编号阶数非零元个数性质114106933413对称正定2155581050618对称正定3176431201668对称不定4184571280016对称不定5193891331962对称不定617445617381对称不定表3第二组测试矩阵问题描述名称阶数非零元个数性质高温压力容器结构分析bcsstk1710974219812对称正定最优化问题cvxbqp150000199984对称正定最优化问题dixmaanl60000179999对称不定PDEdubcova265025547625对称正定结构分析oiplan737521835470对称正定非线性优化c-7176638468096对称不定壳体分析s3dkq4m2904493686223对称正定CFDcfd21234401605669对称不定汽车结构分析bmwcra_11487705396386对称不定结构分析shipsec51798605146478对称正定汽车结构分析bmw3_22273625757996对称不定结构分析msdoor41586310328099对称正定批注[wumx1]:不要边框，图标、轴标要更清晰效果a.第一组测试矩阵b.第二组测试矩阵图1三个稀疏直接求解器求解结果令两组测试矩阵的所有变量的值为1，计算得到方程组的右端项，然后使用三种稀疏直接求解器进行求解，求解时间见图1。从图中可以看出，第一组测试矩阵求解时间都在15s以下，第二组测试矩阵的求解时间都在150s以下。不管线性方程组是对称正定，还是对称不定，稀疏直接求解器都可以快速的求解，且求解时间并不象常规直接解法那样随线性方程组规模的增加而快速增加。在本文所求解的方程组中，最大的含有41万个变量，一千万个非零元。可见，稀疏直接求解器具有稳定可靠的优点，对于各种问题都能够求解，而且求解速度也非常快。研究人员可以根据自己的操作系统、编程语言、编译器，灵活的选择合适的稀疏直接求解器。本课题组开发的力科渗流软件seepage和力科应力变形计算软件3d_stress，都配置了稀疏直接求解器PARDISO。图2是采用3d_stress软件计算某坝体应力变形的有限元模型。该模型采用二次单元，单元总数3528个，节点总数16678个，计算总共分21个加载级。分别采用传统的传统的乔列斯基分解求解器和稀疏直接求解器进行求解。总的求解时间，使用传统直接解法为94472s，约26小时，其中求解方程组花费了82925s，占用了89%的计算时间；使用稀疏直接解法为5072s，约1小时24分钟，其中方程组求解的总时间不到5分钟，只占6％。可见稀疏直接解法能够大大提高求解速度。图2某坝体三维有限元网格2迭代求解器对于超大规模的问题，迭代解法往往更有优势，有时候甚至是唯一的选择。迭代解法的收敛性往往依赖于系数矩阵的性质，因此，对于不同的问题类型，例如对于渗流有限元计算和应力变形有限元，往往需要选用不同的迭代求解器。2.1 渗流有限元计算 渗流有限元计算形成的线性方程组的系数矩阵是对称正定的，可以采用共轭梯度法求解，预处理方法可以采用无填充的不完全乔列斯基分解（IC（0））预处理[7]或者对称超松弛（SSOR）预处理[8][9]，这两种方法叫做ICCG（0）求解器和SSOR-CG求解器。在本课题组进行的大量实际工程渗流问题的计算中，这两种方法均可以高效的求解。本课题组使用ICCG（0）求解器，目前求解问题的最大规模是二次单元30万节点的一个模型。表4列出了采用课题组渗流计算软件seepage求解某坝体渗流时这两种求解器的迭代次数。在非线性求解过程中，SSOR-CG用了4个迭代步，ICCG（0）用了3个迭代步。在一开始的迭代步，SSOR-CG的迭代次数比ICCG多。随着非线性迭代的进行，SSOR-CG的迭代次数减小比较快，在第三个非线性迭代步其迭代次数已经少于ICCG（0）了。表4求解某坝体渗流问题时迭代解法的迭代次数比较非线性迭代步1234SSOR-PCG7201268540ICCG（0）105101101图3是采用seepage计算某坝三维渗流场的模型。该模型采用二次单元，总共有68723个单元，290998个节点。分别使用ICCG（0）和SSOR-CG求解器进行求解，迭代精度为10-9，非线性迭代步设为10步。总的求解时间，ICCG（0）约为2个小时，SSOR-CG约为6个小时。从渗流问题的求解来看，ICCG（0）的收敛速度比SSOR-CG收敛速度要快。SSOR-CG的主要优势在于预处理时