三维饱和沉积盆地地震动模拟方法

摘要

本发明提供一种三维饱和半空间中沉积场地地震动模拟方法，该方法基于间接边界元方法，推导了全空间饱和斜面圆盘荷载格林函数，解决了求解边界的奇异性问题；对于饱和两相模型模拟固液动力耦合效应、孔压及流量的变化规律；建立三维饱和场地模型，基于Biot饱和多介质理论推导饱和格林函数，根据边界条件建立方程求解，最后进行频时域地震响应模拟分析以期实现未来地震动预测。本发明效果是为饱和不均匀介质工程波动问题求解提供一种高效算法和实用程序，为沿海复杂饱和场地区域地震区划、结构抗震设防提供理论依据，为实际复杂工程的模拟带来显著的经济效益。

主权项

一种三维饱和半空间中沉积场地地震动模拟方法，该方法是在建立三维饱和沉积盆地模型基础上，通过由饱和土中的土体骨架和流体之间的相对关系，推导饱和格林函数，采用间接边界元法解决饱和局部场地的地震动模拟，该方法包括有以下步骤：量化(一)建立三维饱和场地模型：考虑平面波P或S和Rayleigh波的入射，基于单层位势理论，在沉积盆地模型边界表面施加三个正交方向的虚拟荷载以构造散射波场，采用边界元方法实施中，需首先将边界表面离散为三角形或四边形单元，为便于处理所述三角形或四边形单元格林函数，采用斜面圆盘均布荷载近似覆盖在所述边界表面的离散单元上，即分所述离散单元上虚拟荷载对自身单元作用为本单元，所述离散单元上虚拟荷载对自身单元之外的离散单元作用为非本单元，推导出饱和介质中三个方向的力F_x，F_y，F_z集中力和流量元(Q)作用下的格林函数，进而根据所述边界表面的透水和不透水边界条件构造矩阵方程组：[H][φ]＝[B]，求解得到虚拟荷载密度：[φ]＝[H]^‑1[B]其中，[H]为格林函数影响矩阵，[B]为边界条件解最后将边界表面上各个单元上的虚拟荷载的综合作用而得散射波场和不含局部不均匀场地时弹性波入射下的半空间自由波场叠加即得到总波场；由于虚拟荷载直接施加在边界表面上，能够解决多种形状的饱和局部场地；(二)三维饱和格林函数推导：基于Biot饱和两相介质理论，推导饱和介质中三个方向的力F_x，F_y，F_z集中力和流量元(Q)作用下的格林函数，格林函数解如下：1非本单元格林函数解1.1集中力作用类似于三维弹性半空间问题，在均匀各向同性弹性固体介质中，位移格林函数表达为：G_ij(x,ξ)＝[f₂δ_ij+(f₁‑f₂)γ_iγ_j]/4πμr       (1)其中，γ_j＝(x_j‑ξ_j)/r，r表示波源点(ξ₁,ξ₂,ξ₃)与接收点(x₁,x₂,x₃)之间的距离，即r²＝(x₁‑ξ₁)²+(x₂‑ξ₂)²+(x₃‑ξ₃)²，δ_ij表示Kronecker delta，μ表示拉梅常数；k_sv＝ω/β，k_α1＝ω/α₁，k_α2＝ω/α₂分别表示S，P₁(快波)，P₂(慢波)波数；β，α₁，α₂表示其各自对应波速，f₁和f₂定义为：f₁＝A₁(β²/α₁²)[1‑2i/(k_α1r)‑2/(k_α1r)²]exp(‑ik_α1r)+A₂(β²/α₂²)[1‑2i/(k_α2r)‑2/(k_α2r)²]exp(‑ik_α2r)+A₃[2i/(k_βr)+2/(k_βr)²]exp(‑ik_βr)              (2a)f₂＝A₁(β²/α₁²)[i/(k_α1r)+1/(k_α1r)²]exp(‑ik_α1r)+A₂(β²/α₂²)[i/(k_α2r)+1/(k_α2r)²]exp(‑ik_α2r)+A₃[1‑i/(k_βr)‑1/(k_βr)²]exp(‑ik_βr)             (2b)式中系数A₁＝(χ₂‑χ₃)/(χ₂‑χ₁)；A₂＝(χ₁‑χ₃)/(χ₁‑χ₂)；A₃＝1(下同)应力格林函数T_ij表达式如下：$\left\{\begin{array}{c}{T}_{xj}=n_1{\sigma }_{xx}^j+n_2{\sigma }_{yx}^j+n_3{\sigma }_{zx}^j\\ T_{yj}=n_1{\sigma }_{xy}^j+n_2{\sigma }_{yy}^j+n_3{\sigma }_{zy}^j\\ T_{zj}=n_1{\sigma }_{xz}^j+n_2{\sigma }_{yz}^j+n_3{\sigma }_{zz}^j\end{array}\right.---\left(3\right)$其中：σ为单元应力；n为单元法向量同理，另给出流体相对位移W_ij、孔隙水压力P_j格林函数表达式如下：$W_{ij}(x,\xi )=[f_2^w{\delta }_{ij}+(f_1^w-f_2^w){\gamma }_i{\gamma }_j]/4\pi \mu r---\left(4\right)$$P_j(x,\xi )={\mathrm{M\gamma }}_j[k_{\alpha 1}^2(\alpha +{\chi }_1)A_1\exp (-{\mathrm{ik}}_{\alpha 1}r)({\mathrm{ik}}_{\alpha 1}r+1)/r+k_{\alpha 2}^2(\alpha +{\chi }_2)A_2\exp (-{\mathrm{ik}}_{\alpha 2}r)({\mathrm{ik}}_{\alpha 2}r+1)/r]/4{{\mathrm{\pi \mu rk}}_{\beta }}^2---\left(5\right)$其中和为：$\begin{array}{c}{f}_1^w={\chi }_1{A}_1({\beta }^2/{{\alpha }_1}^2)[1-2i/\left(k_{\alpha 1}r\right)-2/{\left(k_{\alpha 1}r\right)}^2]\exp (-{\mathrm{ik}}_{\alpha 1}r)+{\chi }_2{A}_2({\beta }^2/{{\alpha }_2}^2)[1-2i/\left(k_{\alpha 2}r\right)-2/{\left(k_{\alpha 2}r\right)}^2]\exp \left(-{\mathrm{ik}}_{\alpha 2}r\right)\\ +{\chi }_3{A}_3\left[2i/\left(k_{\beta }r\right)+2/{\left(k_{\beta }r\right)}^2\right]\exp \left(-{\mathrm{ik}}_{\beta }r\right)\end{array}---\left(6a\right)$$\begin{array}{c}{f}_2^w={\chi }_1{A}_1\left({\beta }^2/{{\alpha }_1}^2\right)\left[i/\left(k_{\alpha 1}r\right)+1/{\left(k_{\alpha 1}r\right)}^2\right]\exp \left(-{\mathrm{ik}}_{\alpha 1}r\right)\\ +{\chi }_2{A}_2\left({\beta }^2/{{\alpha }_2}^2\right)\left[i/\left(k_{\alpha 2}r\right)+1/{\left(k_{\alpha 2}r\right)}^2\right]\exp \left(-{\mathrm{ik}}_{\alpha 2}r\right)\end{array}$$+{\chi }_3{A}_3[1-i/\left(k_{sv}r\right)-1/{\left(k_{sv}r\right)}^2]\exp (-{\mathrm{ik}}_{sv}r)---\left(6b\right)$1.2流量源作用G_iF(x,ξ)＝‑iBγ_i[exp(‑ik_α1r)(ik_α1r+1)/r‑exp(‑ik_α2r)(ik_α2r+1)/r]/4πωr   (7)W_iF(x,ξ)＝‑iBγ_i[χ₁exp(‑ik_α1r)(ik_α1r+1)/r‑χ₂exp(‑ik_α2r)(ik_α2r+1)/r]/4πωr  (8)$P_F(x,\xi )=iBM[k_{\alpha 1}^2(\alpha +{\chi }_1)\exp (-{\mathrm{ik}}_{\alpha 1}r)/r-k_{\alpha 2}^2(\alpha +{\chi }_2)A_2\exp (-{\mathrm{ik}}_{\alpha 2}r)/r]/4\pi \omega ---\left(9\right)$其中：B＝1/(χ₁‑χ₂),系数2本单元格林函数解2.1集中力作用：(1)固体骨架位移由G_ij(x,ξ)＝[f₂δ_ij+(f₁‑f₂)γ_iγ_j]/4πμr进行推导，其过程以水平方向力F_x作用产生的x方向位移为例，为了方便，以离散单元中心为原点建立局部极坐标系o‑θ，将水平力F_x分解为法向F_z'和切向F_x'，即F_z'＝F_xn_x，在法向分解力F_z'作用下：G_z'z'(x_n,ξ_l)＝[f₂+(f₁‑f₂)γ_z'γ_z']/4πμr＝f₂/4πμr；γ_z'＝0    (10)则法向分解力F_z'作用下产生的x方向位移为：在切向分解力F_x'作用下：G_x'x'(x_n,ξ_l)＝[f₂+(f₁‑f₂)γ_x'γ_x']/4πμr＝[f₁(1‑cos2θ)+f₂(1+cos2θ)]/8πμr；γ_x'＝‑cosθ                             (12)则切向分解力F_x'作用下产生的x方向位移为：由公式(11)、(13)得出，水平方向力F_x作用在整体坐标系下产生的x方向的总位移为：$u_{xx}={\int }_0^R[(f_2+f_1)+(f_2-f_1)n_x{n}_x/4\mu dr---\left(14\right)$则集中力作用(F_x，F_y，F_z)下的位移格林函数积分结果为：$g_{ij}(x_n,{\xi }_l)={\int }_{{\mathrm{\Delta S}}_l}{G}_{ij}(x_n,{\xi }_l){\mathrm{dS}}_{\xi }=[(F_2+F_1){\delta }_{ij}+(F_2-F_1)n_i{n}_j]/4\mu ;i,j=x,y,z---\left(15\right)$其中，ΔS_l为离散单元面，δ_ij为Kronecker delta(当i＝j时，δ_ij＝1；当i≠j时，δ_ij＝0)，n_j为离散单元法向量；式(15)中F₁、F₂分别表示式(2)中f₁、f₂从0到R积分，R为离散单元半径，积分结果表达式如下：$\begin{array}{c}{F}_1=\left({\beta }^2/{{\alpha }_1}^2\right)\exp \left(-{\mathrm{ik}}_{\alpha 1}R-1\right)A_{1}i/q_1+\left({\beta }^2/{{\alpha }_2}^2\right)\exp \left(-{\mathrm{ik}}_{\alpha 2}R-1\right)A_{2}i/k_{\alpha 2}\\ +2\left[A_1\exp \left(-{\mathrm{ik}}_{\alpha 1}R\right)+A_2\exp \left(-{\mathrm{ik}}_{\alpha 2}R\right)-A_3\exp \left(-{\mathrm{ik}}_{sv}R\right)\right]/\left({\mathrm{Rk}}_{sv}^2\right)\\ +2\left({\mathrm{ik}}_{\alpha 1}{A}_1+{\mathrm{ik}}_{\alpha 2}{A}_2-{\mathrm{ik}}_{sv}{A}_3\right)/k_{sv}^2\end{array}---\left(16a\right)$$\begin{array}{c}{F}_2=A_{3}i\left[\exp \left(-{\mathrm{ik}}_{sv}R\right)-1\right]/k_{sv}-[A_1\exp \left(-{\mathrm{ik}}_{\alpha 1}R\right)+A_2\exp \left(-{\mathrm{ik}}_{\alpha 2}R\right)\\ -A_3\exp \left(-{\mathrm{ik}}_{sv}R\right)]/\left({\mathrm{Rk}}_{sv}^2\right)-\left({\mathrm{ik}}_{\alpha 1}{A}_1+{\mathrm{ik}}_{\alpha 2}{A}_2-{\mathrm{ik}}_{sv}{A}_3\right)/k_{sv}^2;\end{array}---\left(16b\right)$(2)流体相对位移$w_{ij}(x_n,{\xi }_l)={\int }_{{\mathrm{\Delta S}}_l}{W}_{ij}(x_n,{\xi }_l){\mathrm{dS}}_{\xi }=[(F_2^w+F_1^w){\delta }_{ij}+(F_2^w-F_1^w)n_i{n}_j]/4\mu ;i,j=x,y,z---\left(17\right)$其中：$\begin{array}{c}{F}_1^w=({\beta }^2/{{\alpha }_1}^2)\exp (-{\mathrm{ik}}_{\alpha 1}R-1){\chi }_1{A}_{1}i/k_{\alpha 1}+({\beta }^2/{{\alpha }_2}^2)\exp \left(-{\mathrm{ik}}_{\alpha 2}R-1\right){\chi }_2{A}_{2}i/k_{\alpha 2}+2\left({\mathrm{ik}}_{\alpha 1}{\chi }_1{A}_1+{\mathrm{ik}}_{\alpha 2}{\chi }_2{A}_2-{\mathrm{ik}}_{\beta }{\chi }_3{A}_3\right)/k_{\beta }^2\\ +2\left[{\chi }_1{A}_1\exp \left(-{\mathrm{ik}}_{\alpha 1}R\right)+{\chi }_2{A}_2\exp \left(-{\mathrm{ik}}_{\alpha 2}R\right)-{\chi }_3{A}_3\exp \left(-{\mathrm{ik}}_{\beta }R\right)\right]/\left(k_{\beta }^{2}R\right)\end{array}---\left(18a\right)$$\begin{array}{c}{F}_2^w={\chi }_3{A}_{3}i\left[\exp \left(-{\mathrm{ik}}_{sv}R\right)-1\right]/k_{sv}-\left[{\chi }_1{A}_1\exp \left(-{\mathrm{ik}}_{\alpha 1}R\right)+{\chi }_2{A}_2\exp \left(-{\mathrm{ik}}_{\alpha 2}R\right)-{\chi }_3{A}_3\exp \left(-{\mathrm{ik}}_{\beta }R\right)\right]/\left(k_{\beta }^{2}R\right)\\ -\left({\mathrm{ik}}_{\alpha 1}{\chi }_1{A}_1+{\mathrm{ik}}_{\alpha 2}{\chi }_2{A}_2-{\mathrm{ik}}_{\beta }{\chi }_3{A}_3\right)/k_{\beta }^2;\end{array}---\left(18b\right)$(3)孔隙水压力由对称性原理可知：$p_j(x_n,{\xi }_l)={\int }_{{\mathrm{\Delta S}}_l}{P}_j(x_n,{\xi }_l){\mathrm{dS}}_{\xi }=0;---\left(19\right)$(4)牵引力$t_{ij}(x_n,{\xi }_l)=0.5{\delta }_{ij}{\delta }_{nl}+{\int }_{{\mathrm{\Delta S}}_l}{T}_{ij}(x_n,{\xi }_l){\mathrm{dS}}_{\xi }---\left(20\right)$当x_n＝ξ_l时，被积函数T_ij具有奇异性其柯西主值为零，即当离散单元为光滑圆面时，牵引力格林函数张量T_ij(x_n,ξ_l)为零，则：t_ij(x_n,ξ_l)＝0.5δ_ij；        (21)2.2流量源作用(1)固体骨架位移：由式(5)、(7)可知，P_j和G_iF之间的关系为：P_j＝iωG_iF由对称性可知，当x_n＝ξ_l时：$g_{iF}(x_n,{\xi }_l)={\int }_{{\mathrm{\Delta S}}_l}{G}_{iF}(x_n,{\xi }_l){\mathrm{dS}}_{\xi }=0---\left(22\right)$(2)流体相对位移$w_{iF}(x_n,{\xi }_l)=\frac{1}{2}\sqrt{\frac{\rho }{\mu }}\frac{i}{\delta }{\delta }_{nl}+{\int }_{{\mathrm{\Delta S}}_l}{W}_{iF}(x_n,{\xi }_l){\mathrm{dS}}_{\xi }---\left(23\right)$其中，当x_n＝ξ_l时，W_iF(x_n,ξ_l)为零，则：$w_{iF}(x_n,{\xi }_l)=\frac{1}{2}\sqrt{\frac{\rho }{\mu }}\frac{i}{\delta }---\left(24\right)$(3)孔隙水压力$p_F(x_n,{\xi }_l)={\int }_{{\mathrm{\Delta S}}_l}{P}_F(x_n,{\xi }_l){\mathrm{dS}}_{\xi }=-MB[k_{\alpha 1}(\alpha +{\chi }_1)\left(\exp \right(-{\mathrm{ik}}_{\alpha 1}R)-1)+k_{\alpha 2}(\alpha +{\chi }_2)\left(\exp \right(-{\mathrm{ik}}_{\alpha 2}R)-1)]/2\omega ---\left(25\right)$(4)牵引力$t_{iF}(x_n,{\xi }_l)={\int }_{{\mathrm{\Delta S}}_l}{T}_{iF}(x_n,{\xi }_l){\mathrm{dS}}_{\xi }={\int }_{{\mathrm{\Delta S}}_l}{n}_i{\sigma }_{iF}{\mathrm{dS}}_{\xi }=F_{iF}---\left(26\right)$其中：所述散射波场通过施加在沉积盆地模型边界上的虚拟波源与格林函数求得，具体表达为：$u_i^{\left(s\right)}\left(x\right)={\int }_S[{\phi }_j\left(\xi \right)G_{ij}(x,\xi )+{\phi }_F\left(\xi \right)G_{iF}(x,\xi )]{\mathrm{dS}}_{\xi }---\left(28\right)$$t_i^{\left(s\right)}\left(x\right)={\int }_S[{\phi }_j\left(\xi \right)T_{ij}(x,\xi )+{\phi }_F\left(\xi \right)T_{iF}(x,\xi )]{\mathrm{dS}}_{\xi }---\left(29\right)$$w_i^{\left(s\right)}\left(x\right)={\int }_S[{\phi }_j\left(\xi \right)W_{ij}(x,\xi )+{\phi }_F\left(\xi \right)W_{iF}(x,\xi )]{\mathrm{dS}}_{\xi }---\left(30\right)$$p^{\left(s\right)}\left(x\right)={\int }_S[{\phi }_j\left(\xi \right)P_j(x,\xi )+{\phi }_F\left(\xi \right)P_F(x,\xi )]{\mathrm{dS}}_{\xi };i,j=1,2,3---\left(31\right)$式中，φ_j,φ_F分别为模型边界固体骨架j方向虚拟波源和流量源；G_ij,W_ij,T_ij和P_j依次为集中力作用下的固体骨架位移函数、流体相对位移格林函数、牵引力格林函数和孔压格林函数，下标j＝1,2,3分别对应虚拟波源在直角坐标系下三个正交方向x,y,z；G_iF,W_iF,T_iF和P_F依次为流量源作用下的固体骨架位移函数、流体相对位移格林函数、牵引力格林函数和孔压格林函数；(三)总波场计算根据场地边界条件，即自由表面零应力和耦合边界位移、应力连续条件构造矩阵方程组：[H][φ]＝[B]，并求解虚拟荷载密度[φ]＝[H]^‑1[B]，将所求得的虚拟波源和流量源用于求解散射波场进而叠加自由波场得总波场，即其中自由波场表达式如下：$\begin{array}{c}{u}_x^{\left(f\right)}=-{\mathrm{ik}}_{\alpha 1}\left({\phi }^{\left(i\right)}+{\phi }_1^{\left(r\right)}+{\phi }_2^{\left(r\right)}\right){\mathrm{sin\theta }}_{\alpha 1}-{\mathrm{ik}}_{\beta }{\psi }^{\left(r\right)}{\mathrm{cos\theta }}_{\beta };u_y^{\left(f\right)}=0;\\ u_z^{\left(f\right)}={\mathrm{ik}}_{\alpha 1}\left({\phi }^{\left(i\right)}-{\phi }_1^{\left(r\right)}\right){\mathrm{cos\theta }}_{\alpha 1}-{\mathrm{ik}}_{\alpha 2}{\phi }_2^{\left(r\right)}{\mathrm{cos\theta }}_{\alpha 2}+{\mathrm{ik}}_{\beta }{\psi }^{\left(r\right)}{\mathrm{sin\theta }}_{\beta }\end{array}---\left(32\right)$式中，φ和ψ为波势函数；k_α1、k_α2和k_β为波数；θ为入射波和反射波在直角坐标系下与垂直方向的夹角；所述进行饱和格林函数的推导，采用间接边界元法解决饱和局部场地的地震动模拟，得出饱和局部场地的地震动模拟量化条件。