地震激励下等横截面无限水库频域和时域响应的FEM-SBFEM计算方法

李上明; 肖世富

doi:10.21656/1000-0887.450138

地震激励下等横截面无限水库频域和时域响应的FEM-SBFEM计算方法

doi: 10.21656/1000-0887.450138

李上明^{1, 2, ,},
肖世富¹

1.
中国工程物理研究院总体工程研究所，四川绵阳 621999
2.
工程材料与结构冲击振动四川省重点实验室，四川绵阳 621999

我刊编委肖世富来稿

基金项目:

国家自然科学基金 12472171

国家自然科学基金 U2430201

详细信息

通讯作者:
李上明(1978—)，男，研究员，博士，博士生导师(通讯作者. E-mail: hustmingsl@126.com)

中图分类号: O39
计量
- 文章访问数: 116
- HTML全文浏览量: 35
- PDF下载量: 38
- 被引次数: 0
出版历程
- 收稿日期: 2024-05-13
- 修回日期: 2025-02-24
- 刊出日期: 2025-04-01

An FEM-SBFEM Coupled Method for Infinite Reservoir Responses With Uniform Cross Sections Under Seismic Excitations in Frequency and Time Domains

LI Shangming^{1, 2
, ,},
XIAO Shifu¹

1.
Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang, Sichuan 621999, P.R.China
2.
Shock and Vibration of Engineering Materials and Structures Key Laboratory of Sichuan Province, Mianyang, Sichuan 621999, P.R.China

Contributed by XIAO Shifu, M.AMM Editorial Board

摘要

摘要: 坝体抗震设计和评估需要准确计算无限水库动力响应. 基于比例边界有限元法(scaled boundary finite element method，SBFEM)力学推导技术，推导了顺河向地震激励下等横截面无限水域频域响应计算公式，利用Fourier逆变换建立了时域响应控制方程，通过线性叠加推导了顺河、横河、竖直三向组合地震激励下的无限水域频域和时域响应的SBFEM计算公式. 结合有限元法，建立了无限水库频域和时域响应的FEM-SBFEM耦合方程. 分析了地震激励下的二维、三维等横截面无限水库频域、时域响应，数值验证了所建立计算公式的正确性. 所发展的FEM-SBFEM公式体系可推广应用于库底库岸具有吸收性的、横截面有任意几何形状的无限水库谐响应及瞬态响应分析.
- 无限水库 /
- 瞬态分析 /
- 谐响应分析 /
- 比例边界有限元 /
- 有限元
Abstract: It is necessary to accurately calculate the dynamic responses of infinite reservoirs in seismic design and evaluation of dams. Based on the mechanical derivation technique of the scaled boundary finite element method (SBFEM), the frequency domain response formulas for the infinite reservoirs with uniform cross sections under upstream seismic excitations were derived, and the time domain response governing equations were obtained through the inverse Fourier transform. Then the SBFEM formulas of frequency domain and time domain responses of infinite reservoirs under upstream, cross-stream and vertical seismic excitations were derived by linear superposition. The FEM-SBFEM coupled equations for frequency and time domain responses of infinite reservoirs were established. Based on these proposed formulations, the time domain and frequency domain responses of the 2D or 3D infinite reservoirs with uniform cross sections under seismic excitations were analyzed. The numerical results validate the accuracy of the proposed formulations. The proposed formulations can be used to analyze harmonic and transient responses of infinite reservoirs with arbitrary geometry and absorptive reservoir bottom under seismic excitations.
- infinite reservoir /
- transient analysis /
- harmonic response analysis /
- SBFEM /
- FEM
edited-by

edited-by
注释:

1) 我刊编委肖世富来稿

HTML全文

0. 引言

中国是世界坝工大国，大坝安全在国民经济发展中具有特别重要的意义. 坝体抗震设计和评估时常采用数值方法. 在数值方法中，通常将无限水库截断为两部分：任意几何形状的近场和具有等横截面的半无限域(即远场). 其中，坝体和近场为有限域，可用有限元法模拟，而远场为半无限域，需采用能正确描述无穷远处无反射条件的模拟方法.

传统的Sommerfeld无反射边界条件^[1-5]常用来模拟半水库. Sommerfeld无反射边界条件假定在近场与远场相互作用面上(近远场耦合面)的水动压力、其梯度等物理量满足某一方程，并将该方程认为是耦合面上的边界条件. 通常，该边界条件只能近似描述半无限远场的无反射特性，并且某些传递边界条件只有在其耦合面离坝库耦合面较远时才能获取合理的结果. 为了避免上述缺陷，研究人员逐渐采用比例边界有限元法(SBFEM)^[6]模拟等横截面的无限远场.

SBFEM是一种新型半解析力学数值模拟方法，径向上无需划分网格，为解析解；环向上随着网格数的增多收敛于有限元意义上的精确解. 它既有边界元法只需离散边界、空间维数下降一维的特点，又保持了有限元法的矩阵对称、无需基本解和不出现奇异积分等优点，它能有效地模拟无限水库^[7-17]，精确满足无穷远处的无反射条件. SBFEM计算公式可分别采用比例边界变换推导(scaled boundary transformation-based derivation)^[7]技术和力学推导(mechanically-based derivation)技术^[8]推导，两者是一致的.

在谐响应频域分析中，针对横截面垂直顺河向的等横截面无限水库，文献[8]基于力学推导技术，建立了基于动态刚度矩阵的SBFEM方程，借助Schur分析和模态分析技术，发展了等横截面无限水库的动态刚度矩阵求解技术，并通过与有限元法的结合分析了任意形状水库在顺河向激励下的频域响应. 其中有限元法用来模拟任意形状的有限水库(近场)，SBFEM用来模拟横截面垂直与顺河向的等横截面无限水库(远场). 数值算例显示其结果好于其他方法，但其公式体系只适用于顺河向激励响应分析，无法模拟横河向和竖直向激励响应分析. 文献[9]则在文献[8]的基础上，通过与相应有限元法公式体系对比分析，改进了该项不足，建立了适用于远场在任意激励方向下响应分析的SBFEM公式体系，但其还需进一步归一化处理，以保证其能适用横截面具有任意几何形状的三维问题. 而基于比例边界变换技术的SBFEM公式^[7]则可准确模拟等横截面无限水库在顺河向、横河向和竖直向等任意激励方向下的频域响应. 相比于基于力学推导技术的公式体系，该技术可适用于横截面垂直和不垂直于顺河向的等横截面无限水库，但需采用Hamilton矩阵的特征值问题求解技术进行谐响应分析，求解过程和公式复杂，并且分析自由度增加了一倍. 文献[10]则利用该公式体系分析了垂直面重力坝的频域响应. 上述两种公式在形式上和求解技术上存在差异，但均能有效模拟等横截面无限水库的谐响应. 文献[7-9]互为补充，建立了等横截面无限水库的SBFEM频域分析理论.

对等横截面无限水库的瞬态响应，文献[11]、文献[12]、文献[13-15]分别基于动态质量、动态刚度矩阵和连分式方法，建立了有限元与SBFEM瞬态分析耦合公式，准确获取了顺河向激励下的等横截面无限水库瞬态响应. 而文献[16]则通过SBFEM建立了不可压无限水域附加质量，实现了坝库耦合分析，但不可压假定往往会低估坝体响应，应用受到限制. 文献[11-16]中SBFEM瞬态公式体系共同的特点是均只适合等横截面无限水库库底及库岸无吸收性与外界激励为顺河向激励的情况.

虽然基于比例边界变换的SBFEM^[7]在Hamilton矩阵的特征值及特征向量的求解基础上，通过Fourier逆变换，推导了瞬态分析所需的三类单位脉冲响应矩阵，建立了含有三个时域卷积积分的SBFEM瞬态分析公式，适用于等横截面无限水库库底及库岸有吸收性和任意向激励的情况，但单位脉冲响应矩阵缺少明确的数学表达式，其值依赖于外界激励的大小和频率分布，计算精度不易控制，求解存在一定困难. 文献[7]对上述公式的验证分析中，只给出了单频率周期激励下的竖直向瞬态响应和含多频率成份的顺河向瞬态响应，并未说明其无限水库库底及库岸吸收性的情况，其普适性还需进一步验证. 此外，其分析效率受到三个时域卷积积分^[7]的严重影响.

