基于FrontISTR的大规模地震数值模拟软件的研发

董 勤 喜

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

工程结构的抗震分析通常采用地震波作为地震荷载输入，如中国1975年的海城波、美国1940年的E1-Centro波及1952年的Taft波等。这些地震波与震源断层的基岩、覆盖土层等当地的地质情况密切相关，事实上不具备通用性。地震是断层的错动所致，因此，可以直接利用断层模型计算地震动，再利用波动理论计算某一地点的地震动。以此为基础，计算得到的地震动波形将更有应用价值。这是近年来许多学者正不断努力的方向。Nishimura等[1]针对2011年日本东北部海域发生的里氏9.0级大地震出现的同震变形，采用基于开源并行计算环境FrontISTR[2]的EduS/FrontISTR软件，直接模拟地壳内的大逆冲型震源断层，成功再现了日本新泻—神户间软弱沉积平原区域内的应变集中现象。全球卫星观测系统和合成孔径雷达干涉测量技术InSAR的综合应用，使得观测同震地震变形的详细分布成为可能。这为大规模地震分析研究提供了良好的实测数据，使得大规模地震数值模拟分析成为可能。严珍珍等[3]基于横向各向同性PREM地球模型，利用CMT(centroid moment tensor)震源机制，考虑地表地形及海洋等地球特性，采用Jeroen Tromp等提供的SPECFEM3D_GOLBE谱元法软件，进行了大规模地震数值分析，较准确地重现了长周期理论频率值。许才军等[4]采用弹性半空间矩形位错模型，结合GPS和InSAR观测结果对日本东北部海域发生的里氏9.0级大地震引起的同震滑动分布进行了反演分析。

目前，可以进行大规模地震数值模拟分析的常用软件[5]见表1所示。

表1 常用大规模地震波数值模拟软件

本文将详细介绍基于并行计算有限元法和断层模型研制的EduS/FrontISTR软件及其在地震数值模拟方面的基本原理，并通过标准地震模型对软件的性能进行数值模拟验证。最后，以2011年日本东北地区太平洋近海地震为例，应用EduS/FrontISTR软件对地震发生后的地震波的传播过程等进行大规模地震数值模拟研究，很好地再现了观测结果。

1 并行计算有限元地震数值模拟软件

笔者基于开源并行有限元程序开发平台FrontISTR[2]，研究开发了具有自主知识产权的EduS/FrontISTR软件。该软件是一款通用的高性能并行非线性有限元计算软件，适用于个人PC机到超级计算机系统。该软件可以求解常规的静态和动态线性问题，也可以求解几何非线性、材料非线性和边界非线性(含接触)等问题，还可以有效模拟从微小滑动到有限滑动问题、从无摩擦接触到有限摩擦接触。EduS/FrontISTR软件内的子模块EduS/FrontISTR/Quake可以进行大规模的地震数值模拟计算，包含以下功能模块：1)断层模块和震源模块；2)边界能量耗散吸收模块；3)考虑地表地形和海底形状的快速网格自动生成及地震波后处理模块；4)通过大数据结构，嵌入非匀质地下结构J-SHIS[13]和松原地质结构模型[14]等。在EduS/FrontISTR高性能并行计算支撑平台下，运行子程序EduS/FrontISTR/Quake，可以快速地实现地震波的传播模拟。

1.1 并行有限元法

大规模地震数值模拟的计算区域很大，需要将计算区域分割成若干子区域。在EduS/FrontISTR软件中，采用节点法对计算区域进行分割。以图1(a)所示的计算区域为例，将计算区域分割为4个子区域，然后，将各子区域重新按有限元法进行整理，如图1(b)所示。

计算区域分割时，图1(a)中各子区域内部的节点称为内点，图1(b)中各子区域重新整理时各单元引入的其他子区域的节点称为外点，定义一个子区域的内点是其他子区域的外点的节点为边界点。节点编号原则是先内点，后外点。以图1(b)中子域0为例，节点1到7为内点，节点8到15为外点，节点1、2、4和6为子域1的边界点，节点1到3为子域2的边界点，节点1为子域3的边界点。在软件中，采用通信表建立边界和外点的关系。各子区域构建完成后，采用分布式数据存储，再按常规的有限元法对各子区域独立进行有限元计算，各子域间仅利用通信表进行数据信息交换。

图1 并行有限元法的计算区域分割

1.2 地震断层模拟

为模拟地震断层，子程序EduS/FrontISTR/Quake采用分裂节点法[15](Split-node)处理正断层和逆断层处的有摩擦平面剪切裂纹，如图2所示。通过分裂节点法，可以对断层的上面和下面指定不同的滑动位移，但两者的相对位移等于断层的位错量。

图2 分裂节点法

1.3 震源模拟

