CN104462814A - 一种近地表地震动模拟的网格分级方法 - Google Patents

Abstract

本发明涉及一种近地表地震动模拟的网格分级方法，属于工程波动理论和地震学领域。该方法先建立被研究块体，然后建立被研究块体的动力平衡方程，进一步计算被研究块体围线上的应力和位移，通过递归方法，实现粗细网格连接处被研究块体的应力场、位移场和加速度场从t时刻至t+Δt时刻的计算，速度场从t-Δt/2时刻至t+Δt/2时刻的计算。本发明方法可采用不同步距的不规则网格对计算模型进行离散，模拟过程高效、灵活，且可大幅度提高计算效率。

Description

一种近地表地震动模拟的网格分级方法

技术领域

本发明属于工程波动理论和地震学领域，主要涉及近地表地震动的数值模拟，特别是针对地表土层和基岩剪切波速差异较大的地质情况。

背景技术

地表以下几十米至几千米深度范围内，通常存在着多个低波速的土(岩)层。研究表明，地表附近的低波速土层会很大程度影响地震动的幅值和频谱特性。近地表多为低剪切波速土(岩)层，数值模拟技术要求离散网格的允许最大空间步长满足数值方法频散特性，即对于同一频率的地震波，介质的剪切波速越低，则离散网格的空间步长越小。以上因素造成了目前的数值模拟方法需在整个计算区域内，将高波速介质的空间步长按照低波速介质的空间步长进行离散，因此大幅度提高了计算所需内存并降低计算速度。

鉴于上述问题，一些研究者用一种“非连续网格(Discontinuous Grids)”技术，实现了针对不同波速采用不同步长网格模拟地震波在多层介质中的传播。Jastram等先后提出了分别模拟声波和弹性波的变网格步长算法，该方法在x方向使用傅立叶变换法或有限差分法，在z方向使用高阶有限差分法，时间积分采用快速展开法(REM)。Robertsson和Holliger基于有限差分法，使用细网格离散起伏表面附近的低波速介质，实现了对复杂场地条件地震动的高效率模拟。Aoi和Fujiwara基于有限差分法，以一个粗网格匹配奇数个细网格，将非连续网格扩展到处理三维问题。Hayashi等将有限差分法的非连续网格技术用于处理粘弹性介质的复杂场地条件问题。Tessmer，Kang和Baag在非连续网格技术的基础上，用不同时间步长匹配不同步长的网格，进一步提高了计算效率。以上工作均在不同程度上节约了计算成本，但对于网格形状要求较高，即需用矩形网格对计算区域进行离散，对于处理起伏地表或内部不规则几何边界问题存在局限性。

地震动模拟过程通常需建立很大的计算模型，因此实际操作中要划分大量的计算网格，降低网格数量是降低计算成本的主要手段。根据近地表附近地球介质的不同，用不同尺寸的网格进行离散，可降低网格数量，进一步降低计算成本，具有可观的实际应用价值。

发明内容

发明目的：

本发明涉及一种近地表地震动模拟的网格分级方法，其目的是为了降低数值模拟的计算成本，节省计算机内存和提高计算效率。

技术方案：

本发明是通过以下技术方案来实现的：

一种近地表地震动模拟的网格分级方法，其特征在于：该方法步骤如下：

(1)建立被研究块体：

将地球介质按照数值模拟要求分成区域I和区域II，区域I由不规则细网格组成，区域II由不规则粗网格组成；区域I内的被研究块体为m-n-o-p-q-r-s-t-m，区域II内的被研究块体为1-2-3-4-a-5-6-7-1；区域I和区域II连接处的被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a，是通过连接粗网格、细网格内线段构成；辅助计算网格是虚线形成的四边形网格；

(2)建立被研究块体的动力平衡方程：

被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a的动力平衡方程为：

式中：ρ为被研究块体的密度；

ü_x和ü_z分别为被研究块体沿x轴和z轴方向的加速度；

σ_xx、σ_zz和σ_xz是被研究块体围线上的应力分量；

将式(1)和式(2)进一步改写为

$M_L{\stackrel{\·\·}{u}}_{\mathrm{Lx}}=\underset{\mathrm{ii}=1}{\overset{2(n_{\mathrm{ra}}+1)}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sii}}({\mathrm{\σ}}_{\mathrm{xx}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zx}}\mathrm{dx})+\underset{\mathrm{jj}=1}{\overset{{2n}_{\mathrm{co}}}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sjj}}({\mathrm{\σ}}_{\mathrm{xx}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zx}}\mathrm{dx})---\left(3\right)$

$M_L{\stackrel{\·\·}{u}}_{\mathrm{Lz}}=\underset{\mathrm{ii}=1}{\overset{2(n_{\mathrm{ra}}+1)}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sii}}({\mathrm{\σ}}_{\mathrm{xz}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zz}}\mathrm{dx})+\underset{\mathrm{jj}=1}{\overset{{2n}_{\mathrm{co}}}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sjj}}({\mathrm{\σ}}_{\mathrm{xz}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zz}}\mathrm{dx})---\left(4\right)$

式中：M_L为被研究块体的质量；

n_co为节点L周围粗网格的数目；

n_fi为第j1个节点周围的细网格总数；

n_ra为被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a内部所有节点的总数；

n_ra同时也代表与一个粗网格匹配的细网格的数目，n_ra＝1,3,5,...；

Sii和Sjj分别代表被研究块体L周围在细网格内和粗网格内的围线段；

(3)计算被研究块体围线上的应力：

被研究块体围线上的应力为：

σ_ij＝λδ_iju_k,k+μ(u_j,i+u_i,j)               (5)

式中：λ和μ为介质的拉梅系数，δ_ij为Kronecker张量，u_k,k＝u_x,x+u_z,z；

