改进的特征量方法及其在震源动力学模拟中的应用_董森.pdf

第 卷 第期 年月地球物理学报 ，董森，于湘伟，章文波 改进的特征量方法及其在震源动力学模拟中的应用地球物理学报，（）：，：，.（），（）：，：改进的特征量方法及其在震源动力学模拟中的应用董森，于湘伟，章文波中国科学院大学地球与行星科学学院，北京 摘要基于曲线网格有限差分方法，针对震源动力学模型中的边界条件的问题，本文提出了一种改进的特征量方法，即将牵引力镜像方法和特征量方法相结合，在保守形式的方程中通过镜像操作完成单边导数的计算，然后有约束地调整边界处的速度和应力值，以此保证满足模型中的边界条件改进的特征量方法可用于处理曲线网格当中一般的边界条件，从而也可用于构建曲线网格有限差分法的分裂节点模型我们将改进的特征量方法分别用于自由表面和断层面，并分别与前人计算结果进行了比较计算结果表明，在相同的网格划分的情况下，改进的特征量方法能够得到与其他方法相似的计算结果，但改进方法的异常震荡更小，从而证明了本文提出的改进的特征量方法用于曲线网格有限差分法分裂节点模型的可行性关键词震源动力学；曲线网格有限差分法；边界条件；牵引力镜像方法；特征量方法 ：中图分类号 收稿日期 ，收修定稿基金项目国家自然科学基金项目（）资助第一作者简介董森，男，硕士研究生，主要从事震源动力学研究 ：通讯作者章文波，男，教授，博士生导师，主要从事强震地震学、震源物理及强地震动数值模拟研究 ：，（），；地 球 物 理 学 报（）卷引言地震的发生一般是大地构造应力场长期积累作用的结果，震源动力学模型是以断层附近应力场为基本参数来探究震源破裂过程，因而，震源动力学模型是从物理本质上来探究地震的孕育、发生和停止的过程，其研究结果对认知地震断层破裂的物理过程有着重要的意义由于震源动力学断层自发破裂涉及到非线性的摩擦准则，问题较为复杂，迄今尚没有一个理论解析解，因而，常用数值方法来研究这一问题 目前，最常用的 有 三 种 数 值 方 法：（）边 界 积 分 方 程 方 法（，；钱峰等，），这是一种半解析半数值方法，能够避免高频成分缺失所带来的影响，并且有能力处理具有复杂几何形状的断层模型，但边界积分方程方法极大地依赖于介质的均匀性，且不能处理复杂地形，因此其应用范围受到了局限；（）有限元方法（，；，），该方法既能处理复杂的介质和地形，也能处理一部分复杂几何形状的断层模型，但计 算 量 较 大；（）有 限 差 分 方 法（，；，），其原理简单，使用广泛，能够处理复杂介质，但传统的有限差分方法在面对复杂地形或与地表斜交的断层时会遇到困难 等（）通过调整不同方向的网格间距，使倾斜断层落在有限差分网格的对角线上，但这种改进方案仍有明显的局限性 直到曲线网格有限差分方法的提出，有限差分方法处理断层自发破裂问题才有了突破原有局限性的可能曲线网格有限差分方法由 和 （）提出，在采用 和 （）提出的辅助差分方程复频移完全匹 配 层吸 收边 界 之 后 发 展 完 善（，）曲线网格有限差分法起初用于地震波场模拟，但很快引起了震源动力学领域的学者的注意，因为这种方法能够处理复杂地形，同时也为处理一部分具有复杂几何形状的断层模型提供了可能，并且还保留了有限差分方法便于处理非均匀介质的优势有限差分方法处理断层自发破裂问题有三种断层模型可供选择：厚断层模型、应力过剩模型和分裂节点模型（，）王铭锋等（）将厚断层模型与曲线网格有限差分方法相结合，并在此 基 础 上 讨 论 了 地 形 对 破 裂 过 程 的 影 响但 和 （）将有限差分方法三种模型的计算结果进行对比，指出厚断层模型计算结果偏差较大，而与分裂节点模型相比，应力过剩模型收敛慢、效率低，而结合厚断层模型的曲线网格有限差分方法必然继承厚断层模型固有的缺点 与王铭锋等（）不同，等（）将分裂节点模型与曲线网格有限差分方法相结合，并随后成功处理了断层面与地表相交产生的奇异点（，），目前其方法已经被用于复杂断层破裂影响因素的研究（，）、复杂地表对简单断层破裂影响的研究（，）和实际地震破裂过程的模拟（，；，）本文尝试了与 等（）不同的另一种方式，从一般的边界条件出发，将牵引力镜像方法和特征量方法融合为改进的特征量方法，用于构建曲线网格中的分裂节点模型，并将计算结果与 等（）相比对，证明本文方法的可行性方法 曲线网格有限差分为处理复杂地形和复杂的断层面几何形态，和 （）发展了任意曲线网格有限差分方法 为了方便阐述改进的特征量方法，本文对任意曲线网格有限差分和牵引力镜像方法进行简单的介绍，方法细节可以参考 （）和 （）在不规则的物理区域和规则的计算区域之间建立某种映射关系，即进行坐标变换，如图所示图曲线网格示意图（，）曲线网格有限差分方法引入曲线网格，将不规则的物理区域（左图）映射到规则的计算区域 （）（）期董森等：改进的特征量方法及其在震源动力学模拟中的应用（，），（，），（，），（）其中，（，）代表物理区域中的某一点在笛卡尔坐标系中的坐标，（，）代表这一点在曲线网格中的坐标 而对于曲线网格当中格点的编号，我们用（，）表示以下我们假定，如果所研究的问题涉及自由表面或断层面，自由表面即地表上的每个点都满足，而断层面上的每个点都满足在曲线网格坐标下，三维弹性波动方程可以写成如下的形式：，（），（），（），（），（），（），（），（），（），（）.（）其中.这一组方程可以写成矩阵形式：，.（）在曲线网格有限差分方法中，采用如下具有四阶空间精度的同位网格 差分格式，可以数值计算方程右侧的应力和速度各变量对曲线网格坐标的导数（，；，；张伟，；，）：?，?，（）其中?和?分别代表变量在处对自变量的向前和向后差分，差分系数为：.，.，.，.，.（）在方程的右侧，实际上同时存在对曲线网格的三个坐标即、的偏导数，在求取偏导数的时候可能向前差分或向后差分，因此共有八种组合 将（）式右侧记作（），这八种组合为：?（）?，?，?，?（）?，?，?，?（）?，?，?，（）?（）?，?，?，?（）?，?，?，?（）?，?，?，?（）?，?，?，?（）?，?，?，.（）利用方程，通过各变量对空间坐标的导数，可以求得各变量对时间的导数，对时间进行积分即可得到各变量在下一时间步的值 曲线网格有限差分方法使用 方法进行数值积分，并且交替地使用（）式中不同的组合：（）?（），（）?（），（）?（），（）?（），（）（）（）（），（）其中：，.，.，（）边界条件的牵引力镜像方法为了使方程组（）（）在任意曲线网格上满足自由表面牵引力为零的边界条件：地 球 物 理 学 报（）卷（）和 （）、等（）提出了如下的牵引力镜像方法公式（）展开为：，（，），（，），（，）.（）曲线网格有限差分方法以如下的牵引力镜像方法完成（）式中在靠近自由表面处的更新，并保证更新之后对边界条件（）式的满足：在自由表面以上的区域，额外设置三层网格点，即、.令这三层网格点处的牵引力值为：，（）其中，而在自由表面处牵引力为零，这是由上一个时间步的计算结果对边界条件（）的满足所决定的利用曲线坐标系中算符的保守形式，（）式可改写为如下的形式：（，），（，），（，），（，），（，），（，），（，），（，），（，），（）其中：，.（）在这一形式中，对和的求导不受自由表面的影响，而对求导的量则与牵引力具有相似的形式，等于牵引力与的乘积，而后者只与网格剖分的方式有关 由于（）式的映射关系只存在于自由表面以下，因此实际在自由表面以上是不能求出的 若令自由表面以上的与自由表面以下对称地相 等，则 可 利 用（）式 求 出 自 由 表 面 以 上的值，从而（）式中对求导项可用（）式求出 由此，自由表面及附近的速度分量的更新得以完成若将（）式对时间求导，并将（）式代入，可得：，（）其中：，（，），（，），（，），（，），（，），（，），（，），（，），（，）.在（）式中，速度各分量对和的求导不受自由表面影响，因此可以利用该方程解出速度各分量在自由表面处对的导数，进而利用（）式完成自由表面上应力各分量的更新同时由于（）式是由自由表面边界条件得出的，更新后的应力分量自动满足自由表面边界条件在自由表面附近的和处，速度各分量对的导数同样无法用（）式直接求出用如下的紧致差分格式代替（）式可以解决这一问题：期董森等：改进的特征量方法及其在震源动力学模拟中的应用?，?.（）这一差分格式具有与（）式相同的四阶精度 求出速度各分量对的导数后，即可利用（）式完成自由表面附近的应力各分量的更新 边界条件的特征量方法上述的牵引力镜像方法解决自由表面边界条件的问题 但对于断层面，由于边界条件（即摩擦准则）形式要复杂得多，具有很强的非线性，因此很难利用（）式或与之类似的形式满足断层面上的边界条件 为此，我们将牵引力镜像方法与特征量方法结合起来，用于解决这一问题 先以非曲线的同位网格的自由表面边界条件为例，对特征量方法进行简要的介绍利用特征量方法从第个时间步的变量值推算第（）个时间步的变量值时，首先利用向外插值，得到自由表面以上的网格点的速度和应力值，那么就可以像处理内部的网格点一样，直接推算出自由表面处（设轴以向上为正，自由表面为）下一个时间步的速度和应力值，但这时还不能保证得到的应力值符合边界条件 我们可以对自由表面处的速度和应力值进行调整，使之符合边界条件 设调整后与调整前的差为和，则有：，.（）式中出现在等式右侧的量可根据自由表面边界条件求出特征量方法的稳定性由 等（）给出了数学上的证明（，）；此外，这种方法还可以在物理上作如下的理解：所需要解决的问题，是对自由表面处的速度和应力值进行调整使之满足给定的边界条件如果在调整之前，边界条件与当前的应力值刚好一致，那么对于这一新的问题，边界条件就自然满足，从而不需要进行调整 因此，根据弹性力学的叠加原理，这两个问题的差，决定了我们在自由表面处所需要进行的调整 两个问题的初始条件是相同的，所以它们的初始条件的差应为零，而两个问题的边界条件的差则为：（），（），（），（）其中（）为单位阶跃函数 在时刻，自由表面处速度和其他应力分量的值即为所求 由于在这里只是瞬时的响应，不受此后边界条件的影响，所以可以将边界条件的形式简化为阶跃函数 基于同样的考虑，自由表面其他位置的边界条件也不会对瞬时响应产生影响，因此可以将边界条件认为是均匀的，亦即不随、坐标变化，同时所研究的问题是发生在半无限均匀空间中，从而使这一问题得到充分的简化利用简化后的问题在、方向上的平移不变性，通过拉普拉斯变换法可以求出这一问题的解：（），（），（）.（）并可进一步得到：（），（）.（）在（）、（）式中，令、，即为待求的、，从而得到（）式特征量方法的更多信息及上述物理解释的推导过程列在附录 中 以上介绍的特征量方法，只是针对非曲线的网格（或称为常规网格）的，还没有考虑到自由表面和断层面的复杂几何形态因此，本文将对特征量方法加以改进，使之能与曲线网格及牵引力镜像方法相融合 改进的特征量方法在断层面的应用断层面和自由表面的最显著的差异，在于同时涉及到正反两个表面我们将坐标的数值较大的一侧的表面称为正表面（在设计网格时一般放在上盘一侧，以保证滑动量和滑动速率的符号与一般的习惯一致），将坐标的数值较小的一侧的表面称为负表面有限差 分 方 法 对 边 界 的 处 理 要 完 成 两 个 任务，一是完成边界及附近的变量更新，二是使更新地 球 物 理 学 报（）卷后的变量满足边界条件 在牵引力镜像方法中，第二个任务是利用（）式，在完成第一个任务的同时完成的；而在特征量方法中，第一个任务是利用向外插值，第二个任务是利用（）式进行调整，两个任务的完成是独立的 对于使用分裂节点法描述的断层面来说，两侧的物理量尽管都存在，但由于断层面产生的不连续性，不能跨越断层面进行差分计算 因此，这两个 任 务 仍 然 都 是 需 要 完 成的，没有因为断层面与自由表面的差异而发生变化 由于边界条件的复杂性，像（）式那样在完成第一个任务的同时也完成第二个任务，是不易实现的，因此我们提出像特征量方法那样，分别完成这两个任务对于第一个任务，由于牵引力镜像方法已经由实践证明，可以很好地应用于曲