陈少林²⁾, 柯小飞, 张洪翔南京航空航天大学土木工程系,南京 210000

A UNIFIED COMPUTATIONAL FRAMEWORK FOR FLUID-SOLID COUPLING IN MARINE EARTHQUAKE ENGINEERING ¹⁾

Chen Shaolin²⁾, Ke Xiaofei, Zhang HongxiangDepartment of Civil Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

收稿日期:2018-10-12接受日期:2018-12-18网络出版日期:2019-03-18

基金资助:

国家自然科学基金资助项目.51178222 国家自然科学基金资助项目.51278260

Received:2018-10-12Accepted:2018-12-18Online:2019-03-18
作者简介 About authors
2)陈少林,教授,主要研究方向：地震工程.E-mail:icmcsl@nuaa.edu.cn



摘要
海底地震动的模拟以及海洋工程结构的地震反应分析中,涉及到海水、饱和海床、弹性基岩、结构之间的相互耦合.传统的方法分别采用声波方程描述理想流体、Biot方程描述饱和海床、弹性波方程描述基岩和结构,分别考虑相互之间的耦合,十分不便.本文基于理想流体、固体分别为饱和多孔介质的特殊情形(孔隙率分别为1和0),由饱和多孔介质的Biot方程可退化得到理想流体的声波方程和固体的弹性波方程.然后,以饱和多孔介质方程为基础,经集中质量有限元离散,考虑不同孔隙率的饱和多孔介质之间耦合的一般情形,建立了该耦合情形的求解方法.进一步论证了该一般情形的耦合计算方法可分别退化到流体与固体、流体与饱和多孔介质、固体与饱和多孔介质之间的耦合计算,从而将流体、固体、饱和多孔介质间的耦合问题纳入到统一计算框架,并编制了相应的三维并行分析程序.以P-SV波垂直入射时,半无限层状海水-饱和海床、海水-弹性基岩、海水-饱和海床-弹性基岩三种情形的动力分析为例,采用统一计算框架结合透射边界条件进行求解,并与传递矩阵方法得到的解进行对比,验证了该统一计算框架的有效性以及并行计算的可行性.
关键词： 流固耦合;饱和多孔介质;海洋地震工程;集中质量显式有限元;透射边界;并行计算

Abstract
The simulation of seismic wavefield at seafloor and seismic response of marine structures involves 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 traction and velocity continuity on the interface between porous media with different porosity, 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, examples for analysis of layered water-saturated seabed, water-bedrock, and water-saturated seabed-bedrock semi-infinite systems subjected to plane P-SV wave are given, and the proposed unified framework is verified through comparison between the results obtained through the proposed unified framework combined with tansmitting boundary condition and those obtained through tansfer matrix method.
Keywords：fluid-solid coupling;saturated porous medium;marine earthquake engineering;explicit lumped-mass finite element method;transmitting boundary;parallel computation


