基于FrontISTR的大規(guī)模地震數(shù)值模擬軟件的研發(fā)

董 勤 喜

(教育科學(xué)研究所，日本 川崎市 214-0033)

工程結(jié)構(gòu)的抗震分析通常采用地震波作為地震荷載輸入，如中國1975年的海城波、美國1940年的E1-Centro波及1952年的Taft波等。這些地震波與震源斷層的基巖、覆蓋土層等當(dāng)?shù)氐牡刭|(zhì)情況密切相關(guān)，事實上不具備通用性。地震是斷層的錯動所致，因此，可以直接利用斷層模型計算地震動，再利用波動理論計算某一地點的地震動。以此為基礎(chǔ)，計算得到的地震動波形將更有應(yīng)用價值。這是近年來許多學(xué)者正不斷努力的方向。Nishimura等[1]針對2011年日本東北部海域發(fā)生的里氏9.0級大地震出現(xiàn)的同震變形，采用基于開源并行計算環(huán)境FrontISTR[2]的EduS/FrontISTR軟件，直接模擬地殼內(nèi)的大逆沖型震源斷層，成功再現(xiàn)了日本新瀉—神戶間軟弱沉積平原區(qū)域內(nèi)的應(yīng)變集中現(xiàn)象。全球衛(wèi)星觀測系統(tǒng)和合成孔徑雷達(dá)干涉測量技術(shù)InSAR的綜合應(yīng)用，使得觀測同震地震變形的詳細(xì)分布成為可能。這為大規(guī)模地震分析研究提供了良好的實測數(shù)據(jù)，使得大規(guī)模地震數(shù)值模擬分析成為可能。嚴(yán)珍珍等[3]基于橫向各向同性PREM地球模型，利用CMT(centroid moment tensor)震源機(jī)制，考慮地表地形及海洋等地球特性，采用Jeroen Tromp等提供的SPECFEM3D_GOLBE譜元法軟件，進(jìn)行了大規(guī)模地震數(shù)值分析，較準(zhǔn)確地重現(xiàn)了長周期理論頻率值。許才軍等[4]采用彈性半空間矩形位錯模型，結(jié)合GPS和InSAR觀測結(jié)果對日本東北部海域發(fā)生的里氏9.0級大地震引起的同震滑動分布進(jìn)行了反演分析。

目前，可以進(jìn)行大規(guī)模地震數(shù)值模擬分析的常用軟件[5]見表1所示。

表1 常用大規(guī)模地震波數(shù)值模擬軟件

本文將詳細(xì)介紹基于并行計算有限元法和斷層模型研制的EduS/FrontISTR軟件及其在地震數(shù)值模擬方面的基本原理，并通過標(biāo)準(zhǔn)地震模型對軟件的性能進(jìn)行數(shù)值模擬驗證。最后，以2011年日本東北地區(qū)太平洋近海地震為例，應(yīng)用EduS/FrontISTR軟件對地震發(fā)生后的地震波的傳播過程等進(jìn)行大規(guī)模地震數(shù)值模擬研究，很好地再現(xiàn)了觀測結(jié)果。

1 并行計算有限元地震數(shù)值模擬軟件

筆者基于開源并行有限元程序開發(fā)平臺FrontISTR[2]，研究開發(fā)了具有自主知識產(chǎn)權(quán)的EduS/FrontISTR軟件。該軟件是一款通用的高性能并行非線性有限元計算軟件，適用于個人PC機(jī)到超級計算機(jī)系統(tǒng)。該軟件可以求解常規(guī)的靜態(tài)和動態(tài)線性問題，也可以求解幾何非線性、材料非線性和邊界非線性(含接觸)等問題，還可以有效模擬從微小滑動到有限滑動問題、從無摩擦接觸到有限摩擦接觸。EduS/FrontISTR軟件內(nèi)的子模塊EduS/FrontISTR/Quake可以進(jìn)行大規(guī)模的地震數(shù)值模擬計算，包含以下功能模塊：1)斷層模塊和震源模塊；2)邊界能量耗散吸收模塊；3)考慮地表地形和海底形狀的快速網(wǎng)格自動生成及地震波后處理模塊；4)通過大數(shù)據(jù)結(jié)構(gòu)，嵌入非勻質(zhì)地下結(jié)構(gòu)J-SHIS[13]和松原地質(zhì)結(jié)構(gòu)模型[14]等。在EduS/FrontISTR高性能并行計算支撐平臺下，運行子程序EduS/FrontISTR/Quake，可以快速地實現(xiàn)地震波的傳播模擬。

