高壩抗震高性能并行計(jì)算研究及其應(yīng)用

（中國水利水電科學(xué)研究院，北京市 100048）

0 引言

高壩是復(fù)雜的空間結(jié)構(gòu)，其地震響應(yīng)分析需要同時(shí)計(jì)入壩體-地基-庫水的動態(tài)相互作用、壩體伸縮橫縫開合的邊界非線性、遠(yuǎn)域壩基的振動能量逸散、近域地基的復(fù)雜地形和地質(zhì)構(gòu)造、壩體和地基巖體材料的非線性、沿壩基地震動的不均勻性，以及在拱座和地基動力耦合作用下，壩肩潛在滑動巖塊的動態(tài)穩(wěn)定性等復(fù)雜問題。整個(gè)體系的有限元數(shù)值計(jì)算需要求解幾十萬、甚至幾百萬個(gè)自由度的方程組。大壩在動荷載下響應(yīng)的計(jì)算量隨著壩體自由度的增加，呈幾何級數(shù)增長。特別是非線性動力學(xué)問題，為保證數(shù)值計(jì)算的收斂性和穩(wěn)定性，時(shí)域計(jì)算中的時(shí)間步長常需取十萬分之幾秒，甚至在一個(gè)時(shí)間步長內(nèi)還要進(jìn)行多次迭代，而輸入地震動的持續(xù)時(shí)間長達(dá)幾十秒。例如具有120萬自由度的沙牌拱壩抗震計(jì)算應(yīng)用單機(jī)串行程序計(jì)算，運(yùn)用單機(jī)串行計(jì)算，需要34天才能完成。在時(shí)域內(nèi)求解如此大規(guī)模的非線性動力學(xué)問題，計(jì)算速度已成為高拱壩地震動分析和抗震安全評價(jià)的瓶頸。

高性能并行計(jì)算已經(jīng)成為解決大規(guī)?？茖W(xué)計(jì)算的主要手段。隨著我國超級計(jì)算機(jī)的發(fā)展，天河一號、天河二號、神威等超級計(jì)算機(jī)計(jì)算資源對社會開放，硬件資源使百萬千萬級自由度的大規(guī)模非線性動力計(jì)算成為可能，為精細(xì)化研究高壩抗震提供了條件。高壩抗震高性能并行計(jì)算軟件研發(fā)成為高壩抗震研究工作者的新課題。在2017年公布的國家重點(diǎn)研發(fā)計(jì)劃“300m級特高壩抗震安全評價(jià)與控制關(guān)鍵技術(shù)”申請指南中，將研發(fā)300m級特高壩大壩-地基-庫水系統(tǒng)非線性動力耦合并行計(jì)算仿真分析軟件系統(tǒng)列為考核指標(biāo)，體現(xiàn)了高性能并行計(jì)算在高壩抗震領(lǐng)域的重要性。

本文將主要介紹由中國水利水電科學(xué)研究院陳厚群院士領(lǐng)導(dǎo)的研究團(tuán)隊(duì)多年來自主研發(fā)的高壩抗震分析并行計(jì)算軟件（PSDAP）有關(guān)并行算法和架構(gòu)、分析模塊等，及其在2010年排名世界上最快計(jì)算速度的超級計(jì)算機(jī)“天河一號”上運(yùn)行的若干工程實(shí)踐。

1 計(jì)算模型與參數(shù)

1.1 區(qū)域分解法

區(qū)域分解法（DDM）的基本思想是采用“分塊”的策略，將一個(gè)復(fù)雜的計(jì)算系統(tǒng)分成若干個(gè)子系統(tǒng)，原系統(tǒng)的求解就轉(zhuǎn)化為子系統(tǒng)上的求解，各個(gè)子系統(tǒng)通過信息傳遞完成數(shù)據(jù)交換。