PDF (2525KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
陈少林, 柯小飞, 张洪翔. 海洋地震工程流固耦合问题统一计算框架 ¹⁾. 力学学报[J], 2019, 51(2): 594-606 DOI:10.6052/0459-1879-18-333
Chen Shaolin, Ke Xiaofei, Zhang Hongxiang. A UNIFIED COMPUTATIONAL FRAMEWORK FOR FLUID-SOLID COUPLING IN MARINE EARTHQUAKE ENGINEERING ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2019, 51(2): 594-606 DOI:10.6052/0459-1879-18-333

引 言

在海洋地震工程领域,需要关注海底地震波的传播^[1-2]及海底地震作用下结构的安全性问题^[3-6],其响应分析涉及流固耦合问题.总体上,流固耦合从机理上可分为两类：一类是耦合作用仅仅发生在流体和固体介质的交界面上,由界面协调条件引入的^[7-16];另一类是流固两相混合在一起,难以明显分开,其耦合效应通过描述问题的微分方程来体现,如描述饱和多孔介质的Biot方程^[17-18].在海洋工程中,海床土一般为饱和土,通过Biot方程进行描述,属于第二类流固耦合,海水与饱和海床之间的耦合、海水与结构(或基岩)间的耦合、以及饱和海床土与基岩(或结构)之间的耦合属于第一类流固耦合.因此,海洋工程中的流固耦合问题十分复杂,计算较为困难.若采用现行方法,需不同分析模块之间进行交互耦合^[19].

地震波作用下,海水的响应一般通过声波方程进行描述.Komatitsch等^[20]考虑海水与基岩界面法向速度连续和应力连续,建立平衡方程的耦合弱积分形式,通过谱元离散和显式Newmark时间积分,得到海底地震波的谱元模拟方法.李伟华等^[21]考虑理想流体和饱和土层的界面连续条件,建立了水与饱和场地耦合的显式有限元动力分析方法.两种方法均对海水、基岩和饱和土采用不同方程,分别离散,再通过界面条件进行耦合求解,十分不便.

理论上,固体和流体介质均为饱和多孔介质的特殊情形,孔隙率分别为0和1,上述耦合均可在饱和多孔介质理论体系进行描述,如图1所示.基于此,本文从饱和多孔介质的一般情形出发,考虑不同孔隙率的饱和多孔介质之间的耦合,从而将上述所有耦合统一在同一计算框架,避免不同模块之间的交互.另外,采用集中质量显式有限元方法,避免求解大型线性方程组,效率较高,且易于实现并行计算,有望用于大型海底地震波的模拟以及大型海工结构的地震反应分析.对于上述问题,需要采用人工边界模拟无限域的影响,一般采用透射边界^[22-23]和黏弹性边界^[24-25],并分别以自由场^[26]和等效力^[27]作为输入.

图1

新窗口打开|下载原图ZIP|生成PPT
图1介质关系示意图

Fig.1Schematic diagram of media relations

1 基本理论

1.1 一般饱和多孔介质情形

饱和多孔介质的固相平衡方程

(1)$ {\bf{L}}_{\rm s}^{\rm T} {{{\sigma }'}}-(1-\beta ){\bf{L}}_{\rm w}^{\rm T} P + b({{\dot {U}}}-{{\dot {u}}}) = (1-\beta )\rho _{\rm s} {{\ddot {u}}}$

液相平衡方程

(2)$ -\beta {\bf{L}}_{\rm w}^{\rm T} P + b({{\dot {u}}}-{{\dot {U}}}) = \beta \rho _{\rm w} {{\ddot {U}}}$

相容方程(考虑初始孔压和初始体应变为零时)

(3)$$\tau =-\beta P = E_{\rm w} [\beta e^{\rm w} + (1-\beta )e^{\rm s}]$$

其中,${\bf{L}}_{\rm s}$和${\bf{L}}_{\rm w}$为微分算子矩阵;${{ {\sigma }'}}$为有效应力矢量,$\tau$为平均孔压,以拉为正;$P$为孔隙水压,以压为正;${{U}}$和${{u}}$分别为液相和固相的位移矢量,${{\dot {U}}}$和${{\dot {u}}}$为速度,${{\ddot {U}}}$和${{ \ddot {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}$分别为固相和液相的体应变,${{e}}$为固相的应变矢量.${\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 / {\partial x_2 ,\partial/{\partial x_3 }}}} \right) $$

Dirichlet边界条件

(5a)$${u}-\bar{u} ={0} $$

(5b)$${{U}}-{{\bar {U}}} ={{{{0}}}} $$

Neumann边界条件

(5c)$$\hat {n}{{\sigma}}-{{\bar {\sigma }}} = {0} $$

(5d)$$(P-\bar {P}){{n}} ={\bf 0} $$

其中,${{\bar {u}}}$和${{\bar {U}}}$分别为在边界上给定的固相和液相位移.${{\bar {\sigma }}}$和$\bar {P}$分别为边界上固相平均应力和真实孔压的给定值;${{n}}$为沿边界外法线的方向矢量,${{\hat {n}}}$为由方向导数组成的矩阵,其形式如下

(6a)$${{n}} = \left( {n_x ,n_y ,n_z } \right)^{\rm T} $$

(6b)$${{\hat {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 {u}}}_i {{M}}_{{\rm s}i} + {{F}}_i^{\rm s} + {{T}}_i^{\rm s}-{{S}}_i^{\rm s} = {\bf 0}$

(8)$ {{\ddot {U}}}_i {{M}}_{{\rm w}i} + {{F}}_i^{\rm w} + {{T}}_i^{\rm w}-{{S}}_i^{\rm w} = {\bf 0}$

其中,${{M}}_{{\rm s}i}$和${{M}}_{{\rm w}i}$分别为集中在i节点上的固相质量和液相质量;${{F}}_i^{\rm s}$和${{F}}_i^{\rm w}$分别为集中在节点i的固、液相本构力;${{ T}}_i^{\rm s}$和${{T}}_i^{\rm w}$分别为集中在节点i的固、液相黏性阻力;${{S}}_i^{\rm s}$和${{S}}_i^{\rm w}$分别作用在节点i的固、液相边界力.在同一介质内,由于所有位移和应力都连续,通过单元界面作用在它们之上的应力大小相等、方向相反,在单元组装的过程中,内部节点的${{S}}_i^{\rm s}$和${{S}}_i^{\rm w}$均为零.它们由下面的式子计算得到

(9a)$${{M}}_{{\rm s}i} = \sum\limits_e {\sum\limits_{j = 1}^J {\int_{\mbox{ }\Omega ^e} {{{N}}_i^{\rm T} } (1-\beta )\rho _{\rm s} {{N}}_j {\rm d}V} } $$

(9b)$${{M}}_{{\rm w}i} = \sum\limits_e {\sum\limits_{j = 1}^J {\int_{\mbox{ }\Omega ^e} {{{N}}_i^{\rm T} } \beta \rho _{\rm w} {{ N}}_j {\rm d}V} } $$

(9c)$${{T}}_i^{\rm s} = \sum\limits_e {\int_{\mbox{ }var\Omega ^e} {{{N}}_i ^{\rm T}b{{N}}_j ({{\dot {u}}}_j-{{\dot {U}}}_j )} {\rm d}V} $$

(9d)$${{T}}_i^{\rm s} = \sum\limits_e {\int_{\mbox{ }\Omega ^e} {{{N}}_i ^{\rm T}b{{N}}_j ({{\dot {u}}}_j-{{\dot {U}}}_j )} {\rm d}V} $$

(9e)$${{T}}_i^{\rm s} = \sum\limits_e {\int_{{ }\Omega ^e} {{{N}}_i ^{\rm T}b{{N}}_j ({{\dot {u}}}_j-{{\dot {U}}}_j )} {\rm d}V} $$

(9f)$${{S}}_i^{\rm s} = \sum\limits_e {\int_{{S}^e} {{{N}}_i ^{\rm T}} {{\hat {n}\sigma }}{\rm d}S} $$

(9g)$${{S}}_i^{\rm s} = \sum\limits_e {\int_{{S}^e} {{{N}}_i ^{\rm T}} {{\hat {n}\sigma }}{\rm d}S} $$

(9h)$${{S}}_i^{\rm w} = \sum\limits_e {\int_{{S}^e} {{{N}}_i ^{\rm T}} {{n}}\beta P{\rm d}S} $$

若节点i为内部节点(非界面点),${{S}}_i^{\rm s}$和${{S}}_i^{\rm w}$均为零,给定本构关系,可对式(7)、式(8)式采用时步积分求解^[28].这里讨论节点i为两种不同介质的界面点情形,如图2所示.采用隔离体概念,则介质一中界面点i的动力方程由式(7)、式(8)描述,介质二中界面点k的动力方程表示如下(除微分算子和形函数外,其他物理量均带上划线,以区别于介质一)

图2

新窗口打开|下载原图ZIP|生成PPT
图2界面受力示意图

Fig.2Schematic diagram of interfacial force



(10)$ {{\bar { \ddot {u}}}}_k {{\bar {M}}}_{{\rm s}k} + {{\bar {F}}}_k^{\rm s} + {{\bar {T}}}_k^{\rm s}-{{\bar {S}}}_k^{\rm s} = 0$

(11)$ {{ \bar{\ddot {U}}}}_k {{\bar {M}}}_{{\rm w}k} + {{\bar { F}}}_k^{\rm w} + {{\bar {T}}}_k^{\rm w}-{{\bar {S}}}_k^{\rm w} = 0$

其中

(12a)$${{\bar {M}}}_{{\rm s}k} = \sum\limits_e {\sum\limits_{j = 1}^J {\int_{\mbox{ }\Omega ^e} {{{N}}_k^{\rm T} } (1-\bar {\beta })\bar {\rho }_{\rm s} {{N}}_j{\rm d}V} } $$

(12b)$${{\bar {M}}}_{{\rm w}k} = \sum\limits_e {\sum\limits_{j = 1}^J {\int_{\mbox{ }\Omega ^e} {{{N}}_k^{\rm T} } \bar {\beta }\bar {\rho }_{\rm w} {{N}}_j {\rm d}V} } $$

(12c)$${{\bar {F}}}_k^{\rm s} = \sum\limits_e {{{\bar {f}}}_k^{{\rm s}e} } = \sum\limits_e {\int_{\mbox{ }\Omega ^e} {({\bf {L}}_{\rm s} {{N}}_k )^{\rm T}{{\bar {\sigma }}}} {\rm d}V} $$

(12d)$${{\bar {F}}}_k^{\rm w} = \sum\limits_e {{{\bar {f}}}_k^{we} } = \sum\limits_e {\int_{\mbox{ }\Omega ^e} {({\bf {L}}_{\rm w} {{N}}_k )^{\rm T}\bar {\beta }\bar {P}} {\rm d}V} $$

(12e)$${{\bar {T}}}_k^{\rm s} = \sum\limits_e {\int_{\mbox{ }\Omega ^e} {{{N}}_k ^{\rm T}\bar {b}{{N}}_j ({{\dot {u}}}_j-{{\dot {U}}}_j )} {\rm d}V} $$

(12f)$${{\bar {T}}}_k^{\rm w} =-{{\bar {T}}}_k^{\rm s} $$

(12g)$${{\bar {S}}}_k^{\rm s} = \sum\limits_e {\int_{{ S}^e} {{{N}}_k ^{\rm T}} {{ \bar{\hat{n}}\bar {\sigma }}}{\rm d}S} $$

(12h)$${{\bar {S}}}_k^{\rm w} = \sum\limits_e {\int_{{ S}^e} {{{N}}_k ^{\rm T}} {{\bar {n}}}\bar {\beta }\bar {P}{\rm d}S} $$

此时,${{S}}_i^{\rm s}$和${{S}}_i^{\rm w}$为介质二作用在界面点i上的界面力,${{\bar {S}}}_k^{\rm s}$和${{\bar {S}}}_k^{\rm w}$为介质一作用在界面点k上的界面力,它们之间的关系由如下界面连续条件^[30]确定

(13a)$$\sigma _{zz} + \tau = \bar {\sigma }_{zz} + \bar {\tau } $$

(13b)$${\sigma _{zx}} = {\bar \sigma _{zx}}, {\sigma _{zy}} = {\bar \sigma _{zy}} $$

(13c)$$P-\bar {P} = {k}'\beta (\dot {U}_z-\dot {u}_z ) $$

(13d)$$\dot {u}_x =\bar { \dot {u}}_x , \dot {u}_y = \dot {\bar {u}}_y $$

(13e)$$\dot {u}_z = \bar {\dot{u}}_z,~ \beta (\dot {U}_z-\dot {u}_z) = \bar {\beta }( \bar{\dot {U}}_z-\bar {\dot {u}}_z ) $$

其中,式(13c)中的${k}'$为阻滞系数,这里假设两种介质的孔隙完全连通,${k}' = 0$.另外,假设系统的初始值为零,则式(13d)、式(13e)中的速度可换成位移和加速度,本文的推导中采用这一假设.将式(7)和式(10)中x方向的方程相加,考虑到式(13d)和式(13b),以及${ {\hat {n}}} =-{{ \bar{\hat {n}}}}$,有

(14)$ \bar {\ddot{u}}_{{\rm k}x} (m_{{\rm s}i} + \bar {m}_{{\rm s}k} ) + F_{ix}^{\rm s} + \bar {F}_{kx}^{\rm s} + T_{ix}^{\rm s} + \bar {T}_{kx}^{\rm s} = 0$

这即是通常的单元组装后,节点k的x方向平衡方程.同理,可得节点k的y方向固相方程如下

(15)$ \bar {\ddot{u}}_{ky} (m_{{\rm s}i} + \bar {m}_{{\rm s}k} ) + F_{iy}^{\rm s} + \bar {F}_{ky}^{\rm s} + T_{iy}^{\rm s} + \bar {T}_{ky}^{\rm s} = 0$

这里以水平成层介质为例,则方向导数$n_x = n_y =-\bar {n}_x =-\bar {n}_y = 0$,$n_z =-\bar {n}_z = 1$,由式(9h)、式(12h)可知,$S_{ix}^{\rm w} = S_{iy}^{\rm w} = 0$ ,$\bar {S}_{ix}^{\rm w} = \bar {S}_{iy}^{\rm w} = 0$.因此,i点和k点在x,y方向的液相方程如下

(16)$ m_{{\rm w}i} \ddot {U}_{iy} + F_{iy}^{\rm w} + T_{iy}^{\rm w} = 0$

(17)$ m_{{\rm w}i} \ddot {U}_{iy} + F_{iy}^{\rm w} + T_{iy}^{\rm w} = 0$

(18)$$\bar {m}_{{\rm w}k} \bar{\ddot {U}}_{kx} + \bar {F}_{kx}^{\rm w} + \bar {T}_{kx}^{\rm w} = 0 $$

(19)$$\bar {m}_{{\rm w}k} \bar{\ddot {U}}_{ky} + \bar {F}_{ky}^{\rm w} + \bar {T}_{ky}^{\rm w} = 0 $$

取式(7)、式(8)、式(10)、式(11)式的z方向方程相加,考虑到$n_x = n_y =-\bar {n}_x =-\bar {n}_y = 0$,$n_z =-\bar {n}_z = 1$,由式(13a)可知,$S_{iz}^{\rm s} + S_{iz}^{\rm w} + \bar {S}_{kz}^{\rm s} + \bar {S}_{kz}^{\rm w} = 0$,故可得

(20)$$\bar {m}_{{\rm s}k} \bar{\ddot {u}}_{kz} + \bar {m}_{{\rm w}k} \bar {\ddot {U}}_{kz} + m_{{\rm s}i} \ddot {u}_{iz} + m_{{\rm w}i} \ddot {U}_{iz} + F_{iz}^{\rm s} + F_{iz}^{\rm w} + \\ \bar {F}_{kz}^{\rm s} + \bar {F}_{kz}^{\rm w} + T_{iz}^{\rm s} + T_{iz}^{\rm w} + \bar {T}_{kz}^{\rm s} + \bar {T}_{kz}^{\rm w} = 0$$

将式(8)乘以$\bar {\beta }$,加上式(11)乘以$\beta$,取z方向的方程.由式(9h)、式(12h)和式(13c)可知,$\bar {\beta }S_{iz}^{\rm w} + \beta \bar {S}_{kz}^{\rm w} = 0$,故可得

(21)$$\bar {\beta }m_{{\rm w}i} \ddot {U}_{iz} + \beta \bar {m}_{{\rm w}k} \bar {\ddot{U}}_{kz} + \bar {\beta }F_{iz}^{\rm w} + \beta \bar {F}_{kz}^{\rm w} +\\ \bar {\beta }T_{iz}^{\rm w} + \beta \bar {T}_{kz}^{\rm w} = 0 $$

将式(13e)代入式(20)和式(21),消去$\ddot {u}_{iz}$、$\ddot {U}_{iz}$,整理得

(22)$$A_{11} \bar{\ddot {u}}_{kz} + A_{12} \bar{\ddot{U}}_{kz} + B_{11} = 0 $$

(23)$$A_{21} \bar{\ddot {u}}_{kz} + A_{22} \bar{\ddot {U}}_{kz} + B_{22} = 0$$

其中

(24a)$$A_{11} = m_{{\rm s}i} + \bar {m}_{{\rm s}k} + (1-\bar {\beta }/\beta )m_{{\rm w}i} $$

(24b)$$A_{12} = \bar {m}_{{\rm w}k} + ({\bar {\beta }} / \beta )m_{{\rm w}i} $$

(24c)$$A_{21} = \bar {\beta }(1-\bar {\beta }/ \beta )m_{{\rm w}i} $$

(24d)$$A_{22} = \beta \bar {m}_{{\rm w}k} + (\bar {\beta }^2 / \beta)m_{{\rm w}i} + \varepsilon $$

(24e)$$B_{11} = F_{iz}^{\rm s} + F_{iz}^{\rm w} + \bar {F}_{kz}^{\rm s} +\bar {F}_{kz}^{\rm w} + T_{iz}^{\rm s}+\\ T_{iz}^{\rm w} + \bar {T}_{kz}^{\rm s} + \bar {T}_{kz}^{\rm w} $$

(24f)$$B_{22} = \bar {\beta }F_{iz}^{\rm w} + \beta \bar {F}_{kz}^{\rm w} + \bar {\beta }T_{iz}^{\rm w} + \beta \bar {T}_{kz}^{\rm w} $$

式(24d)中,附加的$\varepsilon$为一非常小的量,其作用是避免饱和多孔介质退化到流体和固体情形时,分母为零.联立式(22)、式(23),可解得

(25)$$\bar{\ddot {u}}_{kz} = \frac{A_{22} B_{11}-A_{12} B_{22} }{A_{21} A_{12}-A_{11} A_{22} } $$

(26)$$\bar{\ddot {U}}_{kz} = \frac{A_{11} B_{22}-A_{21} B_{11} }{A_{21} A_{12}-A_{11} A_{22} } $$

式(14)~式(19),以及式(25)和式(26)可采用时步积分求解。本文采用如下显式积分格式

(27)$$ \ddot {U}^p = \frac{1}{\Delta t^2}(U^{(p + 1)}-2U^p + U^{(p-1)}) $$

(28)$$\dot {U}^p = \frac{1}{\Delta t}(U^p-U^{(p-1)}) $$

其中,上标$p$表示$t = p\Delta t$时刻.对式(14)~式(19)以及式(25)和式(26)式采用上述时步积分格式,可得

(29a)$$ \bar {u}_{kx}^{(p + 1)} = 2\bar {u}_{kx}^p-\bar {u}_{kx}^{(p-1)}-\frac{\Delta t^2}{(m_{si} + \bar {m}_{sk} )}(F_{ix}^{{\rm s}p} +\\ \bar {F}_{kx}^{{\rm s}p} + T_{ix}^{{\rm s}p} + \bar {T}_{kx}^{{\rm s}p} ) $$

(29b)$$\bar {u}_{ky}^{(p + 1)} = 2\bar {u}_{ky}^p-\bar {u}_{ky}^{(p-1)}-\frac{\Delta t^2}{(m_{si} + \bar {m}_{sk} )}(F_{iy}^{{\rm s}p} +\\ \bar {F}_{ky}^{{\rm s}p} + T_{iy}^{{\rm s}p} + \bar {T}_{ky}^{{\rm s}p} ) $$

(29c)$$\bar {u}_{kz}^{(p + 1)} = 2\bar {u}_{kz}^p-\bar {u}_{kz}^{(p-1)} + \frac{A_{22} B_{11}^p-A_{12} B_{22}^p }{A_{21} A_{12}-A_{11} A_{22} }\Delta t^2 $$

(29d)$$\bar {U}_{kx}^{(p + 1)} = 2\bar {U}_{kx}^p-\bar {U}_{kx}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm w}k} }(\bar {F}_{kx}^{{\rm w}p} + \bar {T}_{kx}^{{\rm w}p} ) $$

(29e)$$\bar {U}_{ky}^{(p + 1)} = 2\bar {U}_{ky}^p-\bar {U}_{ky}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm w}k} }(\bar {F}_{ky}^{{\rm w}p} + \bar {T}_{ky}^{{\rm w}p} ) $$