1.1 并行有限元法

大規(guī)模地震數(shù)值模擬的計算區(qū)域很大，需要將計算區(qū)域分割成若干子區(qū)域。在EduS/FrontISTR軟件中，采用節(jié)點法對計算區(qū)域進(jìn)行分割。以圖1(a)所示的計算區(qū)域為例，將計算區(qū)域分割為4個子區(qū)域，然后，將各子區(qū)域重新按有限元法進(jìn)行整理，如圖1(b)所示。

計算區(qū)域分割時，圖1(a)中各子區(qū)域內(nèi)部的節(jié)點稱為內(nèi)點，圖1(b)中各子區(qū)域重新整理時各單元引入的其他子區(qū)域的節(jié)點稱為外點，定義一個子區(qū)域的內(nèi)點是其他子區(qū)域的外點的節(jié)點為邊界點。節(jié)點編號原則是先內(nèi)點，后外點。以圖1(b)中子域0為例，節(jié)點1到7為內(nèi)點，節(jié)點8到15為外點，節(jié)點1、2、4和6為子域1的邊界點，節(jié)點1到3為子域2的邊界點，節(jié)點1為子域3的邊界點。在軟件中，采用通信表建立邊界和外點的關(guān)系。各子區(qū)域構(gòu)建完成后，采用分布式數(shù)據(jù)存儲，再按常規(guī)的有限元法對各子區(qū)域獨立進(jìn)行有限元計算，各子域間僅利用通信表進(jìn)行數(shù)據(jù)信息交換。

圖1 并行有限元法的計算區(qū)域分割

1.2 地震斷層模擬

為模擬地震斷層，子程序EduS/FrontISTR/Quake采用分裂節(jié)點法[15](Split-node)處理正斷層和逆斷層處的有摩擦平面剪切裂紋，如圖2所示。通過分裂節(jié)點法，可以對斷層的上面和下面指定不同的滑動位移，但兩者的相對位移等于斷層的位錯量。

圖2 分裂節(jié)點法

1.3 震源模擬

在子程序EduS/FrontISTR/Quake中采用點源法將地震位錯表示為一個集中力，為雙力偶源。斷層參數(shù)的定義如圖3所示。

斷層斷裂時，數(shù)值模擬計算采用的矩張量[16]的表達(dá)式為：

式中：M0為地震的標(biāo)量地震矩；Mij為矩張量的分量；其余參數(shù)的定義見圖3。

圖3 斷層參數(shù)的定義

通過單元的體積積分得到節(jié)點的集中力，即雙力偶源，為

(2)

式中：Pi表示節(jié)點力；N表示位移形狀函數(shù)；xj表示坐標(biāo)方向。

1.4 網(wǎng)格模型生成

利用非均質(zhì)的地下結(jié)構(gòu)模型J-SHIS[13]，得到圖4所示的地層結(jié)構(gòu)形狀，再采用EduS/FrontISTR軟件，可以自動生成計算區(qū)域的有限元網(wǎng)格，最后完成如圖5所示的并行有限元計算的自動建模。

圖4 地層形狀(富士山周邊)

圖5 EduS/FrontISTR軟件生成的網(wǎng)格模型

在EduS/FrontISTR軟件中，可以采用黏性吸收邊界[17]或無限元處理計算區(qū)域的邊界。

2 數(shù)值計算結(jié)果的理論驗證