式(6)右侧位移的空间导数为：

${\∂u}_{\mathrm{xi}}/\∂x=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}[(1-\mathrm{\η})z_2+(\mathrm{\η}-\mathrm{\ξ})z_3+(\mathrm{\ξ}-1)z_4]u_{x1}+[(\mathrm{\η}-1)z_1+(1+\mathrm{\ξ})z_3+(-\mathrm{\η}-\mathrm{\ξ})z_4]u_{x2}+\\ [(\mathrm{\ξ}-\mathrm{\η})z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\η})z_4]u_{x3}+[(1-\mathrm{\ξ})z_1+(\mathrm{\ξ}+\mathrm{\η})z_2+(-1-\mathrm{\η})z_3]u_{x4}\end{array}\right\}---\left(6\right)$

${\∂u}_{\mathrm{zi}}/\∂z=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}-[(1-\mathrm{\η})x_2+(\mathrm{\η}-\mathrm{\ξ})x_3+(\mathrm{\ξ}-1)x_4]u_{z1}-[(\mathrm{\η}-1)x_1+(1+\mathrm{\ξ})x_3+(-\mathrm{\η}-\mathrm{\ξ})x_4]u_{z2}-\\ [(\mathrm{\ξ}-\mathrm{\η})x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\η})x_4]u_{z3}-[(1-\mathrm{\ξ})x_1+(\mathrm{\ξ}+\mathrm{\η})x_2+(-1-\mathrm{\η})x_3]u_{z4}\end{array}\right\}---\left(7\right)$

${\∂u}_{\mathrm{xi}}/\∂z=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}-[(1-\mathrm{\η})x_2+(\mathrm{\η}-\mathrm{\ξ})x_3+(\mathrm{\ξ}-1)x_4]u_{x1}-[(\mathrm{\η}-1)x_1+(1+\mathrm{\ξ})x_3+(-\mathrm{\η}-\mathrm{\ξ})x_4]u_{x2}-\\ [(\mathrm{\ξ}-\mathrm{\η})x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\η})x_4]u_{x3}-[(1-\mathrm{\ξ})x_1+(\mathrm{\ξ}+\mathrm{\η})x_2+(-1-\mathrm{\η})x_3]u_{x4}\end{array}\right\}---\left(8\right)$

${\∂u}_{\mathrm{zi}}/\∂x=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}[(1-\mathrm{\η})z_2+(\mathrm{\η}-\mathrm{\ξ})z_3+(\mathrm{\ξ}-1)z_4]u_{z1}+[(\mathrm{\η}-1)z_1+(1+\mathrm{\ξ})z_3+(-\mathrm{\η}-\mathrm{\ξ})z_4]u_{z2}+\\ [(\mathrm{\ξ}-\mathrm{\η})z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\η})z_4]u_{z3}+[(1-\mathrm{\ξ})z_1+(\mathrm{\ξ}+\mathrm{\η})z_2+(-1-\mathrm{\η})z_3]u_{z4}\end{array}\right\}---\left(9\right)$

式中：u_xj和u_zj(j＝1,2,3,4)为四边形网格四个角点的位移分量；

x_j和z_j为四边形网格四个角点的坐标；

ξ和η为映射坐标系下的坐标变量；

(4)位移的线性插值：

①一个粗网格匹配三个细网格的情况(粗网格边上内部节点序号为1,2)：

u_1x＝2u_Kx/3+u_Lx/3；u_1y＝2u_Ky/3+u_Ly/3

u_2x＝u_Kx/3+2u_Lx/3；u_2y＝u_Ky/3+2u_Ly/3               (10)

②一个粗网格匹配五个细网格的情况(粗网格边上内部节点序号为1,2,3,4)：

u_1x＝4u_Kx/5+u_Lx/5；u_1y＝4u_Ky/5+u_Ly/5

u_2x＝3u_Kx/5+2u_Lx/5；u_2y＝3u_Ky/5+2u_Ly/5

u_3x＝2u_Kx/5+3u_Lx/5；u_3y＝2u_Ky/5+3u_Ly/5

u_4x＝u_Kx/5+4u_Lx/5；u_4y＝u_Ky/5+4u_Ly/5            (11)

③一个粗网格匹配七个细网格的情况(粗网格边上内部节点序号为1,2,3,4,5,6)：

u_1x＝6u_Kx/7+u_Lx/7；u_1y＝6u_Ky/7+u_Ly/7

u_2x＝5u_Kx/7+2u_Lx/7；u_2y＝5u_Ky/7+2u_Ly/7

u_3x＝4u_Kx/7+3u_Lx/7；u_3y＝4u_Ky/7+3u_Ly/7

u_4x＝3u_Kx/7+4u_Lx/7；u_4y＝3u_Ky/7+4u_Ly/7

u_5x＝2u_Kx/7+5u_Lx/7；u_5y＝2u_Ky/7+5u_Ly/7

u_6x＝u_Kx/7+6u_Lx/7；u_6y＝u_Ky/7+6u_Ly/7                 (12)

④一个粗网格匹配九个细网格的情况(粗网格边上内部节点序号为1,2,3,4,5,6,7,8)：

u_1x＝8u_Kx/9+u_Lx/9；u_1y＝8u_Ky/9+u_Ly/9

u_2x＝7u_Kx/9+2u_Lx/9；u_2y＝7u_Ky/9+2u_Ly/9

u_3x＝6u_Kx/9+3u_Lx/9；u_3y＝6u_Ky/9+3u_Ly/9

u_4x＝5u_Kx/9+4u_Lx/9；u_4y＝5u_Ky/9+4u_Ly/9

u_5x＝4u_Kx/9+5u_Lx/9；u_5y＝4u_Ky/9+5u_Ly/9