(29f)$$\bar {U}_{kz}^{(p + 1)} = 2\bar {U}_{kz}^p-\bar {U}_{kz}^{(p-1)} + \frac{A_{11} B_{22}^p-A_{21} B_{11}^p }{A_{21} A_{12}-A_{11} A_{22} }\Delta t^2 $$

(29g)$$U_{ix}^{(p + 1)} = 2U_{ix}^p-U_{ix}^{(p-1)}-\frac{\Delta t^2}{m_{{\rm w}i} }(F_{ix}^{{\rm w}p} + T_{ix}^{{\rm w}p} ) $$

(29h)$$U_{iy}^{(p + 1)} = 2U_{iy}^p-U_{iy}^{(p-1)}-\frac{\Delta t^2}{m_{{\rm w}i} }(F_{iy}^{{\rm w}p} + T_{iy}^{{\rm w}p} ) $$

得到$\bar {u}_{kx}^{(p + 1)}$,$\bar {u}_{ky}^{(p + 1)}$,$\bar {u}_{kz}^{(p + 1)}$,$\bar {U}_{kz}^{(p + 1)}$后,由式(13d)和式(13e)可进一步求得$u_{ix}^{(p + 1)}$,$u_{iy}^{(p + 1)}$,$u_{iz}^{(p + 1)}$,$U_{iz}^{(p + 1)}$

(29i)$$u_{ix}^{(p + 1)} = \bar {u}_{kx}^{(p + 1)},~ u_{iy}^{(p + 1)} = \bar {u}_{ky}^{(p + 1)},~ u_{iz}^{(p + 1)} = \bar {u}_{kz}^{(p + 1)} $$

(29j)$$U_{iz}^{(p + 1)} = \bar {\beta }/ \beta (\bar {U}_{kz}^{(p + 1)}-\bar {u}_{kz}^{(p + 1)} ) + \bar {u}_{kz}^{(p + 1)} $$

因此,由式(29)可得界面点i和k的响应.

1.2 特殊情形

1.2.1 流体-饱和土情形

考虑饱和土上覆流体情形.此时,图2中介质二为饱和多孔介质,介质一为流体,$\beta = 1$,${{M}}_{{\rm s}i} = {\bf 0}$.考虑无黏的理想流体,动黏度系数为零,则$b = 0$,${{T}}_i^{\rm s} = {\bf 0}$,${{T}}_i^{\rm w} = {\bf 0}$;不存在固相,固相骨架模量可取为零,则${{F}}_i^{\rm s} = {\bf 0}$,${{S}}_i^{\rm s} = {\bf 0}$;由此,式(7)自动满足,式(8)退化为如下理想流体方程

(30)$ {{\ddot {U}}}_i {{M}}_{{\rm w}i} + {{F}}_i^{\rm w}-{{S}}_i^{\rm w} = {\bf 0}$

因此,饱和多孔介质方程可以用来分析理想流体的动力响应.

流体-饱和土界面连续条件为^[28]

(31a)$$ P = \bar {\sigma }_{zz} + \bar {\tau } $$

(31b)$$P = \bar {P} $$

(31c)$$\dot {U}_z = \bar {\beta }(\bar{\dot {U}}_z-\bar{\dot {u}}_z ) + \dot {\bar {u}}_z $$

由动力方程式(10)、式(11)、式(30)以及界面连续条件式(31),按2.1节的方法,可得界面点的运动方程如下

(32a)$$\bar {u}_{kx}^{(p + 1)} = 2\bar {u}_{kx}^p-\bar {u}_{kx}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm s}k} }(\bar {F}_{kx}^{{\rm s}p} + \bar {T}_{kx}^{{\rm s}p} ) $$

(32b)$$\bar {u}_{ky}^{(p + 1)} = 2\bar {u}_{ky}^p-\bar {u}_{ky}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm s}k} }(\bar {F}_{ky}^{{\rm s}p} + \bar {T}_{ky}^{{\rm s}p} ) $$

(32c)$$\bar {u}_{kz}^{(p + 1)} = 2\bar {u}_{kz}^p-\bar {u}_{kz}^{(p-1)} + \frac{A_{22} B_{11}^p-A_{12} B_{22}^p }{A_{21} A_{12}-A_{11} A_{22} }\Delta t^2 $$

(32d)$$\bar {U}_{kx}^{(p + 1)} = 2\bar {U}_{kx}^p-\bar {U}_{kx}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm w}k} }(\bar {F}_{kx}^{{\rm w}p} + \bar {T}_{kx}^{{\rm w}p} ) $$

(32e)$$\bar {U}_{ky}^{(p + 1)} = 2\bar {U}_{ky}^p-\bar {U}_{ky}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm w}k} }(\bar {F}_{ky}^{{\rm w}p} + \bar {T}_{ky}^{{\rm w}p} ) $$

(32f)$$\bar {U}_{kz}^{(p + 1)} = 2\bar {U}_{kz}^p-\bar {U}_{kz}^{(p-1)} + \frac{A_{11} B_{22}^p-A_{21} B_{11}^p }{A_{21} A_{12}-A_{11} A_{22} }\Delta t^2 $$

(32g)$$U_{ix}^{(p + 1)} = 2U_{ix}^p-U_{ix}^{(p-1)}-\frac{\Delta t^2}{m_{{\rm w}i} }F_{ix}^{{\rm w}p} $$

(32h)$$U_{iy}^{(p + 1)} = 2U_{iy}^p-U_{iy}^{(p-1)}-\frac{\Delta t^2}{m_{{\rm w}i} }F_{iy}^{{\rm w}p} $$

(32i)$$U_{iz}^{(p + 1)} = \bar {\beta }(\bar {U}_{kz}^{(p + 1)}-\bar {u}_{kz}^{(p + 1)} ) + \bar {u}_{kz}^{(p + 1)} $$

其中

(33a)$$A_{11} = \bar {m}_{{\rm s}k} + (1-\bar {\beta })m_{{\rm w}i} $$

(33b)$$A_{12} = \bar {m}_{{\rm w}k} + \bar {\beta }m_{{\rm w}i} $$

(33c)$$A_{21} = \bar {\beta }(1-\bar {\beta })m_{{\rm w}i} $$

(33d)$$A_{22} = \bar {m}_{{\rm w}k} + \bar {\beta }^2m_{{\rm w}i} $$