综上，本文将基于力学推导技术，建立适用任意几何形状的等横截面无限水库在顺河向、横河向和竖直向等激励方向下的频域公式，并发展对应的时域瞬态分析公式，进而建立无限水库频域和时域响应的FEM-SBFEM计算公式，并采用算例验证其正确性.

1. 无限水库控制方程与边界条件

对于图 1所示的坝库耦合系统，在数值分析中无限水库常考虑为无旋无黏的、不考虑体力的小扰动流体，控制方程为如下的波动方程：

$\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}+\frac{\partial^{2}}{\partial z^{2}}\right) p=\nabla^{2} p=\frac{1}{c^{2}} \ddot{p}=\frac{1}{c^{2}} \frac{\partial^{2} p}{\partial t^{2}},$

(1)

图 1 坝库系统

Figure 1. The dam-reservoir system

下载: 全尺寸图片幻灯片

式中p和c为水动压力(后续简称为压力)和水库中的声速，x，y和z为直角坐标， $\nabla^2$ 为Laplace算子. 如考虑体力，其可视为静水压力，不改变后续公式推导. 总水压力为水动压力与静水压力标量之和.

式(1)的边界条件如下：

1) 在坝库耦合面处，满足

$\frac{\partial p}{\partial n}=-\rho a_{n},$

(2)

式中a_n为边界法向加速度，指向水库外部，下同；ρ为水库中水的密度.

2) 在库底和库岸处，满足

$\frac{\partial p}{\partial n}+q \dot{p}=-\rho a_{n},$

(3)

其中q定义为

$q=\frac{1}{c}\left(\frac{1-\alpha}{1+\alpha}\right),$

(4)

式中α为库底和库岸的反射系数.

3) 忽略重力波时，在自由表面处，满足

$p=0 .$

(5)

为获得无限水库的响应，通常无限水库截断成近场和远场两部(见图 1). 近场具有任意几何构型，远场为等横截面无限水库. 远场对近场的影响通过近远场耦合面的边界条件描述. 该边界条件与式(1)—(5)构成了初边值问题，求解该初边值问题便可实现无限水库的响应求解. 本文将基于SBFEM建立近远场耦合面上的边界条件，即远场控制方程.

2. 顺河向激励下远场SBFEM控制方程

对于图 2的等横截面远场，其承受的地震激励方向可能有顺河向、横河向、竖直向等. 图 2中x向为顺河向，y向为横河向，z向为竖直向； $a_n^{\mathrm{CS}}(x, y, z)$ 表示顺河向激励，而 $a_n^{\mathrm{BS}}(x, y, z)$ 表示横河向、竖直向激励，因其载荷形式一样，统称横河向激励. 假定沿x向任意x处横河向激励一致，即 $a_n^{\mathrm{BS}}(x, y, z)=a_n^{\mathrm{BS}}(y, z)$ . 当采用SBFEM模拟远场时，整个无限远场只需离散近远场耦合面. 图 3给出了用于模拟三维等横截面无限远场域的八结点SBFEM网格.

图 2 具有等横截面的远场

Figure 2. The far field with a uniform cross section

下载: 全尺寸图片幻灯片

图 3 远场SBFEM网格即近远场耦合面网格

Figure 3. The SBFEM discretization of the far field

下载: 全尺寸图片幻灯片

2.1 频域控制方程

经SBFEM离散近远场耦合面后，根据SBFEM理论^{[6, 8]}，在顺河向激励下，近远场耦合面 $\mathit{\Gamma}_{\mathrm{w}}^{\mathrm{e}}$ 上压力与加速度关系满足

$\boldsymbol{A}_{\mathrm{n}}(\omega)=\boldsymbol{S}^{\infty}(\omega) \boldsymbol{P}(\omega),$

(6)

式中 $\boldsymbol{S}^{\infty}(\omega)$ 为顺河向激励下无限远场动态刚度矩阵； $\boldsymbol{P}$ 为离散面节点压力列向量； $\boldsymbol{A}_{\mathrm{n}}$ 为离散面外法向节点加速度列向量； $\omega$ 为圆频率. $\boldsymbol{A}_{\mathrm{n}}(\omega)$ 满足

$\boldsymbol{A}_{\mathrm{n}}(\omega)=-\rho \Sigma_{e} \int_{\mathit{\Gamma}_{\mathrm{w}}^{\mathrm{e}}} \boldsymbol{N}_{\mathrm{f}}^{\mathrm{T}} a_{n}^{\mathrm{CS}}(\omega) \mathrm{d} \mathit{\Gamma}_{\mathrm{w}}^{\mathrm{e}},$

(7)

式中 $\boldsymbol{N}_{\mathrm{f}}$ 为近远场耦合面离散采用的有限元单元的形函数； $\Sigma_{e}$ 为单元装配； $\mathit{\Gamma}_{w}^{e}$ 为近远场耦合面；上标T表示矩阵转置. 采用文献[, ]中的力学推导技术，可推导 $\boldsymbol{S}^{\infty}({\omega})$ 满足方程

$\left(\boldsymbol{S}^{\infty}(\omega)+\boldsymbol{E}^{1}\right)\left(\boldsymbol{E}^{0}\right)^{-1}\left(\boldsymbol{S}^{\infty}(\omega)+\left(\boldsymbol{E}^{1}\right)^{\mathrm{T}}\right)-\boldsymbol{E}^{2}-\mathrm{j} \omega \boldsymbol{C}^{0}-(\mathrm{j} \omega)^{2} \boldsymbol{M}^{0}=\mathbf{0},$

(8)

式中 $\boldsymbol{E}^{0}, \boldsymbol{E}^{1}, \boldsymbol{E}^{2}, \boldsymbol{C}^{0}, \boldsymbol{M}^{0}$ 为SBFEM系数矩阵；j为纯虚数单位. 式(8)可采用Hamilton矩阵特征值和特征向量进行求解^[7].

2.2 SBFEM系数矩阵

SBFEM系数矩阵 $\boldsymbol{E}^{0}, \boldsymbol{E}^{1}, \boldsymbol{E}^{2}, \boldsymbol{C}^{0}, \boldsymbol{M}^{0}$ 只依赖于离散面，由相应单元系数矩阵装配而成，其计算方法和装配过程均同有限元单元刚度的计算方法和装配过程. 它们的表达式如下(详细推导方法可参考文献[8]):

$\boldsymbol{E}^{0}=\int_{-1}^{1} \int_{-1}^{1} \boldsymbol{B}^{\mathrm{1T}} \boldsymbol{B}^{1}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{~d} \zeta, \boldsymbol{E}^{1}=\int_{-1}^{1} \int_{-1}^{1} \boldsymbol{B}^{2 \mathrm{~T}} \boldsymbol{B}^{1}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{~d} \zeta, \boldsymbol{E}^{2}=\int_{-1}^{1} \int_{-1}^{1} \boldsymbol{B}^{2 \mathrm{~T}} \boldsymbol{B}^{2}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{~d} \zeta ,$

(9a)

$\boldsymbol{M}^{0}=\int_{-1}^{1} \int_{-1}^{1} \frac{1}{c^{2}} \boldsymbol{N}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{N}_{\mathrm{f}}|\boldsymbol{J}| \mathrm{d} \eta \mathrm{~d} \zeta, \boldsymbol{C}^{0}=q \int_{\mathit{\Gamma}_{\mathrm{b}}} \boldsymbol{N}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{N}_{\mathrm{f}} \mathrm{~d} \mathit{\Gamma}_{\mathrm{b}},$

(9b)

式中 $\eta, \zeta$ 为等参单元坐标； $\mathit{\Gamma}_{\mathrm{b}}$ 为横截面内的库底和库岸边界； $\boldsymbol{J}, \boldsymbol{B}^{1}, \boldsymbol{B}^{2}$ 定义为

$\boldsymbol{J}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{x} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{z} \\ \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{x} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{z} \end{array}\right],$

