一种近地表地震动模拟的网格分级方法
【专利摘要】本发明涉及一种近地表地震动模拟的网格分级方法，属于工程波动理论和地震学领域。该方法先建立被研究块体，然后建立被研究块体的动力平衡方程，进一步计算被研究块体围线上的应力和位移，通过递归方法，实现粗细网格连接处被研究块体的应力场、位移场和加速度场从t时刻至t+Δt时刻的计算，速度场从t-Δt/2时刻至t+Δt/2时刻的计算。本发明方法可采用不同步距的不规则网格对计算模型进行离散，模拟过程高效、灵活，且可大幅度提高计算效率。
【专利说明】一种近地表地震动模拟的网格分级方法
【技术领域】
[0001] 本发明属于工程波动理论和地震学领域，主要涉及近地表地震动的数值模拟，特 别是针对地表土层和基岩剪切波速差异较大的地质情况。
【背景技术】
[0002] 地表以下几十米至几千米深度范围内，通常存在着多个低波速的土（岩）层。研 宄表明，地表附近的低波速土层会很大程度影响地震动的幅值和频谱特性。近地表多为低 剪切波速土（岩）层，数值模拟技术要求离散网格的允许最大空间步长满足数值方法频散 特性，即对于同一频率的地震波，介质的剪切波速越低，则离散网格的空间步长越小。以上 因素造成了目前的数值模拟方法需在整个计算区域内，将高波速介质的空间步长按照低波 速介质的空间步长进行离散，因此大幅度提高了计算所需内存并降低计算速度。
[0003] 鉴于上述问题，一些研宄者用一种"非连续网格（Discontinuous Grids) "技术， 实现了针对不同波速采用不同步长网格模拟地震波在多层介质中的传播。Jastram等先后 提出了分别模拟声波和弹性波的变网格步长算法，该方法在X方向使用傅立叶变换法或有 限差分法，在z方向使用高阶有限差分法，时间积分采用快速展开法（REM)。Robertsson和 Holliger基于有限差分法，使用细网格离散起伏表面附近的低波速介质，实现了对复杂场 地条件地震动的高效率模拟。Aoi和Fujiwara基于有限差分法，以一个粗网格匹配奇数个 细网格，将非连续网格扩展到处理三维问题。Hayashi等将有限差分法的非连续网格技术用 于处理粘弹性介质的复杂场地条件问题。Tessmer，Kang和Baag在非连续网格技术的基础 上，用不同时间步长匹配不同步长的网格，进一步提高了计算效率。以上工作均在不同程度 上节约了计算成本，但对于网格形状要求较高，即需用矩形网格对计算区域进行离散，对于 处理起伏地表或内部不规则几何边界问题存在局限性。
[0004] 地震动模拟过程通常需建立很大的计算模型，因此实际操作中要划分大量的计算 网格，降低网格数量是降低计算成本的主要手段。根据近地表附近地球介质的不同，用不同 尺寸的网格进行离散，可降低网格数量，进一步降低计算成本，具有可观的实际应用价值。
【发明内容】
[0005] 发明目的：
[0006] 本发明涉及一种近地表地震动模拟的网格分级方法，其目的是为了降低数值模拟 的计算成本，节省计算机内存和提高计算效率。
[0007] 技术方案：
[0008] 本发明是通过以下技术方案来实现的：
[0009] 一种近地表地震动模拟的网格分级方法，其特征在于：该方法步骤如下：
[0010] (1)建立被研宄块体：
[0011] 将地球介质按照数值模拟要求分成区域I和区域II，区域I由不规则细网格 组成，区域II由不规则粗网格组成；区域I内的被研宄块体为m-n-0-p-q-r-s-t-m，区 域II内的被研宄块体为l-2-3-4-a-5-6-7-l ;区域I和区域II连接处的被研宄块体 a-b-c-d-e-f-g-h-i-j-k-1-a，是通过连接粗网格、细网格内线段构成；辅助计算网格是虚 线形成的四边形网格；
[0012] (2)建立被研宄块体的动力平衡方程：
[0013] 被研宄块体a-b-c-d-e-f-g-h-i-j-k-1-a的动力平衡方程为：
【权利要求】
1. 一种近地表地震动模拟的网格分级方法，其特征在于：该方法步骤如下： (1) 建立被研宄块体： 将地球介质按照数值模拟要求分成区域I和区域II，区域I由不规则细网格组 成，区域II由不规则粗网格组成；区域I内的被研宄块体为m-n-0-p-q-r-s-t-m，区 域II内的被研宄块体为l-2_3-4-a-5-6-7-l;区域I和区域II连接处的被研宄块体 a-b-c-d-e-f-g-h-i-j-k-1-a，是通过连接粗网格、细网格内线段构成；辅助计算网格是虚 线形成的四边形网格； (2) 建立被研宄块体的动力平衡方程： 被研宄块体a-b-c-d-e-f-g-h-i-j-k-1-a的动力平衡方程为：
(1) (2) 式中：P为被研宄块体的密度； 屯和U2分别为被研宄块体沿X轴和z轴方向的加速度；σχχ、σσxz是被研宄块体围线上的应力分量； 将式（1)和式（2)进一步改写为
(3) (4) 式中軋为被研宄块体的质量；η。。为节点L周围粗网格的数目； nfi为第jl个节点周围的细网格总数； nM为被研宄块体a-b-c-d-e-f-g-h-i-j-k-1-a内部所有节点的总数；nM同时也代表与一个粗网格匹配的细网格的数目，ηM= 1，3, 5,...； Sii和Sjj分别代表被研宄块体L周围在细网格内和粗网格内的围线段； (3) 计算被研宄块体围线上的应力： 被研宄块体围线上的应力为： σJj=λδJjUkjk+μ(uJ；^uijj) (5) 式中：λ和μ为介质的拉梅系数，Sij为Kronecker张量，uk，k=ux，x+uz，z; 式（6)右侧位移的空间导数为：
式中：Uxj和uzj(j= 1，2, 3, 4)为四边形网格四个角点的位移分量； Xj和h为四边形网格四个角点的坐标；ξ和η为映射坐标系下的坐标变量； (4) 位移的线性插值： ① 一个粗网格匹配三个细网格的情况（粗网格边上内部节点序号为1，2): Uix=2uKx/3+u Lx/3；uly= 2uKy/3+uLy/3 U2x=uKx/3+2u Lx/3；u2y=uKy/3+2uLy/3 (10) ② 一个粗网格匹配五个细网格的情况（粗网格边上内部节点序号为1，2, 3, 4): Uix=4uKx/5+u Lx/5；uly= 4uKy/5+uLy/5 U2x= 3uKx/5+2u Lx/5；u2y= 3u Ky/5+2uLy/5 U3x= 2uKx/5+3u Lx/5；u3y= 2u Ky/5+3uLy/5 u4x= ukx/5+4uLx/5；u4y= uKy/5+4uLy/5 (11) ③ 一个粗网格匹配七个细网格的情况（粗网格边上内部节点序号为1，2, 3, 4, 5, 6): Uix=6uKx/7+u Lx/7；uly= 6uKy/7+uLy/7 U2x= 5uKx/7+2u Lx/7；u2y= 5uKy/7+2uLy/7 U3x= 4uKx/7+3u Lx/7；u3y= 4uKy/7+3uLy/7 U4x=3uKx/7+4u Lx/7；u4y= 3uKy/7+4uLy/7 U5x= 2uKx/7+5u Lx/7；u5y= 2uKy/7+5uLy/7 U6x=uKx/7+6u Lx/7；u6y= uKy/7+6uLy/7 (12) ④ 一个粗网格匹配九个细网格的情况（粗网格边上内部节点序号为 1,2,3,4,5,6,7,8)： Uix= 8uKx/9+u Lx/9；uly= 8uKy/9+uLy/9U2x= 7uKx/9+2u Lx/9；u2y= 7uKy/9+2uLy/9 U3x= 6uKx/9+3u Lx/9；u3y= 6uKy/9+3uLy/9 u4x= 5uKx/9+4uLx/9 ;u4y= 5uKy/9+4uLy/9 U5x= 4uKx/9+5u Lx/9；u5y= 4uKy/9+5uLy/9 U6x=3uKx/9+6u Lx/9；u6y= 3uKy/9+6uLy/9 U7x= 2uKx/9+7u Lx/9；u7y= 2uKy/9+7uLy/9 u8x= ukx/9+8uLx/9；u8y=uKy/9+8uLy/9 (13)； (5) 网格分级方法的数值实现： 被研宄块体K和被研宄块体L的合力为： FxK= xk7_fxk6_fxk5+fxkl+fxll+fxl2 FzK= Zk7_fZk6_fZk5 +fZkl+fZll+fZ12 FxL= xl3+fxll+fxkl+fxk2+fxk3+fxk4 FzL= _f Z13+fZll+fZkl +fZk2+fZk3+fZk4 (14) 依据t时刻粗细网格内各条围线段上的应力分量σ(Λ.、和σ_〔_，积分可得t时刻被研 宄块体各条围线段上的内力/丨和/J，再由式（14)可将这些内力分配给相应的被研宄块体， 可求得t时刻粗细网格内各被研宄块体的合力S和<，进而由式（3)和式（4)给出各被研 宄块体在t时刻的加速度分量g和圮； 通过时间积分，可求得t+Λt/2时刻各被研宄块体的速度分量<+Λ〃2和<+Δ"2; 进一步通过时间积分，可求得t+Λt时刻各被研宄块体的位移分量4和圮；针对一个 粗网格匹配三个细网格的情况，由式（10)的插值关系可给出t+At时刻粗网格边上节点的 位移；将所有被研宄块体的位移分量代入式出)_(9)，可求得t+At时刻位移的空间导数， 将其代入式（5)可以得到t+At时刻被研宄块体各条围线段上的应力分量、σ=?和 根据以上递归方法，实现粗细网格连接处被研宄块体的应力场、位移场和加速度场从t时刻至t+Δt时刻的计算，速度场从t-Δt/2时刻至t+Δt/2时刻的计算。
2.根据权利要求1所述的近地表地震动模拟的网格分级方法，其特征在于：用不同步 距的不规则网格对计算模型进行离散，且粗细网格间无须共用节点。
【文档编号】G06F19/00GK104462814SQ201410739967
【公开日】2015年3月25日 申请日期:2014年12月3日 优先权日:2014年12月3日
【发明者】栾宇, 鲍文博, 黄志强, 金生吉 申请人:沈阳工业大学