水工结构数值仿真中的有限元分析研究综述

随着我国水利水电工程建设的向西推进，在建或拟建超大规模大坝数量不断增多，再加上西南山区复杂地质条件、高地应力的影响，复杂力学问题求解分析的深入，水工结构工程科学计算的“瓶颈”日趋严重，传统Von Neumann体系下计算机计算能力已不能满足需求。即使将模型及算法简化，一次计算仍然需要几小时、几天乃至数十天的时间，严重制约着水利学科的发展。

　　因此急需采用已在其他诸如气象预报、分子动力学、新药研制等领域广泛应用的高性能并行计算，为水工结构工程的数值仿真研究提供新的技术支撑。虽然并行计算自20世纪70年代起就已有学者研究，但其应用于水工结构数值仿真中的时间较晚。水工问题主要以有限元为分析方法，因而目前研究主要针对有限元展开，并取得了一定成果。

　　1水工有限元并行计算的3种策略

　　在水工结构工程中，有限元分析已成为不可或缺的方法。针对有限元的计算步骤，常用的有限元并行计算方法可以分为方程组并行求解、区域分解方法和EBE（Element-by-element）方法[1].

　　1.1方程组并行求解

　　线性方程组是许多水工数值仿真问题的核心。理论上，有限元分析最终都归结到求解线性方程组Ax=b,A为整体刚度矩阵，x为待求结点位移向量，b为区域荷载向量。对于线性的结构力学分析，刚度矩阵A是对称正定的，而且往往是稀疏带状的。有限元计算过程中，方程组的求解占据了大部分计算时间。即使是对于最简单的静力分析，代数方程组的求解时间也占整个分析时间的70%以上，动力分析则要占到90%以上[2].

　　因此，将上述方程组求解并行化成为提高计算效率的简单、有效方法。对于方程组Ax=b的并行求解，有两方面工作需要做：一是并行计算刚度矩阵A的分解；二是并行求解相应的三角形方程组，此部分较容易实现。对于矩阵A的分解，可以采用Cholesky分解。

　　张健飞、姜弘道[3]研究了带状对称正定矩阵A在等带宽存储方式下的并行Cholesky分解算法，并采用提前发送（Send Ahead）的策略减少了通信时间。数值试验表明该算法具有较高的加速比和效率。但这种并行化的方法在计算过程中往往会引入大量的非零元素的填充，降低了矩阵的稀疏性，同时对内存需求也很大。当计算量大时，采用分布式并行计算的数据交换量也大，因而限制了数值仿真问题的规模。

　　1.2区域分解

　　由于在迭代过程中需要大量的数据交换，采用方程组并行求解会大大降低计算速度，导致并行效率的降低，尤其是在分布式内存并行系统中。因而对于分布式内存并行系统，采用区域分解算法有利于提升大型水工结构仿真计算速度。

　　区域分解算法（Domain Decomposition Meth-od）作为一种偏微分方程数值解的新技术在20世纪80年代骤然兴起。它是为适应并行计算机工作原理而生成的一种算法。

　　在并行计算中，区域分解算法又被称为“分而治之”策略，即将一个大型的计算问题依据某种方法将其分解成为许多个小问题，再利用对应多个并行机来求解这些小问题，最终将问题解决[4].在水工数值仿真问题中，通常将大坝、地基、水工闸门等规模较大的水工结构建模，将有限元网格分解为若干子结构，再在这些子结构上进行相应的有限元计算，最后汇总得到整个模型的解。水工有限元计算问题中，一个区域的系统方程为

　　同时，相应的荷载向量也按同样原理进行分解。这样，式（1）的求解就可转化到子区域上进行，从而整个有限元计算就可并行求解。

　　1.3 EBE方法

　　进行有限元计算时，采用方程组并行求解实现比较容易，但当问题规模增大时，内存的需求也会相应增大，此外还可能涉及内外存交换技术，导致计算速度的降低；而区域分解方法往往只能在计算区域比较规则且剖分的子结构较好时才能获得很高的并行效率。对于三维问题，区域分解算法仍是一个没有得到很好解决的问题[5].

　　有限元EBE（element-by-element）方法则能很好地解决上述方法所面临的问题。它不进行整体刚度矩阵和整体荷载向量的集成，所有计算均在单元一级上进行，因此可以实现有限元整个计算过程的并行化[5].

　　对于一个水工有限元数值仿真问题，在将求解区域进行有限元网格离散化后，得到如下方程组：

    假设Ae与be为单元e对整个问题的贡献，因而它们具有和整体矩阵相同的维数，从而可将式（4）改写成：

　　其中：

　　Ae为仅在与单元e相关的位置元素才不为0的极度稀疏矩阵。在实际计算时，仅需存储Ae的压缩形式。同时，边界条件经过处理也很容易施加到单元一级上。这样，所有的计算均可在单元一级上相互独立地进行，这也就是EBE的基本思想。