筆者在軟件研制過程中進(jìn)行了大量數(shù)值模擬驗證?，F(xiàn)以地震數(shù)值模擬中典型模型為例，將研制的三維數(shù)值模擬計算結(jié)果與理論值進(jìn)行對比驗證。

2.1 地殼形變分析的理論驗證

數(shù)值計算區(qū)域取為長200 km，寬200 km，深100 km，地質(zhì)材料為勻質(zhì)地層。計算區(qū)域采用六面體單元，網(wǎng)格單元數(shù)約為102.4萬個，網(wǎng)格節(jié)點數(shù)約為106.3萬個。地震斷層的理論模型如圖6所示，相應(yīng)的計算參數(shù)見表2。使用的程序為EduS/FrontISTR/Quake，CPU為Quad Core Intel (R) Xeon X3430 2.40GHz, 采用1個Core進(jìn)行地殼變動數(shù)值模擬的時間約為3 min，用4個Core進(jìn)行數(shù)值模擬的時間約為1.5 min。計算得到的地面位移與岡田模型[18]的理論位移的對比如圖7所示?？梢钥闯?，數(shù)值模擬的計算結(jié)果與理論值吻合得非常好。

圖6 地殼變形的斷層理論模型

斷層長度L/km斷層寬度W/km傾斜角δ/(°)滑動角λ/(°)走向角θ/(°)埋深d/km位錯量D0/m201030909027．55

(a)地面東西向分量位移 (b)地面垂向分量位移

2.2 地震波傳播的理論驗證

地震波傳播的理論模型如圖8所示，計算區(qū)域取為長100 km，寬100 km，深50 km，斷層為水平橫向斷層，地質(zhì)材料常數(shù)為Vs=3 km/s,Vp=5 km/s，ρ=2 500 kg/m3的勻質(zhì)地層，斷層用10×10個小斷層進(jìn)行數(shù)值模擬。斷層破壞從震中(0,0,-10 km)開始，以傳播速度vr=3 km/s同心圓的形式傳播，斷層滑動的速度時間函數(shù)是上升時間為3 s的三角函數(shù)。地質(zhì)材料的阻尼假設(shè)為無阻尼和質(zhì)量阻尼，相應(yīng)的計算參數(shù)見表3。

圖8 地震波傳播的斷層理論模型(1/4模型)

斷層長度L/km斷層寬度W/km傾斜角δ/(°)滑動角λ/(°)走向角θ/(°)埋深d/km位錯量D0/m10100090102

計算區(qū)域采用六面體單元，1/4模型的網(wǎng)格單元數(shù)約為51萬個，網(wǎng)格節(jié)點數(shù)約為53萬個。模型的東西向采用無限元消除邊界波反射的影響。與理論值對比的位置(觀測點)取為東西向A(-2, 0, 0)，B(-1, 0, 0)，C(0, 0, 0)，D(1, 0, 0)和E(2, 0, 0)等5點(單位為km)。

使用的程序為EduS/FrontISTR/Quake，CPU為Quad Core Intel (R) Xeon X3430 2.40GHz,模擬地震波傳播時間為30 s，時間步長為5 ms，采用1個Core進(jìn)行地震波傳播模擬的計算時間約為20 min，用4個Core的計算時間約為10 min。

將5個觀測點的理論值與數(shù)值計算結(jié)果進(jìn)行對比，以驗證數(shù)值模擬結(jié)果的精確性，如圖9所示。數(shù)值計算的速度波形、峰值和位相同理論解非常一致。數(shù)值模擬的速度波形，在11 s之后出現(xiàn)振蕩現(xiàn)象，這是由于計算區(qū)域范圍比較狹窄，出現(xiàn)邊界人工波的反射干涉所致。

(a)無阻尼 (b)質(zhì)量阻尼