按各分區(qū)的重疊與否，區(qū)域分解法分為重疊型區(qū)域分解法和不重疊型區(qū)域分解法[1]。重疊型區(qū)域分解法的理論基礎(chǔ)是Schwarz交替法[2]，各分區(qū)間消息傳遞通過各分區(qū)的相關(guān)重疊區(qū)域來實(shí)現(xiàn)。不重疊型區(qū)域分解法的理論基礎(chǔ)是子結(jié)構(gòu)法，各分區(qū)間消息傳遞通過各分區(qū)交界面區(qū)域來實(shí)現(xiàn)。

本文以重疊型區(qū)域分解法（ODDM）為基礎(chǔ)進(jìn)行并行程序設(shè)計(jì)[3]。Schwarz交替法基于“化整為零”的思想把一個(gè)整體復(fù)雜區(qū)域分解為若干相互重疊的簡單分區(qū)求解，這就為分布式并行計(jì)算提供了數(shù)學(xué)基礎(chǔ)。下面以一個(gè)簡單的平面應(yīng)力靜力學(xué)問題為例來說明Schwarz交替法基本原理。如圖1所示，一個(gè)ABCD區(qū)域的平面應(yīng)力問題分解為ABGH和CDFE兩個(gè)相互重疊區(qū)域的平面應(yīng)力問題，EFHG為重疊區(qū)域。

圖1 平面應(yīng)力靜力學(xué)問題的重疊型區(qū)域分解示意圖Fig.1 The diagram of overlapping domain decomposition for plane static problem

定義ABGH區(qū)域：kuu=fu，u=gu（GH邊界）

定義CDFE區(qū)域：kvv=fv，v=gv（EF邊界）

交替求解方法流程如圖2所示。

圖2 Schwarz交替法基本流程Fig.2 The basic process of Schwarz alternating method

1.2 顯式求解格式波動方程的重疊型區(qū)域分解法

波動方程采用顯式求解，由于質(zhì)量矩陣為對角矩陣，方程是解耦的。顯式求解格式的波動方程求解有著大規(guī)模并行計(jì)算的優(yōu)勢，通過交換相鄰區(qū)域的重疊型區(qū)域數(shù)據(jù)以更新下一時(shí)步分區(qū)邊界的位移條件，然后由新的邊界條件求解。因此，顯式計(jì)算的重疊性區(qū)域分解法無需迭代，只需交換一次邊界信息更新邊界條件。

1.3 基于重疊單元的區(qū)域分解

基于節(jié)點(diǎn)分區(qū)的原則對有限元節(jié)點(diǎn)進(jìn)行分割后，求出每個(gè)分區(qū)的內(nèi)部單元和重疊單元以及按節(jié)點(diǎn)屬性進(jìn)行分類，具體流程如圖3所示。

圖4 重疊型分區(qū)Fig.4 Overlapping partition

圖5 分區(qū)節(jié)點(diǎn)分類Fig.5 Node classification of sub-domain

如圖4所示，每個(gè)分區(qū)的單元分為內(nèi)部單元和重疊型單元。如圖5所示，每個(gè)分區(qū)的節(jié)點(diǎn)分為內(nèi)部節(jié)點(diǎn)、邊界節(jié)點(diǎn)和外部節(jié)點(diǎn)，其中內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)是每個(gè)分區(qū)待求解節(jié)點(diǎn)。邊界節(jié)點(diǎn)和外部節(jié)點(diǎn)是與其他分區(qū)交換的信息，本分區(qū)的邊界節(jié)點(diǎn)則是其他某分區(qū)的外部節(jié)點(diǎn)，外部節(jié)點(diǎn)為本分區(qū)提供邊界條件，不在本分區(qū)求解。

1.4 并行計(jì)算結(jié)構(gòu)

采用能夠適用于大型機(jī)群的分布式并行計(jì)算，采用主從式編程模式，由一個(gè)主進(jìn)程和若干個(gè)從進(jìn)程組成。主進(jìn)程是一個(gè)控制程序，不參與計(jì)算，負(fù)責(zé)把數(shù)據(jù)發(fā)送給從進(jìn)程，并接收整理從進(jìn)程數(shù)據(jù)，所有消息傳遞發(fā)生在主進(jìn)程與從進(jìn)程之間，從進(jìn)程之間相互不發(fā)生消息傳遞。從進(jìn)程只負(fù)責(zé)對應(yīng)子區(qū)域的計(jì)算。消息傳遞采用計(jì)算與通信先后進(jìn)行的阻塞式通信來實(shí)現(xiàn)，如圖6所示。