(33e)$$B_{11}^p = F_{iz}^{{\rm w}p} + \bar {F}_{kz}^{{\rm s}p} + \bar {F}_{kz}^{{\rm w}p} + \bar {T}_{kz}^{{\rm s}p} + \bar {T}_{kz}^{{\rm w}p} $$

(33f)$$B_{22}^p = \bar {\beta }F_{iz}^{{\rm w}p} + \bar {F}_{kz}^{{\rm w}p} + \bar {T}_{kz}^{{\rm w}p} $$

将$\beta = 1$,${{M}}_{{\rm s}i} = {\bf 0}$,${{F}}_i^{\rm s} = {\bf 0}$,${{ T}}_i^{\rm s} = {\bf 0}$,${{T}}_i^{\rm w} = {\bf 0}$,${{F}}_i^{\rm s} = {\bf 0}$,${{ S}}_i^{\rm s} = {\bf 0}$代入式(24)和式(29),可分别得到式(33)和式(32).因此流体-饱和土耦合情形是两种不同多孔介质耦合情形的特例,可直接由1.1节中的理论进行分析.

1.2.2 饱和土-干基岩情形

图2中的介质二若为不透水的基岩,不存在液相,则$\bar {\beta } = 0$,${{\bar {M}}}_{{\rm w}k} = {\bf 0}$;液相体积模量和固液之间的黏性力取为零,则${{\bar {F}}}_k^{\rm w} = {\bf 0}$,${{\bar {T}}}_k^{\rm w} = {\bf 0}$,${{\bar {T}}}_k^{\rm s} = {\bf 0}$,${{\bar {S}}}_k^{\rm w} = {\bf 0}$,方程(11)自动满足,式(10)退化为如下干基岩方程

(34)$ {{\ddot { \bar {u}}}}_k {{\bar {M}}}_{{\rm s}k} + {{\bar {F}}}_k^{\rm s}-{{\bar {S}}}_k^{\rm s} = {\bf 0}$

因此,饱和多孔介质方程可以用来分析干基岩的动力响应.

饱和土-干基岩界面条件为^[30]

(35a)$$\sigma _{zz} + \tau = \bar {\sigma }_{zz} $$

(35b)$${\sigma _{zx}} = {\bar \sigma _{zx}},~ {\sigma _{zy}} = {\bar \sigma _{zy}} $$

(35c)$$\dot {u}_x = \bar {\dot {u}}_x ,~ \dot {u}_y = \bar {\dot {u}}_y $$

(35d)$$\dot {u}_z = \bar {\dot {u}}_z ,~ \beta (\dot {U}_z-\dot {u}_z ) = 0 $$

由动力方程式(7)、式(8)、式(34)以及界面连续条件(35)式,按2.1节的方法,可得界面点的运动方程如下

(36a)$$\bar {u}_{kx}^{(p + 1)} = 2\bar {u}_{kx}^p-\bar {u}_{kx}^{(p-1)}-{\Delta t^2}(F_{ix}^{{\rm s}p}+ \\ \bar {F}_{kx}^{{\rm s}p} + T_{ix}^{{\rm s}p} )/{(m_{{\rm s}i} + \bar {m}_{{\rm s}k} )} $$

(36b)$$\bar {u}_{ky}^{(p + 1)} = 2\bar {u}_{ky}^p-\bar {u}_{ky}^{(p-1)}-{\Delta t^2}(F_{iy}^{{\rm s}p}+\\ \bar {F}_{ky}^{{\rm s}p} + T_{iy}^{{\rm s}p} )/{(m_{{\rm s}i} + \bar {m}_{{\rm s}k} )} $$

(36c)$$\bar {u}_{kz}^{(p + 1)} = 2\bar {u}_{kz}^p-\bar {u}_{kz}^{(p-1)} + {B_{11}^p }\Delta t^2/{A_{11} } $$

(36d)$$U_{ix}^{(p + 1)} = 2U_{ix}^p-U_{ix}^{(p-1)}-{\Delta t^2}(F_{ix}^{{\rm w}p} + T_{ix}^{{\rm w}p} )/{m_{{\rm w}i} } $$

(36e)$$U_{iy}^{(p + 1)} = 2U_{iy}^p-U_{iy}^{(p-1)}-{\Delta t^2}(F_{iy}^{{\rm w}p} + T_{iy}^{{\rm w}p} )/{m_{ik} } $$

(36f)$$u_{ix}^{(p + 1)} = \bar {u}_{ix}^{(p + 1)} $$

(36g)$$u_{iy}^{(p + 1)} = \bar {u}_{iy}^{(p + 1)} $$

(36h)$$u_{iz}^{(p + 1)} = \bar {u}_{iz}^{(p + 1)} $$

(36i)$$U_{iz}^{(p + 1)} = \bar {u}_{iz}^{(p + 1)} $$

其中

(37a)$$A_{11} = m_{{\rm s}i} + \bar {m}_{{\rm s}k} + m_{{\rm w}i} $$

(37b)$$B_{11}^p = F_{iz}^{{\rm s}p} + F_{iz}^{{\rm w}p} + \bar {F}_{kz}^{{\rm s}p} + T_{iz}^{{\rm s}p} + T_{iz}^{{\rm w}p} $$

将$\bar {\beta } = {\bf 0}$,${{\bar {M}}}_{{\rm w}k} = {\bf 0}$,${{\bar {F}}}_k^{\rm w} = {\bf 0}$,${{\bar {T}}}_k^{\rm w} = {\bf 0}$,${{\bar {T}}}_k^{\rm s} = {\bf 0}$,${{\bar {S}}}_k^{\rm w} = {\bf 0}$代入式(24)和式(29),可分别得到式(37)和式(36).因此饱和土-干基岩耦合情形也是两种不同多孔介质耦合情形的特例,可直接由2.1节中的理论进行分析.

1.2.3 流体-干基岩情形

此时,图2中介质一为理想流体,介质二为干基岩.按前述参数取值,运动方程分别退化为式(30)和式(34).流体-干基岩界面连续条件为(不透水)

(38a)$$P = \bar {\sigma }_{zz} $$

(38b)$$\dot {U}_z = \bar {\dot{u}}_z $$

由动力方程(30)和(34)以及界面连续条件式(38),可得界面点的运动方程如下

(39a)$$\bar {u}_{kx}^{(p + 1)} = 2\bar {u}_{kx}^p-\bar {u}_{kx}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm s}k} }\bar {F}_{kx}^{{\rm s}p} $$

(39b)$$\bar {u}_{ky}^{(p + 1)} = 2\bar {u}_{ky}^p-\bar {u}_{ky}^{(p-1)}-\frac{\Delta t^2}{\bar {m}_{{\rm s}k} }\bar {F}_{ky}^{{\rm s}p} $$

(39c)$$\bar {u}_{kz}^{(p + 1)} = 2\bar {u}_{kz}^p-\bar {u}_{kz}^{(p-1)} + \frac{B_{11}^p }{A_{11} }\Delta t^2 $$

(39d)$$U_{ix}^{(p + 1)} = 2U_{ix}^p-U_{ix}^{(p-1)}-\frac{\Delta t^2}{m_{{\rm w}i} }F_{ix}^{{\rm w}p} $$

(39e)$$U_{iy}^{(p + 1)} = 2U_{iy}^p-U_{iy}^{(p-1)}-\frac{\Delta t^2}{m_{{\rm w}i} }F_{iy}^{{\rm w}p} $$

(39f)$$U_{iz}^{(p + 1)} = \bar {u}_{iz}^{(p + 1)} $$

其中

(40a)$$A_{11} = \bar {m}_{{\rm s}k} + m_{{\rm w}i} $$

(40b)$$B_{11}^p = F_{iz}^{{\rm w}p} + \bar {F}_{kz}^{{\rm s}p} $$

同样,将$\beta = 1$,${{M}}_{{\rm s}i} = {\bf 0}$,${{F}}_i^{\rm s} = {\bf 0}$,${{T}}_i^{\rm s} = {\bf 0}$,${{T}}_i^{\rm w} = {\bf 0}$,${{S}}_i^{\rm s} = {\bf 0}$,以及$\bar {\beta } = {\bf 0}$,${{\bar {M}}}_{{\rm w}k} = {\bf 0}$,${{\bar {F}}}_k^{\rm w} = {\bf 0}$,${{\bar {T}}}_k^{\rm w} = {\bf 0}$,${{\bar {T}}}_k^{\rm s} = {\bf 0}$,${{\bar {S}}}_k^{\rm w} = {\bf 0}$代入式(24)和式(29),可分别得到式(40)和式(39).因此,流体-干基岩耦合情形同样是两种不同多孔介质耦合情形的特例.

综上可知,流体、固体、饱和多孔介质之间的耦合均可统一在同一理论框架进行分析.

1.3 实施方法

以图3(a)所示的模型为例,考虑海水-饱和海床-基岩水平成层体系在平面P-SV波垂直入射时的反应.对于该问题,可采用Thomson-Haskell传递矩阵方法(transfer matrix method,TMM)求得解析 解^[31-33],称为自由场,做为有限元数值解的验证依据.

图3

新窗口打开|下载原图ZIP|生成PPT
图3模型示意图和计算模型示意

Fig.3Schematic diagram of model



另外,对于复杂地形,或有散射体(如海工结构)存在的情形,需要通过有限元等数值方法进行求解.因此,我们对图3(a)所示的问题采用本文提出的有限元方法求解,计算模型如图3(b)所示,在侧面和底面采用多次透射人工边界^[22],用于模拟无限域.在边界上采用传递矩阵方法得到的自由场作为波动输入.海水-饱和海床-基岩体系中,内部点(除人工边界点和界面点外的其余节点)及界面点的响应计算,均采用统一的饱和多孔介质方程和不同饱和多孔介质间界面点的计算方法,即2.1节中的理论方法.具体实施时,在界面处设置一层"虚单元",该单元的材料参数为零.这样,对于内部点的响应,按饱和多孔介质显式有限元方法计算;对于界面点的响应,没有界面连续条件约束的某些方向(对于文中所述情形,为x和y方向的液相等),其响应计算与内部点一样.对于有连续条件的方向,仅需每时步在节点对(如图3(b)中i和k点)之间传递节点本构力和位移.基于此,编制了相应的三维有限元计算程序,以实现本文方法.