在子程序EduS/FrontISTR/Quake中采用点源法将地震位错表示为一个集中力，为双力偶源。断层参数的定义如图3所示。

断层断裂时，数值模拟计算采用的矩张量[16]的表达式为：

式中：M0为地震的标量地震矩；Mij为矩张量的分量；其余参数的定义见图3。

图3 断层参数的定义

通过单元的体积积分得到节点的集中力，即双力偶源，为

(2)

式中：Pi表示节点力；N表示位移形状函数；xj表示坐标方向。

1.4 网格模型生成

利用非均质的地下结构模型J-SHIS[13]，得到图4所示的地层结构形状，再采用EduS/FrontISTR软件，可以自动生成计算区域的有限元网格，最后完成如图5所示的并行有限元计算的自动建模。

图4 地层形状(富士山周边)

图5 EduS/FrontISTR软件生成的网格模型

在EduS/FrontISTR软件中，可以采用黏性吸收边界[17]或无限元处理计算区域的边界。

2 数值计算结果的理论验证

笔者在软件研制过程中进行了大量数值模拟验证。现以地震数值模拟中典型模型为例，将研制的三维数值模拟计算结果与理论值进行对比验证。

2.1 地壳形变分析的理论验证

数值计算区域取为长200 km，宽200 km，深100 km，地质材料为匀质地层。计算区域采用六面体单元，网格单元数约为102.4万个，网格节点数约为106.3万个。地震断层的理论模型如图6所示，相应的计算参数见表2。使用的程序为EduS/FrontISTR/Quake，CPU为Quad Core Intel (R) Xeon X3430 2.40GHz, 采用1个Core进行地壳变动数值模拟的时间约为3 min，用4个Core进行数值模拟的时间约为1.5 min。计算得到的地面位移与冈田模型[18]的理论位移的对比如图7所示。可以看出，数值模拟的计算结果与理论值吻合得非常好。

图6 地壳变形的断层理论模型

断层长度L/km断层宽度W/km倾斜角δ/(°)滑动角λ/(°)走向角θ/(°)埋深d/km位错量D0/m201030909027．55

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

2.2 地震波传播的理论验证

地震波传播的理论模型如图8所示，计算区域取为长100 km，宽100 km，深50 km，断层为水平横向断层，地质材料常数为Vs=3 km/s,Vp=5 km/s，ρ=2 500 kg/m3的匀质地层，断层用10×10个小断层进行数值模拟。断层破坏从震中(0,0,-10 km)开始，以传播速度vr=3 km/s同心圆的形式传播，断层滑动的速度时间函数是上升时间为3 s的三角函数。地质材料的阻尼假设为无阻尼和质量阻尼，相应的计算参数见表3。

图8 地震波传播的断层理论模型(1/4模型)

断层长度L/km断层宽度W/km倾斜角δ/(°)滑动角λ/(°)走向角θ/(°)埋深d/km位错量D0/m10100090102

计算区域采用六面体单元，1/4模型的网格单元数约为51万个，网格节点数约为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进行地震波传播模拟的计算时间约为20 min，用4个Core的计算时间约为10 min。

将5个观测点的理论值与数值计算结果进行对比，以验证数值模拟结果的精确性，如图9所示。数值计算的速度波形、峰值和位相同理论解非常一致。数值模拟的速度波形，在11 s之后出现振荡现象，这是由于计算区域范围比较狭窄，出现边界人工波的反射干涉所致。

(a)无阻尼 (b)质量阻尼

对比无阻尼和质量阻尼的速度波形，可以看出两者的波形相同，质量阻尼波形的位相稍微滞后，峰值是无阻尼波形的80%左右。综合以上的结果，开发的地壳变形和地震波传播的数值模拟程序的计算结果是精确和可靠的，同时也验证了阻尼处理方法的正确性。

3 日本东北地区太平洋近海地震数值模拟

3.1 再现地壳变形

数值模拟计算范围取为北纬 35°到42°，东经138°到145°，深度80 km。再现地壳变形使用4片断层构成的计算模型[19]，采用匀质模型计算得到的水平位移，与冈田模型及观测值的对比如图10所示。由图知，匀质模型计算结果与冈田模型及观测结果有较好的一致性。进一步，基于非均质的地下结构模型J-SHIS[13]和松原地震波速度模型[14]相结合，建立了水平方向1 km分辨率的地壳变形和地震波传播的大型数据库。采用非均质的地下结构模型计算得到东西向的应变，如图11所示。非均质地下结构模型的数值模拟结果很好地揭示了日本东北地区太平洋近海地震的断层运动引起越后平原的东西向伸缩，再现了日本新泻—神户地区间的应变集中带，与观测结果一致。这是国际上首次使用有限元和非均质的地下结构模型再现了应变集中带。

图10 数值计算的地面水平位移与冈田模型及观测值的对比

