某城市的近断层强地面运动影响场显式有限元数值模拟:模型与方法

　　摘要：
　　本文利用某城市在《城市活断层探测与地震危险性评价》项目前期工作中得到的文字和数据资料，结合实际概要说明并给出了该城市在设定地震作用下近断层强地面运动影响场显式有限元数值模拟所涉及的具体计算模型和实用计算方法，为下一步预测该城市设定地震作用下的近断层强地面运动影响场完成了必要的前期准备工作。
　　
　　关键词：强地面运动影响场；显式有限元；设定地震； 数值模拟； 模型与方法
　　
　　引言
　　
　　合理确定近断层强地面运动影响场，为城市或地区的整体地震安全性评价、城市发展规划、国土利用、重大工程选址等提供科学参考依据，是国家十五重大工程项目《城市活断层探测与地震危险性评价》的一个主要产出目标。近断层强地面运动影响场显式有限元数值模拟方法是国内外公认的研究近断层强地面运动的最有效方法之一，同时也是实现谢礼立院士提出的建设数字地震减灾系统构想的有力工具〔1〕。我国现有20 多个城市的《城市活断层探测与地震危险性评价》项目正在进行之中，其中某些城市已经取得了初步的探测和分析结果。
　　
　　按照中国地震局震防司和徐锡伟首席科学家要求：中国地震局工程力学研究所“近断层强地面运动影响场的确定”项目组（以下简称项目组）的主要任务是：根据上述《城市活断层探测与地震危险性评价》的相关结果，采用显式有限元数值模拟方法研究并近似确定（预测）相关城市的近断层强地面运动影响场。项目组于2004 年底陆续收到中国地震局地质研究所（以下简称地质所）曲国胜研究员提供的有关城市《城市活断层探测与地震危险性评价》前期工作中得到的与建立有限断层震源模型和3 维地壳结构模型相关的部分文字和数据资料。
　　
　　本文利用收到的文字和数据资料、结合实际说明城市在设定地震作用下近断层强地面运动影响场显式有限元数值模拟涉及的模型和方法。
　　
　　1 基于设定地震和相关资料的有限断层运动学模型
　　
　　1.1 设定地震与全局震源参数
　　
　　根据地质所提供的数据资料和对设定地震的规定：某城市的初步近断层强地面运动影响场预测仅考虑第x 号断层（Fx）的影响，设定地震考虑为矩震级6.75 w M = ；断层的产状为NE20° / SED70 -80°，错动方式为正倾滑；断层长度为22千米，震源（初始破裂点）坐标为X = 20605691,Y = 2704476, Z =10千米。根据文献〔2〕，以上数据不能满足构建有限断层运动学模型对全局震源参数的要求，还需要补充断层沿下倾方向的宽度、上界埋深、平均滑移量和平均破裂传播速度，这些补充参数可以根据文献〔3〕给出的相关经验统计公式并结合作者的数值模拟计算经验做出选择。据以上说明，对于设定地震为6.75 w M = 之情况，本文取描述该断层的全局震源参数分别为：走向角20 s f = °，倾向角d = 75°，滑动角l = -90°（正倾滑），断层沿走向的长度L = 22千米，断层沿下倾方向的宽度W =18千米，。此外，取断层上下盘破裂面之间的平均位错D =1.1米，取断层平均破裂传播速度r V为相应介质剪切波速s V 的0.8倍，断层的上界埋深稍后讨论。离散化有限断层的平面图如图1 所示。
　　
　　1.2 局部震源参数与震源离散时空过程
　　
　　本文采用文献〔2〕提出的震源离散时空过程模型（断层的空间离散化在3.4小节讨论），这一模型可以表达为式（1）右端中：D(t)代表离散断层任一节点的标准化滑移时间函数； 0D (l,w)代表离散断层各节点的空间滑动分布函数； 3( ) g D t 代表离散断层任一节点的标准化振动时间函数；0 <a <1是待定常数。本文取a = 0.12，采用的具体震源时空函数示于图2。关于震源时空函数和其中包含的局部震源参数的确定方法可以参考文献〔2-5〕，限于篇幅本文不做展开叙述。需要指出：式（1）代表断层上盘破裂面的运动，断层下盘破裂面的运动与之大小相等、
　　
　　2 有限计算区模型和3 维地壳结构模型
　　
　　2.1 有限计算区模型
　　
　　地震断层是有限断层，由其产生的近断层强地面运动是一个以该断层为几何中心的3维有限区域，定性地说，适当外延这一区域就可以得到显式有限元数值模拟的3 维有限计算区。然而，何为适当在定量上并不容易。一般地，选择有限计算区域的尺度除应考虑地震断层本身的空间展布和预先估计的该地震断层的近断层强地面运动影响场展布范围外，还应考虑所采用的人工边界条件的精度和稳定性对计算区尺度的要求，应保证有限计算区域外的介质可以简化为均匀或水平成层介质，应考虑当前计算机的求解能力〔2,6,7〕。综合考虑以上因素并结合经验，对于图1所示有限断层，本文取有限计算区为一个如图3所示的长、宽、厚分别是X ′Y ′ Z = 60′60′ 40立方千米的长方体模型。该长方体的上表面为自由地表，其余的5 个表面均为人工边界面。另外，为计算方便（不失一般性），取有限计算区域3 维坐标系的原点在自由地表的左下角，x轴平行于断层走向。图3 中断层的上界埋深 留待2.3小节讨论给出。
　　
　　2.2 地壳结构模型
　　
　　进行近断层强地面运动影响场有限元数值模拟需要了解3 维有限计算区域内地壳介质的几何和物理结构，包括各地层介质的空间展布状态、各层介质的P波、S波传播速度以及密度、品质因子和泊桑比等参数，以求得能够正确反映传播途径对近断层强地面运动的影响。但不仅建立3维地壳结构模型本身存在诸多困难（如成本太高、误差很大），3 维地壳结构模型的使用还会大大增加计算机的存储需求和计算时间，有时将导致当前计算机的存储和计算能力无法承受。因此，作者认为，使用简化的均匀地壳结构模型或水平成层地壳结构模型是当前比较切合实际的选择。根据地质所提供的数据资料，某城市的简化水平成层地壳结构模型参数如表1所示。
　　
　　表 1 中各层介质的厚度是作者根据地质所提供的地层3 维勘探数据求得的平均值，其余皆为地质所提供的原始数据。
　　
　　2.3 有限计算区的空间离散化
　　
　　本文采用边长为Dx的正方体单元离散化图3 所示的有限计算区。根据离散化准则〔6〕，Dx应满足条件：
　　
　　p Dx = aV Dt （2）为方便起见并考虑到计算机的存储和计算能力，本文取Dt = 0.04秒，通过对各介质层取不同的a 值使各层的10 p aV = 千米/秒，相应的正方体单元的边长为Dx = 0.4千米。离散化后，有限计算区的节点总数为No =150′150′100 = 2250000 （3）每个节点有3 个自由度，总自由度为6750000个。根据作者的计算经验，对于采用显式有限元方法求解开放系统动力响应而言，这样一个大规模的自由度几乎是当前单个高端PC 计算机的内存和计算能力可以胜任的极限。需要指出：一般地，选择时间步距Dt = 0.04对于位移影响场的求解是适当的，但对于速度特别是加速度影响场的求解是过大和不适当的，不能以足够精度反映高频地震动的影响，可能使计算结果偏于不安全；通常，采用Dt = 0.01（或更小）解出的加速度影响场能以较好精度反映高频地震动的影响，但与此Dt 相应的Dx?0.1千米，采用以此Dx 为边长的正方体单元进行离散化，有限计算区的自由度将增加64 倍从而达到超4亿之多。不仅当前单个高端PC计算机的内存和计算能力无法胜任这一超大规模自由度系统的动力响应计算，就是基于单个高端PC 计算机组成的并行计算系统是否能够胜任也是一个需要研究的课题；本文对上述计算模型和有关参数的选择是根据实际情况和经验，考虑了高低频计算精度、稳定性、计算机的存储和计算能力、地震断层的规模和影响场的可能范围之间的综合平衡，做出的一个相对合理的选择。目前作者正在购置组建并行计算系统，预期至少可以使上述计算模型中的时空离散步距下降到Dt = 0.02秒（对应的理论上限频率为25 赫兹），Dx = 0.2千米，大大提高计算结果的参考和使用价值。另外，综合考虑断层的可能上界埋深和计算方便，在数值算例中将取图3 中断层的上界埋深为H = kDx，其中k =1,5,10为正整数。
　　
　　2.4 有限断层的空间离散化
　　
　　本文采用长方形单元离散化2维有限地震断层（如图1 所示），并且取断层长方形单元沿走向的边长Dl 等于有限计算区正方体单元的边长Dx，取断层长方形单元沿下倾方向的边长Dw = Dl / sin 75°。这即有助于保证使震源运动学模型可以激发出的高频和短波分量的阈值与可以在有限计算区离散网格内传播的高频和短波分量的阈值相一致，也便于将作用在断层平面单元节点上的运动（通过插值）转化为相应空间单元节点上的节点运动。
　　
　　3 内部节点阻尼动力平衡方程与解序列的显式计算方法
　　
　　离散化有限计算区中所有位于人工边界面内的节点通称人工边界节点，此外的所有节点（包括位于自由地表的节点）通称内部节点。离散化后，控制有限计算区连续介质时空状态的偏微分方程被转化为内部节点需要满足的关于时间变量的一组常微分方程，并且通常以矩阵形式表达为：
　　
　　 bb b b Ma +Cv + Kd = P - C v + K d （4）式（4）中：a、v、d 分别为内部节点的加速度、速度、位移列向量；M、C、K 分别为集中质量矩阵（对角阵）、阻尼矩阵和刚度矩阵（对称方阵），P为节点荷载列向量。b b v、d分别为人工边界节点的速度、位移列向量； b b C、K 分别为与人工边界节点相应的阻尼矩阵和刚度矩阵（非方阵），P 为节点荷载列向量。若采用与刚度矩阵成正比的瑞雷阻尼，设比例系数为b （该值可以根据品质因子确定），则式（4）可以改写为：　　
　　
　　本文称式（5）为开放系统内部节点的阻尼动力平衡方程。记所有内部节点的自由度（刚度矩阵K 的行数）为h N ，全部节点（内部与边界之和）的自由度为N ，所有人工边界节点的自由度为l N （刚度矩阵b K 的列数）。则式（5）中仅包含有h N 个微分方程，却包含有h l N = N + N 个待求函数。因此式（5）只可以用于求解内部节点的运动状态，人工边界节点的运动状态必须通过建立补充方程即人工边界条件进行求解。设t ￡ iDt 情况下全部节点的运动状态都是已知的，本文以式（5）作为基本方程采用文献〔8〕提出的3 阶显式方法求解所有内部节点在t = (i +1)Dt 时刻的运动状态，该方法的主要优点是具有良好的阻尼稳定性，限于篇幅本文对此不作展开介绍。需要强调指出：对于震源动力学模型，在每一步逐步积分中，式（5）右端的节点荷载列向量P 必须为已知，才可以求解所有内部节点在t = (i +1)Dt 时刻的运动状态；但震源运动学模型不要求已知P值，它通过将由式（1）规定的t = (i +1)Dt 时刻的断层节点运动状态赋值（通过插值）给相邻的有限计算区内部节点来间接实现地震荷载的输入。
　　
　　4 显式人工边界条件
　　
　　建立显式人工边界条件是近断层强地面运动影响场显式有限元数值模拟的重要环节。记式（5）右端中与人工边界节点相关的项为b Q ，则有Qb = -Kb (bvb + db ) （6）显然b Q 的物理意义是（近似）无限外域通过人工边界节点对有限计算区内部节点施加的节点作用力，并且式（6）本身就构成一组包含未知函数的隐式吸收人工边界条件（具体地说是一组粘弹性人工边界条件）〔6〕。式（6）中的b K 可以由相关单元刚度阵的集成获得，显式求解b Q 完全归结为如何建立求解人工边界节点位移b d 和速度b v 的显式方程，亦即建立显式透射人工边界条件〔6〕。一般地，高精度求解人工边界节点的位移和速度依赖基于平面波假定的单向外行波动运动学模型和由此导出的傍轴近似公式，文献〔6，9-11〕对此有较深入的讨论。本文采用文献〔4〕给出的标准离散形式的2 阶傍轴近似公式作为一般人工边界节点的人工边界条件；对于位于两个人工边界面交线上的人工边界节点特别是三个人工边界面的交点，则采用考虑优化和非平面波情形的2 阶多项式外推公式作为人工边界条件〔10-11〕，限于篇幅本文对此不做展开介绍。
　　
　　5 结语
　　
　　本文利用已经收到的有关城市《城市活断层探测与地震危险性评价》前期工作中得到的与建立有限断层震源模型和3维地壳结构模型相关的部分文字和数据资料、结合实际概要说明城市在设定地震作用下近断层强地面运动影响场显式有限元数值模拟涉及的模型和方法。