對比無阻尼和質(zhì)量阻尼的速度波形，可以看出兩者的波形相同，質(zhì)量阻尼波形的位相稍微滯后，峰值是無阻尼波形的80%左右。綜合以上的結(jié)果，開發(fā)的地殼變形和地震波傳播的數(shù)值模擬程序的計算結(jié)果是精確和可靠的，同時也驗證了阻尼處理方法的正確性。

3 日本東北地區(qū)太平洋近海地震數(shù)值模擬

3.1 再現(xiàn)地殼變形

數(shù)值模擬計算范圍取為北緯 35°到42°，東經(jīng)138°到145°，深度80 km。再現(xiàn)地殼變形使用4片斷層構(gòu)成的計算模型[19]，采用勻質(zhì)模型計算得到的水平位移，與岡田模型及觀測值的對比如圖10所示。由圖知，勻質(zhì)模型計算結(jié)果與岡田模型及觀測結(jié)果有較好的一致性。進(jìn)一步，基于非均質(zhì)的地下結(jié)構(gòu)模型J-SHIS[13]和松原地震波速度模型[14]相結(jié)合，建立了水平方向1 km分辨率的地殼變形和地震波傳播的大型數(shù)據(jù)庫。采用非均質(zhì)的地下結(jié)構(gòu)模型計算得到東西向的應(yīng)變，如圖11所示。非均質(zhì)地下結(jié)構(gòu)模型的數(shù)值模擬結(jié)果很好地揭示了日本東北地區(qū)太平洋近海地震的斷層運動引起越后平原的東西向伸縮，再現(xiàn)了日本新瀉—神戶地區(qū)間的應(yīng)變集中帶，與觀測結(jié)果一致。這是國際上首次使用有限元和非均質(zhì)的地下結(jié)構(gòu)模型再現(xiàn)了應(yīng)變集中帶。

圖10 數(shù)值計算的地面水平位移與岡田模型及觀測值的對比

圖11 數(shù)值計算得到的越后平原東西向應(yīng)變分布[1]

3.2 再現(xiàn)地震波傳播過程

數(shù)值模擬計算范圍和非均質(zhì)的地下結(jié)構(gòu)模型同3.1節(jié)相同，地震斷層位錯分布采用日本氣象廳提供的原始數(shù)據(jù)[20]，該斷層的位錯分布是利用震源附近觀測點的時序擬合曲線得到，最大位錯為37 m，斷層的范圍為475 km×175 km，由走行方向19個，傾斜方向7個，共133個小斷層構(gòu)成，如圖12所示。

圖12 斷層的位錯分布

計算區(qū)域采用六面體和四面體單元，深度方向根據(jù)非均質(zhì)地下結(jié)構(gòu)地震波速結(jié)構(gòu)，通過四面體過渡單元(遷移層)，網(wǎng)格尺寸逐漸增大，以便節(jié)省內(nèi)存和節(jié)約計算時間。整個模擬計算區(qū)域的網(wǎng)格節(jié)點數(shù)約為5 000萬個，約1.5億自由度。計算區(qū)域的網(wǎng)格如圖13所示(由于網(wǎng)格太多，未能詳細(xì)顯示)。網(wǎng)格分割過程中，每個小斷層再進(jìn)一步分割為400個更小的斷層，總共使用5萬3 200個斷層來模擬斷層滑動，斷層滑動的速度時間函數(shù)采用8 s上升的三角函數(shù)來模擬。

圖13 數(shù)值模擬計算區(qū)域網(wǎng)格示意圖

模擬地震波傳播時間為300 s，時間步長為25 ms，時間步為12 000，采用24個Core進(jìn)行地震波傳播模擬的計算時間約為3 h。圖14所示的計算結(jié)果清晰地顯示了地震發(fā)生后，地震波向日本全土傳播的過程。

圖14 地震波傳播過程