u_6x＝3u_Kx/9+6u_Lx/9；u_6y＝3u_Ky/9+6u_Ly/9

u_7x＝2u_Kx/9+7u_Lx/9；u_7y＝2u_Ky/9+7u_Ly/9

u_8x＝u_Kx/9+8u_Lx/9；u_8y＝u_Ky/9+8u_Ly/9               (13)；

(5)网格分级方法的数值实现：

被研究块体K和被研究块体L的合力为：

F_xK＝-f_xk7-f_xk6-f_xk5+f_xk1+f_xl1+f_xl2

F_zK＝-f_zk7-f_zk6-f_zk5+f_zk1+f_zl1+f_zl2

F_xL＝-f_xl3+f_xl1+f_xk1+f_xk2+f_xk3+f_xk4

F_zL＝-f_zl3+f_zl1+f_zk1+f_zk2+f_zk3+f_zk4            (14)

依据t时刻粗细网格内各条围线段上的应力分量和积分可得t时刻被研究块体各条围线段上的内力和再由式(14)可将这些内力分配给相应的被研究块体，可求得t时刻粗细网格内各被研究块体的合力和进而由式(3)和式(4)给出各被研究块体在t时刻的加速度分量和

通过时间积分，可求得t+Δt/2时刻各被研究块体的速度分量和

进一步通过时间积分，可求得t+Δt时刻各被研究块体的位移分量和针对一个粗网格匹配三个细网格的情况，由式(10)的插值关系可给出t+Δt时刻粗网格边上节点的位移；将所有被研究块体的位移分量代入式(6)-(9)，可求得t+Δt时刻位移的空间导数，将其代入式(5)可以得到t+Δt时刻被研究块体各条围线段上的应力分量和

根据以上递归方法，实现粗细网格连接处被研究块体的应力场、位移场和加速度场从t时刻至t+Δt时刻的计算，速度场从t-Δt/2时刻至t+Δt/2时刻的计算。

可用不同步距的不规则网格对计算模型进行离散，且粗细网格间无须共用节点。

优点及效果：

本发明涉及一种近地表地震动模拟的网格分级方法，具有如下优点：

(1)计算模型可采用不规则的粗细网格对介质进行离散，便于建立不规则地质结构的盆地模型，模拟过程高效、灵活，且可大幅度提高计算效率；

(3)可使用1个粗网格匹配奇数个(3,5,7,9)细网格,大幅度的提高了计算效率；

(2)不同波速的介质可采用不同空间步长的四边形网格进行离散(即盆地低波速介质采用较细网格，高波速介质采用较粗的网格)，有效地解决了一般由低波速介质决定整个计算区域空间步长的问题。

附图说明

图1为被研究块体示意图；

图2为粗细网格内被研究块体示意图；

图3为与Lamb解析解对比的计算模型；

图4为数值解与Lamb解析解的对比结果(观测点r1)；

图5为数值解与Lamb解析解的对比结果(观测点r2)；

图6为竖直位移的波场图；

图7为沉积盆地和起伏地表的计算模型；

图8为局部粗细网格的划分；

图9为有无盆地情况下观测点的水平峰值加速度，双力偶点源的倾角为10°；

图10为有无盆地情况下观测点的水平峰值加速度，双力偶点源的倾角为45°；

图11为盆地和无盆地情况下观测点的水平峰值加速度，双力偶点源的倾角为80°。

具体实施方式

下面结合附图和具体的实施例对本发明做进一步的说明，但本发明的保护范围不受实施例的限制。

被研究块体方法用于地震动的模拟，该方法构造的被研究块体物理意义明确，便于处理地表起伏和内部不规则几何边界问题。本发明构造了新的被研究块体以模拟波动在不同波速介质之间相互传播，并使该方法能用不规则网格(或不规则被研究块体)对计算区域进行离散，高效、灵活的进行近地表复杂场地条件地震动的模拟。

本发明提出了一种近地表地震动模拟的网格分级方法，具体步骤如下：

步骤(1)：被研究块体的建立

图1是区域Ⅰ内部、区域Ⅱ内部和两区域连接处被研究块体示意图。区域Ⅰ、Ⅱ均由不规则细、粗网格组成。实线m-n-o-p-q-r-s-t-m包围的封闭区域为区域Ⅰ内部的被研究块体；实线1-2-3-4-a-5-6-7-1包围的封闭区域为区域Ⅱ内部的被研究块体；实线a-b-c-d-e-f-g-h-i-j-k-l-a包围的封闭区域为区域Ⅰ和区域Ⅱ连接处的被研究块体，它是通过连接粗、细网格内线段构成。图中虚线围成的四边形为构造被研究块体的辅助网格，小空心圆代表粗细网格的节点，小实心圆代表网格几何中心或网格四边上的中点。

步骤(2)：动力平衡方程的建立

针对被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a，根据牛顿第二定律可得到被研究块体的动力平衡方程为

以上两式中Ω为围线a-b-c-d-e-f-g-h-i-j-k-l-a所封闭的区域，ρ为被研究块体的密度，ü_x和ü_z分别为被研究块体沿x轴和z轴方向的加速度，σ_xx、σ_zz和σ_xz是被研究块体围线上的应力分量。若先不考虑作用在被研究块体上的外力，则式(1)和式(2)的右侧实际上是来自周围其它被研究块体沿x轴和z轴方向的相互作用力。

假定被研究块体区域Ω内各点的加速度相等，则被研究块体的加速度分量定义在图1中的节点L上。同理被研究块体Ω的质量也可以定义在节点L上，则式(1)和式(2)的左侧可改写为M_Lü_Lx和M_Lü_Lz