　　在求解过程中，需先对矩阵A进行预处理。适用于EBE方法的预处理技术可分为两类：一类是将对整体刚度矩阵的预处理技术应用到共轭梯度法（CG）中，包括Jacobi预处理和多项式预处理；一类是基于矩阵的分裂获取一个A的近似来作为预处理矩阵。EBE并行求解时，具体的计算过程可归纳如下：

　　1）对每一个单元，计算单元刚度矩阵和相应的右端荷载向量，并施加边界条件。

　　2）采用基于EBE策略的预处理共轭梯度法（EBE-PCG）求解方程组，得到单元节点位移。

　　3）利用单元节点位移计算单元应变和应力。

　　4）单元节点位移集成总体节点位移并输出。

　　其中，在采用EBE-PCG方法求解方程组时，各个处理器之间需要交换与边界节点（各处理器间有关联的节点）有关的数据，实现时可采用压缩形式进行传输。

　　当计算规模增大时，可以保证处理器间需交换的数据量不会随之增大。

　　2并行实现的方法

　　2.1 PC网络机群

　　PC网络机群主要是利用高速通讯网络将一组高性能PC机按某种结构连接起来，进行统一调度、协调处理，从而实现高效并行计算。PC网络机群具有可扩展性、高性价比、用户编程方便等特点。机群系统采用通用网络连接，因而当计算需要时易于扩大机群规模。

　　其次，PC机群采用廉价的PC机搭建而成，相对价格高昂的巨型机经济适用。

　　此外，PC机群多采用WINDOWS操作系统，因而科研人员在熟悉的编程环境下易于开发调试并行计算程序，并进行相应问题的求解。

　　基于上述特点，在土木水利工程方面，许多科研人员搭建了PC网络机群并进行了并行计算在该领域的应用研究。

　　张友良等[6]利用100Mbps高速以太网络联结若干台个人计算机搭建了并行系统平台，在此平台上基于消息传递机制的MPI系统进行并行有限元程序开发，对水布娅地下厂房三维模型进行了模拟。

　　茹忠亮等[7]搭建了PC机群，采用C++语言与MPI环境开发了并行有限元计算程序，利用6台计算机对247 871个单元的有限元模型进行了并行求解，并行加速比达到5.26,验证了程序的可靠性和高效性。

　　2.2高性能计算集群

　　高性能计算集群又称分布存储式并行计算系统，是由多个互相连接的独立计算机组成的系统。

　　有别于PC网络机群，高性能计算集群主要是指由商业公司开发的大型系统，如曙光集群、HP集群等。它 主 要 采 用Linux操 作 系 统，常 见 的 有SuSeLinux、Redhat Enterprise Linux等。

　　在并行计算技术的发展过程中，基于X86架构的服务器机群系统因其具备良好的可扩展性、更高的可管理性、更低的使用维护成本、更好的可使用性、更高的系统鲁棒性和更多的应用支持成为目前并行计算领域的主流技术。基于上述优点，高性能集群在数值仿真领域有着广泛的运用。

　　刘耀儒等[5]利用清华大学研制的高性能集群系统，采用C/C++语言和MPI通信标准编制了有限元程序，成功对水布垭大型地下洞室进行了大规模数值计算。
 

2.3 GPU通用计算

　　GPU是一种专门在个人计算机、工作站、游戏机和一些移动设备上进行图像运算工作的微处理器[9],是计算机显卡的核心。GPU通用计算（GPG-PU,General-purpose computing on graphics pro-cessing units）是指利用GPU的计算性能来实现一般意义上的并行计算，区别于在图形图像处理方面的应用。

　　GPU通用计算通常采用CPU+GPU异构模式，由CPU负责执行复杂逻辑处理和事务处理等不适合数据并行的计算，由GPU负责计算密集型的大规模数据并行计算。这种利用GPU强大处理能力和高带宽弥补CPU性能不足的计算方式在发掘计算机潜在性能、成本和性价比方面有显着的优势。

　　自从SIGGRAPH 2003大会首次提出GPU通用计算概念以来，GPGPU经历了飞速的发展。CU-DA架构、AMD APP架构等的先后出现大大方便了编程人员的工作，对GPGPU在各个领域的成功应用起到了很大的作用。

　　在工程方面，刘明贵等[10]综述了GPU各类计算模式的发展特点与成果，展示了其在坝区渗透特性中随机微分方程加速求解过程中的优越性，探讨了采用GPU进行大规模岩土工程数值计算的应用前景；李 路 涛[11]采 用CPU-GPU混 合 平 台 已 将GPU成功应用于大规模岩土工程有限元计算中的迭代加速求解；韩博等[12]采用CPU-GPU协同计算技术，建立城区震害的计算模拟模型，成功完成了城区建筑群震害模拟；董廷星等[13]将GPU应用于流体力学中Navier-Stokes方程、Euler方程的加速求解上，经测试获得了33.2倍的加速比。