圖6 并行計(jì)算結(jié)構(gòu)圖Fig.6 The diagram of parallel computing architecture

2 分析模塊

2.1 粘彈性人工邊界

有限元方法模擬無限域的波動問題中，應(yīng)盡量減小底邊界和側(cè)邊界的地震波反射。Lysmer and Kuhlemeyer（1969）[4]提出黏性邊界的方法來吸收反射到邊界上的地震波。對于黏性邊界可能引起相對較大的誤差和低頻失穩(wěn)問題，研究人員[5，6]提出了黏彈性人工邊界。在有限元方法中，底邊界和側(cè)邊界設(shè)為黏彈性人工邊界，底邊界和側(cè)邊界的節(jié)點(diǎn)上施加彈簧和阻尼器，如圖7和圖8所示。

在這項(xiàng)研究中，在數(shù)學(xué)上實(shí)現(xiàn)這些彈簧和減震器，可在邊界相關(guān)單元矩陣的對角項(xiàng)上增加彈簧和阻尼項(xiàng)，因此將在邊界節(jié)點(diǎn)x、y和z三方向上施加與位移和速度相關(guān)的力，相應(yīng)的彈簧和阻尼系數(shù)為：

以上式中E為彈性模型；G為剪切模量；ρ為密度；A為人工邊界節(jié)點(diǎn)影響面積；rb表示從邊界底部到頂部的距離；cp和cs分別為有限元模型外側(cè)介質(zhì)的壓縮波波速和剪切波波速。

圖7 黏彈性人工邊界模型Fig.7 Viscoelastic artificial boundary model

圖8 三維黏彈性人工邊界示意圖Fig.8 The diagram of three-dimensional viscoelastic artificial boundary

2.2 地震動輸入

在本項(xiàng)研究中，通過在人工邊界節(jié)點(diǎn)上施加等效荷載的方式來實(shí)現(xiàn)地震波動輸入的方法，即將地震波動的位移和速度時(shí)程轉(zhuǎn)換為等效節(jié)點(diǎn)荷載施加于人工邊界上，完成地震動的輸入。人工邊界上的地震動輸入為從底部入射的地震波和從地表反射的地震波兩者疊加形成的自由場。

如圖9所示，對于底部和側(cè)邊界的每個(gè)節(jié)點(diǎn)，自由場位移和速度是從底部入射波和地表反射波的位移和速度疊加。

圖9 箱型模型Fig.9 Box-shape model

因此，等效的邊界節(jié)點(diǎn)荷載時(shí)程可表示為：

其中，F(xiàn)B為邊界荷載向量；σ為自由場應(yīng)力張量；n為邊界外法線向量；CB為邊界節(jié)點(diǎn)阻尼系數(shù)矩陣；KB為邊界節(jié)點(diǎn)彈簧系數(shù)矩陣；u為自由場位移向量；˙u為自由場速度向量。

2.3 接觸非線性模型

接觸問題物理概念明確，接觸界面約束條件容易確定，其關(guān)鍵是對于這種不連續(xù)非線性問題的數(shù)值求解方法。本文基于邊界條件非線性的接觸力模型，以Lagrange乘子法為代表，Lagrange乘子代表界面上的未知接觸力。接觸力模型以在界面上施加作為未知量的接觸力來滿足界面約束條件，不存在人為假定剛度的問題，對界面的處理更符合實(shí)際。

離散后的波動方程可寫為：

其中，A為剛度矩陣；U為位移向量；F為荷載向量。

包含接觸力的平衡方程為：

其中，B為接觸約束矩陣；λ表示接觸力向量。

接觸面位移約束方程為：

其中，γ為位移約束條件向量。

聯(lián)立式（3）和式（4）可得接觸力方程，接觸力方程柔度矩陣一般為滿陣，法向接觸力和局部接觸力相互耦合，采用迭代算法求解，本文不再詳述。該模型嚴(yán)格滿足位移約束條件，可模擬壩體中橫縫縱縫以及壩肩山體內(nèi)斷層的黏結(jié)、張開、滑移等。