由圖14可以看出，數(shù)值模擬在150 s時刻地震波的峰值剛好到達(dá)福島核發(fā)電廠附近。實際上，地震發(fā)生后，第2個波峰到達(dá)日本福島核發(fā)電廠的時間約為150 s，正是該波峰引起海底地殼的變動，產(chǎn)生的海嘯將核發(fā)電廠摧毀。數(shù)值模擬與實際發(fā)生海嘯的時間吻合的很好。這種利用地震斷層位錯模型計算得到的地震波更接近實際情況，有利于工程結(jié)構(gòu)的災(zāi)害評估。

4 結(jié)論

筆者基于開源并行有限元計算環(huán)境FrontISTR，研究開發(fā)了高性能的地殼變形和地震波傳播并行計算軟件EduS/FrontISTR，同時構(gòu)建了非均質(zhì)的地下結(jié)構(gòu)數(shù)據(jù)庫。該軟件可以高效地對大規(guī)模地震數(shù)值計算進(jìn)行建模和并行計算，數(shù)值模擬結(jié)果很好地再現(xiàn)了地震后的地表觀測結(jié)果或理論結(jié)果，能夠很好地再現(xiàn)地震發(fā)生后地震波的持續(xù)傳播過程?；贓duS/FrontISTR軟件計算得到的地震波將更有實際應(yīng)用價值，將有助于地震后工程結(jié)構(gòu)的災(zāi)害評估。

日本京都大學(xué)西村卓也教授(原日本國土地理院主任研究員)為本文繪制了圖10，加工了圖11，日本計算力學(xué)研究中心吉見顯一朗為本文繪制了地震波傳播圖片等，在此一并致謝。

[1]NISHIMURA T, SUITO H, KOBAYASHI T, et al. Excess strain in the Echigo Plain sedimentary basin, NE Japan: evidence from coseismic deformation of the 2011 Tohoku-Oki earthquake[J].Geophysical Journal Inertnation, 2016, 205(3): 1613.

[2]日本計算工學(xué)會GreenCAE研究會共催ユーザ會資料[EB/OL].http://www.multi.k.u-tokyo.ac.jp/FrontISTR/.

[3]嚴(yán)珍珍，張懷，范湘濤，等. 日本Mw9.0 大地震激發(fā)全球地震波傳播特性的數(shù)值模擬研究與對比分析[J]. 中國科學(xué)：地球科學(xué)，2014，44(2)：259.

[4] 許才軍，何平，溫?fù)P茂，等. 日本2011 Tohoku-Oki Mw 9.0級地震的同震形變及其滑動分布反演：GPS和InSAR約束[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版)，2012，37(12)：1387.

[5]張林波，鄭偉英，盧本卓，等. 并行自適應(yīng)有限元軟件平臺PHG及其應(yīng)用[J]. 中國科學(xué)：信息科學(xué)，2016，46(10)：1442.

[6]董勤喜，柴山恭.基于開源并行計算平臺FrontISTR的非均質(zhì)地下結(jié)構(gòu)的地殼變形和地震波傳播的數(shù)值模擬[C]//日本文部科學(xué)省次世代革新數(shù)值模擬軟件的研究開發(fā)最終成果報告會.東京都:東京大學(xué)生產(chǎn)技術(shù)研究所，2013.

[7]KOMATITSCH D, TSUBOI S, CHEN J, et al. A 14.6 billion degree of freedom, 5teraflops, 2.5 terabyte earthquake simulation on the Earth simulator[C]//Proceedings of the 2003 ACM/IEEE International Conference on Supercomputing (SC03). New York: ACM, 2003:4-11.

[8]CUI Y, POYRAZ E, OLSEN K B, et al. Physics-based seismic hazard analysis on petascale heterogenous supercomputers[C]// Proceeding of 2013 ACM/IEEE International Conference on High Performance Computing, Networking, Storage and Analysis (SC13). New York: ACM, 2013:1-12.

[9]HEINECKE A, BREUER A, RETTENBERGER S, et al. Petascale high order dynamic rupture earthquake simulations on heterogenous supercomputers[C]//Proceeding of the ACM/IEEE International Conference on High Performance Computing, Networking, Storage and Analysis (SC14).Piscataway: IEEE Press, 2014:3-14.