　　由上可以看到，GPGPU已经在岩土力学、流体力学方面有了相应的研究应用，并取得了良好的效果，预计随着研究的深入、GPU处理能力的提升与编程环境的优化，GPGPU也将会在水工结构仿真中得到广泛的应用。

　　目前，GPU高性能计算在NVIDIA、AMD等公司的推动下正快速发展。2013年11月由国际TOP500组织发布的超级计算机TOP 500榜单[14]中，前10名中已有No.2Titan和No.6Piz Daint采用NVIDIA GPUs来加速计算。

　　此外由我国自主研发的天 河 一 号A超 级 计 算 机 也 采 用 了7168块NVIDIA Tesla M2050计算卡加速计算，并已在石油勘探数据处理、生物医药研究、航空航天装备研制等[15]领域得到广泛运用。由此可见，采用GPU来加速计算正逐步得到业界的认可，相信随着应用环境的改 善，GPU将 成 为 水 工 数 值 仿 真 计 算 的 硬 件基础。

　　3高性能计算在水工计算仿真中应用的现状

　　3.1高拱坝地震损伤计算

　　随着我国西部山区一系列高坝的建设，特别是2008年汶川地震以来，高拱坝的抗震性能与在地震中的损伤破坏受到越来越多的关注。我国在建或拟建的200~300m级高坝大部分采用混凝土技术方案，从混凝土细观结构入手研究大坝损伤开裂非常热门且卓有成效。

　　由于混凝土具有复杂的细观结构，因此存在着计算规模急剧增大的问题。针对该问题，钟红等[16]利用并行计算技术研究了高拱坝在地震中的损伤破坏。针对坝高210m的大岗山拱坝，拱坝-坝基系统全部采用8结点等参单元进行剖分，并将坝体部分的尺寸控制在3 m以内。整个模型剖分为358 956个单元，386 276个结点，自由度数在百万以上。

　　采用区域分解方法，将问题“分而治之”,从而将原问题分解为若干可独立求解的小问题。在4节点组成的集群上进行并行计算，每节点配置相同，CPU为PIV3.0,内存2GB,硬盘160G.计算结果显示，随着参与计算节点数的增多，加速比相应增大。在4节点计算时，加速比大于3,可认为已达到较高的加速比。最终的非线性动力计算结果与实验结果吻合性良好，说明了并行计算在拱坝损伤计算中的有效性。

　　3.2拱坝-地基系统的有限元并行计算

　　针对三维的拱坝稳定和破坏过程的有限元分析，刘耀儒等[17]基于Jacobi预处理共轭梯度法，采用EBE方法，通过并行计算进行了二滩拱坝-地基系统的有限元分析。

　　对二滩拱坝-地基系统有限元模型，采用EBE方法，将并行计算建立在单元一级的子结构上，并依据单元编号次序将单元依次分配给计算节点处理，实现了并行任务的自动分配，有效解决了区域分解方法在处理三维复杂拱坝-地基系统问题上的难题。计算结果显示，随着参与计算的节点数的增多，并行效率逐渐降低，但当CPU数为16时，并行效率依旧可达到43%,说明了EBE方法对三维有限元计算的有效性。

　　3.3闸门结构分析

　　平面钢闸门是一个复杂的空间结构。在大型水利工程中，由于需要精确的计算才能保证结构的安全，因而需要将闸门结构进行三维分析，扩大了计算的规模。

　　付朝江[18]针对大型工程闸门结构计算规模巨大的情况，将闸门面板、主横梁腹板、纵梁腹板、边梁腹板等离散为板壳单元；次梁、主横梁和纵梁翼缘等离散为梁单元。并以各梁系的中心线为分界线将钢闸门进行分块，继而在各梁格区间将单元细分，从而建立钢闸门有限元计算模型。并行计算时，将单元刚度的形成与方程组的求解进行并行处理。

　　其中单元刚度的并行处理采用隐式区域分解技术，而方程组的并 行 则 采 用 预 处 理 技 术。最 终 在 由2台 双 核DELL工作站用以太网连接而成的机群上实现并行计算，采用的是消息传递标准（MPI）。计算结果显示，双 核 计 算 时，加 速 比 为1.41,并 行 效 率 为70.5%;4核计算时，加速比为2.45,并行效率可达61.2%.

　　3.4蒙特卡罗随机有限元的并行计算