$M_L=\underset{i1=1}{\overset{n_{\mathrm{co}}}{\mathrm{\Σ}}}{M}_{i1}/4+\underset{j1=1}{\overset{n_{\mathrm{ra}}}{\mathrm{\Σ}}}\underset{i2=1}{\overset{n_{\mathrm{fi}}}{\mathrm{\Σ}}}{M}_{i2}/4,---\left(4\right)$

上式中：M_i1为节点L周围第i1个粗网格的质量，n_co为节点L周围粗网格的数目；M_i2为环绕在第j1个节点周围的第i1个细网格的质量，n_fi为第j1个节点周围的细网格总数，n_ra为被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a内部所有节点的总数，n_ra同时也代表与一个粗网格匹配的细网格的数目，n_ra＝1,3,5,...。

因此，可将式(1)和式(2)进一步改写为

$M_L{\stackrel{\·\·}{u}}_{\mathrm{Lx}}=\underset{\mathrm{ii}=1}{\overset{2(n_{\mathrm{ra}}+1)}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sii}}({\mathrm{\σ}}_{\mathrm{xx}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zx}}\mathrm{dx})+\underset{\mathrm{jj}=1}{\overset{{2n}_{\mathrm{co}}}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sjj}}({\mathrm{\σ}}_{\mathrm{xx}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zx}}\mathrm{dx})---\left(4\right)$

$M_L{\stackrel{\·\·}{u}}_{\mathrm{Lz}}=\underset{\mathrm{ii}=1}{\overset{2(n_{\mathrm{ra}}+1)}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sii}}({\mathrm{\σ}}_{\mathrm{xz}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zz}}\mathrm{dx})+\underset{\mathrm{jj}=1}{\overset{{2n}_{\mathrm{co}}}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sjj}}({\mathrm{\σ}}_{\mathrm{xz}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zz}}\mathrm{dx})---\left(5\right)$

上式中Sii和Sjj分别代表被研究块体L周围在细网格内和粗网格内的围线段，例如，S1代表线段l-k，S2代表线段a-b(见图1)。

步骤(3)：被研究块体围线上的应力

依据弹性介质中的应力-应变关系，式(4)和式(5)右侧的应力分量可由位移表示为

σ_ij＝λδ_iju_k,k+μ(u_j,i+u_i,j)               (6)式中λ和μ为介质的拉梅系数，δ_ij为Kronecker张量，u_k,k＝u_x,x+u_z,z。

根据等参元原理，可得式(6)右侧位移的空间导数如下

${\∂u}_{\mathrm{xi}}/\∂x=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}[(1-\mathrm{\η})z_2+(\mathrm{\η}-\mathrm{\ξ})z_3+(\mathrm{\ξ}-1)z_4]u_{x1}+[(\mathrm{\η}-1)z_1+(1+\mathrm{\ξ})z_3+(-\mathrm{\η}-\mathrm{\ξ})z_4]u_{x2}+\\ [(\mathrm{\ξ}-\mathrm{\η})z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\η})z_4]u_{x3}+[(1-\mathrm{\ξ})z_1+(\mathrm{\ξ}+\mathrm{\η})z_2+(-1-\mathrm{\η})z_3]u_{x4}\end{array}\right\}---\left(7\right)$

${\∂u}_{\mathrm{zi}}/\∂z=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}-[(1-\mathrm{\η})x_2+(\mathrm{\η}-\mathrm{\ξ})x_3+(\mathrm{\ξ}-1)x_4]u_{z1}-[(\mathrm{\η}-1)x_1+(1+\mathrm{\ξ})x_3+(-\mathrm{\η}-\mathrm{\ξ})x_4]u_{z2}-\\ [(\mathrm{\ξ}-\mathrm{\η})x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\η})x_4]u_{z3}-[(1-\mathrm{\ξ})x_1+(\mathrm{\ξ}+\mathrm{\η})x_2+(-1-\mathrm{\η})x_3]u_{z4}\end{array}\right\}---\left(8\right)$

${\∂u}_{\mathrm{xi}}/\∂z=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}-[(1-\mathrm{\η})x_2+(\mathrm{\η}-\mathrm{\ξ})x_3+(\mathrm{\ξ}-1)x_4]u_{x1}-[(\mathrm{\η}-1)x_1+(1+\mathrm{\ξ})x_3+(-\mathrm{\η}-\mathrm{\ξ})x_4]u_{x2}-\\ [(\mathrm{\ξ}-\mathrm{\η})x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\η})x_4]u_{x3}-[(1-\mathrm{\ξ})x_1+(\mathrm{\ξ}+\mathrm{\η})x_2+(-1-\mathrm{\η})x_3]u_{x4}\end{array}\right\}---\left(9\right)$

${\∂u}_{\mathrm{zi}}/\∂x=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}[(1-\mathrm{\η})z_2+(\mathrm{\η}-\mathrm{\ξ})z_3+(\mathrm{\ξ}-1)z_4]u_{z1}+[(\mathrm{\η}-1)z_1+(1+\mathrm{\ξ})z_3+(-\mathrm{\η}-\mathrm{\ξ})z_4]u_{z2}+\\ [(\mathrm{\ξ}-\mathrm{\η})z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\η})z_4]u_{z3}+[(1-\mathrm{\ξ})z_1+(\mathrm{\ξ}+\mathrm{\η})z_2+(-1-\mathrm{\η})z_3]u_{z4}\end{array}\right\}---\left(10\right)$

式中u_xj和u_zj(j＝1,2,3,4)为四边形网格四个角点的位移分量，x_j和z_j为四边形网格四个角点的坐标，ξ和η为映射坐标系下的坐标变量， $\left|J\right|=(\∂x/\∂\mathrm{\ξ})(\∂z/\∂\mathrm{\η})-(\∂z/\∂\mathrm{\ξ})(\∂x/\∂\mathrm{\η})$ 为雅克比行列式。