2.4 材料非線性

本軟件包括多種混凝土損傷分析模型、基巖損傷模型、塑性模型等，限于篇幅有限本文不再詳細(xì)介紹。

3 在“天河一號”上運(yùn)行的若干工程實(shí)踐

3.1 并行性能測試

本文采用互聯(lián)網(wǎng)遠(yuǎn)程登錄的方式利用“天河一號”提供的硬件資源，所有計(jì)算工作都在天河一號超級計(jì)算機(jī)上完成。

如圖10所示，沙牌拱壩-地基有限元模型節(jié)點(diǎn)數(shù)為425568，單元數(shù)目為404090，自由度為1276704。整體有限元模型分成1～191個(gè)子區(qū)域，對并行性能進(jìn)行測試，圖11為分成11個(gè)子區(qū)域的信息。本算例計(jì)算時(shí)間40s，一共200萬步，對并行程序性能進(jìn)行測試。

圖10 整體模型Fig.10 The overall model

圖11 11個(gè)子區(qū)域在總體區(qū)域的中的位置（不同顏色代表不同子區(qū)域）Fig.11 11 sub-domains in the overall domain（different colors represents different sub-domains）

并行程序性能分析：

（1）表1表明，本文開發(fā)的并行程序具有較高的計(jì)算性能，運(yùn)行時(shí)間從串行計(jì)算的949.4h（39d13.4h），到191個(gè)進(jìn)程的18.9h，加速比50.23，為高壩抗震研究提供了強(qiáng)有力的計(jì)算手段。

表1 并行計(jì)算性能測試值Tab. 1 Parallel computation test value

（2）圖12和圖13說明，隨著進(jìn)程個(gè)數(shù)的增加，運(yùn)行時(shí)間減少，加速比增大，表明本課題的并行程序具備良好的可擴(kuò)展性。

圖12 進(jìn)程個(gè)數(shù)與運(yùn)行時(shí)間的關(guān)系Fig.12 The relationship between process number and running time

圖13 并行計(jì)算加速比Fig.13 Parallel computation speed-up

（3）圖14說明，隨著進(jìn)程個(gè)數(shù)的增加，雖然加速比在增加，并行效率卻在減小，這也符合并行程序的一般規(guī)律。隨著進(jìn)程個(gè)數(shù)的增加，各個(gè)進(jìn)程計(jì)算量在減小，進(jìn)程之間數(shù)據(jù)通信的比重在增加，雖然并行的效益在提高，并行效率卻在下降。當(dāng)進(jìn)程之間通信的消耗超過各個(gè)進(jìn)程并行計(jì)算帶來的效益時(shí)，并行加速比不會再提高。

圖14 并行計(jì)算效率Fig.14 Parallel computation efficiency

3.2 西部強(qiáng)震區(qū)某特高拱壩整體系統(tǒng)抗震分析

西部強(qiáng)震區(qū)某特高拱壩壩高289.0m，屬300m級的超高拱壩，整個(gè)分析系統(tǒng)全部由三維塊體單元離散，節(jié)點(diǎn)總數(shù)約6.9萬，自由度總數(shù)約20.5萬。大壩壩體-壩肩-地基系統(tǒng)有限元網(wǎng)格見圖15。圖16所示為壩體-壩肩巖體模型。

圖15 大壩壩體-壩肩-地基系統(tǒng)有限元網(wǎng)格Fig.15 Dam-abutment rock-foundation model

圖16 壩體-壩肩巖體模型Fig.16 Dam-abutment rock mode

圖17 設(shè)計(jì)地震靜動綜合最大主應(yīng)力極值（MPa）Fig.17 The maximum principle stress distribution of dam（MPa）

圖18 設(shè)計(jì)地震頂拱橫縫張開度Fig.18 The maximum joint opening value of dam crest upstream and downstream under design earthquake

