Numerical Simulation of Dynamic Interaction Analysis of Cap Pile-soil under Earthquake Action
-
摘要: 基于已开展的非液化场地-群桩基础-结构体系动力响应大型振动台模型试验，进行三维全时程动力数值模拟分析。采用修正的Davidenkov模型反映土体在地震反应过程中的模量衰减，通过“捆绑边界”模拟模型箱的层状剪切运动。通过对比试验中土-结构体系加速度响应时程、土体位移和桩基内力等，验证数值模型的有效性。利用已验证的数值模型，开展承台尺寸对桩-土-上部结构动力响应影响研究。结果表明，承台厚度的增大会导致上部结构和桩顶惯性效应减小；地震作用下沿激振方向前桩大于后桩，随着承台厚度的增大，前桩与后桩峰值弯矩差值率为16.1%～32.1%，群桩效应影响增大；随着承台厚度的增大，承台-土动土压力增大了3～6倍，承台与桩基水平荷载分担比增大，桩基弯矩反弯点位置上移了0.50 m；承台-土的相互摩擦作用会降低结构整体动力响应。
-
关键词:
- 非液化场地 /
- 承台-桩-土动力相互作用 /
- 地震响应 /
- 数值模拟 /
- 试验研究
Abstract: Based on the large-scale shaking table test, the dynamic response of a non-liquefied site pile group foundation system is analyzed using three-dimensional full-time dynamic simulation. In order to simulate the scenario, the modified Davidenkov model is used to represent soil modulus attenuation under seismic response analysis. Additionally, the "binding boundary" technique is employed to describe the layered shear movement of the model box. The impact of pile cap size on the dynamic response of the pile-soil-superstructure is investigated. The numerical model is validated by comparing the acceleration response time history, soil displacement, and pile internal force. The findings demonstrate that increased pile cap thickness can decrease inertial effect of superstructure and top pile, the response of front pile is greater than that of back pile with the difference bing 16.1% to 32.1%, and the group pile effect is increased. With increaing the thickness of cap, the soil dynamic pressure increases 3~6 times, and pile cap- pile lateral load shearing ratio increases, and the bending moment reverse bending point of pile body moves upward about 0.5m along the buried depth. The friction interaction between cap and soil can decrease the overal structure response. -
引言
礼乐盆地有着巨大的油气资源潜力（孙龙涛等，2010；李颂等，2012），但其位于南海南部，处于环太平洋地震带上，地震活动频繁。为保障我国对海洋资源的开发，研究海底盆地的地震动规律具有重要意义。目前用来研究地震波在复杂介质中传播与散射的主要方法有：边界元法（Lee，2013）、有限元法（苏波等，2019）、有限差分法（Graves等，1998；姚铭等，2017）、伪谱法（谭文卓等，2020）以及谱元法（刘少林等，2021）。
谱元法由 Patera（1984）提出，它结合有限元法和伪谱法的思想，先在每个单元上进行谱展开，然后按照有限元的思想进行集成运算，因此该方法不但可以模拟任意复杂介质，而且具有较高的精度。相较边界元法，谱元法的使用条件更为宽松，无需求出问题的基本解；与有限元法相比（Antonietti等，2018），谱元法具有伪谱法的优势，即采用快速傅里叶变换，把波场中的空间域转化到波数域，避免了对空间求导的差分运算，可以有效减少数值频散和伪波；与有限差分法不同的是（彭浩天，2018），谱元法可以在模拟任意复杂形状的模型时，避免产生较多的散射波；比伪谱法更方便的是，由于具有有限元法的特点，所以谱元法可以自动满足自由边界条件。因此本文采用谱元法进行海底盆地地震响应的研究。
在流固耦合方面，李伟华等（2003）和赵成刚等（2008）给出了不同场地模型的波动解析解，但一般只能应用于简单规则地形。Komatitsch等（2000）将谱元法引入场地地震动的流固耦合分析中，开发了相应的显式并行计算软件Specfem 2D 与Specfem 3D。孔曦骏等（2021）使用软件Specfem 2D研究了几何地形、水深、地震波卓越频率以及土体性质对江水地震动水压力的影响规律。以上学者对流固耦合问题进行了大量研究，但并未对考虑流固耦合作用的海底盆地进行深入研究。在边界条件的处理上，采用开源软件Specfem 2D的内置完美匹配层PML（Perfectly Matched Layer）边界（董阳，2016），PML的思想是在模型区域外部加1层吸收层，将模型区域内部和吸收层内的波动方程耦合起来，使得在接触面处没有反射，并在吸收层上使用指数形式的衰减层，当波场传播至吸收层时会被衰减，不会发生反射。
大量研究表明，震源入射角对场地地震响应具有显著影响，丁海平等（2017）采用有限元法，对凹陷地形的地震动峰值放大系数进行了研究，发现放大系数受震源入射角的影响较大，但其并未考虑到盆地内软土覆盖层的影响。朱重洋等（2016）采用有限元法，在考虑覆盖层影响的情况下，研究了震源入射角度对放大系数的影响，但其仅考虑了1层覆盖层，对于覆盖层的研究依然不够深入。魏成前（2021）在考虑盆地内多层覆盖层的情况下，采用二维有限元法研究了震源入射角度对盆地地震动放大特征的影响，发现这些因素的影响较为显著，但其盆地模型为规则的梯形，并不能完全符合真实情况。以上学者在研究震源入射角对地震响应的影响时，建立的模型都比较理想化，不能足够准确地反映真实情况，因此本文使用Fortran语言编写了针对开源软件Specfem 2D的数据处理程序，以南海礼乐盆地为背景，建立了考虑海水作用和盆地内不规则地形以及多层介质影响的盆地模型，研究了震源入射角对盆地地震响应的影响，本文结论可为实际工程分析与设计提供有效的参考。
1. 计算方法以及边界条件
1.1 谱元法理论
谱元法的基本思想是先将计算区域划分为有限个单元，然后在每个单元上使用伪谱法，通过在每个单元中配置不均匀的分布节点，把单元的近似解表示为截断的正交多项式展开式，最后采用Galerkin方法求解得到全局的近似解，从而完成对区域内波传播过程的数值模拟（万子轩，2020）。
波动方程可以表示为：
$$ \rho {\partial }_{t}^{2}\boldsymbol{u}={\boldsymbol{\nabla}} \cdot \boldsymbol{\sigma }+\boldsymbol{f} $$ (1) 式中，
$ \rho $ 是密度；$ \boldsymbol{u} $ 是待求解的位移向量；${\boldsymbol{\nabla}}$ 是梯度算子；$ \boldsymbol{\sigma } $ 是应力张量，$ \boldsymbol{\sigma } $ 与位移的梯度呈线性关系；$ \boldsymbol{f} $ 为震源项。在数值求解时，首先需要对式（1）中的微分形式进行积分。积分的具体方式为将式（1）乘以权函数
$ \boldsymbol{w} $ ，并在整个求解区域上进行分步积分，从而得到相应的积分方程：$$ \iint_{\varOmega }\rho \boldsymbol{w}\cdot {\partial }_{t}^{2}\boldsymbol{u}{\rm{d}}\Omega =-\iint _{\varOmega }\nabla \boldsymbol{w}:\boldsymbol{\sigma }{\rm{d}}\varOmega +\iint _{\varOmega }\boldsymbol{w}\cdot \boldsymbol{f}{\rm{d}}\varOmega +\int _{\varGamma }\boldsymbol{w}\cdot \boldsymbol{t}{\rm{d}}\varGamma $$ (2) 式中，
$ \varOmega $ 是计算区域；$ \varGamma $ 是其边界；$ \boldsymbol{t} $ 为边界上的力。与有限元法一样，谱元法在求解时，需要将计算区域离散为一个个单元，对于单元内待求的未知量u，使用一系列在基点上的
$ {u}_{i}=u\left({\xi }_{i},{\eta }_{i}\right) $ ，$ i $ =1, …, n和一系列的形函数$ {N}_{i}\left(\xi ,\eta \right)，i $ =1, …, n来表示。映射单元有2套坐标系。1套是局部坐标系
$ \left(\xi ,\eta \right) $ ，另1套是全局的笛卡尔坐标系$ \left(x,y\right) $ ，坐标系$ \left(x,y\right) $ 和坐标系$ \left(\xi ,\eta \right) $ 通过坐标映射关系进行转换。$$ \left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}{f}_{x}\left(\xi ,\eta \right)\\ {f}_{y}\left(\xi ,\eta \right)\end{array}\right] $$ (3) 局部坐标系和整体坐标系的转换关系为：
$$ \left[\begin{array}{c}\dfrac{\partial {N}_{i}}{\partial \xi }\\ \dfrac{\partial {N}_{i}}{\partial \eta }\end{array}\right]=\left[\begin{array}{c}\begin{array}{cc}\dfrac{\partial x}{\partial \xi }& \dfrac{\partial y}{\partial \xi }\end{array}\\ \begin{array}{cc}\dfrac{\partial x}{\partial \eta }& \dfrac{\partial y}{\partial \eta }\end{array}\end{array}\right]\left[\begin{array}{c}\dfrac{\partial {N}_{i}}{\partial x}\\ \dfrac{\partial {N}_{i}}{\partial y}\end{array}\right]=\boldsymbol{J}\left[\begin{array}{c}\dfrac{\partial {N}_{i}}{\partial x}\\ \dfrac{\partial {N}_{i}}{\partial y}\end{array}\right] $$ (4) 式中，
$ \boldsymbol{J} $ 为雅可比矩阵。面积元的计算：
$$ \mathrm{d}x\mathrm{d}y=\left|\boldsymbol{J}\right|\mathrm{d}\xi \mathrm{d}\eta \text{，} \boldsymbol{J} = \frac{{{\partial}} \left(x,y\right)}{{{\partial}} \left(\xi ,\eta \right)} $$ (5) 由于谱元法基于对微分方程的积分形式来处理，因此刚度矩阵、阻尼矩阵、质量矩阵形式均为：
$$ \iint_{S}\boldsymbol{G}{\rm{d}}S=\int_{-1}^{1}\int _{-1}^{1}\boldsymbol{G}\left(\boldsymbol{x}\left(\xi ,\eta \right),\boldsymbol{y}\left(\xi ,\eta \right)\right)\left|\boldsymbol{J}\right|\mathrm{d}\xi \mathrm{d}\eta $$ (6) 式中，S为计算区域；
$ \boldsymbol{G} $ 表示刚度矩阵、阻尼矩阵或质量矩阵。积分区间选取[−1，1]内的点
$ {\xi }_{1} $ ,$ {\xi }_{2} $ , …,$ {\xi }_{n} $ 和权重值$ {H}_{1} $ ,$ {H}_{2} $ , …,$ {H}_{n} $ ，使得：$$ I=\int _{-1}^{1}\boldsymbol{G}\left(\xi \right){\rm{d}}\xi \approx \sum _{i=1}^{n}{H}_{i}\boldsymbol{G}\left({\xi }_{i}\right) $$ (7) 谱元法的取值点
$ {\xi }_{i} $ （i=0, …, n）是Lagrange插值多项式中所选取的节点，同时满足n+1阶Legendre多项式的根：$$ \left(1-{\xi }^{2}\right){P}'_{n}\left(x\right)=0 $$ (8) 该取值点被称为Gassus-Lobatto-Legendre点（简称“GLL点”）。
$ {P}'_{n} $ 是n阶Legendre多项式的导数。GLL点具有正交的特点，所以质量矩阵为：$$ {\boldsymbol{M}}_{\bf{e}} = \iint_{{\boldsymbol{\varOmega }}_{\bf{e}}}\rho \boldsymbol{w}\mathrm{d}x\mathrm{d}y = \int_{-1}^{1}\int _{-1}^{1}\rho \boldsymbol{w}\left(\boldsymbol{x}\left(\xi ,\eta \right),\boldsymbol{y}\left(\xi ,\eta \right)\right)\left|\boldsymbol{J}\right|\mathrm{d}\xi \mathrm{d}\eta $$ (9) 式中，
${\boldsymbol{\varOmega }}_{\bf{e}}$ 为计算区域。最终得到：
$$ \boldsymbol{M}\ddot{\boldsymbol{U}}+\boldsymbol{C}\dot{\boldsymbol{U}}+\boldsymbol{K}\boldsymbol{U}=\boldsymbol{F} $$ (10) 式中，
$ \boldsymbol{M} $ 代表质量矩阵；$ \boldsymbol{K} $ 代表刚度矩阵；$ \boldsymbol{F} $ 代表外力。1.2 考虑流-固耦合影响的运动方程
由于Legendre谱元法形成的质量矩阵为对角阵形式，便于建立显式动力算法，故在大型波动分析中应用广泛。采用处理流体动力方程的速度势格式，基于谱元法的流-固整体质量矩阵仍然保持对角阵形式，运动方程可表示为（孔曦骏等，2021）：
$$ \left[\begin{array}{cc}{\boldsymbol{M}}_{\mathbf{s}}& 0\\ 0& {\boldsymbol{M}}_{\mathbf{f}}\end{array}\right]\left[\begin{array}{c}\ddot{\boldsymbol{U}}\\ \ddot{\varPhi }\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{D}}_{\mathbf{s}}& \boldsymbol{A}\\ {\boldsymbol{A}}^{{\rm{T}}}& {\boldsymbol{D}}_{\mathbf{f}}\end{array}\right]\left[\begin{array}{c}\dot{\boldsymbol{U}}\\ \dot{\varPhi }\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{K}}_{\mathbf{s}}& 0\\ 0& {-\boldsymbol{K}}_{\mathbf{f}}\end{array}\right]\left[\begin{array}{c}\boldsymbol{U}\\ \varPhi \end{array}\right]=\left[\begin{array}{c}{\boldsymbol{F}}_{\mathbf{w}\mathbf{s}}\\ 0\end{array}\right] $$ (11) 式中，
$ \boldsymbol{U} $ 为固体位移矢量；$ \varPhi $ 为流体速度势；速度矢量$ \boldsymbol{V}=\nabla \varPhi $ ；$ {\boldsymbol{M}}_{\mathbf{s}} $ 为固体质量矩阵；$ {\boldsymbol{M}}_{\mathbf{f}} $ 为流体质量矩阵；$ {\boldsymbol{D}}_{\mathbf{s}} $ 为固体吸收边界相关的矩阵；$ {\boldsymbol{D}}_{\mathbf{f}} $ 为流体吸收边界相关的矩阵；$ \boldsymbol{A} $ 为流固2种介质的耦合矩阵，这是由交界面的位移和应力的连续条件产生；$ {\boldsymbol{K}}_{\mathbf{s}} $ 为固体刚度矩阵；$ {\boldsymbol{K}}_{\mathbf{f}} $ 为流体刚度矩阵；$ {\boldsymbol{F}}_{\mathbf{w}\mathbf{s}} $ 为固体或流体中波源的等效荷载矢量。流体压力
$ p $ 与速度势$ \varPhi $ 的关系为：$$ p=-\rho \dot{\varPhi } $$ (12) 式中，
$ \rho $ 为流体密度。自由液面处的压力为0，该处的边界条件为：
$$ \varPhi =0 $$ (13) 流固交界面的法向速度和法向力连续：
$$ \boldsymbol{n}\cdot \dot{\boldsymbol{U}}=\boldsymbol{n}\cdot \nabla \varPhi $$ (14) $$ \boldsymbol{n}\cdot {\boldsymbol{\sigma }}_{\mathbf{s}}=\rho \dot{\varPhi } $$ (15) 式中，
$ \boldsymbol{n} $ 为流固交界面处法向单位向量；$ {\boldsymbol{\sigma }}_{\mathbf{s}} $ 为固体介质中的应力矢量。1.3 PML边界条件
为了避免内部入射波在界面发生反射，本文引入了PML（完全匹配层）吸收边界条件。如图1所示，根据PML原理沿水平方向x引入衰减因子
$ {d}^{x} $ ；沿竖直方向y引入衰减因子$ {d}^{y} $ 。$ {d}^{x} $ 和$ {d}^{y} $ 一般选择为（彭浩天，2018）：$$ {d}^{x}=\left\{\begin{array}{c}\begin{array}{l}-\dfrac{{v}_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{ln}\alpha }{L}\left[a\dfrac{{x}_{i}}{L}+b{\left(\dfrac{{x}_{i}}{L}\right)}^{2}\right]，\mathrm{匹}\mathrm{配}\mathrm{层}\mathrm{区}\mathrm{域}\\ 0，\mathrm{非}\mathrm{匹}\mathrm{配}\mathrm{层}\mathrm{区}\mathrm{域}\end{array}\end{array}\right. $$ (16) $$ {d}^{y}=\left\{\begin{array}{l}l-\dfrac{{v}_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{ln}\alpha }{L}\left[a\dfrac{{y}_{i}}{L}+b{\left(\dfrac{{y}_{i}}{L}\right)}^{2}\right]，\mathrm{匹}\mathrm{配}\mathrm{层}\mathrm{区}\mathrm{域}\\ 0，\mathrm{非}\mathrm{匹}\mathrm{配}\mathrm{层}\mathrm{区}\mathrm{域}\end{array}\right. $$ (17) 式中，
$ {x}_{i} $ 、$ {y}_{i} $ 分别为匹配层区域与内部区域界面的横向距离、纵向距离。$ {v}_{\mathrm{m}\mathrm{a}\mathrm{x}} $ 为最大纵波速度值；L为匹配层厚度，在不考虑流固耦合的情况下，通常PML厚度为半个波长时，则可达到满意的吸收效果。一般$ \alpha $ =1e-6，$ a= $ 0.25，$ b= $ 0.75（彭浩天，2018）。1.4 方法验证
本文选取Cristini等（2012）中基于谱元法的水声传播数值模型进行验证，并研究了网格大小对计算精度的影响。计算简图如图2所示，模型长8000 m，海底深2400 m，海水深400 m，声源作用在点M（400，395）处，接收器位于点N（800，380）处。假设海水为无黏、无旋的理想流体，海底介质为均匀和各向同性材料。海水材料参数为：密度
$ {\rho }_{1}= $ 1000$ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ ，压缩波声速$ {V}_{{\rm{p}}1}=1\,500\;\mathrm{ }\mathrm{m}/\mathrm{s} $ 。海底材料参数为：密度$ {\rho }_{2}= $ 2000$ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ ，压缩波声速$ {V}_{{\rm{p}}2}=2\,400\;\mathrm{m}/\mathrm{s} $ ，剪切波声速$ {V}_{{\rm{s}}2}=1\,200\;\mathrm{m}/\mathrm{s} $ 。模型上侧采用自由边界，左侧、右侧和下侧边界采用PML吸收单元。声源选择主频率为30 Hz的Ricker波，波长λ为40 m，波源时程曲线如图3所示。模型1、2、3网格大小分别为5 m、10 m和20 m，边界均为PML边界。计算得到N点的压强曲线，将其峰值取为1，并与Cristini等（2012）结果进行对比。由图4可以看出，当网格尺寸为λ/2（20 m）时（模型3），计算结果存在误差，误差主要为曲线的不规则波动；当网格尺寸为λ/4波长（10 m）时（模型2），计算结果与Cristini等（2012）结果比较接近，也与网格尺寸为λ/8（5 m）时（模型1）的结果一致，说明本文所用研究方法是正确的，且考虑流固耦合问题时，网格尺寸为λ/4时即可取得较高精度。1.5 PML边界参数分析
为研究PML边界吸收效率和厚度的关系，采用图5所示计算模型，该模型除了震源、观测点位置以及材料参数外，其余参数均与图2一样。震源作用在点M（4000，1400）处，以集中力的方式施加，观测点为点N（5000，1400），海水材料参数为：密度
$ {\rho }_{{\rm{w}}}= $ 1000$ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ ，压缩波速$ {V}_{{\rm{pw}}}=1\,500\;\mathrm{ }\mathrm{m}/\mathrm{s} $ 。海底材料参数为：密度$ {\rho }_{{\rm{s}}}= $ 2000$ \mathrm{k}\mathrm{g}/{\mathrm{m}}^{3} $ ，压缩波速$ {V}_{{\rm{ps}}}=1\,000\;\mathrm{m}/\mathrm{s} $ ，剪切波速$ {V}_{{\rm{ss}}}=200\;\mathrm{m}/\mathrm{s} $ 。输入主频率为2 、2.5 、3 Hz的Ricker波，波长分别为100、80、67 m，网格大小为10 m，时程曲线如图6所示。PML边界厚度h分别取60、90、120、150 m，得到观测点的加速度时程曲线如图7所示。从图中可以看到，当h为60 m时，会存在较多反射波，导致加速度不断增大，随着h增大，反射波逐渐减少，直至h为120 m时反射波消失。这说明在考虑流固耦合的情况下，PML边界厚度为2个波长时方可达到理想吸能效果。2. 礼乐盆地模型的建立
礼乐盆地模型如图8所示，地形采用张田升等（2019）对礼乐盆地的勘探数据，并将模型等比例缩小10倍。介质采用Chen等（2021）的数据，如表1所示。假定海水为无黏、无旋的理想流体，海底介质为均匀性和各向同性材料。模型长5000 m，海底介质分别为黏土、砂土和岩石，厚度均为60 m，海水深度为0 ～143.75 m，震源坐标分别为I（2500，500）和G（0，500），以集中力的方式施加，震源从盆地中部入射时，震源为I点，震源从盆地左侧入射时，震源为G点，震源入射角为α，分别取15°、30°、45°和60°。模型最大网格尺寸为10 m，时间步长为0. 2 ms。模型上侧采用自由边界，左侧、右侧、下侧均为PML吸收边界。流固耦合面上布置122个观测点，用以观测流固耦合面上不同位置的地震动峰值放大系数，取A、B、C、D和E 5个特征点，用以观测加速度时程以及频谱。震源分别采用Friuli波、Imperial_Valley波和Ricker子波，输入P_SV型地震波，时程以及频谱曲线如图9所示。
表 1 模型介质参数Table 1. Model media parameters介质 $ \rho $/（kg·m−3） Vp/（m·s−1） Vs/（m·s−1） 海水 1000 1500 — 黏土 1650 1650 218 砂土 1800 1697 264 岩石 2100 2135 485 3. 震源入射角对海底盆地地震响应的影响
取水平地震动峰值放大系数
$ \beta = $ |$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x}} $ /$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}} $ |，即观测点地表水平地震动峰值$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x}} $ 与输入点水平地震动峰值$ {{A}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}} $ 比值的绝对值（朱重洋等，2016）。本研究编写了提取观测点地震响应峰值的Fortran程序，该程序可以批量读取结果文件，并按照观测点顺序提取观测点地震响应峰值，有效减少了工作量。3.1 地震波由盆地中部入射情形
3.1.1 时域放大特征分析
图10为不同震源入射情况下观测点地震动峰值放大系数，可以发现：（1）震源的不同仅影响放大系数的大小，而不会影响其变化规律。不论采用何种震源，盆地中部的放大系数均为最大，盆地两侧的放大系数则最小，二者最大相差约18倍。盆地内的观测点放大系数分布曲线变化剧烈，这可能是由于模型考虑了流固耦合作用、盆地内大量凹凸地形以及多层介质的影响，导致地震波在传播过程中进行了多次反射和折射，从而使放大系数分布曲线变化剧烈，这是以往研究中没有发现的问题。图10（a）中最大放大系数峰值为0.68，图10（b）中最大放大系数峰值为0.62，图10（c）中最大放大系数峰值则为1.8，图10（c）中的放大系数峰值分别是图10（a）、图10（b）中放大系数峰值的2.65倍、2.9倍，这是因为图10（a）和图10（b）工况下的地震波作用时间较长，达到30 s，因此由盆地表面反射回来的反射波，会与后续震源部位的入射波产生相互作用，导致输入点位置的加速度峰值变大，根据本文对放大系数的定义可知，这种影响会使图10（a）和图10（b）的放大系数较图10（c）小。（2）随着震源入射角的增大，放大系数逐渐减小，该结论与丁海平等（2017）的结论一致。且入射角每增加15°，放大系数减少20%～30%。
图 10 观测点放大系数分布Figure 10. Amplification coefficient distribution of observation points3.1.2 特征点的频谱分析
图11为特征点加速度时程以及频谱的对比图，可以看出：（1）在不同震源入射角下，盆地不同位置的加速度时程差异较大。在盆地左侧（A点），入射角为60°时，加速度峰值具有滞后性，且滞后3 s左右；而在盆地其他部位，不同入射角下的加速度则在同一时间达到峰值，且越靠近盆地中部，加速度到达峰值的时间就越短。（2）盆地不同位置的频谱特征具有明显差异。在盆地边缘部位（A点与E点），谱曲线的波动比较强烈，震源入射角的改变不仅会改变谱峰的大小，也会改变谱曲线的变化规律，且当震源为Ricker波时，这种改变更为明显。这与禹乐等（2020）使用Ricker波进行研究而得出的结论一致。而在盆地中部（C点），谱曲线的波动则更具规律性，震源入射角的改变不会使谱峰位置发生改变，仅会影响谱峰大小，且谱峰随震源入射角的增大而增大。（3）震源对不同震源入射角下的频谱特征具有较大影响。当震源为Friuli 波时，在盆地边缘部位（A点与E点），有效幅值在25 Hz 之内，而在盆地中部（C点），则30 Hz 以上的幅值才可忽略；当震源为Imperial_Valley 波时，仅有A点处25 Hz 以上的幅值可以忽略，其余盆地部位的幅值均在33 Hz 左右才会显著减小；当震源为Ricker 波时，有效幅值的频率区间则无明显变化。（4）结合禹乐等（2020）的相关研究可以发现，由于考虑了海水作用、盆地内不规则地形以及多层介质的影响，特征点的加速度时程存在剧烈变化，这是由地震波的多次反射和折射导致的。
图 11 特征点加速度时程以及频谱Figure 11. Acceleration time histories and spectrum of characteristic points3.2 地震波由盆地左侧入射情形
3.2.1 时域放大特征分析
图12为不同震源入射情况下观测点地震动峰值放大系数，与地震波从中部入射的区别如下：（1）当地震波由左侧入射时，放大系数最大值均位于盆地左侧靠近盆地边缘部位，且随着入射角的增大向左偏移，而放大系数最小值则位于盆地右侧，二者最大相差14倍；图12（a）、图12（b）、图12（c）中的最大放大系数峰值分别为0.72、0.88、1.4，图12（c）中的放大系数峰值分别是图12（a）和图12（b）放大系数峰值的1.94倍、1.75倍。（2）随着震源入射角的增大，放大系数逐渐增大。与地震波由中部入射情形不同的是，入射角每增加15°，放大系数减少17%～30%。
图 12 观测点放大系数分布Figure 12. Amplification coefficient distribution of observation points3.2.2 特征点的频谱分析
图13为特征点加速度时程以及频谱的对比图，与地震波从中部入射的区别如下：（1）无论在盆地何位置，加速度时程及谱曲线均具有规律性，震源入射角的改变不会使二者峰值位置发生改变，仅会影响峰值大小，且加速度峰值和谱峰均随着震源入射角的增大而减小。越靠近盆地左侧，加速度峰值到来的时间越短。（2）震源对不同震源入射角下的频谱特征具有较大影响。当震源为Friuli 波时，在盆地左侧（A点），有效幅值在32 Hz之内，越靠近右侧的部位，有效幅值区间越小，在盆地右侧（E点），有效幅值区间在20 Hz之内；当震源为Imperial_Valley 波时，在盆地左侧（A点），有效幅值在33 Hz之内，而在盆地右侧（E点），则在20 Hz之内；当震源为Ricker 波时，有效幅值的频率区间无明显变化。
图 13 特征点加速度时程以及频谱Figure 13. Acceleration time histories and spectrum of characteristic pointsts4. 结论
本文运用谱元法理论建立了礼乐盆地地震动分析二维模型，研究了震源入射角对盆地地震响应的影响。得到以下结论：
（1）在考虑流固耦合的情况下，网格尺寸为四分之一波长时即可取得较高精度，PML边界厚度至少为2个波长时方可达到理想吸能效果。
（2）当地震波从盆地中部入射时，随着震源入射角的增大，放大系数逐渐减小，且入射角每增加15°，放大系数就减少20%～30%；盆地中部的放大系数最大，盆地两侧的最小，二者最大相差18倍左右；震源的不同会影响放大系数的大小，而不会影响其变化规律，但是对频谱特征具有较大影响；在不同的震源入射角下，盆地不同位置的加速度时程以及频谱特征差异较大。
（3）当地震波由左侧入射时，随着震源入射角的增大，放大系数逐渐减小，且入射角每增加15°，放大系数减少17%～30%；放大系数最大值均位于盆地左侧，靠近盆地边缘部位，随着入射角的增大而向左偏移，而放大系数最小值则位于盆地右侧，二者最大相差14倍左右；与地震波从中部入射相比，放大系数最大值变化了6%～30%；震源的不同会影响放大系数的大小以及频谱特征；在不同的震源入射角下，盆地不同位置的加速度时程以及频谱特征差异较小。
（4）由于考虑了海水、盆地内大量不规则地形以及多层介质的影响，地震波在传播过程中进行了多次反射和折射，致使盆地内的观测点放大系数分布曲线以及特征点的加速度时程曲线均存在剧烈变化。因此，在研究盆地地震相关问题时，这些因素的作用是不可忽视的。在后续的工作中，会对相关因素进行更加全面深入的研究。
-
图 1 非液化非自由场振动台试验传感器布置
Figure 1. Non-liquefaction non-free field shaking table test sensor layout
图 2 汶川地震卧龙台地震记录加速度时程曲线与傅里叶谱
Figure 2. Acceleration time history and Fourier spectra of Wenchuan earthquake recorded by Wolong station
图 3 不同时刻土体对桩的侧向约束力分布示意
Figure 3. Distribution diagram of lateral binding force of soil on pile at different time
图 4 非液化非自由场三维计算模型
Figure 4. Non-liquefaction free-field three-dimensional calculation model
图 5 三维计算模型输入加速度时程曲线
Figure 5. Input acceleration time-history for three-dimensional calculation model
图 6 Davidenkov模型应力-应变滞回曲线
Figure 6. Dynamic shear stress-strain hysteresis curves of Davidenkov model
图 7 振动台试验砂土Davidenkov拟合曲线
Figure 7. Davidenkov fitting curve of sandy soil for shaking table test
图 8 桩内土体加速度时程曲线对比
Figure 8. Comparison of acceleration time-history curves of soil in piles
图 9 远场土体加速度时程曲线对比
Figure 9. Comparison of acceleration time-history curve of soil in far field
图 10 远场土体加速度放大系数对比
Figure 10. Comparison of acceleration amplification factor of soil in far field
图 11 桩基-结构加速度时程曲线对比
Figure 11. Comparison of acceleration time-history curves of pile-structure
图 12 不同测点加速度反应谱对比
Figure 12. Comparison of acceleration response spectra with different test points
图 16 不同承台厚度结构加速度时程曲线对比
Figure 16. Comparison of acceleration time-history curves of structure with different thickness of cap
图 17 不同承台厚度结构位移时程曲线对比
Figure 17. Comparison of structural displacement time history with different thickness of cap
图 18 不同承台厚度桩基-结构最大时刻位移
Figure 18. Maximum time displacement of pile foundation and structure with different thickness of cap
图 19 不同承台厚度桩身最大时刻剪力
Figure 19. Maximum instantaneous shear force of pile body with different thickness of cap
图 21 不同承台厚度桩基最大时刻弯矩
Figure 21. Maximum moment bending moment of pile with different thickness of cap
图 22 不同承台厚度土体对桩的侧向约束力分布示意
Figure 22. Distribution diagram of lateral binding force of soil with different thickness of cap on pile
表 1 模型材料参数
Table 1. Physical parameters of model
土层编号
土层厚度/m
密度ρ/（kg·m−3）Davidenkov模型参数 最大剪切模量Gmax/MPa 泊松比ν A B γ0 γult 剪切波速vs/（m·s−1） 1 0.2 1 400 17 0.35 1.02 0.35 0.000 40 0.003 110 2 1.2 1 460 25 0.30 1.10 0.45 0.000 45 0.003 130 3 0.5 1 600 40 0.30 1.10 0.45 0.000 45 0.003 158 
下载: 导出CSV
表 2 结构模型材料参数
Table 2. Physical parameters of model
材料 弹性模量/GPa 密度ρ 泊松比 阻尼比 峰值抗压强度/MPa 峰值抗拉强度/MPa 混凝土 14 2 400 0.20 0.05 29.6 2.95 钢筋/H型钢 200 7 800 0.18 — 240.0 240.00
下载: 导出CSV
表 3 模型结构动力响应对比
Table 3. Comparison of dynamic responses of model structures
动力响应计算值 0.25 m厚承台 0.50 m厚承台 0.75 m厚承台 有摩擦 无摩擦 有摩擦 无摩擦 有摩擦 无摩擦 承台加速度峰值/g 3.50 3.46 2.91 3.06 2.46 2.51 承台位移峰值/mm 16.40 16.61 11.30 12.52 11.26 11.3 上部结构顶部加速度峰值/g 1.64 1.74 1.48 1.67 1.36 1.58 上部结构顶部位移峰值/mm 7.710 7.850 5.690 6.940 3.682 4.470 桩1弯矩峰值/（N·m） 191.30 201.60 198.30 206.80 211.80 213.30 桩1剪力峰值/N 393.30 411.90 582.60 604.70 783.40 788.70
下载: 导出CSV
-
陈国兴, 庄海洋, 2005. 基于Davidenkov骨架曲线的土体动力本构关系及其参数研究. 岩土工程学报, 27（8）: 860—864 doi: 10.3321/j.issn:1000-4548.2005.08.002Chen G. X. , Zhuang H. Y. , 2005. Developed nonlinear dynamic constitutive relations of soils based on Davidenkov skeleton curve. Chinese Journal of Geotechnical Engineering, 27(8): 860—864. (in Chinese) doi: 10.3321/j.issn:1000-4548.2005.08.002 陈跃庆, 吕西林, 李培振等, 2001. 分层土-基础-高层框架结构相互作用体系振动台模型试验研究. 地震工程与工程振动, 21（3）: 104—112 doi: 10.3969/j.issn.1000-1301.2001.03.019Chen Y. Q. , Lü X. L. , Li P. Z. , et al. , 2001. Shaking table testing for layered soil-foundation-structure interaction system. Earthquake Engineering and Engineering Vibration, 21(3): 104—112. (in Chinese) doi: 10.3969/j.issn.1000-1301.2001.03.019 邸博, 2017. 考虑有限基础刚度的高层建筑结构性能研究与应用. 哈尔滨: 哈尔滨工业大学.Di B., 2017. Structural performance study and application of high-rise buildings considering foundation stiffness. Harbin: Harbin Institute of Technology. (in Chinese) 范重, 刘涛, 陈巍等, 2017. 基础刚度对高层建筑抗震性能影响研究. 工程力学, 34（7）: 203—213 doi: 10.6052/j.issn.1000-4750.2016.08.0650Fan Z. , Liu T. , Chen W. , et al. , 2017. Study on the influence of foundation stiffness on the seismic performance of high-rise buildings. Engineering Mechanics, 34(7): 203—213. (in Chinese) doi: 10.6052/j.issn.1000-4750.2016.08.0650 尚守平, 卢华喜, 王海东等, 2006. 大比例模型结构-桩-土动力相互作用试验研究与理论分析. 工程力学, 23（S2）: 155—166Shang S. P. , Lu H. X. , Wang H. D. , et al. , 2006. Test investigation and theoretical analysis of large-scale model on dynamic soil-pile-structure interaction. Engineering Mechanics, 23(S2): 155—166. (in Chinese) 唐亮, 凌贤长, 徐鹏举等, 2010. 承台型式对可液化场地桥梁桩-柱墩地震响应影响振动台试验. 地震工程与工程振动, 30（1）: 155—160Tang L. , Ling X. Z. , Xu P. J. , et al. , 2010. Effect of cap type on seismic responses of bridge pile-column pier in liquefiable ground by shaking table test. Journal of Earthquake Engineering and Engineering Vibration, 30(1): 155—160. (in Chinese) 肖晓春, 2003. 地震作用下土-桩-结构动力相互作用的数值模拟. 大连: 大连理工大学.Xiao X. C., 2003. Numerical modeling of soil-pile-structure dynamic interaction under seismic excitation. Dalian: Dalian University of Technology. (in Chinese) 许成顺, 豆鹏飞, 杜修力等, 2022. 非液化土-群桩基础-结构体系相互作用动力响应振动台试验研究. 建筑结构学报, 43（5）: 185—194, 204 doi: 10.14006/j.jzjgxb.2020.0260Xu C. X. , Dou P. F. , Du X. L. , et al. , 2022. Dynamic interaction and seismic response of non-liquefiable soil-pile group foundation-structure system from shaking table test. Journal of Building Structures, 43(5): 185—194, 204. (in Chinese) doi: 10.14006/j.jzjgxb.2020.0260 赵晓光, 2020. 地震作用下建筑高低承台群桩基础响应规律试验研究. 北京: 中国建筑科学研究院.Zhao X. G., 2020. Experimental research on response law of pile group foundation with high and low cap under earthquake action. Beijing: China Academy of Building Research. (in Chinese) 中华人民共和国交通运输部, 2012. JTS 146—2012 水运工程抗震设计规范. 北京: 人民交通出版社.Ministry of Communications of the People’s Republic of China, 2012. JTS 146—2012 Code for seismic design of water transport engineering. Beijing: China Communications Press. (in Chinese) 中华人民共和国住房和城乡建设部, 2008. JGJ 94—2008 建筑桩基技术规范. 北京: 中国建筑工业出版社.Ministry of Housing and Urban-Rural Development of the People's Republic of China, 2008. JGJ 94—2008 Technical code for building pile foundations. Beijing: China Architecture & Building Press. (in Chinese) 中华人民共和国住房和城乡建设部, 中华人民共和国国家质量监督检验检疫总局, 2010. GB 50011—2010 建筑抗震设计规范. 北京: 中国建筑工业出版社.Ministry of Housing and Urban-Rural Development of the People's Republic of China, General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China, 2010. GB 50011—2010 Code for seismic design of buildings. Beijing: China Architecture & Building Press. (in Chinese) 中华人民共和国住房和城乡建设部, 中华人民共和国国家质量监督检验检疫总局, 2012. GB 50191—2012 构筑物抗震设计规范. 北京: 中国计划出版社.Ministry of Housing and Urban-Rural Development of the People's Republic of China, General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China, 2012. GB 50191—2012 Design code for antiseismic of special structures. Beijing: China Planning Press. (in Chinese) 朱斌, 熊根, 刘晋超等, 2013. 砂土中大直径单桩水平受荷离心模型试验. 岩土工程学报, 35（10）: 1807—1815Zhu B. , Xiong G. , Liu J. C. , et al. , 2013. Centrifuge modelling of a large-diameter single pile under lateral loads in sand. Chinese Journal of Geotechnical Engineering, 35(10): 1807—1815. (in Chinese) 朱志辉, 2006. 土-箱基-框架结构动力相互作用的试验研究与理论分析. 长沙: 湖南大学.Zhu Z. H. , 2006. Theoretical analysis and test research on soil-box foundation-structure dynamic interaction system. Changsha: Hunan University. (in Chinese) 庄海洋, 陈国兴, 梁艳仙等, 2007. 土体动非线性黏弹性模型及其ABAQUS软件的实现. 岩土力学, 28（3）: 436—442 doi: 10.3969/j.issn.1000-7598.2007.03.002Zhuang H. Y. , Chen G. X. , Liang Y. X. , et al. , 2007. A developed dynamic viscoelastic constitutive relations of soil and implemented by ABAQUS software. Rock and Soil Mechanics, 28(3): 436—442. (in Chinese) doi: 10.3969/j.issn.1000-7598.2007.03.002 Chandler A. M. , 1986. Building damage in Mexico City earthquake. Nature, 320(6062): 497—501. doi: 10.1038/320497a0 Deb P. , Pal S. K. , 2019. Numerical analysis of piled raft foundation under combined vertical and lateral loading. Ocean Engineering, 190: 106431. doi: 10.1016/j.oceaneng.2019.106431 Elahi H., Moradi M., Poulos H. G., et al., 2010. Pseudostatic approach for seismic analysis of pile group. Computers and Geotechnics, 37(1—2): 25—39. Flores J. , Novaro O. , Seligman T. H. , 1987. Possible resonance effect in the distribution of earthquake damage in Mexico City. Nature, 326(6115): 783—785. doi: 10.1038/326783a0 Kim B. T. , Yoon G. L. , 2011. Laboratory modeling of laterally loaded pile groups in sand. KSCE Journal of Civil Engineering, 15(1): 65—75. doi: 10.1007/s12205-011-0924-3 Li L. C. , Liu H. , Wu W. B. , et al. , 2021. Investigation on the behavior of hybrid pile foundation and its surrounding soil during cyclic lateral loading. Ocean Engineering, 240: 110006. doi: 10.1016/j.oceaneng.2021.110006 Maheshwari B. K. , Truman K. Z. , El Naggar M. H. , et al. , 2004. Three-dimensional finite element nonlinear dynamic analysis of pile groups for lateral transient and seismic excitations. Canadian Geotechnical Journal, 41(1): 118—133. doi: 10.1139/t03-073 Niemann C. , O'Loughlin C. , Tian Y. H. , et al. , 2019. Response of pile groups in sand due to lateral cyclic loading. International Journal of Physical Modelling in Geotechnics, 19(6): 318—330. doi: 10.1680/jphmg.18.00027 Talebi A. , Derakhshani A. , 2019. Estimation of P-multipliers for laterally loaded pile groups in clay and sand. Ships and Offshore Structures, 14(3): 229—237. doi: 10.1080/17445302.2018.1495542 Xiang B. S. , Huang B. , Yang Z. Y. , et al. , 2017. Influences of pile group effects on wave forces on an offshore bridge pile-cap foundation. Challenges, 8(2): 30. doi: 10.3390/challe8020030 Xu C. S. , Zhang Z. H. , Li Y. , et al. , 2020. Validation of a numerical model based on dynamic centrifuge tests and studies on the earthquake damage mechanism of underground frame structures. Tunnelling and Underground Space Technology, 104: 103538. doi: 10.1016/j.tust.2020.103538 期刊类型引用(0)
其他类型引用(1)
-
下载:
点击查看大图
计量
- 文章访问数: 168
- HTML全文浏览量: 27
- PDF下载量: 16
- 被引次数: 1