式(4)和式(5)右侧的dx和dy，可对ξ和η求导数得到它们的全微分形式如下

$\begin{array}{c}\mathrm{dz}=\frac{\∂z}{\∂\mathrm{\ξ}}\mathrm{d\ξ}+\frac{\∂z}{\∂\mathrm{\η}}\mathrm{d\η}=(\mathrm{\η}-1)z_1+(1-\mathrm{\η})z_2+(1+\mathrm{\η})z_3+(-1-\mathrm{\η})z_4]\mathrm{d\ξ}/4+\\ [(\mathrm{\ξ}-1)z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\ξ})z_3+(1-\mathrm{\ξ})z_4]\mathrm{d\η}/4\end{array}---\left(11\right)$

$\begin{array}{c}\mathrm{dx}=\frac{\∂x}{\∂\mathrm{\ξ}}\mathrm{d\ξ}+\frac{\∂x}{\∂\mathrm{\η}}\mathrm{d\η}=[(\mathrm{\η}-1)x_1+(1-\mathrm{\η})x_2+(1+\mathrm{\η})x_3+(-1-\mathrm{\η})x_4]\mathrm{d\ξ}/4+\\ [(\mathrm{\ξ}-1)x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\ξ})x_3+(1-\mathrm{\ξ})x_4]\mathrm{d\η}/4\end{array}---\left(12\right)$

步骤(4)：位移的线性插值

从式(7)-(10)可以看到，为得到位移的空间导数，应已知域内所有网格的节点位移。在粗细网格连接处，粗网格边上的节点位移可通过粗网格的节点位移的线性插值给出，例如，图1中节点L3和节点L1的位移分量可由节点K和节点L的位移插值得到，粗细网格的共用节点K和L具有相同的位移。

(1)一个粗网格匹配三个细网格的情况(粗网格边上内部节点序号为1,2)：

u_1x＝2u_Kx/3+u_Lx/3；u_1y＝2u_Ky/3+u_Ly/3

u_2x＝u_Kx/3+2u_Lx/3；u_2y＝u_Ky/3+2u_Ly/3           (13)

(2)一个粗网格匹配五个细网格的情况(粗网格边上内部节点序号为1,2,3,4)：

u_1x＝4u_Kx/5+u_Lx/5；u_1y＝4u_Ky/5+u_Ly/5

u_2x＝3u_Kx/5+2u_Lx/5；u_2y＝3u_Ky/5+2u_Ly/5

u_3x＝2u_Kx/5+3u_Lx/5；u_3y＝2u_Ky/5+3u_Ly/5

u_4x＝u_Kx/5+4u_Lx/5；u_4y＝u_Ky/5+4u_Ly/5            (14)

(3)一个粗网格匹配七个细网格的情况(粗网格边上内部节点序号为1,2,3,4,5,6)：

u_1x＝6u_Kx/7+u_Lx/7；u_1y＝6u_Ky/7+u_Ly/7

u_2x＝5u_Kx/7+2u_Lx/7；u_2y＝5u_Ky/7+2u_Ly/7

u_3x＝4u_Kx/7+3u_Lx/7；u_3y＝4u_Ky/7+3u_Ly/7

u_4x＝3u_Kx/7+4u_Lx/7；u_4y＝3u_Ky/7+4u_Ly/7

u_5x＝2u_Kx/7+5u_Lx/7；u_5y＝2u_Ky/7+5u_Ly/7

u_6x＝u_Kx/7+6u_Lx/7；u_6y＝u_Ky/7+6u_Ly/7               (15)

(4)一个粗网格匹配九个细网格的情况(粗网格边上内部节点序号为1,2,3,4,5,6,7,8)：

u_1x＝8u_Kx/9+u_Lx/9；u_1y＝8u_Ky/9+u_Ly/9

u_2x＝7u_Kx/9+2u_Lx/9；u_2y＝7u_Ky/9+2u_Ly/9

u_3x＝6u_Kx/9+3u_Lx/9；u_3y＝6u_Ky/9+3u_Ly/9

u_4x＝5u_Kx/9+4u_Lx/9；u_4y＝5u_Ky/9+4u_Ly/9

u_5x＝4u_Kx/9+5u_Lx/9；u_5y＝4u_Ky/9+5u_Ly/9

u_6x＝3u_Kx/9+6u_Lx/9；u_6y＝3u_Ky/9+6u_Ly/9

u_7x＝2u_Kx/9+7u_Lx/9；u_7y＝2u_Ky/9+7u_Ly/9

u_8x＝u_Kx/9+8u_Lx/9；u_8y＝u_Ky/9+8u_Ly/9            (16)

由以上公式，可得到各种情况下粗网格边上各节点的位移分量。

步骤(5)：算法的实现

以粗网格K1-K2-L-K和与其相匹配的三个细网格为例(图2)，被研究块体围线上的内力f可由以下各式分配给被研究块体K和被研究块体L：

F_xK＝-f_xk7-f_xk6-f_xk5+f_xk1+f_xl1+f_xl2

F_zK＝-f_zk7-f_zk6-f_zk5+f_zk1+f_zl1+f_zl2

F_xL＝-f_xl3+f_xl1+f_xk1+f_xk2+f_xk3+f_xk4

F_zL＝-f_zl3+f_zl1+f_zk1+f_zk2+f_zk3+f_zk4             (17)

式中F_xK和F_zK为被研究块体K在x轴和z轴方向的合力，F_xL和F_zL为被研究块体L在x轴和z轴方向的合力。将式(17)在所有粗细网格上循环计算，可得到粗细网格连接处所有被研究块体的合力。

依据t时刻粗细网格内各条围线段上的应力分量和积分可得t时刻被研究块体各条围线段上的内力和再由式(17)可将这些内力分配给相应的被研究块体，可求得t时刻粗细网格内各被研究块体的合力和进而由式(4)和式(5)给出各被研究块体在t时刻的加速度分量和