(10)

式中 $\boldsymbol{x}, \boldsymbol{y}, \boldsymbol{z}$ 为单元节点坐标列向量. 式(10)中 $\boldsymbol{J}$ 的逆 $\boldsymbol{J}^{-1}$ 定义为

$\boldsymbol{J}^{-1}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ j_{21} & j_{22} & j_{23} \\ j_{31} & j_{32} & j_{33} \end{array}\right].$

(11)

在此基础上， $\boldsymbol{B}^{1}$ 满足

$\boldsymbol{B}^{1}=\left[\begin{array}{c} 1 \\ j_{21} \\ j_{31} \end{array}\right] \boldsymbol{N}_{\mathrm{f}},$

(12)

$\boldsymbol{B}^{2}$ 满足

$\boldsymbol{B}^{2}=\left[\begin{array}{c} 0 \\ j_{22} \\ j_{32} \end{array}\right] \frac{\mathrm{d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta}+\left[\begin{array}{c} 0 \\ j_{23} \\ j_{33} \end{array}\right] \frac{\mathrm{d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta}.$

(13)

式(9)—(13)适用于横截面具有任意形状的远场，并且远场横截面可垂直于顺河向，也可不垂直于顺河向.

2.3 横截面垂直于顺河向的SBFEM系数矩阵

对于垂直x轴(顺河向)的近远场耦合面，式(10)变为

$\boldsymbol{J}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{z} \\ 0 & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{z} \end{array}\right],$

(14)

其模 $|\boldsymbol{J}|$ 为

$|\boldsymbol{J}|=\left|\begin{array}{ccc} 1 & 0 & 0 \\ 0 & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{z} \\ 0 & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{z} \end{array}\right|=\left|\begin{array}{ll} \frac{\mathrm{d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \eta} \boldsymbol{z} \\ \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{y} & \frac{\mathrm{~d} \boldsymbol{N}_{\mathrm{f}}}{\mathrm{~d} \zeta} \boldsymbol{z} \end{array}\right|,$

(15)

$\boldsymbol{J}^{-1}$ 为

$\boldsymbol{J}^{-1}=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & j_{22} & j_{23} \\ 0 & j_{32} & j_{33} \end{array}\right],$

(16)

$\boldsymbol{B}^{1}$ 为

$\boldsymbol{B}^{1}=\left[\begin{array}{l} 1 \\ 0 \\ 0 \end{array}\right] \boldsymbol{N}_{\mathrm{f}} .$

(17)

将式(14)—(17)代入式(9a)得

$\boldsymbol{E}^1=\mathbf{0} .$

(18)

将式(18)代入式(8)，变形后可得横截面垂直于x轴的远场动态刚度矩阵 $\boldsymbol{S}^{\infty}(\omega)$ 为

$\boldsymbol{S}^{\infty}(\omega)=\sqrt{\left(\boldsymbol{E}^{2}+\mathrm{j} \omega \boldsymbol{C}^{0}-\omega^{2} \boldsymbol{M}^{0}\right)\left(\boldsymbol{E}^{0}\right)^{-1}} \boldsymbol{E}^{0} .$

(19)

式(19)可利用Schur分解进行矩阵平方根求解.

当无限水库远场域库底和库岸无吸收性，即α=1时，

$\boldsymbol{C}^{0}=\mathbf{0} .$

(20)

将式(20)代入式(19)得

$\boldsymbol{S}^{\infty}(\omega)=\sqrt{\left(\boldsymbol{E}^{2}-\omega^{2} \boldsymbol{M}^{0}\right)\left(\boldsymbol{E}^{0}\right)^{-1}} \boldsymbol{E}^{0}.$

(21)

利用文献[17]的对角化技术，式(21)变为

$\boldsymbol{S}^{\infty}(\omega)=\left(\boldsymbol{X}^{\mathrm{T}}\right)^{-1} \sqrt{\boldsymbol{\mathit{\Lambda}}-\left(\frac{\omega}{c}\right)^{2}} \boldsymbol{I X}^{-1},$

(22)

其中

$\boldsymbol{E}^{2} \boldsymbol{X}=\boldsymbol{E}^{0} \boldsymbol{X} \boldsymbol{\mathit{\Lambda}}, \boldsymbol{X}^{\mathrm{T}} \boldsymbol{E}^{2} \boldsymbol{X}=\boldsymbol{\mathit{\Lambda}}, \boldsymbol{X}^{\mathrm{T}} \boldsymbol{E}^{0} \boldsymbol{X}=\boldsymbol{I},$

(23)

式中 $\boldsymbol{\mathit{\Lambda}}, \boldsymbol{X}$ 分别为 $\boldsymbol{E}^{2}, \boldsymbol{E}^{0}$ 构成的特征值问题的特征值和特征向量矩阵； $\boldsymbol{I}$ 为单位矩阵.

式(22)主要特点为对所有频率ω只需进行一次特征值问题的求解，避免了式(19)和式(8)因频率改变而需重新进行矩阵特征值的求解，极大地提高了分析效率. 当 $\boldsymbol{S}^{\infty}(\omega)$ 已知后，利用式(6)、(7)可获取无限水库远场谐响应.

2.4 时域控制方程

用 $\boldsymbol{S}^{\infty}(\omega)$ 方程式(8)除以 $(\mathrm{j} \omega)^{4}$ 变形后得

$\begin{gather*} \boldsymbol{M}^{\infty}(\omega)\left(\boldsymbol{E}^{0}\right)^{-1} \boldsymbol{M}^{\infty}(\omega)+\boldsymbol{E}^{1}\left(\boldsymbol{E}^{0}\right)^{-1} \frac{\boldsymbol{M}^{\infty}(\omega)}{(\mathrm{j} \omega)^{2}}+\frac{\boldsymbol{M}^{\infty}(\omega)}{(\mathrm{j} \omega)^{2}}\left(\boldsymbol{E}^{0}\right)^{-1}\left(\boldsymbol{E}^{1}\right)^{\mathrm{T}}- \\ \frac{1}{(\mathrm{j} \omega)^{4}}\left(\boldsymbol{E}^{2}-\boldsymbol{E}^{1}\left(\boldsymbol{E}^{0}\right)^{-1}\left(\boldsymbol{E}^{1}\right)^{\mathrm{T}}\right)-\frac{1}{(\mathrm{j} \omega)^{3}} \boldsymbol{C}^{0}-\frac{1}{(\mathrm{j} \omega)^{2}} \boldsymbol{M}^{0}=0 , \end{gather*}$

(24)

式中 $\boldsymbol{M}^{\infty}(\omega)$ 为顺河向激励下无限远场动态质量矩阵，等于 $\boldsymbol{S}^{\infty}(\omega) /(\mathrm{j} \omega)^{2}$ .对式(24)进行Fourier逆变换得

$\begin{aligned} & \int_{0}^{t} \boldsymbol{M}^{\infty}(t-\boldsymbol{\tau})\left(\boldsymbol{E}^{0}\right)^{-1} \boldsymbol{M}^{\infty}(\boldsymbol{\tau}) \mathrm{d} \boldsymbol{\tau}+\boldsymbol{E}^{1}\left(\boldsymbol{E}^{0}\right)^{-1} \int_{0}^{t} \int_{0}^{\tau} \boldsymbol{M}^{\infty}\left(\boldsymbol{\tau}^{\prime}\right) \mathrm{d} \boldsymbol{\tau}^{\prime} \mathrm{d} \boldsymbol{\tau}+ \\ & \quad \int_{0}^{t} \int_{0}^{\tau} \boldsymbol{M}^{\infty}\left(\boldsymbol{\tau}^{\prime}\right) \mathrm{d} \tau^{\prime} \mathrm{d} \boldsymbol{\tau}\left(\boldsymbol{E}^{0}\right)^{-1}\left(\boldsymbol{E}^{1}\right)^{\mathrm{T}}-\frac{t^{3}}{6}\left(\boldsymbol{E}^{2}-\boldsymbol{E}^{1}\left(\boldsymbol{E}^{0}\right)^{-1}\left(\boldsymbol{E}^{1}\right)^{\mathrm{T}}\right)-\frac{t^{2}}{2} \boldsymbol{C}^{0}-t \boldsymbol{M}^{0}=\mathbf{0} , \end{aligned}$

(25)

式中 $t$ 为时间； $\tau, \tau^{\prime}$ 为积分变量； $\boldsymbol{M}^{\infty}(t)$ 与 $\boldsymbol{M}^{\infty}(\omega)$ 构成Fourier变换对.简化式（25）为

$\begin{aligned} & \int_{0}^{t} \boldsymbol{m}^{\infty}(t-\tau) \boldsymbol{m}^{\infty}(\tau) \mathrm{d} \tau+\boldsymbol{e}^{1} \int_{0}^{t} \int_{0}^{\tau} \boldsymbol{m}^{\infty}\left(\tau^{\prime}\right) \mathrm{d} \tau^{\prime} \mathrm{d} \tau+ \\ & \quad \int_{0}^{t} \int_{0}^{\tau} \boldsymbol{m}^{\infty}\left(\tau^{\prime}\right) \mathrm{d} \tau^{\prime} \mathrm{d} \tau\left(\boldsymbol{e}^{1}\right)^{\mathrm{T}}-\frac{t^{3}}{6} \boldsymbol{e}^{2}-\frac{t^{2}}{2} \boldsymbol{c}^{0}-t \boldsymbol{m}^{0}=\mathbf{0} , \end{aligned}$

(26)

式中 $\boldsymbol{m}^{\infty}(t)$ 为顺河向激励下远场等效动态质量矩阵； $\boldsymbol{e}^{1}, \boldsymbol{e}^{2}, \boldsymbol{c}^{0}, \boldsymbol{m}^{0}$ 为系数矩阵，满足

$\left\{\begin{array}{l} \boldsymbol{m}^{\infty}(t)=\left(\boldsymbol{U}^{-1}\right)^{\mathrm{T}} \boldsymbol{M}^{\infty}(t) \boldsymbol{U}^{-1}, \boldsymbol{e}^{1}=\left(\boldsymbol{U}^{-1}\right)^{\mathrm{T}} \boldsymbol{E}^{1} \boldsymbol{U}^{-1}, \boldsymbol{e}^{2}=\left(\boldsymbol{U}^{-1}\right)^{\mathrm{T}}\left(\boldsymbol{E}^{2}-\boldsymbol{E}^{1}\left(\boldsymbol{E}^{0}\right)^{-1}\left(\boldsymbol{E}^{1}\right)^{\mathrm{T}}\right) \boldsymbol{U}^{-1} , \\ \boldsymbol{c}^{0}=\left(\boldsymbol{U}^{-1}\right)^{\mathrm{T}} \boldsymbol{C}^{0} \boldsymbol{U}^{-1}, \boldsymbol{m}^{0}=\left(\boldsymbol{U}^{-1}\right)^{\mathrm{T}} \boldsymbol{M}^{0} \boldsymbol{U}^{-1}, \boldsymbol{E}^{0}=\boldsymbol{U}^{\mathrm{T}} \boldsymbol{U}. \end{array}\right.$

(27)

对于横截面垂直x轴(顺河向)的远场，有 $\boldsymbol{E}^{1}=\mathbf{0}$ ，则

$\boldsymbol{e}^{1}=\mathbf{0}, \boldsymbol{e}^{2}=\left(\boldsymbol{U}^{-1}\right)^{\mathrm{T}} \boldsymbol{E}^{2} \boldsymbol{U}^{-1} .$

(28)

利用式(28)，式(26)简化后为

$\int_{0}^{t} \boldsymbol{m}^{\infty}(t-\tau) \boldsymbol{m}^{\infty}(\tau) \mathrm{d} \tau-\frac{t^{3}}{6} \boldsymbol{e}^{2}-\frac{t^{2}}{2} \boldsymbol{c}^{0}-t \boldsymbol{m}^{0}=\mathbf{0}.$

(29)

式(29)可采用下述方法进行近似求解.

假定 $\boldsymbol{m}^{\infty}(t)$ 在每个时间步 $\Delta t$ 内为分段常数，并用符号 $\boldsymbol{m}_{n}^{\infty}$ 表示第 $n$ 个时间步 $\boldsymbol{m}^{\infty}(t=n \Delta t)$ 的值.式（29）时间离散后，在 $t=\Delta t$ 时，由式（29）得

$\boldsymbol{m}_{1}^{\infty}=\sqrt{\frac{\Delta t^{2}}{6} \boldsymbol{e}^{2}+\frac{\Delta t}{2} \boldsymbol{c}^{0}+\boldsymbol{m}^{0}} ,$

(30)

在 $t=n \Delta t$ 时，由式(29)得

$\boldsymbol{m}_{1}^{\infty} \boldsymbol{m}_{n}^{\infty}+\boldsymbol{m}_{n}^{\infty} \boldsymbol{m}_{1}^{\infty}=-\sum\limits_{j=2}^{n-1} \boldsymbol{m}_{n-j+1}^{\infty} \boldsymbol{m}_{j}^{\infty}+\frac{n^{3} \Delta t^{2}}{6} \boldsymbol{e}^{2}+\frac{n^{2} \Delta t}{2} \boldsymbol{c}^{0}+n \boldsymbol{m}^{0} .$

(31)

式(31)满足Lyapunov方程形式，可利用相应算法进行求解，即利用式(26)—(31)可计算 $\boldsymbol{M}^{\infty}(t)$ 的值. 根据式(26)—(31)可知：顺河向激励下等横截面无限远场动态质量矩阵 $\boldsymbol{M}^{\infty}(t)$ 与其响应无关，只与无限远场横截面形状、库底和库岸吸收性、无限远场域声速等有关.

根据SBFEM理论^{[6, 8]}，等横截面无限远场在顺河向激励下的瞬态响应即与式(6)对应的时域控制方程满足

$\boldsymbol{a}_{\mathrm{n}}(t)=\int_{0}^{t} \boldsymbol{M}^{\infty}(t-\tau) \ddot{\boldsymbol{p}}(\tau) \mathrm{d} \tau,$

(32)

式中 $\boldsymbol{a}_{\mathrm{n}}(t)$ 与 $\boldsymbol{A}_{\mathrm{n}}(\omega)$ 为Fourier变换对； $\boldsymbol{p}(t)$ 与 $\boldsymbol{P}(\omega)$ 为Fourier变换对； $\ddot{\boldsymbol{p}}(t)$ 是 $\boldsymbol{p}(t)$ 对时间的二阶导数.

对式(32)卷积积分进行时间离散后，满足

$\boldsymbol{M}_{1}^{\infty} \dot{\boldsymbol{p}}^{n}=\boldsymbol{a}_{\mathrm{n}}^{n}-\sum\limits_{j=1}^{n-1}\left(\boldsymbol{M}_{n-j+1}^{\infty}-\boldsymbol{M}_{n-j}^{\infty}\right) \dot{\boldsymbol{p}}^{j}+\boldsymbol{M}_{n}^{\infty} \dot{\boldsymbol{p}}^{0} ,$

(33)

式中 $\boldsymbol{a}_{\mathrm{n}}^{n}=\boldsymbol{a}_{\mathrm{n}}(n \Delta t); \dot{\boldsymbol{p}}^{j}=\dot{\boldsymbol{p}}(j \Delta t)$ 为 $\boldsymbol{p}(t)$ 对时间的一阶导数， $j=1, 2, \cdots, n; \boldsymbol{M}_{j}^{\infty}=\boldsymbol{M}^{\infty}(j \Delta t)$ .

式(32)的初始条件为无限水库远场静止，即满足

$\boldsymbol{p}(0)=\mathbf{0}, \dot{\boldsymbol{p}}(0)=\dot{\boldsymbol{p}}^{0}=\mathbf{0}.$

(34)

结合上述初始条件，式(33)变为

(35)

采用式(35)可获取等横截面无限远场在顺河向激励下的瞬态响应.

3. 横河向、顺河向、竖直向等激励方向下远场控制方程

对于横截面垂直x轴(顺河向)的等横截面无限水库远场，远场整体在横河向与竖直向激励下，在近远场耦合面(离散面)上，频域控制方程满足^[9]

$\boldsymbol{A}_{\mathrm{I}}(\omega)=\left(\boldsymbol{E}^{2}+\mathrm{j} \omega \boldsymbol{C}^{0}-\omega^{2} \boldsymbol{M}^{0}\right) \boldsymbol{P}_{\mathrm{I}}=\boldsymbol{S}_{\mathrm{I}}^{\infty}(\omega) \boldsymbol{P}_{\mathrm{I}},$

(36)

$\boldsymbol{A}_{\mathrm{I}}(\omega)=-\rho \Sigma_{e} \int_{\mathit{\Gamma}_{\mathrm{b}}} \boldsymbol{N}_{\mathrm{f}}^{\mathrm{T}} a_{n}^{\mathrm{BS}}(\omega) \mathrm{d} \mathit{\Gamma}_{\mathrm{b}},$

(37)

式中 $S_{1}^{\infty}(\omega)$ 为横河向及竖直向激励下远场动态刚度矩阵； $\boldsymbol{P}_{\mathrm{I}}(\omega)$ 为横河向及竖直向激励下离散面节点压力列向量； $\boldsymbol{A}_{\mathrm{I}}(\omega)$ 为横河向及竖直向激励下离散面外法向节点加速度列向量； $\mathit{\Gamma}_{\mathrm{b}}$ 为离散面的库底和库岸边界； $a_{n}^{\mathrm{BS}}(\omega)$ 为离散面边界 $\mathit{\Gamma}_{\mathrm{b}}$ 的外法向加速度，即横河向及竖直向激励，见图 2.

式(36)对应的时域控制方程满足

$\boldsymbol{a}_{\mathrm{I}}(t)=\boldsymbol{E}^{2} \boldsymbol{p}_{\mathrm{I}}(t)+\boldsymbol{C}^{0} \dot{\boldsymbol{p}}_{\mathrm{I}}(t)+\boldsymbol{M}^{0} \ddot{\boldsymbol{p}}_{\mathrm{I}}(t),$

(38)

式中 $\boldsymbol{a}_{\mathrm{I}}(t)$ 与 $\boldsymbol{a}_{\mathrm{I}}(\omega)$ 为Fourier变换对； $\boldsymbol{p}_{\mathrm{I}}(t)$ 与 $\boldsymbol{p}_{\mathrm{I}}(\omega)$ 为Fourier变换对； $\ddot{p}_{\mathrm{I}}(t), \dot{\boldsymbol{p}}_{\mathrm{I}}(t)$ 为 $\boldsymbol{p}_{\mathrm{I}}(t)$ 对时间的二阶和一阶导数.

由线性叠加原理知，等横截面无限水库远场在顺河向、横河向及竖直向激励共同作用下的节点压力列向量 $\boldsymbol{P}_{\mathrm{T}}(\omega), \boldsymbol{p}_{\mathrm{T}}(t)$ (它们互为Fourier变换对)，满足

$\boldsymbol{P}_{\mathrm{T}}(\omega)=\boldsymbol{P}(\omega)+\boldsymbol{P}_{\mathrm{I}}(\omega), \boldsymbol{p}_{\mathrm{T}}(t)=\boldsymbol{p}(t)+\boldsymbol{p}_{\mathrm{I}}(t) .$

(39)

将式(6)、(36)、(33)代入式(39)得

$\boldsymbol{P}_{\mathrm{T}}(\omega)=\left(\boldsymbol{S}^{\infty}(\omega)\right)^{-1} \boldsymbol{A}_{\mathrm{n}}(\omega)+\left(\boldsymbol{S}_{\mathrm{I}}^{\infty}(\omega)\right)^{-1} \boldsymbol{A}_{\mathrm{I}}(\omega),$

(40)

$\dot{\boldsymbol{p}}_{\mathrm{T}}^{n}=\left(\boldsymbol{M}_{1}^{\infty}\right)^{-1}\left(\boldsymbol{a}_{\mathrm{n}}^{n}-\sum\limits_{j=1}^{n-1}\left(\boldsymbol{M}_{n-j+1}^{\infty}-\boldsymbol{M}_{n-j}^{\infty}\right)\left(\dot{\boldsymbol{p}}_{\mathrm{T}}^{j}-\dot{\boldsymbol{p}}_{\mathrm{I}}^{j}\right)\right)+\dot{\boldsymbol{p}}_{1}^{\mathrm{n}} ,$

(41)

式中 $\dot{\boldsymbol{p}}_{\mathrm{T}}^{j}=\dot{\boldsymbol{p}}_{\mathrm{T}}(j \Delta t)$ 为 $\boldsymbol{p}_{\mathrm{T}}(t)$ 对时间的一阶导数， $\dot{\boldsymbol{p}}_{\mathrm{I}}^{j}=\dot{\boldsymbol{p}}_{\mathrm{I}}(j \Delta t)$ 为 $\boldsymbol{p}_{\mathrm{I}}(t)$ 对时间的一阶导数， $j=1, 2, \cdots, n$ .

式(38)、(40)和(41)给出了顺河向、横河向及竖直向激励下横截面垂直x轴(顺河向)的等横截面无限水库远场的谐响应和瞬态响应计算公式.

4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式

对于图 1坝库耦合系统中的无限水库，采用FEM离散近场(任意构型)，SBFEM离散远场(等横截面无限水库). 其中，近场有限元方程满足

$\left[\begin{array}{cc} \boldsymbol{m}_{11} & \boldsymbol{m}_{10} \\ \boldsymbol{m}_{01} & \boldsymbol{m}_{00} \end{array}\right]\left\{\begin{array}{l} \ddot{\boldsymbol{p}}_{1} \\ \ddot{\boldsymbol{p}}_{0} \end{array}\right\}+\left[\begin{array}{cc} \boldsymbol{c}_{11} & \boldsymbol{c}_{10} \\ \boldsymbol{c}_{01} & \boldsymbol{c}_{00} \end{array}\right]\left\{\begin{array}{l} \dot{\boldsymbol{p}}_{1} \\ \dot{\boldsymbol{p}}_{0} \end{array}\right\}+\left[\begin{array}{cc} \boldsymbol{k}_{11} & \boldsymbol{k}_{10} \\ \boldsymbol{k}_{01} & \boldsymbol{k}_{00} \end{array}\right]\left\{\begin{array}{l} \boldsymbol{p}_{1} \\ \boldsymbol{p}_{0} \end{array}\right\}=\left\{\begin{array}{l} \boldsymbol{a}_{1} \\ \boldsymbol{a}_{0} \end{array}\right\},$

(42)

式中下标1，0分别表示近远场耦合面及近场其他位置处的物理量； $\boldsymbol{m}, \boldsymbol{c}, \boldsymbol{k}$ 分别表示近场整体质量、阻尼、刚度矩阵，均可采用标准有限元程序计算获取； $\boldsymbol{a}$ 为近场边界外法向节点加速度列向量； $\boldsymbol{p}$ 为近场节点压力列向量； $\ddot{p}, \dot{p}$ 为p对时间的二阶和一阶导数.

在近远场流固耦合面上，法向加速度相反而压力相同，即

$\boldsymbol{a}_{1}=-\boldsymbol{a}_{\mathrm{n}}, \boldsymbol{p}_{1}=\boldsymbol{p}_{\mathrm{T}} .$

(43)

在 $t=n \Delta t$ 时刻和零初始条件下，联合式(41)—(43)得

$\boldsymbol{a}_{1}=-\boldsymbol{M}_{1}^{\infty}\left(\dot{\boldsymbol{p}}_{1}^{n}-\dot{\boldsymbol{p}}_{1}^{n}\right)-\sum\limits_{j=1}^{n-1}\left(\boldsymbol{M}_{n-j+1}^{\infty}-\boldsymbol{M}_{n-j}^{\infty}\right)\left(\dot{\boldsymbol{p}}_{1}^{j}-\dot{\boldsymbol{p}}_{1}^{j}\right),$

(44)

$\begin{align*} & {\left[\begin{array}{ll} \boldsymbol{m}_{11} & \boldsymbol{m}_{10} \\ \boldsymbol{m}_{01} & \boldsymbol{m}_{00} \end{array}\right]\left\{\begin{array}{l} \ddot{\boldsymbol{p}}_{1} \\ \ddot{\boldsymbol{p}}_{0} \end{array}\right\}+\left[\begin{array}{cc} \boldsymbol{c}_{11}+\boldsymbol{M}_{1}^{\infty} & \boldsymbol{c}_{10} \\ \boldsymbol{c}_{01} & \boldsymbol{c}_{00} \end{array}\right]\left\{\begin{array}{l} \dot{\boldsymbol{p}}_{1} \\ \dot{\boldsymbol{p}}_{0} \end{array}\right\}+\left[\begin{array}{ll} \boldsymbol{k}_{11} & \boldsymbol{k}_{10} \\ \boldsymbol{k}_{01} & \boldsymbol{k}_{00} \end{array}\right]\left\{\begin{array}{l} \boldsymbol{p}_{1} \\ \boldsymbol{p}_{0} \end{array}\right\}=} \\ & \;\;\;\left\{\begin{array}{c} \boldsymbol{M}_{1}^{\infty} \dot{\boldsymbol{p}}_{1}^{n}-\sum\limits_{j=1}^{n-1}\left(\boldsymbol{M}_{n-j+1}^{\infty}-\boldsymbol{M}_{n-j}^{\infty}\right)\left(\dot{\boldsymbol{p}}_{1}^{j}-\dot{\boldsymbol{p}}_{1}^{j}\right) \\ \boldsymbol{a}_{0} \end{array}\right\} \cdot \end{align*}$

(45)

式(42)对应的频域公式为

$(\mathrm{j} \omega)^{2}\left[\begin{array}{ll} \boldsymbol{m}_{11} & \boldsymbol{m}_{10} \\ \boldsymbol{m}_{01} & \boldsymbol{m}_{00} \end{array}\right]\left\{\begin{array}{l} \boldsymbol{P}_{1} \\ \boldsymbol{P}_{0} \end{array}\right\}+\mathrm{j} \omega\left[\begin{array}{ll} \boldsymbol{c}_{11} & \boldsymbol{c}_{10} \\ \boldsymbol{c}_{01} & \boldsymbol{c}_{00} \end{array}\right]\left\{\begin{array}{l} \boldsymbol{P}_{1} \\ \boldsymbol{P}_{0} \end{array}\right\}+\left[\begin{array}{ll} \boldsymbol{k}_{11} & \boldsymbol{k}_{10} \\ \boldsymbol{k}_{01} & \boldsymbol{k}_{00} \end{array}\right]\left\{\begin{array}{l} \boldsymbol{P}_{1} \\ \boldsymbol{P}_{0} \end{array}\right\}=\left\{\begin{array}{l} \boldsymbol{A}_{1} \\ \boldsymbol{A}_{0} \end{array}\right\} .$

(46)

式(46)与式(42)中的 $\boldsymbol{P}$ 与 $\boldsymbol{p}, \boldsymbol{A}$ 与 $\boldsymbol{a}$ 互为Fourier变换对. 联合式(40)、(43)、(46)得

$\begin{align*} & \left((\mathrm{j} \omega)^{2}\left[\begin{array}{ll} \boldsymbol{m}_{11} & \boldsymbol{m}_{10} \\ \boldsymbol{m}_{01} & \boldsymbol{m}_{00} \end{array}\right]+\mathrm{j} \omega\left[\begin{array}{ll} \boldsymbol{c}_{11} & \boldsymbol{c}_{10} \\ \boldsymbol{c}_{01} & \boldsymbol{c}_{00} \end{array}\right]+\left[\begin{array}{cc} \boldsymbol{k}_{11}+\boldsymbol{S}^{\infty}(\omega) & \boldsymbol{k}_{10} \\ \boldsymbol{k}_{01} & \boldsymbol{k}_{00} \end{array}\right]\right)\left\{\begin{array}{l} \boldsymbol{P}_{1} \\ \boldsymbol{P}_{0} \end{array}\right\}= \\ & \;\;\;\;\left\{\begin{array}{c} \boldsymbol{S}^{\infty}(\omega)\left(\boldsymbol{S}_{1}^{\infty}(\omega)\right)^{-1} \boldsymbol{A}_{1}(\omega) \\ \boldsymbol{A}_{0} \end{array}\right\} \cdot \end{align*}$

(47)

式(45)和式(47)为含横截面垂直x轴(顺河向)的等横截面无限水库远场与任意几何形状近场的无限水库的FEM-SBFEM耦合公式. 该公式可模拟远场具有等横截面特性的、具有任意几何构型的无限水库地震激励响应.

5. 算例

5.1 频域分析

考虑图 4等水深矩形截面坝库系统，水库两侧库岸无吸收性，库底具有吸收性(α=0.8)，坝体和水库两侧面为刚性，忽略自由液面波. 由于整个水库具有相同横截面，故整个水库用SBFEM离散. 离散面为坝库耦合面. 网格为10个八节点单元，即y向1个单元，z向10个单元. 坝体受顺河向地震加速度 $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 激励时，坝面z=0.6H处的水动压力p如所示；而整个坝库系统受竖直向地震加速度 $a_{z} \mathrm{e}^{\mathrm{i} \omega t}$ 激励时，相同位置处的水动压力p如图 5(b)所示. 图中ρ为库水密度，H为水深. 图中analytical、matrix square root、Hamilton的结果分别由解析解^[18]、式(40)和Hamilton方法^[7]求解获得. 三种方法的结果几乎完全一致.

图 4 三维矩形截面坝库系统

Figure 4. The 3D dam-reservoir system with a rectangular channel

下载: 全尺寸图片幻灯片

图 5 顺河和竖直向激励下坝体z=0.6H处水动压力

Figure 5. Dam's pressures at z=0.6H induced by upstream and vertical excitations

下载: 全尺寸图片幻灯片

为了验证式(40)模拟非矩形截面的可行性，考虑半圆形等横截面无限水库整体在顺河向 $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 、横河向 $a_{z} \mathrm{e}^{\mathrm{i} \omega t}$ 与竖直向 $a_{y} \mathrm{e}^{\mathrm{i} \omega t}$ 激励下的响应. $r=H$ 边界的反射系数为 $\alpha, y=0$ 边界为自由液面. 整个半圆形无限水库用SBFEM离散，离散网格见. $r=0.6 H, \theta=0$ 点处在顺河向 $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 、竖直向 $a_{y} \mathrm{e}^{\mathrm{i} \omega t}$ 激励下的水动压力p见图 8、图 9. r=0.6H，θ=45°点处在横河向 $a_{z} \mathrm{e}^{\mathrm{i} \omega t}$ 激励下的水动压力p见. 图中 $T_{x}=p /\left(\rho H a_{x}\right), T_{y}=p /\left(\rho H a_{y}\right), T_{z}=p /\left(\rho H a_{z}\right)$ . 图中analytical、matrix square root、Hamilton的结果分别由解析解^[18]、式(40)和Hamilton方法^[7]求解获得. 三种方法的结果差异小，这表明所建立公式能正确模拟任意等横截面无限水库在任意方向激励下的响应.

图 6 半圆形等横截面无限水库

Figure 6. The semi-circular infinite reservoir with a uniform cross section

下载: 全尺寸图片幻灯片

图 7 半圆形横截面离散网格

Figure 7. The semi-circular infinite reservoir cross section mesh

下载: 全尺寸图片幻灯片

图 8 顺河向激励

$a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 下

$r=0.6 H, \theta=0$ 的压力

Figure 8. Pressures of point

$r=0.6 H, \theta=0$ under upstream excitations

$a_{x} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

图 9 坚直向激励

$a_{y} \mathrm{e}^{\mathrm{i} \omega t}$ 下

$r=0.6 H, \theta=0$ 的压力

Figure 9. Pressures of point

$r=0.6 H, \theta=0$ under vertical excitations

$a_{y} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

图 10 横河向激励

$a_{z} \mathrm{e}^{\mathrm{i} \omega t}$ 下

$r=0.6 H, \theta=45^{\circ}$ 的压力

Figure 10. Pressures of point

$r=0.6 H, \theta=45^{\circ}$ under cross－stream excitations

$a_{z} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

考虑图 11二维重力坝，其目的为检验FEM-SBFEM耦合公式(即式(47))在复杂坝库系统中应用的可行性和准确性. 坝体考虑为刚性，忽略自由液面波和库底吸收性. 无限水库分解成近场与远场. 近远场耦合面选择在离坝底6.1 m处，以确保近场能相对容易地用四边形网格离散. 无限水库离散网格如图 11所示. 由FEM-SBFEM耦合公式(即式(47))获取的、顺河向地震激励下的坝面压力分布如图 12所示. 基于模态分析技术(式(22))和基于Schur分解的(即矩阵平方根(式(19))的解一致，除了坝底附近的值与Tsai的解^[19]有偏差外，其他位置处吻合较好. 其中，Ω=ωH/c=1.5处的差异较大，这是因为Tsai的解无法准确地获取Ω接近库水一阶固有频率Ω＝π/２、坝面斜面角度θ较大时坝底附近的解^[19].

图 11 二维重力坝及其网格

Figure 11. The 2D gravity dam and its mesh

下载: 全尺寸图片幻灯片

图 12 顺河向激励

$a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 下重力坝坝面压力分布

Figure 12. Pressure distributions along the gravity dam surface under upstream exciations

$a_{x} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

5.2 时域分析

考虑图 13重力坝体在图 14顺河向地震激励下的坝库耦合瞬态响应. 分析中，忽略水库自由表面波及坝体的材料阻尼，不考虑材料和几何非线性，库底刚性并固定. 文献[11]论述了库底无吸收性时SBFEM时域公式的正确性，本文则进一步论述了库底有吸收性时时域控制公式的正确性，即式(26)—(35)和式(45)的正确性. 由于地震激励下、整个水库库底有吸收性的坝库耦合系统的瞬态响应相应研究较少，故文中利用FEM分析结果来校核式(26)—(35)和式(45)的正确性. 图 15给出了库底无吸收性(α=1.0)和库底有吸收性(α=0.8)的、利用FEM与式(26)—(35)和式(45)分析的坝角压力P，图中H=120 m. FEM与SBFEM的结果几乎相似，即说明了式(26)—(35)和式(45)的正确性. 分析中采用的时间增量为0.001 s，网格见图 13. 在有限元分析中，为了获取较准确的结果，近远场耦合面离坝角1 926 m，单元6 681个，单元边长约6 m；近远场耦合面、库底的吸收性通过定义材料的吸收性来模拟实现. 相比于有限元分析技术，文中技术具有离散区域小、自由度少的特点.

图 13 重力坝库系统及其网格

Figure 13. The gravity dam-reservoir system and its mesh

下载: 全尺寸图片幻灯片

图 14 顺河向斜面加速度

Figure 14. The upsteam ramp acceleration

下载: 全尺寸图片幻灯片

图 15 α=1.0，0.8重力坝坝角压力

Figure 15. Dam heel's pressures of the gravity dam for α=1.0, 0.8

下载: 全尺寸图片幻灯片

6. 结论

本文推导了横河向激励下等横截面无限水库响应的SBFEM频域公式，建立了库底和库边具有吸收性的等横截面无限水库的时域响应公式及任意方向激励下等横截面无限水库谐响应和瞬态响应分析的FEM-SBFEM耦合公式，可为坝体地震激励下开裂、损伤等模拟提供水库水动压力计算方法. 本文的计算公式具有以下特点：

1) 其可以模拟任意几何形状的等横截面无限水库的谐响应和瞬态响应；

2) 其可以模拟等横截面无限水库库底和库边的能量吸收性；

3) 其满足无限水库无穷远处的无反射条件，为半解析解，具有较高的计算精度，能与有限元法实现无缝耦合.