2 算例验证

采用如下线弹性本构

(41a)$${\sigma }'_{ij}-(1-\alpha )\delta _{ij} P = 2\mu e_{ij} + \delta _{ij} \lambda e $$

(41b)$$P =-\alpha Me + M\zeta $$

其中,根据文献[17,18],$\lambda$和$\mu$为固相骨架在排水情形时的拉梅常数,$\alpha$和M则为表征饱和多孔介质压缩性的常数,以压缩模量表示时为

(42)$$\alpha = 1-\frac{E_{\rm b} }{E_{\rm u} }$$

(43)$$M = \frac{E_{\rm u}^2 }{E_{\rm u} \left[n\left(\frac{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×40 m,单元尺寸$\Delta x = 1$ m,时间步距为$\Delta t = 0.000~2$ s,输入图4所示的单位脉冲,$\Delta x \le \lambda _{\min } / 10$,满足精度要求,其中$\lambda _{\min }$为所需模拟的最小波长.所有材料参数见表1.

图4

新窗口打开|下载原图ZIP|生成PPT
图4输入脉冲

Fig.4Input pulse



Table 1
表1
表1材料参数表
Table 1Parameters of material

Media	Prosity yS	仲	Ps/(kg.m_3)	Pw/(kg.m_3)	V	G/MPa	Ew/GPa	M/GPa	a	k0 / ^m2
seawater	1	0	0	1000	0.49	0	2.25	2.25	1	1
saturated soil	0.26	0.001	2000	1000	0.49	83.2	2.25	4.78	0.697	10_7
bedrock	0	0	2500	0	0.2	480	0	—	0	0

新窗口打开|下载CSV

2.1 海水-干土情形

海水-干基岩模型如图5所示.图6和图7分别为P波和SV波入射时的响应.从图中可以看出,传递矩阵方法得到的解析解^[33]与有限元计算的数值解结果完全吻合.对于SV波垂直入射情形,由于没考虑流体黏性,流体中不能传播SV波,A点没有响应,该情形与基岩半空间相同,图中结果也验证了这点.

图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
图7SV波入射时海水-基岩系统的位移响应

Fig.7Displacement of seawater-bedrock system for SV wave incidence

2.2 海水-饱和土情形

将海床考虑为饱和土情形,海水-饱和海床模型如图8所示.图9和图10分别为P波和SV波入射时的响应.从图中可以看出,数值解与解析解完全重合.

图8

新窗口打开|下载原图ZIP|生成PPT
图8海水-饱和土模型

Fig.8Seawater-saturated soil model

图9

新窗口打开|下载原图ZIP|生成PPT
图9P波入射时海水-饱和土系统的位移响应

Fig.9Displacement of seawater-saturated soil system for P wave incidence

图10

新窗口打开|下载原图ZIP|生成PPT
图10SV波入射时海水-饱和土系统的位移响应

Fig.10Displacement of seawater-saturated soil system for SV wave incidence

2.3 海水-饱和海床-基岩情形(并行算例)

大规模计算通常采用并行技术以提高效率.这里,对于海水-饱和海床-基岩模型,我们采用串行和并行(三进程,如图11所示)两种计算方法.两进程的计算区域重叠一层单元(如图11中阴影部分),i点为进程1的内部点,为进程2的边界点,因此在进程1中计算i点的响应,通过MPI通讯协议传递给进程2,同样$j$点的响应在进程2计算,然后传递给进程1.将计算结果与解析解进行对比,见图13和图14.由图可知,除了SV波垂直入射时,饱和土液相位移略有误差外,其余响应均与解析解重合.串行和并行计算时间分别为91 min和34 min,并行计算时间大致为串行的1/3,大大提高了计算效率.

图11

新窗口打开|下载原图ZIP|生成PPT
图11并行计算示意图

Fig.11Schematic diagram of parallel computation

图12

新窗口打开|下载原图ZIP|生成PPT
图12海水-饱和土-基岩模型

Fig.12Seawater-saturated soil-bedrock

图13

新窗口打开|下载原图ZIP|生成PPT
图13P波入射时海水-饱和土-基岩体系的位移响应

Fig.13Displacement of seawater-saturated soil-bedrock system for P wave incidence

图14

新窗口打开|下载原图ZIP|生成PPT
图14SV波入射时海水-饱和土-基岩体系的位移响应

Fig.14Displacement of seawater-saturated soil-bedrock system for SV wave incidence



通过上述三个算例,验证了本文统一计算框架的有效性,以及并行计算的可行性.

3 结 论

本文将理想流体、固体、饱和多孔介质之间的耦合分析纳入到统一计算框架,建立了集中质量显式有限元求解方法,并编制了相应的三维并行分析程序.分析了半无限海水-饱和海床、海水-弹性基岩、海水-饱和海床-弹性基岩三种情形在P-SV波垂直入射时的动力响应,通过与传递矩阵方法得到的结果进行对比,验证了该统一计算框架的有效性以及并行计算的可行性.

文中的算例采用线弹性本构,但公式推导只涉及本构力,可以适用于非线性情形.界面点的精度与内域点一致,取决于有限元空间离散精度和时步积分格式的精度,可以采用其他格式满足所需精度要求.文中算例未出现失稳现象,但耦合计算格式的稳定性需另做讨论.另外,本文只考虑了水平成层情形,通过考虑界面的方向导数,该方法可以推广到复杂海底地形情形,结合并行技术,可以模拟大规模海底地震动问题以及海水-海床-结构体系的地震响应问题,也可为海洋声学及海洋地震勘探提供正演分析方法.

The authors have declared that no competing interests exist.

作者已声明无竞争性利益关系。

参考文献 View Option 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[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 layer
Bulletin of Seismological Society of America, 2012,102(6):2420-2435
DOIURL
ABSTRACT Feasibility studies on simulating seismic-wave propagation in media for a suboceanic earthquake, including both land and ocean-bottom topographies and a seawater layer, are scarce. Some of the conventional staggered-grid finite-difference method (FDM) simulations use a simplified structure model without a seawater layer and with a flat ocean-bottom topography. In this study, we apply our heterogeneity, oceanic layer, and topography scheme (HOT)-FDM and a 3D structure model including land and ocean-bottom topographies, a seawater layer, and a fluid-solid boundary condition to an aftershock (M-w 5.8) of the 2009 Suruga Bay earthquake. We attempt to reproduce observations at seismic stations near the coast and then simulate waveforms in the ocean-bottom stations. Our results show that a large difference between the cases with and without topographies can be seen in the coda part after the S wave in the simulated waveforms in terms of amplitudes and elongations. The synthetic waveforms in the model with topographies are in agreement with the observed waveforms. Our results also show that a significant difference in the amplification of the coda part between the cases with and without the seawater layer can be found at the ocean-bottom stations. This coda part is the S-and Rayleigh-wave propagation associated with the ocean and the underlying sediment layers. Our results show that a realistic model with topographies and a seawater layer is needed in FDM simulations in order to precisely reproduce observed waveforms or predict seismic motion for a suboceanic earthquake.

[2] Petukhin A, Iwata T, Kagawa T . Study on the effect of the oceanic water layer on strong ground motion simulations
Earth Planets Space, 2010,62:621-630
DOIURL
We studied the effect of the oceanic water layer on strong ground motion simulations. Source faults of subduction zone earthquakes, such as the Nankai-Tonankai earthquake, West Japan, are situated in the offshore area, under a thick water layer. The necessity of including the oceanic water layer in the velocity model for simulations employing the finite difference method is debated by many researchers, and consideration given to the possibility of neglecting by this layer, which would reduce the computation time and stabilize the calculations. Although the oceanic water layer has a low velocity and density, it can affect surface wave generation. In this study, for demonstration purposes, we calculated and compared strong ground motions from three source fault models, placed into the boundary between the crust and subducting plate, where source rupture of the Tonankai earthquake is expected. Simulations were made for two realistic three-dimensional velocity models: without and with the oceanic water layer. The model without the oceanic water layer was constructed simply by subtracting the depth of the oceanic layer from the depth of all velocity interfaces under the ocean. This procedure keeps the thickness of the layers (oceanic sediments, surface low-velocity layer, upper crust, and lower crust) the same as in the model with the oceanic layer and reduces simulation errors. Simulations were made for the set of sites on a line across the subduction zone and directed to the Osaka basin. The results show that the water layer has a strong effect on the fundamental mode of the Rayleigh wave, which can be generated by the shallow (approx. 5 km) source. Considering that all asperities of the expected Nankai-Tonankai earthquake are deep (> 10 km), we conclude that the effect of the water layer can be neglected for ground motions at land sites.

[3] Jianhong Y , Seismic response of poroelastic seabed and composite breakwater under strong earthquake loading
Bulletin of Earthquake Engineering, 2012,10:1609-1633
DOIURL
Abstract 61 ” as the governing equation for porous seabed foundation, the seismic response of a composite breakwater and its porous seabed foundation under the seismic wave recorded in the Japan 311 off the pacific coast of Tohoku earthquake (M L =029.0) is investigated using a FEM numerical model. The numerical results indicate that the seismic response of composite breakwater is very strong in the earthquake process. The amplification of the input seismic wave occurs both in seabed foundation and composite breakwater; and this amplification is positively related to the buried depth of points. The horizontal seismic response is much strong than the vertical seismic response. The seismic wave induced excess pore pressure and effective stresses in seabed foundation vibrates; the vibration amplitude is also positively related to the buried depth of points. Under strong seismic loading, the surface region of seabed foundation could liquefy. The parametric study shows that the young’s modulus of seabed foundation has significant effect on the seismic response of composite breakwater.

[4] Cheng XS, Xu WW, Yue CQ , et al. Seismic response of fluid-structure interaction of undersea tunnel during bidirction earthquake
Ocean Engineering, 2014,75:64-70


[5] Jin HL, Seo SI, Mun HS . Seismic behaviors of a floating submerged tunnel with a rectangular cross-section
Ocean Engineering, 2016,127:32-47
DOIURL
61Two-dimensional seismic behaviors of a submerged floating tunnel (SFT) are examined.61The effects of seawater upon the seismic behaviors of an SFT system are investigated.61The depth of the ocean water is one of the important factors for the SFT system.61Behaviors of the SFT system are affected by the location of the tunnel in the water.61The energy absorption by seabed influences the dynamic responses of the SFT system.

[6] Farhat C, Lesoinne M , LeTallec P. Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: Momentum and energy conservation, optimal discretization and application to aeroelasticity
Computer Methods in Applied Mechanics & Engineering, 1998,157(1):95-114


[7] Farhat C, Lesoinne M . Two effcient staggered algorithms for serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems
Computer Methods in Applied Mechanics & Engineering, 2000,182:499-515
DOIURL
Partitioned procedures and staggered algorithms are often adopted for the solution of coupled fluid/structure interaction problems in the time domain. In this paper, we overview two sequential and parallel partitioned procedures that are popular in computational nonlinear aeroelasticity, and address their limitation in terms of accuracy and numerical stability. We propose two alternative serial and parallel staggered algorithms for the solution of coupled transient aeroelastic problems, and demonstrate their superior accuracy and computational efficiency with the flutter analysis of the AGARD Wing 445.6. We contrast our results with those computed by other investigators and validate them with experimental data.

[8] Farhat C, Zee KGVD, Geuzaine P . Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity
Computer Methods in Applied Mechanics & Engineering, 2006,195(17):1973-2001
DOIURL
A methodology for designing formally second-order time-accurate and yet loosely-coupled partitioned procedures for the solution of nonlinear fluid tructure interaction (FSI) problems on moving grids is presented. Its key components are a fluid time-integrator that is provably second-order time-accurate on moving grids, the midpoint rule for advancing in time the solution of the structural dynamics equations of motion, a second-order structure predictor for bypassing the inner-iterations encountered in strongly-coupled solution procedures, and a carefully designed algorithm for time-integrating the motion of the fluid-mesh. Following this methodology, two different loosely-coupled schemes are constructed for the solution of transient nonlinear FSI problems and proved to be second-order time-accurate. Three-dimensional numerical results pertaining to the simulation of the aeroelastic response to a gravity excitation of a complete F-16 configuration are also presented. In addition to confirming the theoretical results discussed in this paper, these numerical results highlight a very stable behavior of the designed loosely-coupled partitioned procedures.

[9] Bathe KJ, Zhang H . Finite element developments for general fluid flows with structural interactions
International Journal for Numerical Methods in Engineering, 2010,60(1):213-232
DOIURL
The objective in this paper is to present some developments for the analysis of Navier-Stokes incompressible and compressible fluid flows with structural interactions. The incompressible fluid is discretized with a new solution approach, a flow-condition-based interpolation finite element scheme. The high-speed compressible fluids are solved using standard finite volume methods. The fluids are fully coupled to general structures that can undergo highly non-linear response due to large deformations, inelasticity, contact and temperature. Particular focus is given on the scheme used to couple the fluid media with the structures. The fluids can also be modelled as low-speed compressible or slightly compressible media, which are important models in engineering practice. Some solutions obtained using ADINA are presented to indicate the analyses that can be performed

[10] Degroote J, Haelterman R, Annerel S , et al. Performance of partitioned procedures in fluid-structure interaction
Computers & Structures, 2010,88(7):446-457
DOIURL
Partitioned simulations of fluid–structure interaction can be solved for the interface’s position with Newton–Raphson iterations but obtaining the exact Jacobian is impossible if the solvers are “black boxes”. It is demonstrated that only an approximate Jacobian is needed, as long as it describes the reaction to certain components of the error on the interface’s position. Based on this insight, a quasi-Newton coupling algorithm with an approximation for the inverse of the Jacobian (IQN-ILS) has been developed and compared with a monolithic solver in previous work. Here, IQN-ILS is compared with other partitioned schemes such as IBQN-LS, Aitken relaxation and Interface-GMRES(R).

[11] Hou G, Wang J, Layton A . Numerical methods for fluid-structure Interaction-A review
Communications in Computational Physics, 2012,12(2):337-377
DOIURL
The interactions between incompressible fluid flows and immersed structures are nonlinear multi-physics phenomena that have applications to a wide range of scientific and engineering disciplines. In this article, we review representative numerical methods based on conforming and non-conforming meshes that are currently available for computing fluid-structure interaction problems, with an emphasis on some of the recent developments in the field. A goal is to categorize the selected methods and assess their accuracy and efficiency. We discuss challenges faced by researchers in this field, and we emphasize the importance of interdisciplinary effort for advancing the study in fluid-structure interactions.

[12] Habchi C, Russeil S, Bougeard D , et al. Partitioned solver for strongly coupled fluid-structure interaction
Computers & Fluids, 2013,71(1):306-319
DOI
In this work a fluid-structure interaction solver is developed in a partitioned approach using block Gauss-Seidel implicit scheme. Finite volume method is used to discretize the fluid flow problem on a moving mesh in an arbitrary Lagrangian-Eulerian formulation and by using an adaptive time step. The pressure-velocity coupling is performed by using the PIMPLE algorithm, a combination of both SIMPLE and PISO algorithms, which permits the use of larger time steps in a moving mesh. The structural elastic deformation is analyzed in a Lagrangian formulation using the St. Venant-Kirchhoff constitutive law, for non-linear large deformations. The solid structure is discretized by the finite volume method in an iterative segregated approach. The automatic mesh motion solver is based on Laplace smoothing equation with variable mesh diffusion. The strong coupling between the different solvers and the equilibrium on the fluid-structure interface are achieved by using an iterative implicit fixed-point algorithm with dynamic Aitken's relaxation method. The solver, which is called vorflexFoam, is developed using the open source C++ library OpenFOAM. The solver is validated on two different benchmarks largely used in the open literature. In the first one the structural deformation is induced by incompressibility. The second benchmark consists on a vortex excited elastic flap in a Von Karman vortex street. Finally, a more complex case is studied including two elastic flaps immersed in a pulsatile flow. The present solver detects accurately the interaction between the complex flow structures generated by the flaps and the effect of the flaps oscillations between each other. (c) 2012 Elsevier Ltd. All rights reserved.

[13] Mehl M, Uekermann B, Bijl H , et al. Parallel coupling numerics for partitioned fluid-structure interaction simulations
Computers & Mathematics with Applications, 2016,71(4):869-891
DOIURL
Within the last decade, very sophisticated numerical methods for the iterative and partitioned solution of fluid tructure interaction problems have been developed that allow for high accuracy and very complex scenarios. The combination of these two aspects ccuracy and complexity emands very high computational grid resolutions and, thus, high performance computing methods designed for massively parallel hardware architectures. For those architectures, currently used coupling methods, which mainly work with a staggered execution of the fluid and the structure solver, i.e., the execution of one solver after the other in every outer iteration, lead to severe load imbalances: if the flow solver, e.g., scales on a very large number of processors but the structural solver does not due to its limited amount of data and required operations, almost all processors assigned to the coupled simulations are idle during the execution of the structure solver. We propose two new iterative coupling methods that allow for the simultaneous execution of flow and structure solvers. In both cases, we show that pure fixed-point iterations based on the parallel execution of the solvers do not lead to good results, but the combination of parallel solver execution and so-called quasi-Newton methods yields very efficient and robust methods. Those methods are known to be very efficient also for the stabilization of critical scenarios solved with the standard staggered solver execution. We demonstrate the competitive convergence of our methods for various established benchmark scenarios. Both methods are perfectly suited for use with black-box solvers because the quasi-Newton approach uses solely input and output information of the solvers to approximate the effect of the unknown Jacobians that would be required in a standard Newton solver.

[14] Bungartz HJ, Lindner F, Gatzhammer B , et al. preCICE-A fully parallel library for multi-physics surface coupling
Computers & Fluids, 2016,141:250-258
DOIURL
In the emerging field of multi-physics simulations, we often face the challenge to establish new connections between physical fields, to add additional aspects to existing models, or to exchange a solver for one of the involved physical fields. If in such cases a fast prototyping of a coupled simulation environment is required, a partitioned setup using existing codes for each physical field is the optimal choice. As accurate models require also accurate numerics, multi-physics simulations typically use very high grid resolutions and, accordingly, are run on massively parallel computers. Here, we face the challenge to combine flexibility with parallel scalability and hardware efficiency. In this paper, we present the coupling tool preCICE which offers the complete coupling functionality required for a fast development of a multi-physics environment using existing, possibly black-box solvers. We hereby restrict ourselves to bidirectional surface coupling which is too expensive to be done via file communication, but in contrast to volume coupling still a candidate for distributed memory parallelism between the involved solvers. The paper gives an overview of the numerical functionalities implemented in preCICE as well as the user interfaces, i.e., the application programming interface and configuration options. Our numerical examples and the list of different open-source and commercial codes that have already been used with preCICE in coupled simulations show the high flexibility, the correctness, and the high performance and parallel scalability of coupled simulations with preCICE as the coupling unit.

[15] Banks JW, Henshaw WD, Kapila AK , et al. An added-mass partition algorithm for fluid-structure interactions of compressible fluids and nonlinear solids
Journal of Computational Physics, 2016,305(C):1037-1064
DOIURL
We describe an added-mass partitioned (AMP) algorithm for solving fluid–structure interaction (FSI) problems involving inviscid compressible fluids interacting with nonlinear solids that undergo large rotations and displacements. The computational approach is a mixed Eulerian–Lagrangian scheme that makes use of deforming composite grids (DCG) to treat large changes in the geometry in an accurate, flexible, and robust manner. The current work extends the AMP algorithm developed in Banks et al. [1] for linearly elasticity to the case of nonlinear solids. To ensure stability for the case oflightsolids, the new AMP algorithm embeds an approximate solution of a nonlinear fluid–solid Riemann (FSR) problem into the interface treatment. The solution to the FSR problem is derived and shown to be of a similar form to that derived for linear solids: the state on the interface being fundamentally an impedance-weighted average of the fluid and solid states. Numerical simulations demonstrate that the AMP algorithm is stable even for light solids when added-mass effects are large. The accuracy and stability of the AMP scheme is verified by comparison to an exact solution using the method of analytical solutions and to a semi-analytical solution that is obtained for a rotating solid disk immersed in a fluid. The scheme is applied to the simulation of a planar shock impacting a light elliptical-shaped solid, and comparisons are made between solutions of the FSI problem for a neo-Hookean solid, a linearly elastic solid, and a rigid solid. The ability of the approach to handle large deformations is demonstrated for a problem of a high-speed flow past a light, thin, and flexible solid beam.

[16] Basting S, Quaini A, Glowinski R . Extended ALE Method for fluid-structure interaction problems with large structural displacements
Journal of Computational Physics, 2016,331(C):312-336
DOIURL
Standard Arbitrary Lagrangian-Eulerian (ALE) methods for the simulation of fluid-structure interaction (FSI) problems fail due to excessive mesh deformations when the structural displacement is large. We propose a method that successfully deals with this problem, keeping the same mesh connectivity while enforcing mesh alignment with the structure. The proposed Extended ALE Method relies on a variational mesh optimization technique, where mesh alignment with the structure is achieved via a constraint. This gives rise to a constrained optimization problem for mesh optimization, which is solved whenever the mesh quality deteriorates. The performance of the proposed Extended ALE Method is demonstrated on a series of numerical examples involving 2D FSI problems with large displacements. Two way coupling between the fluid and structure is considered in all the examples. The FSI problems are solved using either a Dirichlet-Neumann algorithm, or a Robin-Neumann algorithm. The Dirichlet-Neumann algorithm is enhanced by an adaptive relaxation procedure based on Aitken's acceleration. We show that the proposed method has excellent performance in problems with large displacements, and that it agrees well with a standard ALE method in problems with mild displacement.

[17] Biot MA . Theory of propagation of elastic waves in a fluid-saturated porous solid
Acoust Soc Am, 1956,28:168-191
[本文引用: 1]

[18] Biot MA . Mechanics of deformation and acoustic propagation in porous media
Journal of Applied Physics, 1962,33(4):1482-1498
DOIURL [本文引用: 1]
A unified treatment of the mechanics of deformation and acoustic propagation in porous media is presented, and some new results and generalizations are derived. The writer's earlier theory of deformation of porous media derived from general principles of nonequilibrium thermodynamics is applied. The fluid﹕olid medium is treated as a complex physicalヽhemical system with resultant relaxation and viscoelastic properties of a very general nature. Specific relaxation models are discussed, and the general applicability of a correspondence principle is further emphasized. The theory of acoustic propagation is extended to include anisotropic media, solid dissipation, and other relaxation effects. Some typical examples of sources of dissipation other than fluid viscosity are considered.

[19] Komatitsch D, Barnes C, Tromp J . Wave propagation near a fluid-solid interface: A spectral-element approach
Geophysics, 2000,65(2):623-631
[本文引用: 1]

[20] Link G, Kaltenbacher M, Breuer M , et al. A 2D finite-element scheme for fluid-solid-acoustic interactions and its application to human phonation
Computer Methods in Applied Mechanics & Engineering, 2009,198(41):3321-3334
DOIURL [本文引用: 1]
We present a recently developed approach for the modeling of fluid–solid–acoustic interaction problems. For the efficient numerical solution of the coupled three-field problem we apply the finite-element method. The mechanical and the acoustic fields are approximated by a standard Galerkin scheme. A residual-based stabilization method is chosen for the fluid field. The interaction of the Eulerian fluid field with the Lagrangian mechanical field is based on the Arbitrary-Lagrangian–Eulerian (ALE) method and is iteratively coupled in a strong sense. The solid–acoustic interaction is based on continuum mechanics, and the fluid–acoustic coupling on Lighthill’s analogy. The new steps of our scheme are verified through validation examples. Finally, a fluid–solid–acoustic simulation of the human phonation process is presented based on a realistic model.

[21] 李伟华 . 考虑水-饱和土场地-结构耦合时的沉管隧道地震反应分析
防灾减灾工程学报, 2010,30(6):607-613
DOIURL [本文引用: 1]
用理想流体介质模拟水层、流体饱和多孔介质模拟饱和土地层,在理 想流体介质与流体饱和多孔介质相连接边界的连续条件基础上,结合已有的对理想流体介质、流体饱和多孔介质进行动力反应分析的显式有限元方法,并考虑地层与 结构的动力相互作用,建立了进行水与场地、结构耦合动力分析的方法.利用该方法对沉管隧道的地震响应进行了研究,重点分析了水深、地质条件等因素对沉管隧 道地震反应的影响,并从中得出了一些可供相关人员进行沉管隧道抗震分析时参考的结论.
( Li Weihua . Seismic Response analysis of immersed tube tunnels considering water saturated soil site structure coupling
Journal of Disaster Prevention and Mitigation Engineering, 2010,30(6):607-613(in Chinese))
DOIURL [本文引用: 1]
用理想流体介质模拟水层、流体饱和多孔介质模拟饱和土地层,在理 想流体介质与流体饱和多孔介质相连接边界的连续条件基础上,结合已有的对理想流体介质、流体饱和多孔介质进行动力反应分析的显式有限元方法,并考虑地层与 结构的动力相互作用,建立了进行水与场地、结构耦合动力分析的方法.利用该方法对沉管隧道的地震响应进行了研究,重点分析了水深、地质条件等因素对沉管隧 道地震反应的影响,并从中得出了一些可供相关人员进行沉管隧道抗震分析时参考的结论.

[22] 廖振鹏 . 工程波动理论导论. 第2版. 北京: 科学出版社, 2002: 136-285
[本文引用: 1]

( Liao Zhenpeng. Introduction to Wave Motion Theories in Engineering(2nd edn). Beijing: Science Press, 2002: 136-285(in Chinese))
[本文引用: 1]

[23] 邢浩洁, 李鸿晶 . 透射边界条件在波动谱元模拟中的实现:二维波动
力学学报, 2017,49(4):894-906
URL
将邢浩洁和李鸿晶提出的多次透射公式(multi-transmitting formula,MTF)的谱元格式应用于均匀介质中线弹性SH波动问题的谱元模拟.假定紧邻人工边界的一层谱单元为具有直线边界的四边形单元,以保证每个人工边界节点都唯一对应一条指向内域的离散网格线.人工边界节点在某时刻的位移由该离散网格线上的节点在前若干时刻的位移确定,按照MTF谱元格式进行计算.通过平面波以一定角度传播的外源问题算例和点源脉冲自由扩散的内源问题算例,验证了方法的可行性以及对实际复杂波动问题的适用性.通过不同类型初值问题算例,在时域内分析了插值多项式阶次、人工波速和透射阶次三个参数对反射误差的影响.结果表明:插值多项式阶次较高的格式会表现出更好的精度,但总体上对反射误差的影响较小;人工波速对反射误差具有显著影响,当人工波速小于介质物理波速时反射误差较大,而当人工波速等于或稍大于介质物理波速时反射误差处于较低水平;透射阶次对反射误差具有决定性影响,表现在不失稳的情形下提高透射阶次能够迅速降低反射误差,但内源问题从三阶MTF开始出现飘移失稳,外源问题从二阶MTF开始出现轻微的飘移失稳.
( Xing Haojie, Li Hongjing . Implementation of multi-transmitting boundary condition for wave motion simulation by spectral element method: Two dimension case
Chinese Journal of Theoretical and Applied Mechanics, 2017,49(4):894-906 (in Chinese))
URL
将邢浩洁和李鸿晶提出的多次透射公式(multi-transmitting formula,MTF)的谱元格式应用于均匀介质中线弹性SH波动问题的谱元模拟.假定紧邻人工边界的一层谱单元为具有直线边界的四边形单元,以保证每个人工边界节点都唯一对应一条指向内域的离散网格线.人工边界节点在某时刻的位移由该离散网格线上的节点在前若干时刻的位移确定,按照MTF谱元格式进行计算.通过平面波以一定角度传播的外源问题算例和点源脉冲自由扩散的内源问题算例,验证了方法的可行性以及对实际复杂波动问题的适用性.通过不同类型初值问题算例,在时域内分析了插值多项式阶次、人工波速和透射阶次三个参数对反射误差的影响.结果表明:插值多项式阶次较高的格式会表现出更好的精度,但总体上对反射误差的影响较小;人工波速对反射误差具有显著影响,当人工波速小于介质物理波速时反射误差较大,而当人工波速等于或稍大于介质物理波速时反射误差处于较低水平;透射阶次对反射误差具有决定性影响,表现在不失稳的情形下提高透射阶次能够迅速降低反射误差,但内源问题从三阶MTF开始出现飘移失稳,外源问题从二阶MTF开始出现轻微的飘移失稳.

[24] 谷音, 刘晶波, 杜修力 . 三维一致粘弹性人工边界及等效粘弹性边界单元
工程力学, 2007,24(12):31-37
DOIURL
基于粘弹性人工边界推导了三维一致粘弹性人工边界单元的刚度及阻尼矩阵，利用单元矩阵等效原理采用普通有限单元构造了等效粘弹性边界单元来模拟三维粘弹性边界。均匀半空间算例与成层半空间算例证明三维粘弹性边界单元具有与集中粘弹性人工边界相近的精度，并且施加更为简便。
( Gu Lin, Liu Jingbo, Du Xiuli . Three-dimensional uniform viscoelastic artificial boundary and equivalent viscoelastic boundary element
Journal of Engineering Mechanics, 2007,24(12):31-37(in Cinese))
DOIURL
基于粘弹性人工边界推导了三维一致粘弹性人工边界单元的刚度及阻尼矩阵，利用单元矩阵等效原理采用普通有限单元构造了等效粘弹性边界单元来模拟三维粘弹性边界。均匀半空间算例与成层半空间算例证明三维粘弹性边界单元具有与集中粘弹性人工边界相近的精度，并且施加更为简便。

[25] 刘晶波, 宝鑫, 谭辉 等. 波动问题中流体介质的动力人工边界
力学学报, 2017,49(6):1418-1427
DOIURL
无限域流体介质的波动辐射效应是影响海域工程动力反应的重要因素,人工边界是实现此类开放系统近场波动问题数值分析的有效方法.基于位移格式的流体波动理论推导开放域流体介质的人工边界,分别给出一维、二维和三维空间中平面波、柱面波和球面波的流体介质动力人工边界条件,其中一维平面波动人工边界为经典的黏性边界,二维柱面波、三维球面波的人工边界处节点应力与节点速度和加速度成正比,可等效为由阻尼与质量系统构成的人工边界条件.讨论相应的数值模拟技术,给出流体介质动力人工边界在ANSYS软件平台的具体实现方法.近场流体介质动力反应问题的算例表明,所发展的流体动力人工边界对于轴对称波动与非轴对称波动在近场有限域截断处的透射吸收作用的模拟计算精度均较为良好,说明此流体介质人工边界具有较高的可靠性与实用性.所发展的流体介质动力人工边界可较为方便地与大型商用有限元软件结合,可为包括海域地形和海岛在内的海域工程的动力分析提供一定的方法借鉴.
( Liu Jingbo , BaoXin, Tan Hui, et al. Dynamical artificial boundary for fluid medium in wave motion problems
Chinese Journal of Theoretical and Applied Mechanics, 2017,49(6):1418-1427 (in Chinese))
DOIURL
无限域流体介质的波动辐射效应是影响海域工程动力反应的重要因素,人工边界是实现此类开放系统近场波动问题数值分析的有效方法.基于位移格式的流体波动理论推导开放域流体介质的人工边界,分别给出一维、二维和三维空间中平面波、柱面波和球面波的流体介质动力人工边界条件,其中一维平面波动人工边界为经典的黏性边界,二维柱面波、三维球面波的人工边界处节点应力与节点速度和加速度成正比,可等效为由阻尼与质量系统构成的人工边界条件.讨论相应的数值模拟技术,给出流体介质动力人工边界在ANSYS软件平台的具体实现方法.近场流体介质动力反应问题的算例表明,所发展的流体动力人工边界对于轴对称波动与非轴对称波动在近场有限域截断处的透射吸收作用的模拟计算精度均较为良好,说明此流体介质人工边界具有较高的可靠性与实用性.所发展的流体介质动力人工边界可较为方便地与大型商用有限元软件结合,可为包括海域地形和海岛在内的海域工程的动力分析提供一定的方法借鉴.

[26] 赵宇昕, 陈少林 . 关于传递矩阵法分析饱和成层介质响应问题的讨论
力学学报, 2016,48(5):1145-1158
DOIMagsci [本文引用: 1]
<p>水平成层土体的地震响应分析（自由场分析）是地震工程领域地震波散射问题的前提基础，由于饱和多孔方程的复杂性，以往的研究大多集中于干土情形，对于饱和土情形的研究相对较少.而实际工程中，地下水位以下，土体孔隙中充满流体，应考虑饱和多孔介质模型.基于Biot多孔介质模型，考虑饱和土中固液相对运动引起的衰减，采用Thomson-Haskell传递矩阵方法得到了饱和成层土体在地震波入射情形时的稳态反应，经傅里叶反变换，可得到时域暂态反应.通过SV波从基岩入射至上覆饱和土层的数值算例，验证了该方法的有效性.发现和初步阐明了计算中出现的两类违背因果律（即响应先于输入）的现象：（1）当SV波入射角度大于导致基岩中反射P波为非均匀波的临界角时，会使得计算结果违背因果律.因此，当入射角超过临界角时，非均匀波的表示尚需进一步完善；（2）由于P2波的衰减，当与稳态波衰减有关的渗透率、土层厚度、入射波频率等参数导致衰减系数超过计算机表示精度时，会出现结果违背因果律现象，并据此得到了满足因果律的参数范围，该范围可作为实际计算时的一个上界.该工作为采用传递矩阵法分析水平饱和土层自由场响应提供了指导依据，且地下水位以上可采用干土模型，水位以下采用饱和土模型，更符合实际情形.</p>
( Zhao Yuxin, Chen Shaolin . Discussion on the matrix propagator method to analyze the response od saturated layered media
Chinese Journal of Theoretical and Applied Mechanics, 2016,48(5):1145-1158 (in Chinese))
DOIMagsci [本文引用: 1]
<p>水平成层土体的地震响应分析（自由场分析）是地震工程领域地震波散射问题的前提基础，由于饱和多孔方程的复杂性，以往的研究大多集中于干土情形，对于饱和土情形的研究相对较少.而实际工程中，地下水位以下，土体孔隙中充满流体，应考虑饱和多孔介质模型.基于Biot多孔介质模型，考虑饱和土中固液相对运动引起的衰减，采用Thomson-Haskell传递矩阵方法得到了饱和成层土体在地震波入射情形时的稳态反应，经傅里叶反变换，可得到时域暂态反应.通过SV波从基岩入射至上覆饱和土层的数值算例，验证了该方法的有效性.发现和初步阐明了计算中出现的两类违背因果律（即响应先于输入）的现象：（1）当SV波入射角度大于导致基岩中反射P波为非均匀波的临界角时，会使得计算结果违背因果律.因此，当入射角超过临界角时，非均匀波的表示尚需进一步完善；（2）由于P2波的衰减，当与稳态波衰减有关的渗透率、土层厚度、入射波频率等参数导致衰减系数超过计算机表示精度时，会出现结果违背因果律现象，并据此得到了满足因果律的参数范围，该范围可作为实际计算时的一个上界.该工作为采用传递矩阵法分析水平饱和土层自由场响应提供了指导依据，且地下水位以上可采用干土模型，水位以下采用饱和土模型，更符合实际情形.</p>

[27] 刘晶波, 谭辉, 宝鑫 等. 土-结构动力相互作用分析中基于人工边界子结构的地震波动输入方法
力学学报, 2018,50(1):32-43
DOIURL [本文引用: 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 boundaries
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(1):32-43 (in Chinese))
DOIURL [本文引用: 1]
数值模拟是解决土-结构动力相互作用问题的重要手段,而合理地实现地震波动输入直接影响地震作用下土-结构动力相互作用问题数值模拟的精度。波动法是目前常用的地震动输入方法之一,该方法将输入地震动转化为人工边界上的等效荷载,相较于其他地震动输入方法,波动法模拟精度高,但实施上相对复杂。从有限元模型入手,推导了采用波动法确定等效输入地震荷载的另一种形式,以此为基础,提出了一种在人工边界上实现地震动输入的新方法。新方法通过对土-结构有限元模型中由包含人工边界节点的单元组成的子结构施加自由场位移时程并进行动力分析,直接获得可实现地震波动有效输入的等效荷载,然后将等效输入地震荷载施加在土-结构模型的人工边界节点上,从而完成土-结构动力相互作用问题分析中地震动输入和地震反应计算。与原有波动法相比,新方法避免了需分别计算人工边界上自由场应力和由引入人工边界条件带来的附加力,以及需要根据不同人工边界面确定荷载的作用方向等较为复杂的处理过程,具有等效地震荷载计算简便、地震动输入过程更易于实施的特点。采用竖直入射和斜入射地震波动作用下的弹性半空间和成层半空间地震反应算例验证了新方法的有效性。

[28] 陈少林, 廖振鹏, 陈进 . 两相介质近场波动模拟的解耦方法,
地球物理学报, 2005,48(4):909-917
Magsci [本文引用: 3]
本文将求解近场波动问题的一种解耦技术推广到两相介质，得到了一种求解两相介质近场波动问题的直接解耦方法，包括集中质量有限元模型、时域显式积分格式和局部人工边界条件. 首先应用加权残数法，并依据波动模拟的精度要求，得到了两相介质集中质量有限元模型. 然后，结合两相介质中波动的衰减特性，实现了透射边界在两相介质近场波动中的运用. 最后，通过数值实验，并与解析解对比，验证了本文方法的有效性.
( Chen Shaolin, Liao Zhenpeng, Chen Jin . Decoupling method for near-field wave simulation of two-phase media
Journal of Geophysics, 2005,48(4):909-917 (in Chinese))
Magsci [本文引用: 3]
本文将求解近场波动问题的一种解耦技术推广到两相介质，得到了一种求解两相介质近场波动问题的直接解耦方法，包括集中质量有限元模型、时域显式积分格式和局部人工边界条件. 首先应用加权残数法，并依据波动模拟的精度要求，得到了两相介质集中质量有限元模型. 然后，结合两相介质中波动的衰减特性，实现了透射边界在两相介质近场波动中的运用. 最后，通过数值实验，并与解析解对比，验证了本文方法的有效性.

[29] Deresiewicz H, Rice JT . The effect of boundaries on wave propagation in a liquid-filled porous solid: V. Transmission across a plane interface
Bull Seis Soc Am, 1964,54(1):409-416


[30] Deresiewicz H . The effect of boundaries on wave propagation in a liquid-filled porous solid: VII. Surface waves in a half-space in the presence of a liquid layer
Bull Seis Soc Am, 1964,54(1):425-430
[本文引用: 2]

[31] Thomson WT . Transmission of elastic waves through a stratified solid media
Journal of Applied Physics, 1950,21:89-93
DOIURL
The transmission of a plane elastic wave at oblique incidence through a stratified solid medium consisting of any number of parallel plates of different material and thickness is studied theoretically. The matrix method is used to systematize the analysis and to present the equations in a form suitable for computation.

[32] Haskell NA . The dispersion of surface waves on multilayered media
Bull Seismol Soc Am, 1953,43:17-34
DOIURL
A matrix formalism developed by W. T. Thomson is used to obtain the phase velocity dispersion equations for elastic surface waves of Rayleigh and Love type on multi-layered solid media. The method is used to compute phase and group velocities of Rayleigh waves for two assumed three-layer and one two-layer model of the earth's crust in the continents. The computed group velocity curves are compared with published values of the group velocities at various frequencies of Rayleigh waves over continental paths. The scatter of the observed values is larger than the difference between the three computed curves. It is believed that not all of this scatter is due to observational errors but probably represents a real horizontal heterogeneity of the continental crusts. (Author)

[33] 柯小飞, 陈少林, 张洪翔 .P-SV波入射时海水-层状海床体系的自由场分析
. 振动工程学报, 2018, 录用
[本文引用: 1]

( Ke Xiaofei, Chen Shaolin, Zhang Hongxiang . Freefield analysis of seawater-layered seabed system at P-SV wave incident
Journal of Vibration Engineering, 2018, Accepted (in Chinese))
[本文引用: 1]