圖17表示設(shè)計(jì)地震作用下大壩上下游面最大主應(yīng)力分布，圖18表示設(shè)計(jì)地震作用下頂拱橫縫張開度，圖19表示地震結(jié)束后橫縫處于閉合狀態(tài)，圖20表示設(shè)計(jì)地震作用下左岸滑塊特征點(diǎn)錯(cuò)動時(shí)程，圖21表示設(shè)計(jì)地震作用下右岸滑塊特征點(diǎn)錯(cuò)動時(shí)程。

圖19 設(shè)計(jì)地震作用下3縫壩頂上游縫張開度時(shí)程Fig.19 Opening time history of dam crest upstream of joint No. 3 under design earthquake

圖20 設(shè)計(jì)地震作用左岸底滑面特征點(diǎn)順河向錯(cuò)動時(shí)程Fig.20 Sliding distance time history of bottom slip surface of left bank abutment rock block under design earthquake

圖21 設(shè)計(jì)地震作用右岸底滑面特征點(diǎn)順河向錯(cuò)動時(shí)程Fig.21 Sliding distance time history of bottom slip surface of right bank abutment rock block under design earthquake

3.3 國外某重力壩全壩段抗震分析

建立壩體-地基體系三維有限元模型。模型模擬了全部壩段及壩段間橫縫、壩基交界面接縫以及壩體上、下游折坡處對應(yīng)的475m高程層面、485m高程層面的水平層間縫，模擬了沿橫河向、順河向及豎向壩底高程以下2倍最大壩高的地基范圍。有限元模型節(jié)點(diǎn)總數(shù)603230，單元總數(shù)558201，總自由度數(shù)約1800000個(gè)。圖22為壩體-地基體系有限元網(wǎng)格，圖23所示為大壩震后變形圖，圖24所示為地震作用下各壩段沿壩基交界面順河向滑移時(shí)程。

圖22 壩體-地基體系有限元網(wǎng)格Fig.22 Dam-foundation model

圖23 大壩震后變形圖Fig.23 Deformation diagram of dam after earthquake

圖24 各壩段沿壩基交界面順河向滑移時(shí)程Fig.24 Time history of the sliding displacement along the dam-foundation interface of each monolith

4 結(jié)束語

高性能并行計(jì)算在高壩抗震分析中正發(fā)揮著越來越重要的作用。限于篇幅和時(shí)間限制，本文僅列出了中國水利水電科學(xué)研究院自主開發(fā)的高壩地震動響應(yīng)并行計(jì)算軟件（PSDAP）部分應(yīng)用成果，同時(shí)也在積極與高性能計(jì)算領(lǐng)域的專家團(tuán)隊(duì)展開合作，進(jìn)一步提高大規(guī)模計(jì)算性能。

[1] 金先龍，李淵印. 結(jié)構(gòu)動力學(xué)并行計(jì)算方法及應(yīng)用. 北京：國防工業(yè)出版社，2008.JIN Xianlong，LI Yuanyin. Parallel computing method of structural dynamics and its application. Beijing ：National Defend Industry Press，2008.

[2] Schwarz，H. A. Gesammelete mathematische abhandlungen.Volume 2，Springer，Berlin，1890. First published in Vierteljahrsschrift der Naturforschenden Gesellschaft in Zurich，volume 15，1870，272-186.

[3] 呂濤，石濟(jì)民，林振寶. 區(qū)域分解算法. 北京：科學(xué)出版社，1992.LV Tao，Shi Jimin，LIN Zhenbao，Decomposition algorithm.Beijing ：Science Press，1992.

[4] LYSMER J，KUHLEMEYER R L. Finite dynamic model for infinite media[J]. Journal of the Engineering Mechanics Division（ASCE），1969，95 ：759-877.

[5] DEEKS A J，RANDOLPH M F. Axisymmetric time-domain transmitting boundary[J]. Journal of Engineering Mechanics，1994，120（1）：25-42.

[6] LIU Jingbo，LU Yandong. A direct method for analysis of dynamic soil-structure interaction based on interface idea[M]//Dynamic Soil-structure Interaction. Beijing ：International Academic Publishers，1997，258-273.