[发明专利]一种近地表地震动模拟的网格分级方法

摘要

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

主权项

一种近地表地震动模拟的网格分级方法，其特征在于：该方法步骤如下：(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/3u_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/5u_2x＝3u_Kx/5+2u_Lx/5；u_2y＝3u_Ky/5+2u_Ly/5u_3x＝2u_Kx/5+3u_Lx/5；u_3y＝2u_Ky/5+3u_Ly/5u_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/7u_2x＝5u_Kx/7+2u_Lx/7；u_2y＝5u_Ky/7+2u_Ly/7u_3x＝4u_Kx/7+3u_Lx/7；u_3y＝4u_Ky/7+3u_Ly/7u_4x＝3u_Kx/7+4u_Lx/7；u_4y＝3u_Ky/7+4u_Ly/7u_5x＝2u_Kx/7+5u_Lx/7；u_5y＝2u_Ky/7+5u_Ly/7u_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/9u_2x＝7u_Kx/9+2u_Lx/9；u_2y＝7u_Ky/9+2u_Ly/9u_3x＝6u_Kx/9+3u_Lx/9；u_3y＝6u_Ky/9+3u_Ly/9u_4x＝5u_Kx/9+4u_Lx/9；u_4y＝5u_Ky/9+4u_Ly/9u_5x＝4u_Kx/9+5u_Lx/9；u_5y＝4u_Ky/9+5u_Ly/9u_6x＝3u_Kx/9+6u_Lx/9；u_6y＝3u_Ky/9+6u_Ly/9u_7x＝2u_Kx/9+7u_Lx/9；u_7y＝2u_Ky/9+7u_Ly/9u_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_xl2F_zK＝‑f_zk7‑f_zk6‑f_zk5+f_zk1+f_zl1+f_zl2F_xL＝‑f_xl3+f_xl1+f_xk1+f_xk2+f_xk3+f_xk4F_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时刻的计算。