　　蒙特卡罗（Monte Carle）方法是指利用生成的随机变量的样本作为输入获得功能函数的样本，再统计失效区样本的数量从而获得失效概率的一种方法。它具有概念明确、模拟次数与精度可确定的优点。

　　但为精确估算失效概率，所采用的样本数必须足够大，特别是在功能函数没有解析式和失效概率比较小的情况下，因而蒙特卡罗模拟的计算量将变得非常大。

　　此时，采用并行计算进行蒙特卡罗模拟则可大大缩短计算时间。工程结构要求具有一定的可靠性，因为工程结构在设计、施工、使用中存在种种影响其安全的不确定性。在边坡工程中，评价边坡稳定问题往往采用定值分析的方法，如简化毕肖普法等，它们均采用极限平衡理论，是一种粗糙的综合性的近似解。在实际过程中，边坡稳定受到各种不确定、随机因素的影响，如土体性质就具有非均质性。

　　为此，许多学者致力于随机有限元理论在边坡工程中应用的研究。

　　而将蒙特卡罗方法和有限元技术结合起来，把输入的各土工参数均看作具有一定统计特征的随机变量进行有限元分析，能较好地模拟边坡的实际工作状态[19].

　　4水工高性能计算展望

　　随着国家水利建设的不断推进，水利工程规模越来越大，其复杂程度越来越高，对应的数值仿真问题面临着非线性、高复杂边界条件等难题；同时，水利工程的信息化施工、坝体位移的实时监测和反馈分析、地震作用下大坝安全的实时分析等均对水工高性能并行分析提出了越来越高的要求。

　　近年来，随着信息技术的高速发展，信息化已成为各行各业新的发展方向。自20世纪90年代起，我国水利行业在信息化建设方面进行了大量探索尝试。

　　进入新世纪以来，以“国家防汛指挥系统”和“数字流域”的建设为代表，水利信息化进入了快速发展阶段。现阶段，云计算逐渐成为信息化产业发展新方向，其在水利中的应用也逐渐引起了相关部门、研究人员的重视[20-22].

　　所谓云计算（Cloud Computing），即指基于网络的计算，它旨在利用网络连接将分散的、远程的计算资源统一管理和调度，构成一个具有强大计算能力的计算池，从而向用户提供服务。作为并行计算、网络计算的升级发展，云计算技术具有以下几大特点：

　　1）服务性；2）资源整合促使效率提高，虚拟化使得整合多个IT资源成为可能；3）可靠性；4）高可用性；5）高层次的编程模型；6）经济性。

　　在水工计算方面，云计算可在水利工程施工、大坝位移检测分析等方面发挥极大的作用。水利施工过程中，实时监测回馈分析是保证水利工程建筑安全的重要因素。大坝的浇筑高程、混凝土温度、地基应力变形等均对水工建筑物的安全施工产生影响。因而采用云计算，实时远程了解掌控施工过程中各项指标的变化，并作出实时相应分析，对保障施工过程的安全与施工质量有着重要的作用。

　　其次，已建水工建筑物的监测也是云计算重要的应用。近年来，各地地震灾害频发，许多高坝、水库建在地震多发的西部高山峡谷区，针对此类大坝的监测关系到大坝的安全运行与当地居民的生命财产安全。

　　基于云计算，可实时监测地震波影响下坝体位移等指标，从而作出反馈分析，对保障大坝安全有着重要作用。目前，随着2009年《水利信息化顶层设计》的完成，水利信息化正逐步深入推进。

　　基于云计算的水利信息化发展已逐步受到相关部门的重视，研究云计算在水工数值仿真方面的应用也是水利信息化的重要组成部分。

　　加大对该问题的投入力度，发展水利高性能计算，可以促进水工结构工程的发展，为水工结构工程的数值仿真研究提供新的技术支撑。

　　参考文献：

　　[1] 刘耀儒，周维垣，黄岩松，等。水利工程中的大规模数值计算[C]//成都：中国岩石力学与工程学会第8次学术大会，2004:392-395.

　　[2] 林志祥。混凝土大坝温度应力数值仿真分析关键技术研究[D].南京：河海大学，2005:56-74.

　　[3] 张健飞，姜弘道。带状对称正定矩阵的并行Cholesky分解及其实现[C]//香港：第7届全国并行计算学术交流会，2003:262-265.

　　[4] 张友良，冯夏庭，茹忠亮。基于区域分解算法的岩土大规模高性能并行有限元系统研究[J].岩石力学与工程学报，2004,23（21）：3636-3641.

　　[5] 刘耀儒，周维垣，杨强，等。三维有限元并行计算及其工程应 用 [J].岩 石 力 学 与 工 程，2005,24（14）：2434-2438.

　　[6] 张友良。岩土工程大规模高性能并行有限元系统研究[D].武汉：中国科学院武汉岩土力学研究所，2005:58-93.