通过时间积分，可求得t+Δt/2时刻各被研究块体的速度分量和

进一步通过时间积分，可求得t+Δt时刻各被研究块体的位移分量和针对一个粗网格匹配三个细网格的情况，由式(13)的插值关系可给出t+Δt时刻粗网格边上节点的位移。将所有被研究块体的位移分量代入式(7)-(10)，可求得t+Δt时刻位移的空间导数，将其代入式(6)可以得到t+Δt时刻被研究块体各条围线段上的应力分量和

根据以上递归方法，可实现粗细网格连接处被研究块体的应力场、位移场和加速度场从t时刻至t+Δt时刻的计算，速度场从t-Δt/2时刻至t+Δt/2时刻的计算。

实施例1：与Lamb解析解结果对比

图3为与Lamb解析解对比的计算模型。图中作用于弹性半空间表面竖直向下的力源F(t)＝exp[-α(t-t₀)]，α＝90和t₀＝0.3s。计算模型的较小区域采用10m的正方形网格进行离散，较大的区域分别采用30m、50m、70m和90m的正方形网格进行离散。两观测点r₁和r₂距加载点S1分别为2000m和5310m，观测点r₁位于细网格表面，观测点r₂位于粗网格表面且与粗细网格分界线(图3中的虚线)相距360m。弹性半空间中介质的P波波速为5196m/s，SV波波速为3000m/s，密度为2300Kg/m³。计算采用的时间步长为0.001s，计算时间为4s。

图4、5为观测点r₁、r₂处Lamb解析解与本文提出方法得到数值解的对比结果。结果显示，粗网格采用不同步长(30m、50m、70m和90m)，观测点r₁、r₂处曲线峰值与解析解的误差分别为1.4％和1.0％。由此可得出，本文提出的方法在处理Lamb问题时具有良好的计算精度，同时粗网格步长的改变并不影响计算结果。

图6为竖直位移的波场图。图中的位移波场清晰显示了当波通过粗细网格时，无论是水平方向还是竖直方向，波动都能够很好的从细网格一侧向粗网格一侧传播，没有发生任何的异常现象(如反射现象)。

实施例1证实了本发明所开发的网格分级方法的正确性和可行性。

实施例2：模拟复杂场地条件近地表地震动

图7为沉积盆地和起伏地表的计算模型。图中观测点Ra位于左侧山峰的顶点，观测点Rc(位于盆地左侧的边缘)与Ra的相对水平距离和竖直距离分别为2.4km和1.2km；观测点Rb位于右侧山峰的顶点，观测点Rd(位于盆地右侧的边缘)与Rb的相对水平距离和竖直距离分别为1.9km和0.8km。位于左侧的山腰上，观测点Ra与Rc之间共设置19个观测点；位于沉积盆地表面，观测点Rc与Rd之间共设置28个观测点；位于右侧的山腰上，观测点Rd与Rb之间共设置12个观测点。震源为双力偶点源并位于盆地中心表面下方6km处，震源的时间函数为D(t)＝D(∞)[1-(1+t/τ)e^-t/τ]，式中D(∞)为断层上质点运动从零到静止的滑动位移，τ为上升时间。盆地深度约为650m，其内部介质的P波波速为1100m/s，SV波波速为600m/s，密度为2000Kg/m³；基岩介质的P波波速为3500m/s，SV波波速为2000m/s，密度为2600Kg/m³。计算采用的时间步长为0.005s，计算时间为15s。通过改变双力偶震源模型的倾角，研究沉积盆地和起伏地表区域表面地震动的特性。

图8展示了沉积盆地和起伏地表的计算模型的网格划分情况。根据盆地内介质与基岩介质波速的关系，盆地内部较低波速介质采用约50m的四边形网格进行离散，基岩内部较高波速介质采用约150m的四边形网格进行离散。图8(a)和(c)展示了盆地两侧边缘粗细网格的离散情况，图8(b)展示了盆地下侧粗细网格的离散情况。从离散化模型可以看到，整个计算区域均采用非规则的四边形网格进行划分，对不同波速的介质，可分别采用适应其波速(或波长)的不同空间步距的四边形网格进行离散，粗细网格在连接区域能够良好的匹配。以下为不同倾角的双力偶源在地表各观测点的数值模拟结果：

(1)倾角为10°双力偶点源

图9为双力偶点源倾角10°，有无盆地情况观测点的水平峰值加速度。从以上结果可以看到，山谷表面的水平峰值地面加速度(PGA)多数大于两侧的山峰；当盆地存在时，山谷表面的水平PGA由盆地两侧边缘(观测点Rc和Rd)向中央呈现急剧放大的趋势。由于盆地的存在，使山谷表面各点PGA的平均值较无盆地情况放大1.37倍。

(2)倾角为45°双力偶点源

图10为双力偶点源倾角45°，有无盆地情况观测点的水平峰值加速度。从以上结果可以看到，山谷表面的水平峰值加速度(PGA)多数大于两侧的山峰；当盆地存在时，山谷表面的水平PGA由盆地左侧边缘(观测点Rc)向中心方向呈现急剧放大的趋势，到达最大值后，较缓向盆地右侧边缘(观测点Rd)衰减。由于盆地的存在，使山谷表面各点PGA的平均值较无盆地情况放大1.66倍。

(3)倾角为80°双力偶点源

图11为双力偶点源倾角80°，有无盆地情况观测点的水平峰值加速度。从以上结果可以看到，山谷表面的水平峰值加速度(PGA)多数大于两侧的山峰；当盆地存在时，山谷表面的水平PGA由盆地左侧边缘(观测点Rc)向中心方向呈现急剧放大的趋势，到达最大值后，较缓向盆地右侧边缘(观测点Rd)衰减。由于盆地的存在，使山谷表面各点PGA的平均值较无盆地情况放大2.25倍。