图 1 坝库系统

Figure 1. The dam-reservoir system

下载: 全尺寸图片幻灯片

图 2 具有等横截面的远场

Figure 2. The far field with a uniform cross section

下载: 全尺寸图片幻灯片

图 3 远场SBFEM网格即近远场耦合面网格

Figure 3. The SBFEM discretization of the far field

下载: 全尺寸图片幻灯片

图 4 三维矩形截面坝库系统

Figure 4. The 3D dam-reservoir system with a rectangular channel

下载: 全尺寸图片幻灯片

图 5 顺河和竖直向激励下坝体z=0.6H处水动压力

Figure 5. Dam's pressures at z=0.6H induced by upstream and vertical excitations

下载: 全尺寸图片幻灯片

图 6 半圆形等横截面无限水库

Figure 6. The semi-circular infinite reservoir with a uniform cross section

下载: 全尺寸图片幻灯片

图 7 半圆形横截面离散网格

Figure 7. The semi-circular infinite reservoir cross section mesh

下载: 全尺寸图片幻灯片

顺河向激励 $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 下 $r=0.6 H, \theta=0$ 的压力

Pressures of point $r=0.6 H, \theta=0$ under upstream excitations $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

坚直向激励 $a_{y} \mathrm{e}^{\mathrm{i} \omega t}$ 下 $r=0.6 H, \theta=0$ 的压力

