陈少林
, 2) , 程书林, 柯小飞南京航空航天大学土木工程系, 南京 210016
A UNIFIED COMPUTATIONAL FRAMEWORK FOR FLUID-SOLID COUPLING IN MARINE EARTHQUAKE ENGINEERING: IRREGULAR INTERFACE CASE1) Chen Shaolin
, 2) , Cheng Shulin, Ke XiaofeiDepartment of Civil Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China
通讯作者: 2) 陈少林, 教授, 主要研究方向：地震工程, E-mail:
icmcsl@nuaa.edu.cn 收稿日期: 2019-06-17
网络出版日期: 2019-09-18
基金资助: 1) 国家自然科学基金资助项目 .51278260
Received: 2019-06-17
Online: 2019-09-18
作者简介 About authors
摘要 海底地震动场及海洋声场的模拟中,需要考虑复杂海床介质及海底地形的影响,涉及到海水、饱和海床、弹性基岩之间的相互耦合.传统的方法分别采用声波方程描述理想流体、Biot方程描述饱和海床、弹性波方程描述基岩,分别进行空间离散和界面耦合, 十分不便.本文基于理想流体、固体分别为饱和多孔介质的特殊情形(孔隙率分别为1和0),由饱和多孔介质的Biot方程可退化得到理想流体的声波方程和固体的弹性波方程.然后, 以饱和多孔介质方程为基础, 经集中质量有限元离散,严格考虑不同孔隙率的饱和多孔介质在不规则界面的耦合条件,通过求解法向和切向界面力的途径,建立了不同孔隙率的饱和多孔介质耦合情形的求解方法,将流体、固体、饱和多孔介质间的耦合问题纳入到统一计算框架,并编制了相应的三维并行分析程序.考虑海水--弹性基岩、海水--饱和海床--弹性基岩体系中凹陷地形情形,采用本文提出的统一计算框架, 结合透射边界条件,分析了P波入射时的动力反应, 并通过结果是否满足界面条件,验证了该统一计算框架的有效性以及并行计算的可行性.
关键词： 流固耦合 ;
饱和多孔介质 ;
海洋地震工程 ;
不规则海底地形 ;
透射边界 ;
并行计算 Abstract The simulation of seismic wavefield at seafloor and ocean acoustic field, the influence of complex seabed media and seabed topography needs to be considered, involving the coupling between seawater, saturated seabed, elastic bedrock and structure. That means, we target simulation where several types of equations are involved such as fluid, solid and saturated porous media equation. The conventional method for this fluid-solid-saturated porous media interaction problem is to use exsisting solvers of different equations and coupling method, which needs data mapping, communication and coupling algorithm between different solvers. Here, we present an alternative method, in which the coulping between different solvers is avoided. In fact, when porosity equals to one and zero,the saturated porous media is reduced to fluid and solid respectively, so we can use the porous media equation to describe the ideal fluid and solid, and the coupling between porous media,solid and fluid turns to the coupling between porous media with different porosity. Based on this idea, firstly the Biot's equations are approximated by Galerkin scheme and the explicit lumped-mass FEM is chosen, that are well suited to parallel computation. Then considering the conditions of coupling on the irregular interface between porous media with different porosity,by solving the normal and tangential interface forces, the coupled algorithm is derived, which is proved to be suitable for the coupling between fluid,solid and saturated porous media. Thus, the coupling problem between fluid, solid and saturated porous media can be brought into a unified framework, in which only the solver of saturated porous media is used. The three-dimensional parallel code for this proposed method is programed. Considering the situation of sag terrain in water-bedrock, water-saturated seabed-bedrock system, the unified framework proposed in this paper is combined with the transmission boundary conditions to analyze the dynamic response of P wave incident, and the unified framework are verified by the results meeting the interface conditions.
Keywords： fluid-solid coupling ;
saturated porous medium ;
marine earthquake engineering ;
explicit lumped-mass finite element method ;
transmitting boundary ;
parallel computation PDF (8787KB) 元数据 多维度评价 相关文章 导出 EndNote |
Ris |
Bibtex 收藏本文 本文引用格式 陈少林, 程书林, 柯小飞. 海洋地震工程流固耦合问题统一计算框架------不规则界面情形
1) .
力学学报 [J], 2019, 51(5): 1517-1529 DOI:
10.6052/0459-1879-19-156 Chen Shaolin, Cheng Shulin, Ke Xiaofei.
A UNIFIED COMPUTATIONAL FRAMEWORK FOR FLUID-SOLID COUPLING IN MARINE EARTHQUAKE ENGINEERING: IRREGULAR INTERFACE CASE1) .
Chinese Journal of Theoretical and Applied Mechanics [J], 2019, 51(5): 1517-1529 DOI:
10.6052/0459-1879-19-156 引言 海底地震动模拟中, 需要考虑震声(seismoacoustic)散射,若不适当计入流固边界条件, 将导致沿界面的面波误差较大. Nakamura等
[1 ] 通过有限差分方法, 考虑海底地形和海水层的影响,模拟了2009年Suruga Bay地震的地震波传播; 为了满足流固界面条件,在界面附近的方程采用二阶近似, 其余部分采用四阶近似,两者精度不一致. Okamoto和Takenaka
[2 ] 采用反射/透射矩阵方法,考虑流固界面条件,模拟了P-SV波入射情形二维不规则流固界面情形的震声散射,但该方法仅适用于坡度较缓地形, 对于陡峭地形, 容易导致计算失稳. Qian和Yamanaka
[3 ] 利用边界元法适于陡峭地形的特点,以及全局矩阵传播子在模拟多层介质中波动对内存需求较少的特点,对全局矩阵传播子进行扩展, 考虑流固界面条件,将边界元法与全局矩阵传播子方法结合,模拟了P-SV波入射情形二维不规则流固界面情形的震声散射.Komatitsch等
[4 ] 考虑海水与基岩界面法向速度连续和应力连续,建立平衡方程的耦合弱积分形式, 通过谱元离散和显式Newmark时间积分,得到海底地震波的谱元模拟方法.
海洋声学中, 也需考虑海底地形对声波传播的影响
[5 -6 ] ,一般采用抛物方程(PE)方法, 该方法可快速建立远距离声场,在海洋声学领域得到广泛应用, 但考虑海底不规则复杂地形时, 较为困难.Murphy等
[7 ] 采用有限元方法模拟海洋声场, 相比于抛物方程法,可以计算全波场及海底复杂地形的影响, 但效率较低,文中忽略了海底弹性介质的剪切效应.
在海洋物理勘探中, 对于含油气的页岩, 需要考虑为多孔介质.海洋工程结构地基液化失效问题,也需要考虑海床为多孔介质
[8 -11 ] .李伟华等
[11 ] 考虑理想流体和饱和土层的界面连续条件,建立了水与饱和场地耦合的显式有限元动力分析方法.
以上研究中, 海水、基岩、饱和多孔介质分别采用不同波动方程进行描述,再通过界面条件进行耦合求解, 十分不便
[12 -21 ] . 理论上,固体和流体介质均为饱和多孔介质的特殊情形, 孔隙率分别为0和1,上述耦合均可在饱和多孔介质理论体系进行描述, 如
图1 所示. 基于此,可从饱和多孔介质的一般情形出发,考虑不同孔隙率的饱和多孔介质之间的耦合,从而将上述所有耦合统一在同一计算框架, 避免不同模块之间的交互.在文献
[22 ] 中考虑了海域水平成层场地情形.本文将考虑不规则海底地形情形, 需要考虑界面的法向方向导数,若按文献
[22 ] 中的方法, 变得较为困难.这里采用求解法向和切向界面力的方法, 可较为方便地解决此问题.
图1 新窗口打开 |
下载原图ZIP |
生成PPT 图1介质关系示意图 Fig. 1Schematic diagram of media relations 采用有限元法进行海底地震动和海洋声场的模拟,需要采用人工边界模拟无限域的影响,一般采用透射边界
[23 -25 ] 和黏弹性边界
[26 -27 ] ,并分别以自由场
[28 ] 和等效力
[29 ] 作为输入.本文采用集中质量显式有限元方法结合透射边界条件,可避免求解大型线性方程组, 效率较高, 且易于实现并行计算,有望用于大型海底地震波以及海洋声场的模拟中. 通过两凹陷地形的算例,验证了本文统一计算框架的有效性.
1 基本理论 1.1 一般饱和多孔介质情形 饱和多孔介质的固相平衡方程
[30 ] (1) \begin{equation} \label{eq1} {\bf L}_{\rm s}^{\rm T} \boldsymbol \sigma'-(1- \beta ){\bf L}_{\rm w}^{\rm T} P + b(\dot{\boldsymbol U}- \dot{\boldsymbol u}) = (1-\beta )\rho _{\rm s} \ddot{\boldsymbol u} \end{equation}液相平衡方程
(2) \begin{equation} \label{eq2} -\beta {\bf L}_{\rm w}^{\rm T} P + b(\dot{\boldsymbol u}-\dot{\boldsymbol U}) = \beta \rho _{\rm w} \ddot{\boldsymbol U} \end{equation}相容方程(考虑初始孔压和初始体应变为零时)
(3) $$\tau =-\beta P = E_{\rm w} [\beta e^{\rm w} + (1-\beta )e^{\rm s}] $$其中, ${\bf L}_{\rm s}$和${\bf L}_{\rm w} $为微分算子矩阵
(4a) $${\bf L}_{\rm s} = \left[ {{\begin{array}{*{20}c} \partial/\partial x_1 & 0 & 0 \\ 0 & \partial/\partial x_2 & 0 \\ 0 & 0 & \partial/\partial x_3\\ \partial/\partial x_2 & \partial/\partial x_1 & 0 \\ 0 & \partial/\partial x_3 & \partial/\partial x_2 \\ \partial/\partial x_3 & 0 & \partial/\partial x_1\\ \end{array} }} \right] $$(4b) $${\bf L}_{\rm w} = \left(\partial/\partial x_1,\partial/x_2, \partial/\partial x_3\right) $$式中, $\boldsymbol \sigma'$为有效应力矢量, $\tau $为平均孔压,以拉为正; $P$为孔隙水压, 以压为正; $\boldsymbol U$, $\boldsymbol u$分别为液相和固相的位移矢量, $\dot{\boldsymbol U}$,$\dot{\boldsymbol u}$为速度, $\ddot{\boldsymbol U}$,$\ddot{\boldsymbol u}$为加速度; $\rho _{\rm s}$, $\rho _{\rm w}$分别为固相和液相的密度, $\beta $为孔隙率, $b = \beta ^2\mu_0/k_0 $, $k_0 $为流体渗透系数, $\mu _0 $为动黏度系数; $E_{\rm w}$为流体的体变模量, $e^{\rm s}$, $e^{\rm w}$分别为固相和液相的体 应变.Dirichlet边界条件
(5a) $$\boldsymbol u-\bar{\boldsymbol u} = 0 $$(5b) $$\boldsymbol U-\bar{\boldsymbol U} = 0 $$Neumann边界条件
(5c) $$\hat{\boldsymbol n}\boldsymbol \sigma-\bar{\boldsymbol \sigma} = 0$$(5d) $$(P-\bar{P})\boldsymbol n = 0 $$其中, $\bar{\boldsymbol u}$, $\bar{\boldsymbol U}$分别为在边界上给定的固相和液相位移. $\bar{\boldsymbol \sigma}$,$\bar{P}$分别为边界上固相平均应力和真实孔压的给定值; $\boldsymbol n$为沿边界外法线的方向矢量, $\hat{\boldsymbol n}$为由方向导数组成的矩阵, 其形式如下
(6a) $$\boldsymbol n = \left( n_x, n_y, n_z\right)^{\rm T} $$(6b) $$\hat{\boldsymbol n} = \left[ {{\begin{array}{*{20}c} {n_x } & 0 & 0 & {n_y } & 0 & {n_z } \\ 0 & {n_y } & 0 & {n_x } & {n_z } & 0 \\ 0 & 0 & {n_z } & 0 & {n_y } & {n_x } \\ \end{array} }} \right] $$对方程(1)和(2)采用伽辽金方法离散, 考虑边界条件,可得任一结点$i$的解耦运动平衡方程为
[28 ] (7) $\ddot{\boldsymbol u}_i \boldsymbol M_{{\rm s}i} + \boldsymbol F_i^{\rm s} + \boldsymbol T_i^{\rm s}-\boldsymbol S_i^{\rm s} = 0$(8) \end{equation} \begin{equation} \label{eq6} \ddot{\boldsymbol U}_i \boldsymbol M_{{\rm w}i} + \boldsymbol F_i^{\rm w} + \boldsymbol T_i^{\rm w}-\boldsymbol S_i^{\rm w} = 0 \end{equation}其中, $\boldsymbol M_{{\rm s}i}$, $\boldsymbol M_{{\rm w}i} $分别为集中在$i$结点上的固相质量和液相质量; $\boldsymbol F_i^{\rm s} $, $\boldsymbol F_i^{\rm w} $分别为集中在节点$i$的固、液相本构力; $\boldsymbol T_i^{\rm s} $,$\boldsymbol T_i^{\rm w} $分别为集中在节点$i$的固、液相黏性阻力; $\boldsymbol S_i^{\rm s} $, $\boldsymbol S_i^{\rm w}$分别作用在节点$i$的固、液相边界力. 在同一介质内,由于所有位移和应力都连续,通过单元界面作用在它们之上的应力大小相等、方向相反,在单元组装的过程中, 内部节点的$\boldsymbol S_i^{\rm s} $和$\boldsymbol S_i^{\rm w}$均为零.
(9a) $$\boldsymbol M_{{\rm s}i} = \sum\limits_e \sum\limits_{j = 1}^J \int_{\varOmega^e} \boldsymbol N_i^{\rm T}(1-\beta )\rho _{\rm s} \boldsymbol N_j {\rm d}V $$(9b) $$\boldsymbol M_{{\rm w}i} = \sum\limits_e \sum\limits_{j = 1}^J \int_{\varOmega ^e} \boldsymbol N_i^{\rm T} \beta \rho _{\rm w} \boldsymbol N_j {\rm d}V $$(9c) $$\boldsymbol F_i^{\rm s} = \sum\limits_e { f}_i^{\rm se} = \sum\limits_e \int_{\varOmega ^e} (\boldsymbol L_{\rm s} \boldsymbol N_i )^{\rm T} \boldsymbol \sigma{\rm d}V $$(9d) $$\boldsymbol F_i^{\rm w} = \sum\limits_e { f}_i^{\rm we} = \sum\limits_e \int_{\varOmega^e} (\boldsymbol L_{\rm w} \boldsymbol N_i )^{\rm T} \beta P {\rm d}V $$(9e) $$\boldsymbol T_i^{\rm s} = \sum\limits_e \int_{\varOmega ^e} \boldsymbol N_i^{\rm T}b\boldsymbol N_j (\dot{\boldsymbol u}_j- \dot{\boldsymbol U}_j) {\rm d}V $$(9f) $$\boldsymbol T_i^{\rm w} =-\boldsymbol T_i^{\rm s} $$(9g) $$\boldsymbol S_i^{\rm s} = \sum\limits_e \int_{S^e} \boldsymbol N_i ^{\rm T} \hat{\boldsymbol n}\boldsymbol \sigma{\rm d}S $$(9h) $$\boldsymbol S_i^{\rm w} = \sum\limits_e \int_{S^e} \boldsymbol N_i^{\rm T} \boldsymbol n\beta P{\rm d}S $$若节点$i$为内部节点(非界面点), $\boldsymbol S_i^{\rm s}$和$\boldsymbol S_i^{\rm w} $均为零, 给定本构关系,可对式(7)、式(8)采用时步积分求解
[30 ] .这里讨论节点$i$为两种不同介质的界面点情形, 如
图2 所示.采用隔离体概念, 则介质一中界面点$i$的动力方程由式(7)、式(8)描述,介质二中界面点$k$的动力方程表示如下(除微分算子和形函数外,其他物理量均带上划线, 以区别于介质一)
图2 新窗口打开 |
下载原图ZIP |
生成PPT 图2界面受力示意图 Fig. 2Schematic diagram of interfacial force (10) $$\ddot{\overline{\boldsymbol u}}_k \bar{\boldsymbol M}_{{\rm s}k} + \bar{\boldsymbol F}_k^{\rm s} + \bar{\boldsymbol T}_k^{\rm s}- \bar{\boldsymbol S}_k^{\rm s} = 0 $$(11) $$\ddot{\overline{\boldsymbol U}}_k \bar{\boldsymbol M}_{{\rm w}k} + \bar{\boldsymbol F}_k^{\rm w} + \bar{\boldsymbol T}_k^{\rm w}- \bar{\boldsymbol S}_k^{\rm w} = 0 $$其中
(12a) $$\bar{\boldsymbol M}_{{\rm s}k} = \sum\limits_e \sum\limits_{j = 1}^J \int_{\varOmega ^e} \boldsymbol N_k^{\rm T} (1-\bar{\beta })\bar{\rho}_{\rm s} \boldsymbol N_j {\rm d}V $$(12b) $$\bar{\boldsymbol M}_{{\rm w}k} = \sum\limits_e \sum\limits_{j = 1}^J \int_{\varOmega ^e} \boldsymbol N_k^{\rm T} \bar{\beta }\bar {\rho}_{\rm w}\boldsymbol N_j {\rm d}V $$(12c) $$\bar{\boldsymbol F}_k^{\rm s} = \sum\limits_e \bar{\boldsymbol f}_k^{{\rm s}e} = \sum\limits_e \int_{\varOmega ^e} (\boldsymbol L_{\rm s} \boldsymbol N_k )^{\rm T} \bar{\boldsymbol \sigma}{\rm d}V $$(12d) $$\bar{\boldsymbol F}_k^{\rm w} = \sum\limits_e \bar{\boldsymbol f}_k^{{\rm w}e} = \sum\limits_e \int_{\varOmega ^e} (\boldsymbol L_{\rm w} \boldsymbol N_k )^{\rm T}\bar{\beta}\bar {P} {\rm d}V $$(12e) $$\bar{\boldsymbol T}_k^{\rm s} = \sum\limits_e \int_{\varOmega ^e} \boldsymbol N_k^{\rm T}\bar {b}\boldsymbol N_j (\dot{\overline{\boldsymbol u}}_j-\dot{\overline{\boldsymbol U}}_j ) {\rm d}V $$(12f) $$\bar{\boldsymbol T}_k^{\rm w} =-\bar{\boldsymbol T}_k^{\rm s} $$(12g) $$\bar{\boldsymbol S}_k^{\rm s} = \sum\limits_e \int_{S^e} \boldsymbol N_k ^{\rm T} \hat{\overline{\boldsymbol n}}\bar{\boldsymbol \sigma}{\rm d}S $$(12h) $$\bar{\boldsymbol S}_k^{\rm w} = \sum\limits_e \int_{S^e} \boldsymbol N_k ^{\rm T} \bar{\boldsymbol n}\bar{\beta}\bar{P}{\rm d}S $$此时, $\boldsymbol S_i^s $和$\boldsymbol S_i^w$为介质二作用在界面点$i$上的界面力, $\bar{\boldsymbol S}_k^s$和$\bar{\boldsymbol S}_k^w $为介质一作用在界面点$k$上的界面力,它们之间的关系由如下界面连续条件
[31 ] 确定
(13a) $$\boldsymbol \sigma_{\rm N} + \boldsymbol \tau = \bar{\boldsymbol \sigma}_{\rm N} + \bar{\boldsymbol \tau} $$(13b) $$\boldsymbol \sigma_{\rm T} = \bar{\boldsymbol \sigma}_{\rm T} $$(13c) $$P-\bar {P} = 0 $$(13d) $$\boldsymbol u_{\rm N} = \bar{\boldsymbol u}_{\rm N},\quad \boldsymbol u_{\rm T} = \bar{\boldsymbol u}_{\rm T} $$(13e) $$\beta (\boldsymbol U_{\rm N}-\boldsymbol u_{\rm N} ) = \bar {\beta }(\bar{\boldsymbol U}_{\rm N}-\bar{u}_{\rm N} ) $$其中, 下标N表示矢量的法向分量, 下标T表示矢量得切向分量.在文献
[22 ] 中, 已讨论水平界面情形, 即水平成层海域场地情形.对于非水平成层情形, 采用文献
[22 ] 中的方法, 将变得十分繁琐.本文采用求解界面力的方法.
采用如下显式积分格式
(14) $$\ddot{U}^p = \frac{1}{\Delta t^2}(U^{(p + 1)}-2U^p + U^{(p-1)}) $$(15) $$\dot {U}^p = \frac{1}{\Delta t}(U^p-U^{(p-1)}) $$其中, 上标$p$表示$t = p\Delta t$时刻, 下同.
对式(7)、式(8)~$p$时刻的平衡方程进行时步积分得
(16) $$\boldsymbol u_i^{(p + 1)} = \hat{\boldsymbol u}_i^{(p + 1)} +\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} + \Delta \boldsymbol u_{{\rm T}i}^{(p + 1)} $$(17) $$\boldsymbol U_i^{(p + 1)} = \hat{\boldsymbol U}_i^{(p + 1)} +\Delta \boldsymbol U_{{\rm N}i}^{(p + 1)} + \Delta \boldsymbol U_{{\rm T}i}^{(p + 1)} $$其中
(18a) $$\hat{\boldsymbol u}_i^{(p + 1)} = 2\boldsymbol u_i^p- \boldsymbol u_i^{(p-1)}-\frac{(\Delta t)^2}{m_i^{\rm s} }(\boldsymbol F_i^{{\rm s}p} + \boldsymbol T_i^{{\rm s}p} ) $$(18b) $$\hat{\boldsymbol U}_i^{(p + 1)} = 2 \hat{\boldsymbol U}_i^p- \hat{\boldsymbol U}_i^{(p-1)}-\frac{(\Delta t)^2}{m_i^w }(\boldsymbol F_i^{{\rm w}p} + \boldsymbol T_i^{{\rm w}p} ) $$(18c) $$\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} = \frac{(\Delta t)^2}{m_i^s }\boldsymbol S_{{\rm N}i}^{{\rm s}p} $$(18d) $$\Delta \boldsymbol u_{{\rm T}i}^{(p + 1)} = \frac{(\Delta t)^2}{m_i^s }\boldsymbol S_{{\rm T}i}^{{\rm s}p} $$(18e) $$\Delta \boldsymbol U_{{\rm N}i}^{(p + 1)} = \frac{(\Delta t)^2}{m_i^w }\boldsymbol S_{{\rm N}i}^{{\rm w}p} $$(18f) $$\Delta \boldsymbol U_{{\rm T}i}^{(p + 1)} = \frac{(\Delta t)^2}{m_i^w }\boldsymbol S_{{\rm T}i}^{{\rm w}p} $$其中, $\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)}$是由法向界面力$\boldsymbol S_{{\rm N}i}^{{\rm s}p}$引起的固相位移, $\Delta\boldsymbol u_{{\rm T}i}^{(p + 1)}$是由切向界面力$\boldsymbol S_{{\rm T}i}^{{\rm s}p}$引起的固相位移, $\Delta \boldsymbol U_{{\rm N}i}^{(p + 1)}$是由法向界面力$\boldsymbol S_{{\rm N}i}^{{\rm w}p}$引起的液相位移, $\Delta \boldsymbol U_{{\rm T}i}^{(p + 1)}$是由切向界面力$\boldsymbol S_{{\rm T}i}^{{\rm w}p}$引起的液相位移.
同样可得
(19) $$\bar{\boldsymbol u}_k^{(p + 1)} = \hat{\overline{\boldsymbol u}}_k^{(p + 1)} + \Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)} + \Delta \bar{\boldsymbol u}_{{\rm T}k}^{(p + 1)} $$(20) $$\bar{\boldsymbol U}_k^{(p + 1)} = \hat{\overline{\boldsymbol U}}_k^{(p + 1)} + \Delta \bar{\boldsymbol U}_{{\rm N}k}^{(p + 1)} + \Delta \bar{\boldsymbol U}_{{\rm T}k}^{(p + 1)} $$右端各项如同式(18).
由界面连续条件可知(其中式(21c)利用了式(13c)、式(9h)、式(12h))
(21a) $$\boldsymbol S_{{\rm N}i}^{\rm s} + \boldsymbol S_{{\rm N}i}^{\rm w} =-(\bar{\boldsymbol S}_{{\rm N}k}^{\rm s} + \bar{\boldsymbol S}_{{\rm N}k}^{\rm w}) $$(21b) $$\boldsymbol S_{{\rm T}i}^{\rm s} =-\bar{\boldsymbol S}_{{\rm T}k}^{\rm s} $$(21c) $$\bar{\beta }\boldsymbol S_{{\rm N}i}^{\rm w} =-\beta \bar{\boldsymbol S}_{{\rm N}k}^{\rm w} $$(21d) $$\boldsymbol S_{{\rm T}i}^{\rm w} = \bar{\boldsymbol S}_{{\rm T}k}^{\rm w} = 0 $$(21e) $$\boldsymbol u_{{\rm N}i} = \bar{\boldsymbol u}_{{\rm N}k} $$(21f) $$\boldsymbol u_{{\rm T}i} = \bar{\boldsymbol u}_{{\rm T}k} $$(21g) $$\beta \boldsymbol U_{{\rm N}i} = \bar {\beta }(\bar{\boldsymbol U}_{{\rm N}k}-\bar{\boldsymbol u}_{{\rm N}k} ) + \beta \bar{\boldsymbol u}_{{\rm N}k} $$将式(16)、式(19)代入界面处固相法向位移连续条件式(21e),可得
(22) \begin{equation} \label{eq47} \boldsymbol n_i \left( \boldsymbol n_i \cdot \hat{\boldsymbol u}_i^{(p + 1)} \right) + \Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} = \boldsymbol n_i \left( \boldsymbol n_i \cdot \hat{\overline{\boldsymbol u}}_k^{(p + 1)} \right) + \Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)} \end{equation}利用式(18c)、式(21a)、式(21c), 得
(23) $$\boldsymbol n_i \left( \boldsymbol n_i \cdot \hat{\boldsymbol u}_i^{(p + 1)} \right) + \frac{(\Delta t)^2}{m_i^{\rm s} }\boldsymbol S_{{\rm N}i}^{{\rm s}p} = \boldsymbol n_i \left( {\boldsymbol n_i \cdot \hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right)+ \\ \frac{(\Delta t)^2}{m_k^{\rm s}}\left[-\boldsymbol S_{{\rm N}i}^{{\rm s}p}-\boldsymbol S_{{\rm N}i}^{{\rm w}p} + (\bar{\beta}/\beta )\boldsymbol S_{{\rm N}i}^{{\rm w}p} \right] $$由界面处法向连续条件(21g), 考虑到式(16)、式(19)、式(20)、式(21a)、式(21c), 可得
(24) $$\boldsymbol n_i \left( {\beta \boldsymbol n_i \cdot \hat{\boldsymbol U}_i^{(p + 1)} } \right) + \beta \left( {\frac{(\Delta t)^2}{m_i^{\rm w}}\boldsymbol S_{{\rm N}i}^{{\rm w}p} } \right) = \\ \boldsymbol n_i \left\{ {\boldsymbol n_i \cdot \left[ {\bar {\beta }\left( {\hat{\overline{\boldsymbol U}}_k^{(p + 1)}- \hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right) + \beta \hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right]} \right\} + \bar {\beta }\cdot\\ \bigg[-\left( {\bar {\beta }/\beta } \right)\frac{(\Delta t)^2}{m_k^{\rm w}}\boldsymbol S_{{\rm N}i}^{{\rm w}p}-\frac{(\Delta t)^2}{m_k^{\rm s}}\Big( - \boldsymbol S_{{\rm N}i}^{{\rm s}p}-\boldsymbol S_{{\rm N}i}^{{\rm w}p} + \\ \left( {\bar {\beta }/\beta } \right)\boldsymbol S_{{\rm N}i}^{{\rm w}p} \Big) \bigg]+ \beta \left[ {\frac{(\Delta t)^2}{m_k^{\rm s} }\left( {- \boldsymbol S_{{\rm N}i}^{{\rm s}p}-\boldsymbol S_{{\rm N}i}^{{\rm w}p} + \left( {\bar {\beta }/\beta } \right)\boldsymbol S_{{\rm N}i}^{{\rm w}p} } \right)} \right] $$由式(23)和式(24), 可解得
(25) $$\boldsymbol S_{{\rm N}i}^{{\rm s}p} = \frac{A_{22} B_1-A_{12} B_2 }{A_{22} A_{11}-A_{12} A_{21} } $$(26) $$\boldsymbol S_{{\rm N}i}^{{\rm w}p} = \frac{A_{21} B_1-A_{11} B_2 }{A_{21} A_{12}-A_{11} A_{22} } $$其中
(27a) $$A_{11} = (\Delta t)^2\left( {\frac{1}{m_i^{\rm s}} + \frac{1}{m_k^{\rm s}}} \right) $$(27b) $$A_{12} = \left( {1-\frac{\bar {\beta }}{\beta }} \right)\frac{(\Delta t)^2}{m_k^{\rm s}} $$(27c) $$B_1 = \boldsymbol n_i \left[ {\boldsymbol n_i \cdot \left( {\hat{\overline{\boldsymbol u}}_k^{(p + 1)}-\hat{\boldsymbol u}_i^{(p + 1)} } \right)} \right] $$(27d) $$A_{21} = (\Delta t)^2\left( {-\frac{\bar {\beta }}{m_k^{\rm s} } + \frac{\beta }{m_k^{\rm s}}} \right) $$(27e) $$A_{22} = (\Delta t)^2\left[ {\frac{\beta }{m_i^{\rm w}} + \frac{\bar {\beta }^2}{\beta m_k^{\rm w}} + \frac{\left( {\bar {\beta }-\beta } \right)\left( {\bar {\beta }-\beta } \right)}{\beta m_k^{\rm s}}} \right] $$(27f) $$B_2 = \boldsymbol n_i \left\{ {\boldsymbol n_i \cdot \left[ {\bar {\beta }\left( {\hat{\overline{\boldsymbol U}}_k^{(p + 1)}- \hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right)-\beta \left( {\hat{\boldsymbol U}_i^{(p + 1)}-\hat{\boldsymbol u}_k^{(p + 1)} } \right)} \right]} \right\} $$由式(25)、式(26)求得$\boldsymbol S_{{\rm N}i}^{{\rm s}p}$和$\boldsymbol S_{{\rm N}i}^{{\rm w}p} $后,可进一步由式(21a)、式(21c)解得$\bar{\boldsymbol S}_{{\rm N}k}^{{\rm s}p} $和$\bar{\boldsymbol S}_{{\rm N}k}^{{\rm w}p} $.
切向方向导数不确定, 所以固相切向位移连续条件不能直接用.由于界面处固相法向和切向位移连续, 因此界面处固相位移连续
(28) \begin{equation} \label{eq58} \hat{\boldsymbol u}_i^{(p + 1)} + \Delta \boldsymbol u_{Ni}^{(p + 1)} + \Delta \boldsymbol u_{Ti}^{(p + 1)} = \hat{\overline{\boldsymbol u}}_k^{(p + 1)} + \Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)} + \Delta \bar{\boldsymbol u}_{{\rm T}k}^{(p + 1)} \end{equation}将式(18d)、式(21b)代入式(28), 整理得
(29) $$\left( {\frac{(\Delta t)^2}{m_i^{\rm s}} + \frac{(\Delta t)^2}{m_k^{\rm s}}} \right)\boldsymbol S_{{\rm T}i}^{{\rm s}p} = \\ \left( {\hat{\overline{\boldsymbol u}}_k^{(p + 1)} + \Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)}-\hat{\boldsymbol u}_i^{(p + 1)}-\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} } \right) $$由此可解得
(30) \begin{equation} \label{eq60} \boldsymbol S_{{\rm T}i}^{{\rm s}p} = \frac{\left( {\hat{\overline{\boldsymbol u}}_k^{(p + 1)} + \Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)}-\hat{\boldsymbol u}_i^{(p + 1)}-\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} } \right)m_i^{\rm s} m_k^{\rm s} }{(\Delta t)^2(m_i^{\rm s} + m_k^{\rm s} )} \end{equation}式中, $\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} $可由式(18c)求得,同样可求得$\Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)} $.求得$\boldsymbol S_{{\rm T}i}^{{\rm s}p} $后,由式(21b)可得$\boldsymbol S_{{\rm T}k}^{{\rm s}p} $. 有了界面力,由式(16)、式(17)可求得$\boldsymbol u_i^{(p + 1)} $, $\boldsymbol U_i^{(p + 1)} $. 同样, 由式(19)、式(20)可得$\bar{\boldsymbol u}_k^{(p + 1)} $, $\bar{\boldsymbol U}_k^{(p + 1)} $.
1.2 特殊情形 (1)流体--饱和土情形
考虑饱和土上覆流体情形. 此时,
图2 中介质二为饱和多孔介质,介质一为流体, 则$\beta = 1$, $\boldsymbol M_{{\rm s}i} = 0$.考虑无黏的理想流体, 动黏度系数为零, 则$b = 0$, $\boldsymbol T_i^{\rm s} = 0$, $\boldsymbol T_i^{\rm w} = 0$;不存在固相, 固相骨架模量可取为零, 则$\boldsymbol F_i^{\rm s} = 0$, $\boldsymbol S_i^{\rm s} = 0$; 由此, 式(7)自动满足,式(8)退化为如下理想流体方程
(31) \begin{equation}\label{eq61} \ddot{\boldsymbol U}_i \boldsymbol M_{{\rm w}i} +\boldsymbol F_i^{\rm w}-\boldsymbol S_i^{\rm w} = 0 \end{equation}因此, 饱和多孔介质方程可以用来分析理想流体的动力响应.
流体--饱和土界面连续条件为
[32 ] (32a) $$\boldsymbol S_{{\rm N}i}^{\rm w} =-(\bar{\boldsymbol S}_{{\rm N}k}^{\rm s} + \bar{\boldsymbol S}_{{\rm N}k}^{\rm w}) $$(32b) $$\bar {\beta }\boldsymbol S_{{\rm N}i}^{\rm w} =-\bar{\bf S}_{{\rm N}k}^{\rm w} $$(32c) $$\boldsymbol U_{{\rm N}i} = \bar {\beta }(\bar{\boldsymbol U}_{{\rm N}k}-\bar{\boldsymbol u}_{{\rm N}k} ) + \bar{\boldsymbol u}_{{\rm N}k} $$由动力方程(10)、(11)和(31)以及界面连续条件式(32), 按2.1节的方法,可得界面点的运动方程如下
(33a) $$\boldsymbol n_i \left( {\boldsymbol n_i \cdot \hat{\boldsymbol U}_i^{(p + 1)} } \right) + \frac{(\Delta t)^2}{m_i^{\rm w}}\boldsymbol S_{{\rm N}i}^{{\rm w}p} = \boldsymbol n_i \\ \left[ {\boldsymbol n_i \cdot \left( {\bar {\beta }\hat{\overline{\boldsymbol U}}_k^{(p + 1)} + \left( {1-\bar {\beta }} \right)\hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right)} \right] + \\ \left( {-\frac{(\Delta t)^2}{m_k^{\rm w}}\bar {\beta }^2\boldsymbol S_{{\rm N}i}^{{\rm w}p}-\frac{(\Delta t)^2}{m_k^{\rm s} }\left( {\bar {\beta }-1} \right)^2\boldsymbol S_{{\rm N}i}^{{\rm w}p} } \right) $$由式(33a)可解得
(33b) $$\label{eq66} \boldsymbol S_{{\rm N}i}^{{\rm w}p} = \frac{ B_2 }{A_{22} } $$其中
(33c) $$A_{22} = (\Delta t)^2\left( {\frac{1}{m_i^{\rm w} } + \frac{\bar {\beta }^2}{m_k^{\rm w}} + \frac{\left( {\bar {\beta }-1} \right)^2}{m_k^{\rm s} }} \right) $$(33d) $$B_2 = \boldsymbol n_i \left\{ {\boldsymbol n_i \cdot \left[ {\left( {\bar {\beta }\hat{\overline{\boldsymbol U}}_k^{(p + 1)} + \left( {1-\bar {\beta }} \right)\hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right)-\hat{\boldsymbol U}_i^{(p + 1)} } \right]} \right\}\qquad\quad $$求得$\boldsymbol S_{{\rm N}i}^{{\rm w}p} $后,可由式(32a)和式(32b)可得$\bar{\boldsymbol S}_{{\rm N}k}^{{\rm s}p}$和$\bar{\boldsymbol S}_{{\rm N}k}^{{\rm w}p} $.
由于流体为无黏理想流体, 切向界面力为零. 有了界面力,则可由式(17)、式(19)、式(20)求得位移及其他响应.
(2)饱和土--干基岩情形
图2 中的介质二若为不透水的基岩, 不存在液相, 则$\bar {\beta } = 0$,$\bar{\boldsymbol M}_{{\rm w}k} = 0$;液相体积模量和固液之间的黏性力取为零, 则$\bar{\boldsymbol F}_k^{\rm w} = 0$, $\bar{\boldsymbol T}_k^{\rm w} = 0$,$\bar{\boldsymbol T}_k^{\rm s} = 0$, $\bar{\boldsymbol S}_k^{\rm w} = 0$, 方程(11)自动满足,式(10)退化为如下干基岩方程
(34) \begin{equation} \label{eq69} \ddot{\overline{\boldsymbol u}}_k \bar{\boldsymbol M}_{{\rm s}k} + \bar{\boldsymbol F}_k^{\rm s}-\bar{\boldsymbol S}_k^{\rm s} = \boldsymbol 0 \end{equation}因此, 饱和多孔介质方程可以用来分析干基岩的动力响应.
饱和土--干基岩界面条件为
[31 ] (35a) $$\boldsymbol S_{{\rm N}i}^{\rm s} + \boldsymbol S_{{\rm N}i}^{\rm w} =-\bar{\boldsymbol S}_{{\rm N}k}^{\rm s} $$(35b) $$\boldsymbol S_{{\rm T}i}^{\rm s} =-\bar{\boldsymbol S}_{{\rm T}k}^{\rm s} $$(35c) $$\boldsymbol u_{{\rm N}i} = \bar{\bf u}_{{\rm N}k} $$(35d) $$\boldsymbol u_{{\rm T}i} = \bar{\boldsymbol u}_{{\rm T}k} $$(35e) $$\beta (\boldsymbol U_{{\rm N}i}-\boldsymbol u_{{\rm N}i} ) = 0 $$由界面处固相法向位移连续,将式(16)、式(19)代入式(35c)可得
(36) $$\boldsymbol n_i \left( {\boldsymbol n_i \cdot \hat{\boldsymbol u}_i^{(p + 1)} } \right) + \Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} = \boldsymbol n_i \left( {\boldsymbol n_i \cdot \hat{\bar{\boldsymbol u}}_k^{(p + 1)} } \right) + \Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)} $$整理得
(37a) $$\bar {u}_{ky}^{(p + 1)} = 2\bar {u}_{ky}^p-\bar {u}_{ky}^{(p- 1)}- \frac{\Delta t^2}{(m_{{\rm s}i} + \bar {m}_{{\rm s}k} )} \cdot \\ (F_{iy}^{{\rm s}p} + \bar {F}_{ky}^{{\rm s}p} + T_{iy}^{{\rm s}p} ) $$(37b) $$\boldsymbol n_i \left( {\boldsymbol n_i \cdot \hat{\boldsymbol u}_i^{(p + 1)} } \right) + \frac{(\Delta t)^2}{m_i^{\rm s} }\boldsymbol S_{{\rm N}i}^{{\rm s}p} = \boldsymbol n_i \left( {\boldsymbol n_i \cdot \hat{\bar{\boldsymbol u}}_k^{(p + 1)} } \right)+ \\ \frac{(\Delta t)^2}{m_k^{\rm s}}\left( {-\boldsymbol S_{{\rm N}i}^{{\rm s}p}-\boldsymbol S_{{\rm N}i}^{{\rm w}p} } \right) $$由界面处法向连续条件(35e), 并将式(16)、式(17)代入可得
(38) $$\boldsymbol n_i \left[ {\beta \boldsymbol n_i \cdot \left( {\hat{\boldsymbol U}_i^{(p + 1)}-\hat{\boldsymbol u}_i^{(p + 1)} } \right)} \right] + \\ \beta (\Delta \boldsymbol U_{{\rm N}i}^{(p + 1)}-\Delta \boldsymbol u_{{\rm N}i}^{(p + 1)} ) = 0 $$整理得
(39) $$\boldsymbol n_i \left[ {\beta \boldsymbol n_i \cdot \left( {\hat{\boldsymbol U}_i^{(p + 1)}-\hat{\boldsymbol u}_i^{(p + 1)} } \right)} \right] + \\ \beta \left( {\frac{(\Delta t)^2}{m_i^{\rm w}}\boldsymbol S_{{\rm N}i}^{{\rm w}p}-\frac{(\Delta t)^2}{m_i^{\rm s}}\boldsymbol S_{{\rm N}i}^{{\rm s}p} } \right) = 0 $$由式(37)和式(39), 可解得
(40a) $$\boldsymbol S_{{\rm N}i}^{{\rm s}p} = \frac{A_{22} B_1- A_{12} B_2 }{A_{22} A_{11}-A_{12} A_{21} } $$(40b) $$\boldsymbol S_{{\rm N}i}^{{\rm w}p} = \frac{A_{21} B_1- A_{11} B_2 }{A_{21} A_{12}-A_{11} A_{22} } $$其中
(41a) $$A_{11} = (\Delta t)^2\left( {\frac{1}{m_i^{\rm s}} + \frac{1}{m_k^{\rm s}}} \right) $$(41b) $$A_{12} = \frac{(\Delta t)^2}{m_k^{\rm s}} $$(41c) $$B_1 = \boldsymbol n_i \left[ {\boldsymbol n_i \cdot \left( {\hat{\overline{\boldsymbol u}}_k^{(p + 1)}-\hat{\boldsymbol u}_i^{(p + 1)} } \right)} \right] $$(41d) $$A_{21} = \frac{-\beta (\Delta t)^2}{m_i^{\rm s}} $$(41e) $$A_{22} = \frac{(\Delta t)^2\beta }{m_i^{\rm w}} $$(41f) $$B_2 = \boldsymbol n_i \left\{ {\boldsymbol n_i \cdot \left[ {\beta \left( {\hat{\boldsymbol u}_i^{(p + 1)}-\hat{\boldsymbol U}_i^{(p + 1)} } \right)} \right]} \right\} $$由式(35a)可求得$\bar{\boldsymbol S}_{{\rm N}k}^{{\rm s}p} $.与1.1节中相同, 可求得$\boldsymbol S_{{\rm T}i}^{{\rm s}p}$和$\bar{\boldsymbol S}_{{\rm T}k}^{{\rm s}p} $,进而可求得界面点的位移响应.
(3)流体--干基岩情形
此时,
图2 中介质一为理想流体, 介质二为干基岩. 按前述参数取值,运动方程分别退化为式(31)、式(34).流体--干基岩界面连续条件为(不透水)
(42a) $$\boldsymbol S_{{\rm N}i}^{\rm w} =-\bar{\boldsymbol S}_{{\rm N}k}^{\rm s} $$(42b) $$\boldsymbol U_{{\rm N}i} = \bar{\boldsymbol u}_{{\rm N}k} $$将式(17)、式(19)代入式(42b), 可得
(42c) $$\label{eq89} \boldsymbol n_i \left[ {\boldsymbol n_i \cdot \left( {\hat{\boldsymbol U}_i^{(p + 1)}-\hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right)} \right] + {\Delta \boldsymbol U_{{\rm N}i}^{(p + 1)}-\Delta \bar{\boldsymbol u}_{{\rm N}k}^{(p + 1)} } = 0 $$利用式(18), 整理得
(43a) $$\boldsymbol n_i \left[ {\boldsymbol n_i \cdot \left( {\hat{\boldsymbol U}_i^{(p + 1)}-\hat{\overline{\boldsymbol u}}_k^{(p + 1)} } \right)} \right] + \\ \left( {\frac{(\Delta t)^2}{m_i^{\rm w} }\boldsymbol S_{{\rm N}i}^{{\rm w}p} + \frac{(\Delta t)^2}{m_k^{\rm s}}\boldsymbol S_{{\rm N}i}^{{\rm w}p} } \right)= 0 $$可解得
(43b) $$\label{eq91} \boldsymbol S_{{\rm N}i}^{{\rm w}p} = \frac{\boldsymbol B_2 }{A_{22} } $$其中
(43c) $$A_{22} = (\Delta t)^2\left( {\frac{1}{m_i^{\rm w}} + \frac{1}{m_k^{\rm s}}} \right) $$(43d) $$\boldsymbol B_2 = \boldsymbol n_i \left[ {\boldsymbol n_i \cdot \left( {\hat{\overline{\boldsymbol u}}_k^{(p + 1)}- \hat{\boldsymbol U}_i^{(p + 1)} } \right)} \right] $$由式(42a)可得$\bar{\boldsymbol S}_{{\rm N}k}^{{\rm s}p} $,切向界面力为零, 进而可得界面点的位移响应.
综上可知,流体、固体、饱和多孔介质之间的耦合均可统一在饱和多孔介质理论框架进行分析.
1.3 实施方法 以
图3 所示的模型为例,考虑饱和海床表面一凹陷地形情形的海水--饱和海床--基岩非水平成层体系,在平面P波垂直入射时的反应.
图3 新窗口打开 |
下载原图ZIP |
生成PPT 图3模型示意图 Fig. 3Schematic diagram of model 对于该复杂地形, 需要通过有限元等数值方法进行求解.因此对
图3 (a)所示的问题采用本文提出的方法求解, 计算模型如
图3 (b),在侧面和底面采用多次透射人工边界
[23 -24 ] ,用于模拟无限域.在人工边界上采用传递矩阵方法得到的自由场
[35 ] 作为波动输入.海水--饱和海床--基岩体系中,内部点(除人工边界点和界面点外的其余节点)及界面点的响应计算,均采用统一的饱和多孔介质方程和界面点的计算方法,即1.1节和1.2节中的理论方法; 人工边界点的响应通过多次透射公式计算.基于此, 编制了相应的三维有限元计算程序, 以实现本文方法.
2 算例验证 采用如下线弹性本构
(44a) $$\sigma '_{ij}-(1-\alpha )\delta _{ij} P = 2\mu e_{ij} +\delta _{ij} \lambda e $$(44b) $$P =-\alpha Me + M\zeta $$其中, 根据文献
[33 ,34 ] , $\lambda $和$\mu$为固相骨架在排水情形时的拉梅常数, $\alpha$和$M$则为表征饱和多孔介质压缩性的常数, 以压缩模量表示时, 为
(45) $$\alpha = 1-\frac{E_{\rm b}}{E_{\rm u}} $$(46) $$M = \frac{E_u^2 }{E_{\rm u} \left[n\left(\dfrac{E_{\rm u}}{E_{\rm w}}-1\right) + 1\right]-E_{\rm b}} $$其中, $E_{\rm u}$和$E_{\rm b}$分别为不排水和排水时饱和多孔介质的压缩模量, $E_{\rm w}$为孔隙流体的压缩模量, $k_0 $为渗透系数.
下面所有算例模型, $x$和$y$方向尺寸均为40 m,$z$方向的尺寸如模型图所示. 将计算域划分为八节点六面体单元,单元尺寸为1 m $\times$1 m $\times$1 m. 输入
图4 所示的单位脉冲,时间步距为$\Delta t = 0.000~1$s. 单元尺寸$\varDelta \le \lambda_{\min } / 10$, 满足精度要求, 其中$\lambda _{\min }$为所需模拟的最小波长. 所有材料参数见
表1 .
图4 新窗口打开 |
下载原图ZIP |
生成PPT 图4输入脉冲 Fig. 4Input pulse Table 1 表1 表1 材料参数表
Table 1
Parameters of material 新窗口打开 |
下载CSV 2.1 海水--基岩情形 考虑基岩海床表面一三维凹陷(如
图5 (b))地形的海水--基岩非水平成层模型,其剖面如
图5 (a)所示.
图6 为P波垂直入射时各点的响应, 图例中,位移$U$的下标表示方向, 上标w表示海水, 上标b表示基岩,上标s表示饱和土, 如$u_x^{\rm b}$为基岩固相$x$方向的位移, $U_{{\rm N}}^{\rm w}$为界面处海水法向位移. 由图中可以看出, 由于凹陷的存在,即使P波垂直入射时, 也会引起海水层水平方向的响应. $G$点由于对称性,其水平向位移为零. 由于三维海域非水平成层的波动散射问题较为复杂,目前还没有相关解析解. 因此,
图7 给出界面点法向方向海水和基岩的位移,检查其是否满足界面连续条件(42b). 从图中可以看出,界面上海水与基岩法向位移完全吻合, 满足海水--干基岩界面连续条件.
图5 新窗口打开 |
下载原图ZIP |
生成PPT 图5海水--基岩模型 Fig. 5Seawater-bedrock model 图6 新窗口打开 |
下载原图ZIP |
生成PPT 图6P波入射时海水--基岩系统的位移响应 Fig. 6Displacement of seawater-bedrock system for P wave incidence 图7 新窗口打开 |
下载原图ZIP |
生成PPT 图7P波入射时海水--基岩系统界面点法向位移响应 Fig. 7Normal displacement of interface points in seawater-bedrock system for P wave incidence 2.2 海水--饱和海床--基岩情形 海水--饱和海床--基岩模型如
图8 所示.
图9 为P波垂直入射时的响应.
图8 新窗口打开 |
下载原图ZIP |
生成PPT 图8海水--饱和土--基岩模型 Fig. 8Seawater-saturated soil-bedrock model 图9 新窗口打开 |
下载原图ZIP |
生成PPT 图9P波入射时海水--饱和土--基岩体系的位移响应 Fig. 9Displacement of seawater-saturated soil-bedrock system for P wave incidence 同样, 由于凹陷地形的存在, 即使P波垂直入射时也会引起海水$A$,$B$点处的水平位移.
界面点$B$的海水位移和饱和土固、液相位移在法向(该点处为$z$方向)满足位移连续(式32(c)),因此该点处海水$z$方向位移与固相$z$方向的位移接近,与液相$z$方向位移差异较大, 但满足(式32(c))的连续条件(如
图10 (a));
切向位移不连续(此处为$x$, $y$方向),因此该点处海水水平方向位移与饱和土水平向位移差异较大. 界面点$C$,界面为倾斜面, 位移仅在法向满足连续条件(式32(c)), 在$x$, $y$,$z$方向不连续,因此该点处饱和土位移和海水位移差异较大(
图9 (e)$\sim$
图9 (g)),但满足法向位移连续(如
图10 (b)). $D$点为角点, 其法向方向不存在,具体实施计算时, 取周围界面在该点处的单位法向矢量叠加,再归一化为单位矢量, 作为该点的单位法向导数; 同样,该点处饱和土位移和海水位移差异较大(
图9 (h)$\sim$
图9 (i)),但满足法向位移连续(如
图10 (b)).
图10 新窗口打开 |
下载原图ZIP |
生成PPT 图10P波入射时海水--饱和土--基岩体系的法向位移响应 Fig. 10Normal displacement of interface points in seawater-saturated soil-bedrock system for P wave incidence $E$和$F$分别为饱和土与基岩界面点和基岩底 面点.
图10 为P波入射时的法向位移. 从图中可以看出,海水与饱和土海床接触面上同一位置点的$\boldsymbol U_{{\rm N}1} $,$\boldsymbol U_{{\rm N}2} $完全吻合($\boldsymbol U_{{\rm N}1} $,$\boldsymbol U_{{\rm N}2} $如式(47)),满足海水--饱和土海床界面连续条件,且饱和土与基岩接触面也满足界面连续条件.
(47a) $$\boldsymbol U_{{\rm N}1} = \beta \boldsymbol U_{{\rm N}i} $$(47b) $$\boldsymbol U_{{\rm N}2} = \bar {\beta }(\bar{\boldsymbol U}_{{\rm N}k}-\bar{\boldsymbol u}_{{\rm N}k} ) + \beta \bar{\boldsymbol u}_{{\rm N}k} $$通过上述两个算列验证了本文非水平成层海域场地有限元方法的有效性.
3 结论 本文考虑海底不规则地形,将理想流体、固体、饱和多孔介质之间的耦合分析纳入到饱和多孔介质统一计算框架,建立了集中质量显式有限元求解方法, 并编制了相应的三维并行分析程序.通过分析半无限海水--弹性基岩、海水--饱和海床--弹性基岩体系凹陷地形模型在P波垂直入射时的动力响应,验证了该统一计算框架的有效性以及并行计算的可行性.
本文采用集中质量显式有限元方法结合透射边界条件,可避免求解大型线性方程组, 效率较高, 且易于实现并行计算,有望用于大型海底地震波以及海洋声场的模拟中.
[1] Nakamura T Takenaka H Okamoto T et al . FDM Simulation of seismic-wave propagation for an aftershock of the 2009 Suruga Bay earthquake: Effects of ocean-bottom topography and seawater layerBulletin of Seismological Society of America , 2012 ,102(6 ):2420 -2435 [本文引用: 1] [2] Okamoto T Takenaka H Geophysical Journal International , 1999 ,139:531 -546 [本文引用: 1] [3] Qian ZH Yamanaka H Geophysical Journal International , 2012 ,189:524 -540 [本文引用: 1] [4] Komatitsch D Barnes C Tromp J Geophysics , 2000 ,65(2 ):623 -631 [本文引用: 1] [5] Collins MD Siegmann WL IEEE Journal of Ocean Engineering , 2017 ,42(2 ):488 -493 [本文引用: 1] [6] Tang J Piao SC Zhang HG Chinese Physics B , 2017 ,26(11 ):114301 [本文引用: 1] [7] Murphy JE Chin-Bing SA The Journal of Acoustic Society of America , 1989 ,86:1478 -1483 [本文引用: 1] [8] Zhao HY Jeng DS Liao CC et al . Three-dimensional modeling of wave-induced residual seabed response around a mono-pile foundationCoastal Engineering , 2017 ,128:1 -21 [本文引用: 1] [9] Lin ZB Guo YK Jeng DS et al . An integrated numerical model for wave-soil-pipeline interactionsCoastal Engineering , 2016 ,108:25 -35 [10] Ye JH Bulletin of Earthquake Engineering , 2012 ,10:1609 -1633 [11] 李伟华 . 考虑水--饱和土场地--结构耦合时的沉管隧道地震反应分析防灾减灾工程学报 , 2010 ,30(6 ):607 -613 [本文引用: 2] ( Li Weihua Journal of Disaster Prevention and Mitigation Engineering , 2010 ,30(6 ):607 -613 (in Chinese)) [本文引用: 2] [12] Farhat C Lesoinne M Computer Methods in Applied Mechanics & Engineering , 1998 ,157(1 ):95 -114 [本文引用: 1] [13] Farhat C Lesoinne M Computer Methods in Applied Mechanics & Engineering , 2000 ,182:499 -515 [14] Farhat C Zee KGVD Geuzaine P Computer Methods in Applied Mechanics & Engineering , 2006 ,195(17 ):1973 -2001 [15] Bathe KJ Zhang H International Journal for Numerical Methods in Engineering , 2010 ,60(1 ):213 -232 [16] Degroote J Haelterman R Annerel S et al . Performance of partitioned procedures in fluid-structure interactionComputers & Structures , 2010 ,88(7 ):446 -457 [17] Hou G Wang J Layton A Communications in Computational Physics , 2012 ,12(2 ):337 -377 [18] Habchi C Russeil S Bougeard D et al . Partitioned solver for strongly coupled fluid--structure interactionComputers & Fluids , 2013 ,71(1 ):306 -319 [19] Mehl M Uekermann B Bijl H et al . Parallel coupling numerics for partitioned fluid--structure interaction simulationsComputers & Mathematics with Applications , 2016 ,71(4 ):869 -891 [20] Bungartz HJ Lindner F Gatzhammer B et al . Precice--A fully parallel library for multi-physics surface couplingComputers & Fluids , 2016 ,141:250 -258 [21] Link G Kaltenbacher M Breuer M et al . A 2D finite-element scheme for fluid--solid--acoustic interactions and its application to human phonationComputer Methods in Applied Mechanics & Engineering , 2009 ,198(41 ):3321 -3334 [本文引用: 1] [22] 陈少林 , 柯小飞 , 张洪翔 . 海洋地震工程流固耦合问题统一计算框架力学学报 , 2019 ,51(2 ):1 -13 [本文引用: 4] ( Chen Shaolin Ke Xiaofei Zhang Hongxiang Chinese Journal of Theoretical and Applied Mechanics , 2019 ,51(2 ):1 -13 (in Chinese)) [本文引用: 4] [23] 廖振鹏 . 工程波动理论导论. 第2版 . 北京: 科学出版社, 2002 : 136 -285 [本文引用: 2] ( Liao Zhenpeng Introduction to Wave Motion Theories in Engineering(2nd edition). Beijing: Science Press, 2002 : 136 -285 (in Chinese)) [本文引用: 2] [24] Liao ZP Wong HL Soil Dyn Earthq Eng , 1984 ,3:174 -183 [本文引用: 1] [25] 邢浩洁 , 李鸿晶 . 透射边界条件在波动谱元模拟中的实现:二维波动力学学报 , 2017 ,49(4 ):894 -906 [本文引用: 1] ( Xing Haojie Li Hongjing Chinese Journal of Theoretical and Applied Mechanics , 2017 ,49(4 ):894 -906 (in Chinese)) [本文引用: 1] [26] 谷音 , 刘晶波 , 杜修力 . 三维一致粘弹性人工边界及等效粘弹性边界单元工程力学 , 2007 ,24(12 ):31 -37 [本文引用: 1] ( Gu Lin Liu Jingbo Du Xiuli Journal of Engineering Mechanics , 2007 ,24(12 ):31 -37 (in Chinese)) [本文引用: 1] [27] 刘晶波 , 宝鑫 , 谭辉 等 . 波动问题中流体介质的动力人工边界力学学报 , 2017 ,49(6 ):1418 -1427 [本文引用: 1] ( Liu Jingbo Bao Xin Tan Hui Wang Jianping Guo Dong Chinese Journal of Theoretical and Applied Mechanics , 2017 ,49(6 ):1418 -1427 (in Chinese)) [本文引用: 1] [28] 赵宇昕 , 陈少林 . 关于传递矩阵法分析饱和成层介质响应问题的讨论力学学报 , 2016 ,48(5 ):1145 -1158 [本文引用: 2] ( Zhao Yuxin Chen Shaolin Chinese Journal of Theoretical and Applied Mechanics , 2016 ,48(5 ):1145 -1158 (in Chinese)) [本文引用: 2] [29] 刘晶波 , 谭辉 , 宝鑫 等 . 土--结构动力相互作用分析中基于人工边界子结构的地震波动输入方法力学学报 , 2018 ,50(1 ):32 -43 [本文引用: 1] ( Liu Jingbo Tan Hui Bao Xin et al . The seismic wave input method for soil-structure dynamic interaction analysis based on the substructure of artificial boundariesChinese Journal of Theoretical and Applied Mechanics , 2018 ,50(1 ):32 -43 (in Chinese)) [本文引用: 1] [30] 陈少林 , 廖振鹏 , 陈进 . 两相介质近场波动模拟的解耦方法,地球物理学报 , 2005 ,48(4 ):909 -917 [本文引用: 2] ( Chen Shaolin Liao Zhenpeng Chen Jin Journal of Geophysics , 2005 ,48(4 ):909 -917 (in Chinese)) [本文引用: 2] [31] Deresiewicz H Rice JT Bull Seis Soc Am , 1964 ,54(1 ):409 -416 [本文引用: 2] [32] Deresiewicz H Bull Seis Soc Am , 1964 ,54(1 ):425 -430 [本文引用: 1] [33] Biot MA Acoust Soc Am , 1956 ,28:168 -191 [本文引用: 1] [34] Biot MA Journal of Applied Physics , 1962 ,33(4 ):1482 -1498 [本文引用: 1] [35] 柯小飞 , 陈少林 , 张洪翔 . P-SV波入射时海水--层状海床体系的自由场分析振动工程学报 , 2018 , 录用 [本文引用: 1] ( Ke Xiaofei Chen Shaolin Zhang Hongxiang Journal of Vibration Engineering , 2018 , accepted (in Chinese)) [本文引用: 1] 
,2), 程书林, 柯小飞南京航空航天大学土木工程系, 南京 210016
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