结果表明：当震源倾角接近45°时，水平PGA的值最大，随着倾角的增大或减小，水平PGA的值也逐渐减小；盆地对地震动的放大作用随着震源倾角的增大而增大。

实施例2展示了本发明所开发的网格分级方法模拟复杂场地条件地震动的灵活性。

结论：

本发明方法具有良好的计算精度，且波动能够无障碍通过粗细网格的连接区域，不会发生任何异常现象。网格分级方法可高效、灵活进行复杂场地条件地震动模拟。

Claims (2)

1.一种近地表地震动模拟的网格分级方法，其特征在于：该方法步骤如下：

(1)建立被研究块体：

将地球介质按照数值模拟要求分成区域I和区域II，区域I由不规则细网格组成，区域II由不规则粗网格组成；区域I内的被研究块体为m-n-o-p-q-r-s-t-m，区域II内的被研究块体为1-2-3-4-a-5-6-7-1；区域I和区域II连接处的被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a，是通过连接粗网格、细网格内线段构成；辅助计算网格是虚线形成的四边形网格；

(2)建立被研究块体的动力平衡方程：

被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a的动力平衡方程为：

式中：ρ为被研究块体的密度；

ü_x和ü_z分别为被研究块体沿x轴和z轴方向的加速度；

σ_xx、σ_zz和σ_xz是被研究块体围线上的应力分量；

将式(1)和式(2)进一步改写为

$M_L{\stackrel{..}{u}}_{\mathrm{Lx}}=\underset{\mathrm{ii}=1}{\overset{2(n_{\mathrm{ra}}+1)}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sii}}({\mathrm{\σ}}_{\mathrm{xx}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zx}}\mathrm{dx})+\underset{\mathrm{jj}=1}{\overset{{2n}_{\mathrm{co}}}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sjj}}({\mathrm{\σ}}_{\mathrm{xx}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zx}}\mathrm{dx})---\left(3\right)$

$M_L{\stackrel{..}{u}}_{\mathrm{Lz}}=\underset{\mathrm{ii}=1}{\overset{2(n_{\mathrm{ra}}+1)}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sii}}({\mathrm{\σ}}_{\mathrm{xz}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zz}}\mathrm{dx})+\underset{\mathrm{jj}=1}{\overset{{2n}_{\mathrm{co}}}{\mathrm{\Σ}}}{\∫}_{\mathrm{Sjj}}({\mathrm{\σ}}_{\mathrm{xz}}\mathrm{dz}-{\mathrm{\σ}}_{\mathrm{zz}}\mathrm{dx})---\left(4\right)$

式中：M_L为被研究块体的质量；

n_co为节点L周围粗网格的数目；

n_fi为第j1个节点周围的细网格总数；

n_ra为被研究块体a-b-c-d-e-f-g-h-i-j-k-l-a内部所有节点的总数；

n_ra同时也代表与一个粗网格匹配的细网格的数目，n_ra＝1,3,5,...；

Sii和Sjj分别代表被研究块体L周围在细网格内和粗网格内的围线段；

(3)计算被研究块体围线上的应力：

被研究块体围线上的应力为：

σ_ij＝λδ_iju_k,k+μ(u_j,i+u_i,j)    (5)

式中：λ和μ为介质的拉梅系数，δ_ij为Kronecker张量，u_k,k＝u_x,x+u_z,z；

式(6)右侧位移的空间导数为：

${\∂u}_{\mathrm{xi}}/\∂x=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}[(1-\mathrm{\η})z_2+(\mathrm{\η}-\mathrm{\ξ})z_3+(\mathrm{\ξ}-1)z_4]u_{x1}+[(\mathrm{\η}-1)z_1+(1+\mathrm{\ξ})z_3+(-\mathrm{\η}-\mathrm{\ξ})z_4]u_{x2}+\\ [(\mathrm{\ξ}-\mathrm{\η})z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\η})z_4]u_{x3}+[(1-\mathrm{\ξ})z_1+(\mathrm{\ξ}+\mathrm{\η})z_2+(-1-\mathrm{\η})z_3]u_{x4}\end{array}\right\}---\left(6\right)$

${\∂u}_{\mathrm{zi}}/\∂z=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}-[(1-\mathrm{\η})x_2+(\mathrm{\η}-\mathrm{\ξ})x_3+(\mathrm{\ξ}-1)x_4]u_{z1}-[(\mathrm{\η}-1)x_1+(1+\mathrm{\ξ})x_3+(-\mathrm{\η}-\mathrm{\ξ})x_4]u_{z2}-\\ [(\mathrm{\ξ}-\mathrm{\η})x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\η})x_4]u_{z3}-[(1-\mathrm{\ξ})x_1+(\mathrm{\ξ}+\mathrm{\η})x_2+(-1-\mathrm{\η})x_3]u_{z4}\end{array}\right\}---\left(7\right)$

${\∂u}_{\mathrm{xi}}/\∂z=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}-[(1-\mathrm{\η})x_2+(\mathrm{\η}-\mathrm{\ξ})x_3+(\mathrm{\ξ}-1)x_4]u_{x1}-[(\mathrm{\η}-1)x_1+(1+\mathrm{\ξ})x_3+(-\mathrm{\η}-\mathrm{\ξ})x_4]u_{x2}-\\ [(\mathrm{\ξ}-\mathrm{\η})x_1+(-1-\mathrm{\ξ})x_2+(1+\mathrm{\η})x_4]u_{x3}-[(1-\mathrm{\ξ})x_1+(\mathrm{\ξ}+\mathrm{\η})x_2+(-1-\mathrm{\η})x_3]u_{x4}\end{array}\right\}---\left(8\right)$