图11 数值计算得到的越后平原东西向应变分布[1]

3.2 再现地震波传播过程

数值模拟计算范围和非均质的地下结构模型同3.1节相同，地震断层位错分布采用日本气象厅提供的原始数据[20]，该断层的位错分布是利用震源附近观测点的时序拟合曲线得到，最大位错为37 m，断层的范围为475 km×175 km，由走行方向19个，倾斜方向7个，共133个小断层构成，如图12所示。

图12 断层的位错分布

计算区域采用六面体和四面体单元，深度方向根据非均质地下结构地震波速结构，通过四面体过渡单元(迁移层)，网格尺寸逐渐增大，以便节省内存和节约计算时间。整个模拟计算区域的网格节点数约为5 000万个，约1.5亿自由度。计算区域的网格如图13所示(由于网格太多，未能详细显示)。网格分割过程中，每个小断层再进一步分割为400个更小的断层，总共使用5万3 200个断层来模拟断层滑动，断层滑动的速度时间函数采用8 s上升的三角函数来模拟。

图13 数值模拟计算区域网格示意图

模拟地震波传播时间为300 s，时间步长为25 ms，时间步为12 000，采用24个Core进行地震波传播模拟的计算时间约为3 h。图14所示的计算结果清晰地显示了地震发生后，地震波向日本全土传播的过程。

图14 地震波传播过程

由图14可以看出，数值模拟在150 s时刻地震波的峰值刚好到达福岛核发电厂附近。实际上，地震发生后，第2个波峰到达日本福岛核发电厂的时间约为150 s，正是该波峰引起海底地壳的变动，产生的海啸将核发电厂摧毁。数值模拟与实际发生海啸的时间吻合的很好。这种利用地震断层位错模型计算得到的地震波更接近实际情况，有利于工程结构的灾害评估。

4 结论

笔者基于开源并行有限元计算环境FrontISTR，研究开发了高性能的地壳变形和地震波传播并行计算软件EduS/FrontISTR，同时构建了非均质的地下结构数据库。该软件可以高效地对大规模地震数值计算进行建模和并行计算，数值模拟结果很好地再现了地震后的地表观测结果或理论结果，能够很好地再现地震发生后地震波的持续传播过程。基于EduS/FrontISTR软件计算得到的地震波将更有实际应用价值，将有助于地震后工程结构的灾害评估。

日本京都大学西村卓也教授(原日本国土地理院主任研究员)为本文绘制了图10，加工了图11，日本计算力学研究中心吉见显一朗为本文绘制了地震波传播图片等，在此一并致谢。

[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]日本計算工学会GreenCAE研究会共催ユーザ会資料[EB/OL].http://www.multi.k.u-tokyo.ac.jp/FrontISTR/.

[3]严珍珍，张怀，范湘涛，等. 日本Mw9.0 大地震激发全球地震波传播特性的数值模拟研究与对比分析[J]. 中国科学：地球科学，2014，44(2)：259.

[4] 许才军，何平，温扬茂，等. 日本2011 Tohoku-Oki Mw 9.0级地震的同震形变及其滑动分布反演：GPS和InSAR约束[J]. 武汉大学学报(信息科学版)，2012，37(12)：1387.

[5]张林波，郑伟英，卢本卓，等. 并行自适应有限元软件平台PHG及其应用[J]. 中国科学：信息科学，2016，46(10)：1442.

[6]董勤喜，柴山恭.基于开源并行计算平台FrontISTR的非均质地下结构的地壳变形和地震波传播的数值模拟[C]//日本文部科学省次世代革新数值模拟软件的研究开发最终成果报告会.东京都:东京大学生产技术研究所，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年)东北地方太平洋冲地震调查报告(日文)，第I编:日本气象厅技术报告:第133号[R].东京都:日本气象厅,2012.

特约专家介绍

董勤喜(1962—)，男，汉族，山东兰陵人，中共党员，科学技术博士。1997年博士毕业于奥地利因斯布鲁克大学土木与建筑工程学院。1998年8月至2000年7月日本学术振兴会(JSPS)博士后。主要从事隧道开挖数值模拟，材料参数反演，道路及机场跑道设计修复软件的研发，流固耦合并行计算的研发，高性能并行计算及大规模地震数值模拟的研究，日本下一代超级计算机京[K-computer；E级]上先进制造及热塑性复合材料数值模拟软件等的研究。已在《应用数学和力学》《日本土木学会论文集》《日本道路工程学会论文集》《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》等重要学术刊物上发表论文60余篇。主持和参与日本铁道运输基础研究基金、JST自然科学基金、日本文部科学省火山研究基金、日本文部科学省面向E级(百亿亿次级；旗舰2020计划)超级计算机上复合材料数值模拟研究基金资助的研究课题。