[10]ICHIMURA T, FUJITA K, QUINAY P E B, et al.Implicit nonlinear wave simulation with 1.08T DOF and 0.270T unstructured finite elements to enhance comprehensive earthquake simulation[C]// Proceeding of the ACM/IEEE International Conference on High Performance Computing, Networking, Storage and Analysis (SC15). New York: ACM, 2015:1-12.

[11] MAZZIERI I, STUPAZZINI M, GUIDOTTI R, et al. SPEED: spectral elements in elastodynamics with discontinuous Galerkin: a non-conforming approach for 3D multi-scale problems[J].Int J for Numerical Meth in Eng, 2013, 95(12): 991.

[12]BAO H, BIELAK J, GHATTAS O, et al. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers[J].Computer Methods in Applied Mechanics & Engineering, 1998, 152:85.

[13]J-SHIS 地震ハザードステーション[EB/OL].http://www.j-shis.bosai.go.jp/.

[14]MATSUBARA M, OBARA K,KASAHARA K. Three-dimensional P-and S-wave velocity structures beneath the Japan islands obtained by high-density seismic stations by seismic tomography[J].Tectonophysics, 2008, 454(1): 86.

[15]MELOSH H J, RAEFSKY A. A simple and efficient method or introducing faults into finite element computations[J]. Bulletin of the Seismological Society of America, 1981, 71： 1391.

[16] AKI K, RICHARDS P G. Quantitative seismology theory and methods[M].New York : W H Freeman and Company, 1980:61.

[17] LYSMER J, KUHLEMEYER R L.Finite dynamic model for infinite media[J]. ASCE Journal of Engineering Mechanics Division, 1969, 95： 859.

[18]OKADA Y. Internal deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 1992, 82： 1018.

[19]NISHIMURA T, MUNEKANE H, YARAI H. The 2011 off the Pacific coast of Tohoku earthquake and its aftershocks observed by GEONET[J]. Earth Planets Space, 2011, 63(7):631.

[20]平成23年(2011年)東北地方太平洋沖地震調(diào)查報告(日文)，第I編:日本氣象廳技術(shù)報告:第133號[R].東京都:日本氣象廳,2012.

特約專家介紹

董勤喜(1962—)，男，漢族，山東蘭陵人，中共黨員，科學(xué)技術(shù)博士。1997年博士畢業(yè)于奧地利因斯布魯克大學(xué)土木與建筑工程學(xué)院。1998年8月至2000年7月日本學(xué)術(shù)振興會(JSPS)博士后。主要從事隧道開挖數(shù)值模擬，材料參數(shù)反演，道路及機(jī)場跑道設(shè)計修復(fù)軟件的研發(fā)，流固耦合并行計算的研發(fā)，高性能并行計算及大規(guī)模地震數(shù)值模擬的研究，日本下一代超級計算機(jī)京[K-computer；E級]上先進(jìn)制造及熱塑性復(fù)合材料數(shù)值模擬軟件等的研究。已在《應(yīng)用數(shù)學(xué)和力學(xué)》《日本土木學(xué)會論文集》《日本道路工程學(xué)會論文集》《J. Beijing Inst. of Techn.》《Appl. Math. Mech.》《Nonlinear Analysis-Theory, Methods & Applications》《Phys. Rev. E》《Phys. Lett. A》《J. Trauma-Injury Infection and Critical Care》《Int. J. Geomech., ASCE》《Int. J. Numer. Anal. Meth. Geomech》《Finite Elements in Analysis and Design》《Geophys. J. Int.》及《Int. J. Pavement》等重要學(xué)術(shù)刊物上發(fā)表論文60余篇。主持和參與日本鐵道運輸基礎(chǔ)研究基金、JST自然科學(xué)基金、日本文部科學(xué)省火山研究基金、日本文部科學(xué)省面向E級(百億億次級；旗艦2020計劃)超級計算機(jī)上復(fù)合材料數(shù)值模擬研究基金資助的研究課題。