${\∂u}_{\mathrm{zi}}/\∂x=\frac{1}{8\left|J_i\right|}\left\{\begin{array}{c}[(1-\mathrm{\η})z_2+(\mathrm{\η}-\mathrm{\ξ})z_3+(\mathrm{\ξ}-1)z_4]u_{z1}+[(\mathrm{\η}-1)z_1+(1+\mathrm{\ξ})z_3+(-\mathrm{\η}-\mathrm{\ξ})z_4]u_{z2}+\\ [(\mathrm{\ξ}-\mathrm{\η})z_1+(-1-\mathrm{\ξ})z_2+(1+\mathrm{\η})z_4]u_{z3}+[(1-\mathrm{\ξ})z_1+(\mathrm{\ξ}+\mathrm{\η})z_2+(-1-\mathrm{\η})z_3]u_{z4}\end{array}\right\}---\left(9\right)$

式中：u_xj和u_zj(j＝1,2,3,4)为四边形网格四个角点的位移分量；

x_j和z_j为四边形网格四个角点的坐标；

ξ和η为映射坐标系下的坐标变量；

(4)位移的线性插值：

①一个粗网格匹配三个细网格的情况(粗网格边上内部节点序号为1,2)：

u_1x＝2u_Kx/3+u_Lx/3；u_1y＝2u_Ky/3+u_Ly/3

u_2x＝u_Kx/3+2u_Lx/3；u_2y＝u_Ky/3+2u_Ly/3    (10)

②一个粗网格匹配五个细网格的情况(粗网格边上内部节点序号为1,2,3,4)：

u_1x＝4u_Kx/5+u_Lx/5；u_1y＝4u_Ky/5+u_Ly/5

u_2x＝3u_Kx/5+2u_Lx/5；u_2y＝3u_Ky/5+2u_Ly/5

u_3x＝2u_Kx/5+3u_Lx/5；u_3y＝2u_Ky/5+3u_Ly/5

u_4x＝u_Kx/5+4u_Lx/5；u_4y＝u_Ky/5+4u_Ly/5    (11)

③一个粗网格匹配七个细网格的情况(粗网格边上内部节点序号为1,2,3,4,5,6)：

u_1x＝6u_Kx/7+u_Lx/7；u_1y＝6u_Ky/7+u_Ly/7

u_2x＝5u_Kx/7+2u_Lx/7；u_2y＝5u_Ky/7+2u_Ly/7

u_3x＝4u_Kx/7+3u_Lx/7；u_3y＝4u_Ky/7+3u_Ly/7

u_4x＝3u_Kx/7+4u_Lx/7；u_4y＝3u_Ky/7+4u_Ly/7

u_5x＝2u_Kx/7+5u_Lx/7；u_5y＝2u_Ky/7+5u_Ly/7

u_6x＝u_Kx/7+6u_Lx/7；u_6y＝u_Ky/7+6u_Ly/7    (12)

④一个粗网格匹配九个细网格的情况(粗网格边上内部节点序号为1,2,3,4,5,6,7,8)：

u_1x＝8u_Kx/9+u_Lx/9；u_1y＝8u_Ky/9+u_Ly/9

u_2x＝7u_Kx/9+2u_Lx/9；u_2y＝7u_Ky/9+2u_Ly/9

u_3x＝6u_Kx/9+3u_Lx/9；u_3y＝6u_Ky/9+3u_Ly/9

u_4x＝5u_Kx/9+4u_Lx/9；u_4y＝5u_Ky/9+4u_Ly/9

u_5x＝4u_Kx/9+5u_Lx/9；u_5y＝4u_Ky/9+5u_Ly/9

u_6x＝3u_Kx/9+6u_Lx/9；u_6y＝3u_Ky/9+6u_Ly/9

u_7x＝2u_Kx/9+7u_Lx/9；u_7y＝2u_Ky/9+7u_Ly/9

u_8x＝u_Kx/9+8u_Lx/9；u_8y＝u_Ky/9+8u_Ly/9    (13)；

(5)网格分级方法的数值实现：

被研究块体K和被研究块体L的合力为：

F_xK＝-f_xk7-f_xk6-f_xk5+f_xk1+f_xl1+f_xl2

F_zK＝-f_zk7-f_zk6-f_zk5+f_zk1+f_zl1+f_zl2

F_xL＝-f_xl3+f_xl1+f_xk1+f_xk2+f_xk3+f_xk4

F_zL＝-f_zl3+f_zl1+f_zk1+f_zk2+f_zk3+f_zk4    (14)

依据t时刻粗细网格内各条围线段上的应力分量和积分可得t时刻被研究块体各条围线段上的内力和再由式(14)可将这些内力分配给相应的被研究块体，可求得t时刻粗细网格内各被研究块体的合力和进而由式(3)和式(4)给出各被研究块体在t时刻的加速度分量和

通过时间积分，可求得t+Δt/2时刻各被研究块体的速度分量和

进一步通过时间积分，可求得t+Δt时刻各被研究块体的位移分量和针对一个粗网格匹配三个细网格的情况，由式(10)的插值关系可给出t+Δt时刻粗网格边上节点的位移；将所有被研究块体的位移分量代入式(6)-(9)，可求得t+Δt时刻位移的空间导数，将其代入式(5)可以得到t+Δt时刻被研究块体各条围线段上的应力分量和

根据以上递归方法，实现粗细网格连接处被研究块体的应力场、位移场和加速度场从t时刻至t+Δt时刻的计算，速度场从t-Δt/2时刻至t+Δt/2时刻的计算。

2.根据权利要求1所述的近地表地震动模拟的网格分级方法，其特征在于：用不同步距的不规则网格对计算模型进行离散，且粗细网格间无须共用节点。