Pressures of point $r=0.6 H, \theta=0$ under vertical excitations $a_{y} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

横河向激励 $a_{z} \mathrm{e}^{\mathrm{i} \omega t}$ 下 $r=0.6 H, \theta=45^{\circ}$ 的压力

Pressures of point $r=0.6 H, \theta=45^{\circ}$ under cross－stream excitations $a_{z} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

图 11 二维重力坝及其网格

Figure 11. The 2D gravity dam and its mesh

下载: 全尺寸图片幻灯片

顺河向激励 $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$ 下重力坝坝面压力分布

Pressure distributions along the gravity dam surface under upstream exciations $a_{x} \mathrm{e}^{\mathrm{i} \omega t}$

下载: 全尺寸图片幻灯片

图 13 重力坝库系统及其网格

Figure 13. The gravity dam-reservoir system and its mesh

下载: 全尺寸图片幻灯片

图 14 顺河向斜面加速度

Figure 14. The upsteam ramp acceleration

下载: 全尺寸图片幻灯片

图 15 α=1.0，0.8重力坝坝角压力

Figure 15. Dam heel's pressures of the gravity dam for α=1.0, 0.8

下载: 全尺寸图片幻灯片

参考文献(19)

[1]	BEHROOZI A M, VAGHEFI M. Application of coupled RBFDQ-FEM as a meshless method for time-domain analysis of dam-reservoir systems subjected to earthquake excitation[J]. Engineering Analysis With Boundary Elements, 2023, 153 : 160-171. doi: 10.1016/j.enganabound.2023.05.017
[2]	YE J, ZHOU H, ZHOU X. Hydrodynamic pressure on lateral side of dam excited by harmonic seismic vibration: a novel formulation[J]. Soil Dynamics and Earthquake Engineering, 2023, 164 : 107626. doi: 10.1016/j.soildyn.2022.107626
[3]	KHIAVI M P, FERDOUSI A, KHIAVI A M. A probabilistic model for evaluation of the dynamic behavior of a concrete gravity dam considering the fluid-structure interaction[J]. Advances in Civil Engineering, 2023, 2023 : 9927608.
[4]	YA S, EISENTRÄGER S, QU Y L, et al. Seismic analysis of post-tensioned concrete gravity dams using scaled boundary finite elements implemented as ABAQUS UEL[J]. Soil Dynamics and Earthquake Engineering, 2023, 164 : 107620. doi: 10.1016/j.soildyn.2022.107620
[5]	李艳朋, 林皋, 胡志强, 等. 基于声固耦合法的拱坝-库水-地基相互作用分析[J]. 哈尔滨工程大学学报, 2022, 43 (4): 451-457. LI Yanpeng, LIN Gao, HU Zhiqiang, et al. Arch dam-reservoir-foundation interaction analysis based on the acoustic-structure coupling method[J]. Journal of Harbin Engineering University, 2022, 43 (4): 451-457. (in Chinese)
[6]	WOLF J P, SONG C. Finite-Element Modeling of Unbounded Media[M]. New York, USA: Wiley, 1996.
[7]	LIN G, WANG Y, HU Z Q. An efficient approach for frequency-domain and time-domain hydrodynamic analysis of dam-reservoir systems[J]. Earthquake Engineering & Structural Dynamics, 2012, 41 (13): 1725-1749.
[8]	LI S M, LIANG H, LI A M. A semi-analytical solution for characteristics of a dam-reservoir system with absorptive reservoir bottom[J]. Journal of Hydrodynamics, 2008, 20 (6): 727-734. doi: 10.1016/S1001-6058(09)60008-1
[9]	LI S M. Scaled boundary finite element method for semi-infinite reservoir with uniform cross section[J]. International Journal of Computational Methods, 2012, 9 (1): 1240006. doi: 10.1142/S0219876212400063
[10]	王毅, 胡志强, 郭维东. 考虑坝体柔性的重力坝坝面地震动水压力计算[J]. 水利水电科技进展, 2019, 39 (3): 38-43. WANG Yi, HU Zhiqiang, GUO Weidong. Calculation of earthquake-induced hydrodynamic pressure considering gravity dam flexibility[J]. Advances in Science and Technology of Water Resources, 2019, 39 (3): 38-43. (in Chinese)
[11]	李上明. 基于比例边界有限元法的坝库瞬态耦合分析[J]. 华中科技大学学报(自然科学版), 2011, 39 (9): 108-111. LI Shangming. Dam-reservoir transient coupling analysis based on SBFEM[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2011, 39 (9): 108-111. (in Chinese)
[12]	李上明. 基于比例边界有限元法动态刚度矩阵的坝库耦合分析方法[J]. 工程力学, 2013, 30 (2): 313-317. LI Shangming. Transient analysis method for dam-reservoir interaction based on dynamic stiffness of sbfem[J]. Engineering Mechanics, 2013, 30 (2): 313-317. (in Chinese)
[13]	李上明, 吴连军. 基于连分式与有限元法的坝库耦合瞬态分析方法[J]. 工程力学, 2016, 33 (4): 9-16. LI Shangming, WU Lianjun. Dam-reservoir interaction transient-method based on continued fraction formulation and fem[J]. Engineering Mechanics, 2016, 33 (4): 9-16. (in Chinese)
[14]	WANG X, JIN F, PREMPRAMOTE S, et al. Time-domain analysis of gravity dam-reservoir interaction using high-order doubly asymptotic open boundary[J]. Computers & Structures, 2011, 89 (7/8): 668-680.
[15]	QU Y L, CHEN D H, LIU L, et al. A direct time-domain procedure for the seismic analysis of dam-foundation-reservoir systems using the scaled boundary finite element method[J]. Computers and Geotechnics, 2021, 138 : 104364. doi: 10.1016/j.compgeo.2021.104364
[16]	许贺, 邹德高, 孔宪京. 基于FEM-SBFEM的坝-库水动力耦合简化分析方法[J]. 工程力学, 2019, 36 (12): 37-43. XU He, ZOU Degao, KONG Xianjing. A simplified dam-reservoir dynamic coupling analysis method based on fem-sbfem[J]. Engineering Mechanics, 2019, 36 (12): 37-43. (in Chinese)
[17]	LI S M. Diagonalization procedure for scaled boundary finite element method in modeling semi-infinite reservoir with uniform cross-section[J]. International Journal for Numerical Methods in Engineering, 2009, 80 (5): 596-608. doi: 10.1002/nme.2647
[18]	WEBER B. Rational Transmitting Boundaries for Time-Domain Analysis of Dam-Reservoir Interaction[M]. New York: Springer-Verlag, 1994.
[19]	TSAI C S. Semi-analytical solution for hydrodynamic pressures on dams with arbitrary upstream face considering water compressibility[J]. Computers & Structures, 1992, 42 (4): 497-502.

施引文献

资源附件(0)

访问统计

图(15)

计量

文章访问数: 116
HTML全文浏览量: 35
PDF下载量: 38
被引次数: 0

0. 引言
1. 无限水库控制方程与边界条件
2. 顺河向激励下远场SBFEM控制方程
2.1 频域控制方程
2.2 SBFEM系数矩阵
2.3 横截面垂直于顺河向的SBFEM系数矩阵
2.4 时域控制方程
3. 横河向、顺河向、竖直向等激励方向下远场控制方程
4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式
5. 算例
5.1 频域分析
5.2 时域分析
6. 结论

0. 引言
1. 无限水库控制方程与边界条件
2. 顺河向激励下远场SBFEM控制方程
2.1 频域控制方程
2.2 SBFEM系数矩阵
2.3 横截面垂直于顺河向的SBFEM系数矩阵
2.4 时域控制方程
3. 横河向、顺河向、竖直向等激励方向下远场控制方程
4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式
5. 算例
5.1 频域分析
5.2 时域分析
6. 结论

参考文献(19)

施引文献

资源附件(0)

访问统计

姓名
邮箱
手机号码
标题
留言内容
验证码

留言板

地震激励下等横截面无限水库频域和时域响应的FEM-SBFEM计算方法

doi: 10.21656/1000-0887.450138

通讯作者:
李上明(1978—)，男，研究员，博士，博士生导师(通讯作者. E-mail: hustmingsl@126.com)

计量

An FEM-SBFEM Coupled Method for Infinite Reservoir Responses With Uniform Cross Sections Under Seismic Excitations in Frequency and Time Domains

0. 引言

1. 无限水库控制方程与边界条件

2. 顺河向激励下远场SBFEM控制方程

2.1 频域控制方程

2.2 SBFEM系数矩阵

2.3 横截面垂直于顺河向的SBFEM系数矩阵

2.4 时域控制方程

3. 横河向、顺河向、竖直向等激励方向下远场控制方程

4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式

5. 算例

5.1 频域分析

5.2 时域分析

6. 结论

计量

目录

0. 引言

1. 无限水库控制方程与边界条件

2. 顺河向激励下远场SBFEM控制方程

2.1 频域控制方程

2.2 SBFEM系数矩阵

2.3 横截面垂直于顺河向的SBFEM系数矩阵

2.4 时域控制方程

3. 横河向、顺河向、竖直向等激励方向下远场控制方程

4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式

5. 算例

5.1 频域分析

5.2 时域分析

6. 结论

留言板

地震激励下等横截面无限水库频域和时域响应的FEM-SBFEM计算方法

doi: 10.21656/1000-0887.450138

通讯作者: 李上明(1978—)，男，研究员，博士，博士生导师(通讯作者. E-mail: hustmingsl@126.com)

计量

出版历程

An FEM-SBFEM Coupled Method for Infinite Reservoir Responses With Uniform Cross Sections Under Seismic Excitations in Frequency and Time Domains

0. 引言

1. 无限水库控制方程与边界条件

2. 顺河向激励下远场SBFEM控制方程

2.1 频域控制方程

2.2 SBFEM系数矩阵

2.3 横截面垂直于顺河向的SBFEM系数矩阵

2.4 时域控制方程

3. 横河向、顺河向、竖直向等激励方向下远场控制方程

4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式

5. 算例

5.1 频域分析

5.2 时域分析

6. 结论

计量

出版历程

目录

0. 引言

1. 无限水库控制方程与边界条件

2. 顺河向激励下远场SBFEM控制方程

2.1 频域控制方程

2.2 SBFEM系数矩阵

2.3 横截面垂直于顺河向的SBFEM系数矩阵

2.4 时域控制方程

3. 横河向、顺河向、竖直向等激励方向下远场控制方程

4. 含等横截面无限水域的任意构型无限水库FEM-SBFEM耦合公式

5. 算例

5.1 频域分析

5.2 时域分析

6. 结论

通讯作者:
李上明(1978—)，男，研究员，博士，博士生导师(通讯作者. E-mail: hustmingsl@126.com)