MULTI COUPLING PERIODIC FINITE ELEMENT METHOD FOR CALCULATING THE DYNAMIC RESPONSE OF UNSATURATED SOIL-STRUCTURE SYSTEM
-
摘要: 真实的地基土体-隧道系统中土体及结构性质往往沿线路纵向变化. 为考虑土体与结构沿纵向的变化特性, 提出了一种非饱和土-结构系统动力响应分析的多耦合周期性有限元法. 首先基于非饱和土的实用波动方程, 采用 Galerkin 法推导了单节点5个自由度的非饱和土ub-pl-pg格式有限元表达式, 相比于单节点9个自由度的ub-v-w格式有限元表达式, 节省了计算资源. 其次引入复拉伸函数, 构建完美匹配层边界单元截断无限域. 最后采用多周期性模拟结构沿纵向的变化特性, 引入自由波传播理论, 结合各周期性结构间的连续性条件, 实现纵向各周期性结构间的耦合. 通过与既有的2.5维有限元-完美匹配层法和波函数法的计算结果进行对比, 验证了本文方法的可靠性. 该方法与既有的解析、数值方法相比, 具有可高效考虑结构纵向变化特性的优点. 基于该方法, 以非饱和土-隧道-隔离桩系统为例, 讨论了隔离桩的减(隔)振效果, 结果表明: 与咬合桩相比, 考虑隔离桩桩间距后地表动力响应规律有所改变, 因此在进行地铁隧道系统振动响应预测时, 应考虑系统结构沿纵向的变化特性.Abstract: In the actual foundation-tunnel system, the properties of soil and structures often change longitudinally along the tunnel line. In order to consider the longitudinal variation characteristics of structures, a multi coupling periodic finite element method for dynamic response analysis of unsaturated soil-structure system is proposed in this paper. Firstly, based on the practical wave equation of unsaturated soil, the finite element expression of ub-pl -pg scheme for unsaturated soil with 5 degrees of freedom at a single node is derived by Galerkin method, which can save the computing efficiency compared with the finite element expression of ub-v-w scheme with 9 degrees of freedom at a single node. Then, by introducing the stretching function with complex number, the perfectly matched layer boundary element is constructed to truncate the infinite domain. Finally, multi periodicity is used to simulate the longitudinal variation characteristics of the structure. By introducing the free wave propagation theory and combining the continuity conditions between periodic structures, the coupling between periodic structures can be realized. The results of the proposed method are compared with the results obtained by the existing 2.5 dimensional coupled FE-PML model and analytical method respectively, which verified the reliability of the proposed method. Compared with the existing analytical methods or numerical methods, the proposed method has the advantage of considering the longitudinal variation characteristics of the structure efficiently. Based on the proposed method, taking unsaturated soil-tunnel-isolation pile system as an example, the vibration reduction and isolation effect of isolation pile is discussed. The results show that compared with gnawing piles, the regular of surface dynamic response changes after considering the spacing between isolated piles. Therefore, in order to accurately predict the vibration response of subway tunnel system, the longitudinal variation characteristics of the system structure should be taken into consideration.
-
引 言
近年来, 地铁隧道车致环境振动问题日益突出, 严重影响着线路周边居民的正常生活、精密仪器的正常使用及古建筑与文物的保护等[1-3]. 为采取有效的减(隔)振措施, 迫切需要可靠的车致动力响应预测模型.
常见的地铁隧道动力响应预测模型主要包括解析(半解析)模型以及数值法模型[4]. 解析(半解析)模型主要有欧拉梁−地基模型[5-7]、Pip in Pip(PiP)模型[8-10]及波函数法模型[11-12]等. 解析(半解析)模型计算效率较高, 但难以处理不规则边界. 与解析(半解析)模型相比, 数值模型的优势在于可以进行更为精细化的建模, 便于处理复杂的边界(界面), 可考虑复杂的地基结构. 常见的数值模型有二维有限元模型[13-14]、三维有限元模型[15-17]. 二维有限元忽略了弹性波沿隧道纵向的传播、三维有限元则占用大量的计算机资源, 为兼顾二维有限元计算效率, 同时考虑隧道纵向的动力传播规律, 相关性学者假定模型沿纵向几何和物理性质的不变性, 基于纵向Fourier或Floquet变换, 分别提出了2.5维有限元[18-19]和周期性有限元模型[20-21]. 此外, 为避免人工截断模型边界产生的误差, 相关学者提出了边界元[22-23]、无限元[24-25]、比例边界有限元[26]、多阻尼层[27]等边界处理方法. 然而上述既有的数值模型仅将土体视为单相弹性或两相饱和多孔介质, 鲜见考虑地基土体的三相特性, 同时无法高效考虑模型的纵向变化特性.
为此本文基于非饱和土波动方程以及渗流运动方程, 结合应力、渗流边界条件, 采用 Galerkin 法, 推导了非饱和ub-pl-pg格式有限元方程. 同时引入拉伸函数构建完美匹配层(PML)处理边界, 提出了非饱和地基−结构动力响应分析的有限元-完美匹配层方法. 在此基础上进一步引入自由波传播理论, 将结构沿一个方向的变化处理为多耦合周期性结构, 基于自由波传播常数及传播向量进行各周期性结构间的耦合, 最终提出了多耦合周期性有限元法, 并进行了地基-隧道-排桩系统动力响应的算例分析.
1. 多耦合周期性有限元法
1.1 非饱和ub-pl-pg格式有限元推导
为实现频域内的模型求解, 首先定义关于时间变量t的傅里叶变化形式
$$ \tilde{f}(\omega )={\displaystyle {\int }_{-\infty }^{ + \infty }f(t){{\rm{e}}}^{-\text{i}\omega t}{\rm{d}}t}\text{, }\;\;f(t)=\frac{1}{2\text{π}}{\displaystyle {\int }_{-\infty }^{ + \infty }\tilde{f}(\omega ){{\rm{e}}}^{\text{i}\omega t}{\rm{d}}\omega } $$ (1) 式中, ω代表时间t对应的频率; 顶标“ ~ ”表示频域中的变量.
基于式(1)中傅里叶变化形式, 根据三相介质理论[28], 可得频域内三相介质波动控制方程及渗流方程分别为
$$ \left. \begin{array}{l} \mu {{\tilde u}_{bi,jj}} + (\lambda + \mu ){{\tilde u}_{bi,ij}} - {a_1}\gamma {\delta _{ij}}{{\tilde p}_{l,i}} - {a_1}(1 - \\ \;\;\;\;\;\;\;\;\;\;\gamma ){\delta _{ij}}{{\tilde p}_{g,i}} = - {\omega ^2}{{\overset{\frown} \rho }_s}{{\tilde u}_{bi}} - {\omega ^2}{{\overset{\frown} \rho }_l}{{\tilde v}_i} - {\omega ^2}{{\overset{\frown} \rho }_g}{{\tilde w}_i}\\ {B_{11}}{{\tilde p}_l} + {B_{12}}{{\tilde p}_g} + {B_{13}}{{\tilde u}_{b,j}} + {B_{14}}{{\tilde v}_{,j}} = 0\\ {B_{21}}{{\tilde p}_l} + {B_{22}}{{\tilde p}_g} + {B_{24}}{{\tilde u}_{b,j}} + {B_{25}}{{\tilde w}_{,j}} = 0 \end{array} \right\}$$ (2) $$ \left. \begin{array}{l} {\rm{i}}\omega n{S_r}({{\tilde v}_i} - {{\tilde u}_{b{\text{ }}i}}) = - \dfrac{{{k_l}}}{{{\rho _l}g}}({{\tilde p}_{w,i}} - {\omega ^2}{\rho _l}{{\tilde v}_i}) \hfill \\ {\rm{i}}\omega n(1 - {S_r})({{\tilde w}_i} - {{\tilde u}_{b{\text{ }}i}}) = - \dfrac{{{k_g}}}{{{\rho _g}g}}({{\tilde p}_{g,i}} - {\omega ^2}{\rho _g}{{\tilde w}_i}) \hfill \\ \end{array} \right\} $$ (3) 式中, B11 ~ B25表达式详见附录A. 其余各变量定义可参考文献[28].
整理式(3)可得
$$ \left. \begin{array}{l} {{\tilde v}_i} = ({M_l}{{\tilde u}_{b{\text{ }}i}} - {{\tilde p}_{l,i}})/({M_l} - {\omega ^2}{\rho _l}) \hfill \\ {{\tilde w}_i} = ({M_g}{{\tilde u}_{b{\text{ }}i}} - {{\tilde p}_{g,i}})/({M_g} - {\omega ^2}{\rho _g}) \hfill \\ \end{array} \right\} $$ (4) 式中,
${M}_{l}=\dfrac{{\rm{i}}\omega n{S}_{r}{\rho }_{l}g}{{k}_{l}}, {M}_{g}=\dfrac{{\rm{i}}\omega n(1-{S}_{r}){\rho }_{g}g}{{k}_{g}}$ .引入应力、渗流边界条件
$$ \left. \begin{array}{l} {{\tilde \sigma }_{ij}}{n_j} - {{\tilde f}_i} = 0 \hfill \\ {\rm{i}}\omega n{S_r}{{\tilde u}_l}{n_j} - {{\tilde c}_l} = 0 \hfill \\ {\rm{i}}\omega n(1 - {S_r}){{\tilde u}_g}{n_j} - {{\tilde c}_g} = 0 \hfill \\ \end{array} \right\} $$ (5) 式中, f, cl, cg分别表示边界节点处外力及输入的液体、气体流速.
基于Galerkin法, 在式(2)和式(5)的基础上, 引入虚位移δub以及虚孔压δpl, δpg, 同时利用式(4)在计算过程中消去v, w分量, 可得 ub-pl-pg格式的虚原理表达式为
$$\left. \begin{array}{l} \displaystyle\int\limits_V {\Big[ - \delta {{\tilde \varepsilon }_{ij}}{{\tilde \sigma }_b} + \delta {{\tilde u}_b}{\omega ^2}{c_{11}}{{\tilde u}_{bi}} + \delta {{\tilde \varepsilon }_{ij}}{a_1}\gamma {\delta _{ij}}{{\tilde p}_l} - } \\ \qquad \qquad \delta {{\tilde u}_b}{\omega ^2}{c_{12}}{{\tilde p}_{l,i}} + \delta {{\tilde \varepsilon }_{ij}}{a_1}(1 - \gamma ){\delta _{ij}}{{\tilde p}_g} - \\ \qquad \qquad \delta {{\tilde u}_b}{\omega ^2}{c_{13}}{{\tilde p}_{g,i}}\Big]{\rm{d}}V = \displaystyle\int\limits_S { {\delta {{\tilde u}_b}( - {{\tilde f}_i})} {\rm{d}}S} \\ \displaystyle\int\limits_V {(\delta {{\tilde p}_l}{c_{21}}{{\tilde u}_{b,j}} - \delta {{\tilde p}_{l,j}}{c_{22}}{{\tilde u}_b} + \delta {{\tilde p}_l}{c_{23}}{{\tilde p}_l} + \delta {{\tilde p}_{l,j}}{c_{24}}{{\tilde p}_{l,i}} + } \\ \qquad \qquad \delta {{\tilde p}_l}{c_{25}}{{\tilde p}_g}){\rm{d}}V = \displaystyle\int\limits_S { {\delta {{\tilde p}_l}\left( {\frac{{ - {{\tilde c}_l}}}{{{\rm{i}}\omega }}} \right)} {\rm{d}}S} \\ \displaystyle\int\limits_V {(\delta {{\tilde p}_g}{c_{31}}{{\tilde u}_{b,j}} - \delta {{\tilde p}_{g,j}}{c_{32}}{{\tilde u}_b} + \delta {{\tilde p}_g}{c_{33}}{{\tilde p}_l}} + \delta {{\tilde p}_g}{c_{35}}{{\tilde p}_g} + \\ \qquad \qquad \delta {{\tilde p}_{g,j}}{c_{36}}{{\tilde p}_{g,i}}){\rm{d}}V = \displaystyle\int\limits_S { {\delta {{\tilde p}_g}\left( {\frac{{ - {{\tilde c}_g}}}{{{\rm{i}}\omega }}} \right)} {\rm{d}}S} \end{array} \right\}$$ (6) 式中, 各变量表达式详见附录B.
采用空间8节点等参单元进行离散处理, 可得ub-pl-pg格式的三相介质有限元表达式为
$$ \left. \begin{array}{l} ({K_{11}} + {K_{12}}){{\tilde u}_b} + ({L_{11}} + {L_{12}}){{\tilde p}_l} + ({G_{11}} + {G_{12}}){{\tilde p}_g} = \tilde f \hfill \\ ({K_{21}} + {K_{22}}){{\tilde u}_b} + ({L_{21}} + {L_{22}}){{\tilde p}_l} + ({G_{21}} + {G_{22}}){{\tilde p}_g} = \dfrac{{ - {{\tilde c}_l}}}{{{\rm{i}}\omega }} \hfill \\ ({K_{31}} + {K_{32}}){{\tilde u}_b} + ({L_{31}} + {L_{32}}){{\tilde p}_l} + ({G_{31}} + {G_{32}}){{\tilde p}_g} = \dfrac{{ - {{\tilde c}_l}}}{{{\rm{i}}\omega }} \hfill \\ \end{array} \right\} $$ (7) 式中, 各变量表达式详见附录C.
1.2 完美匹配层单元
采用完美匹配层处理模型边界, 可在有限的单元内实现弹性波无反射的迅速衰减, 进而截断无限域. 因此引入拉伸函数[29], 推导三相非饱和多孔介质完美匹配层单元
$$ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over s} } = {s_0} + \int_{{s_0}}^{{s_p}} {\lambda (s){\rm{d}}s} \;\;({s_0} \leqslant s \leqslant {s_p}) $$ (8) 式中, 顶标“↔”表示拉伸含义.
弹性波在传播过程中可分为两种波, 即传播波(propagating)与隐逝波(evanescent), 因此定义拉伸函数为
$$ \lambda (s) = 1 + f_0^e{\left( {\frac{{s - {s_0}}}{{{L_{PML}}}}} \right)^q} - f_0^p\frac{i}{{{a_0}}}{\left( {\frac{{s - {s_0}}}{{{L_{PML}}}}} \right)^q} $$ (9) 式中, a0为无量纲频率常数, a0=ωLPML/c; LPML为PML区间长度; ω为角频率; c为弹性波波速.
为验证完美匹配层对流体中弹性波衰减的有效性, 以孔隙流体影响的一维压缩波的衰减为例, 假定压缩波传播中位移函数为eik(k为压缩波波数), 计算从20 m处开始坐标拉伸处理时, 其位移幅值随传播距离的变化, 并与无坐标拉伸时计算结果进行对比, 如图1所示. 可以发现, 经坐标拉伸变化后, 压缩波迅速衰减, 且对无拉伸处理区域内弹性波传播无影响. 由此可见, 合理选取坐标拉伸参数构建完美匹配层, 同样可实现气相及液相中弹性波的无反射衰减.
![]()
图 1 一维压缩波位移衰减对比图Figure 1. The comparison diagram of displacement attenuation of one-dimensional compression wave将式(9)代入式(7)的推导过程中, 可得完美匹配层有限元表达式为
$$ \left. \begin{array}{l} ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over K} }_{11}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over K} }_{12}}){{\tilde u}_b} + ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over L} }_{11}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over L} }_{12}}){{\tilde p}_l} + ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over G} }_{11}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over G} }_{12}}){{\tilde p}_g} = \tilde f\\ ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over K} }_{21}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over K} }_{22}}){{\tilde u}_b} + ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over L} }_{21}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over L} }_{22}}){{\tilde p}_l} + ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over G} }_{21}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over G} }_{22}}){{\tilde p}_g} = \dfrac{{ - {{\tilde c}_l}}}{{{\rm{i}}\omega }}\\ ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over K} }_{31}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over K} }_{32}}){{\tilde u}_b} + ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over L} }_{31}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over L} }_{32}}){{\tilde p}_l} + ({{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over G} }_{31}} + {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over G} }_{32}}){{\tilde p}_g} = \dfrac{{ - {{\tilde c}_l}}}{{{\rm{i}}\omega }} \end{array} \right\} $$ (10) 附录C中未拉伸的各变量表达式仅需将
$ \dfrac{\partial }{{\partial x}} $ 改为$ \dfrac{1}{{{\lambda _x}(x)}}\dfrac{\partial }{{\partial x}} $ ,$ \dfrac{\partial }{{\partial y}} $ 改为$ \dfrac{1}{{{\lambda _y}(y)}}\dfrac{\partial }{{\partial y}} $ , dx改为$ {\lambda _x}(x){\rm{d}}x $ , dy改为$ {\lambda _y}(y){\rm{d}}y $ , 即可得到拉伸后的各变量表达式.1.3 多耦合周期性有限元模型
将结构根据纵向性质的差异划分为若干个周期性结构, 如图2中的周期性结构1, 2, j等. 同时采用半无限周期性结构模拟纵向边界. 每个周期性结构由若干相同的单元组成. 为便于模型耦合求解, 将载荷作用结构j划分为有载荷作用区域(jB), 与无载荷作用区域(jA, jC).
整理式(7), 可得控制方程的矩阵表达式为
$$ \tilde L = K\tilde Z $$ (11) 式中
$$ \begin{split} & \tilde{{\boldsymbol{Z}}}={\left[\begin{array}{ccc}{\tilde{u}}_{b}& {\tilde{p}}_{l}& {\tilde{p}}_{g}\end{array}\right]}^{{\rm{T}}},\;\; \tilde{{\boldsymbol{L}}}={\left[\begin{array}{ccc}\tilde{f}& {\tilde{c}}_{l}& {\tilde{c}}_{g}\end{array}\right]}^{{\rm{T}}} \\ & {\boldsymbol{K}}=\left[\begin{array}{ccc}{K}_{11}+{K}_{12}& {L}_{11}+{L}_{12}& {G}_{11}+{G}_{12}\\ {K}_{21}+{K}_{22}& {L}_{21}+{L}_{22}& {G}_{21}+{G}_{22}\\ {K}_{31}+{K}_{32}& {L}_{31}+{L}_{32}& {G}_{31}+{G}_{32}\end{array}\right] \end{split} $$ 对于施加外力的周期性区域(jB), 例如图2中单元jB1, jB2, 其控制方程矩阵表达式可写为
$$ \left[ {\begin{array}{*{20}{c}} {{K_{LL}}}&{{K_{LI}}}&{{K_{LR}}} \\ {{K_{IL}}}&{{K_{II}}}&{{K_{IR}}} \\ {{K_{RL}}}&{{K_{RI}}}&{{K_{RR}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde Z}_L}} \\ {{{\tilde Z}_I}} \\ {{{\tilde Z}_R}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\tilde L}_L}} \\ {{{\tilde L}_I}} \\ {{{\tilde L}_R}} \end{array}} \right] $$ (12) 没有施加外力的周期性单元, 例如图2中单元jAi, 其单元控制方程矩阵表达式可写为
$$ \left[ {\begin{array}{*{20}{c}} {{K_{LL}}}&{{K_{LR}}} \\ {{K_{RL}}}&{{K_{RR}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde Z}_{Lk}}} \\ {{{\tilde Z}_{Rk}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\tilde L}_{Lk}}} \\ {{{\tilde L}_{Rk}}} \end{array}} \right] $$ (13) 根据自由波传播理论[30], 当弹性波经过一个周期性单元后(图3), 其前后的位移、应力、孔压等服从指数函数的传递规律, 即
$$ {\tilde{Z}}_{L{j}_{Ai + 1}}={{\rm{e}}}^{-\zeta }{\tilde{Z}}_{L{j}_{Ai}}\text{, }\;\;{\tilde{L}}_{L{j}_{Ai + 1}}={{\rm{e}}}^{-\zeta }{\tilde{L}}_{L{j}_{Ai}} $$ (14) 式中, ς为自由波传播常数.
由图3可知, 相邻单元间存在如下连续性条件
$$ {\tilde{Z}}_{R{j}_{Ai}}={\tilde{Z}}_{L{j}_{Ai + 1}}\text{, }\;\;-{\tilde{L}}_{R{j}_{Ai}}={\tilde{L}}_{L{j}_{Ai + 1}} $$ (15) 将式(14)和式(15)代入式(13), 可得
$$ \left[ {\begin{array}{*{20}{c}} {{K_{RL}}}&{{K_{RR}}} \\ 0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde Z}_{L{j_{Ai}}}}} \\ {{{\tilde Z}_{R{j_{Ai}}}}} \end{array}} \right] = {{\rm{e}}^{ - \zeta }}\left[ {\begin{array}{*{20}{c}} { - {K_{LL}}}&{ - {K_{LR}}} \\ 1&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde Z}_{L{j_{Ai}}}}} \\ {{{\tilde Z}_{R{j_{Ai}}}}} \end{array}} \right] $$ (16) 进一步整理可得
$$ \left[ {{K_{RL}} + {{\rm{e}}^{ - \zeta }}({K_{RR}}{\text{ + }}{K_{LL}}) + {{\rm{e}}^{ - 2\zeta }}{K_{LR}}} \right]{\tilde Z_{L{j_{Ai}}}} = 0 $$ (17) 观察式(17), 其为一元高次方程组, 对其进行求解可得2m(m为单元左面(或右面)的自由度)个ς. 由自由波传播理论可知, 当ς实数部为正数时, 其为由单元左面向单元右面传播的正向自由波. 当ς实数部为负数时, 其为由单元右面向单元左面传播的负向自由波.
求解式(17)时, 每个ς对应一个特征向量即ϕ, 将自由波传播常数ς及其特征向量ϕ代入式(13), 可得
$$ \left. \begin{array}{l} {\psi _p} = ({K_{LL}} + {{\rm{e}}^{ - {\zeta _p}}}{K_{LR}}){\phi _p} \hfill \\ {\psi _n} = ({K_{RL}}{{\rm{e}}^{{\zeta _n}}} + {K_{RR}}){\phi _n} \hfill \\ \end{array} \right\} $$ (18) 式中, ϕ, ψ分别为Z, L对应的特征向量; p, n分别代表正向传播波与负向传播波分量.
而单元的位移及孔压可表示为各自由波特征向量的线性组合, 即
$$ \left. \begin{array}{l} {\tilde{Z}}_{p}={\displaystyle \sum _{i=1}^{m}{\phi }_{pm}{\alpha }_{pm}}\text{, }\;\;{\tilde{Z}}_{n}={\displaystyle \sum _{i=1}^{m}{\phi }_{nm}{\alpha }_{nm}}\\ {\tilde{L}}_{p}={\displaystyle \sum _{i=1}^{m}{\psi }_{pm}{\alpha }_{pm}}\text{, }\;\;{\tilde{L}}_{n}={\displaystyle \sum _{i=1}^{m}{\psi }_{nm}{\alpha }_{nm}} \end{array} \right\}$$ (19) 式中, α为特征向量的线性组合系数.
根据式(19)推导各部分周期性结构刚度矩阵, 主要分为4类进行讨论, 即左边界、右边界、无载荷作用的周期性结构以及载荷作用的周期性结构.
(1)模型左边界
由于模型左边界左侧为半无限空间, 因此不存在波的反射和折射面, 其只存在负向自由波, 故
$$ {\tilde{Z}}_{n}={\phi }_{n}{\alpha }_{n}\text{, }\;\;{\tilde{L}}_{n}={\psi }_{n}{\alpha }_{n} $$ (20) 由式(20)整理可得
$$ {\tilde L_n} = {\psi _n}{\phi _n}^{ - 1}{\tilde Z_n} = {\overset{\frown} K} _L^1{\tilde Z_n} $$ (21) (2)模型右边界
同理, 由于模型右边界右侧为半无限空间, 故
$$ {\tilde Z_p} = {\phi _p}{\alpha _p},\;\;{\tilde L_p} = {\psi _p}{\alpha _p} $$ (22) 由式(22)整理可得
$$ {\tilde L_p} = {\psi _p}{\phi _p}^{ - 1}{\tilde Z_p} = {\overset{\frown} K} _R^N{\tilde Z_p} $$ (23) (3)无载荷作用的周期性结构
假定周期性结构j上无外载荷作用, 如图4所示, 根据式(19)整理可得
![]()
图 4 周期性结构自由波传播示意图Figure 4. Schematic diagram of free wave propagation in periodic structure$$ \left. \begin{array}{l} \tilde Z_{Rk}^j = \tilde Z_{Lk + 1}^j = \phi _n^j{{\rm{e}}^{(N - k)\mu _n^j}}\alpha _n^j + \phi _p^j{{\rm{e}}^{ - k\mu _p^j}}\alpha _p^j \hfill \\ \tilde L_{Rk}^j = - \tilde L_{Lk + 1}^j = - \psi _n^j{{\rm{e}}^{(N - k)\mu _n^j}}\alpha _n^j + \psi _p^j{{\rm{e}}^{ - k\mu _p^j}}\alpha _p^j \hfill \\ \end{array} \right\} $$ (24) 故周期性结构j左右应力、位移及孔压关系矩阵为
$$ \left. \begin{array}{l} \left\{ {\begin{array}{*{20}{c}} {\tilde L_L^j} \\ {\tilde L_R^j} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} {\psi _p^j}&{ - \psi _n^j{{\rm{e}}^{N\mu _n^j}}} \\ { - \psi _p^j{{\rm{e}}^{ - N\mu _p^j}}}&{\psi _n^j} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} {\phi _p^j}&{\phi _n^j{{\rm{e}}^{N\mu _n^j}}} \\ {\phi _p^j{{\rm{e}}^{ - N\mu _p^j}}}&{\phi _n^j} \end{array}} \right]^{ - 1}}\\ \left\{ {\begin{array}{*{20}{c}} {\tilde Z_L^j} \\ {\tilde Z_R^j} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} {{\overset{\frown} K} _{LL}^j}&{{\overset{\frown} K} _{LR}^j} \\ {{\overset{\frown} K} _{RL}^j}&{{\overset{\frown} K} _{RR}^j} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\tilde Z_L^j} \\ {\tilde Z_R^j} \end{array}} \right\}\\[-15pt] \end{array} \right\} $$ (25) 式中
$$ \begin{split} {\overset{\frown} K} _{LL}^j = & \left[ {\psi _p^j + \psi _n^j{{\rm{e}}^{ + N\mu _n^j}}{{\left( {\phi _n^j} \right)}^{ - 1}}\phi _p^j{{\rm{e}}^{ - N\mu _p^j}}} \right] \cdot \\ & {\left[ {\phi _p^j - \phi _n^j{{\rm{e}}^{ + N\mu _n^j}}{{\left( {\phi _n^j} \right)}^{ - 1}}\phi _p^j{{\rm{e}}^{ - N\mu _p^j}}} \right]^{ - 1}} \hfill \\ {\overset{\frown} K} _{LR}^j = & - \left[ {\psi _p^j{{\left( {\phi _p^j} \right)}^{ - 1}}\phi _n^j{{\rm{e}}^{ + N\mu _n^j}} + \psi _n^j{{\rm{e}}^{ + N\mu _n^j}}} \right] \cdot \\ & {\left[ {\phi _n^j - \phi _p^j{{\rm{e}}^{ - N\mu _p^j}}{{\left( {\phi _p^j} \right)}^{ - 1}}\phi _n^j{{\rm{e}}^{ + N\mu _n^j}}} \right]^{ - 1}} \hfill \\ {\overset{\frown} K} _{RL}^j = & - \left[ {\psi _n^j{{\left( {\phi _n^j} \right)}^{ - 1}}\phi _p^j{{\rm{e}}^{ - N\mu _p^j}} + \psi _p^j{{\rm{e}}^{ - N\mu _p^j}}} \right] \cdot \\ & {\left[ {\phi _p^j - \phi _n^j{{\rm{e}}^{ + N\mu _n^j}}{{\left( {\phi _n^j} \right)}^{ - 1}}\phi _p^j{{\rm{e}}^{ - N\mu _p^j}}} \right]^{ - 1}} \hfill \\ {\overset{\frown} K} _{RR}^j = & \left[ {\psi _n^j + \psi _p^j{{\rm{e}}^{ - N\mu _p^j}}{{\left( {\phi _p^j} \right)}^{ - 1}}\phi _n^j{{\rm{e}}^{ + N\mu _n^j}}} \right] \cdot \\ & {\left[ {\phi _n^j - \phi _p^j{{\rm{e}}^{ - N\mu _p^j}}{{\left( {\phi _p^j} \right)}^{ - 1}}\phi _n^j{{\rm{e}}^{ + N\mu _n^j}}} \right]^{ - 1}} \end{split}$$ (4)载荷作用周期性结构
如图2所示, 载荷作用周期性结构可分为A, B, C, 3个区域, 其中A, C两区参考无载荷作用的周期性结构处理方式即可, 对于载荷作用的B区采用传统3维有限元进行分析, 即结合式(12), 式(21), 式(23)以及式(25)可得模型的总体刚度矩阵为
$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}^L}}& \cdots &0& \cdots &0\\ \vdots & \ddots &{}& {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & \vdots \\ 0& \cdots &{{{\boldsymbol{K}}^I}}& \cdots &0\\ \vdots & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} &0& \ddots & \vdots \\ 0& \cdots &0& \cdots &{{{\boldsymbol{K}}^R}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde {\boldsymbol{Z}}}^L}}\\ \vdots \\ {{{\tilde {\boldsymbol{Z}}}^I}}\\ \vdots \\ {{{\tilde {\boldsymbol{Z}}}^R}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} 0\\ \vdots \\ {{{\tilde {\boldsymbol{L}}}^I}}\\ \vdots \\ 0 \end{array}} \right]$$ (26) 式中,
$$\begin{split} & {{\boldsymbol{K}}^L} = \left[ {\begin{array}{*{20}{c}} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _L^1 + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^1}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LR}^1}&0\\ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RL}^1}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^1{\rm{ + }}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^2}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LR}^2}\\ 0&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RL}^2}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^2{\rm{ + }}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^3} \end{array}} \right] \\ & {{\boldsymbol{K}}^I} = \\ &\left[ {\begin{array}{*{20}{c}} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^{j - 1} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^{jA}}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LR}^{jA}}&0&0&0\\ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RL}^{jA}}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^{jA} + K_{LL}^{jB}}&{K_{LI}^{jB}}&{K_{LL}^{jB}}&0\\ 0&{K_{IL}^{jB}}&{K_{II}^{jB}}&{K_{IR}^{jB}}&0\\ 0&{K_{RL}^{jB}}&{K_{RI}^{jB}}&{K_{RR}^{jB} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^{jC}}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LR}^{jC}}\\ 0&0&0&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RL}^{jC}}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^{jC} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^{j + 1}} \end{array}} \right] \\ & {{\boldsymbol{K}}^R} = \left[ {\begin{array}{*{20}{c}} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^{N - 2} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^{N - 1}}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LR}^{N - 1}}&0\\ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RL}^{N - 1}}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^{N - 1} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LL}^N}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{LR}^N}\\ 0&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RL}^N}&{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _{RR}^N + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over K} _R^N} \end{array}} \right]\\ & {\tilde {\boldsymbol{Z}}^L} = {[\begin{array}{*{20}{c}} {{{\tilde Z}_1}}&{{{\tilde Z}_2}}&{{{\tilde Z}_3}} \end{array}]^{\rm{T}}}\\ & {\tilde {\boldsymbol{Z}}^I} = {\left[ {\begin{array}{*{20}{c}} {{{\tilde Z}_{j - i}}}&{\tilde Z_L^{jB}}&{\tilde Z_I^{jB}}&{\tilde Z_R^{jB}}&{{{\tilde Z}_{j + i}}} \end{array}} \right]^{\rm{T}}}\\ & {\tilde {\boldsymbol{Z}}^R} = {[\begin{array}{*{20}{c}} {{{\tilde Z}_{N - 2}}}&{{{\tilde Z}_{N - 1}}}&{{{\tilde Z}_N}} \end{array}]^{\rm{T}}}\\ & {\tilde {\boldsymbol{L}}^I} = {\left[ {\begin{array}{*{20}{c}} 0&0&{\tilde L_I^{jB}}&0&0 \end{array}} \right]^{\rm{T}}} \end{split}$$ 由Gaussian消元法, 可得
$$ \left. \begin{array}{l} {\overset{\frown} K} _L^l = {\overset{\frown} K} _{RR}^l - {\overset{\frown} K} _{RL}^l{\left( {{\overset{\frown} K} _L^{l - 1} + {\overset{\frown} K} _{LL}^l} \right)^{ - 1}}{\overset{\frown} K} _{LR}^l \hfill \\ {\overset{\frown} K} _R^l = {\overset{\frown} K} _{LL}^l - {\overset{\frown} K} _{LR}^l{\left( {{\overset{\frown} K} _R^{l + 1} + {\overset{\frown} K} _{RR}^l} \right)^{ - 1}}{\overset{\frown} K} _{RL}^l \hfill \\ \end{array} \right\} $$ (27) 将式(27)代入式(26), 化简得
$$ \left[ {\begin{array}{*{20}{c}} {{\overset{\frown} K} _L^{jA} + K_{LL}^{jB}}&{K_{LI}^{jB}}&{K_{LL}^{jB}} \\ {K_{IL}^{jB}}&{K_{II}^{jB}}&{K_{IR}^{jB}} \\ {K_{RL}^{jB}}&{K_{RI}^{jB}}&{K_{RR}^{jB} + {\overset{\frown} K} _R^{jC}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\tilde Z_L^{jB}} \\ {\tilde Z_I^{jB}} \\ {\tilde Z_R^{jB}} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} 0 \\ {\tilde L_I^{jB}} \\ 0 \end{array}} \right] $$ (28) 由式(28)可求得载荷作用的相邻单元位移、孔压值. 将ZLjB, ZRjB代入式(26)可解得
$ {\tilde{Z}}_{j-1},\;{\tilde{Z}}_{j + 1} $ 为$$ \left. \begin{array}{l} {{\tilde Z}_{j - 1}} = - {\left({\overset{\frown} K} _{LL}^{jA} + {\overset{\frown} K} _L^{j - 1}\right)^{ - 1}}{\overset{\frown} K} _{LR}^{jA}\tilde Z_L^{jB} \hfill \\ {{\tilde Z}_{j}} = - {\left({\overset{\frown} K} _{RR}^{jC} + {\overset{\frown} K} _R^{j + 1}\right)^{ - 1}}{\overset{\frown} K} _{RL}^{jC}\tilde Z_R^{jC} \hfill \\ \end{array} \right\} $$ (29) 同理, 可解得任意周期性结构左右界面位移、孔压值, 将之代入式(24)可解的各周期性结构模态坐标值, 最终由式(24)可解得频域内周期性结构内部单元的位移、孔压值. 在此基础上进行Fourier逆变化, 可得时域中的解.
当载荷为移动载荷时
$$ f = \delta (z - {v_0}t){{\rm{e}}^{{\rm{i}}{\omega _0}t}} $$ (30) 式中, v0为载荷移动速度, ω0为载荷激振角频率
通过离散的Fourier变化即可实现移动载荷作用下的模型求解.
2. 模型的验证
为了验证模型的可靠性, 首先采用将本文模型中非饱和土参数退化的方式, 来计算单相弹性地基土-隧道系统的动力响应, 并与文献[31]中单相弹性地基土-隧道2.5 D FE-PML模型的计算结果进行对比. 2.5D FE-PML模型长为14 m, 宽为17 m, PML边界厚度h = 1 m. 多耦合周期性有限元模型横断面尺寸与2.5 D FE-PML模型相同, 纵向一个周期单元长度为0.2 m, 模型简化示意如图5所示. 模型计算参数及参数退化方式参考文献[12].
图6给出了隧道仰拱处作用单位简谐载荷(f = 30 Hz, v0 = 0 m/s)时, 不同埋深处竖向动位移ux沿水平向变化曲线. 由图6可得, 两者计算结果吻合较好, 一定程度上验证了本文提出的多耦合周期性有限元方法的可靠性.
![]()
图 6 多耦合周期性有限元方法与2.5 D FE-PML模型计算结果竖向动位移ux对比图Figure 6. Comparison of vertical dynamic displacement ux of multi coupling periodic finite element model and 2.5 D FE-PML model为了进一步验证本文方法中对于非饱和土考虑的正确性, 将本文方法的计算结果与既有的解析算法[28]的计算结果进行对比, 计算参数参考表1.
表 1 验证模型计算参数Table 1. Calculation parameters of validation modelMedium Parameter Value tunnel E/Pa 5.0 × 1010 υ 0.3 ρt/(kg·m−3) 2500 r1/m 3 h/m 0.25 d/m 5 βt 0.03 air Kg/kPa 145 ρg/(kg·m−3) 1.29 ηg/(Pa∙s) 1.5 × 10−5 water Kl/MPa 2 ρl/(kg·m−3) 1000 ηl/(Pa∙s) 1.0 × 10−3 soil grain Ks/GPa 36 ρs/(kg·m−3) 2650 soil framework Kb/MPa 43.3 μg/GPa 20 n0 0.4 Sr 0.9 Sw0 0.05 γ 0.4 βs 0.04 κ 1.0 × 10−4 φ/(°) 27 [α1, α2, α3] [1 × 10−4, 0.5, 2] 图7给出了两个非饱和地基-隧道模型在隧道仰拱处(x = 0 m, y = −2.75 m, z = 0 m)作用单位简谐载荷(f = 10 Hz, 20 Hz, 30 Hz, 40 Hz)时, 竖向动位移沿y = −3 m的分布规律. 从图7可以看出, 两者结果较为一致, 进一步验证了本文提出的多耦合周期性有限元方法的可靠性.
![]()
图 7 多耦合周期性有限元方法与解析算法计算结果对比图Figure 7. Comparison of calculation results of multi coupling periodic finite element model and analytical algorithm需要指出的是, 本文方法计算时主要耗时为自由波特征值的求解, 以验证算例二为例, 一个周期性结构单个面的自由度为5158, 计算自由波传播常数耗时7065 s, 后续耗时为203 s. 总的来说, 计算效率比解析方法略低, 但与有限元-边界元方法相比, 该方法在计算效率上具有较大的优势.
3. 数值计算与分析
基于上述方法, 以地基、隧道、隔离桩系统为例, 讨论隔离桩对地表振动响应的影响, 在隧道仰拱处作用固定单位简谐载荷(f = 20 Hz), 如图8所示. 地基、隧道模型计算参数如表1所示, 隔离桩计算参数如表2所示.
表 2 隔离桩计算参数Table 2. The calculation parameters of isolation pileelastic modulus
E/PaPoisson's ratio
νdensity
ρ/(kg·m−3)pile diameter
Rp/mpile length
Le/m3.0 × 1010 0.3 2200 1 14 图9给出了无隔离桩、单排咬合隔离桩以及单排不咬合隔离桩(存在桩间距)情况下地表振动响应分布特征. 可以发现, 隔离桩的存在使得地表振动响应规律发生了明显变化, 这主要是由于隔离桩的存在使得地基中产生了新的波反射(折射)面. 靠近隔离桩附近, 地表动力响应有所衰减, 但距离隔离桩较远的局部位置存在振动响应放大的现象. 还可以看到, 与单排不咬合隔离桩相比, 咬合桩的隔离效果更好. 由此可见, 地基中结构的纵向多耦合周期性会影响系统动力响应预测结果, 因此为准确预测系统动力响应, 应考虑地基结构的纵向变化特性.
![]()
图 9 不同工况下地表振动响应分布Figure 9. The distribution characteristics of surface vibration response of double row pile wall4. 结论
(1)基于非饱和土运动微分方程、连续性方程以及渗流运动方程, 结合应力、渗流边界等条件, 采用 Galerkin 法, 推导了ub-pl-pg格式的固、液、气三相非饱和介质有限元表达式. 相比于ub-v-w格式有限元单节点的9个自由度, ub-pl-pg格式的有限元单节点减少至5个自由度, 大大节省了计算资源.
(2)基于自由波传播理论, 求解自由波传播常数及自由波特征向量, 基于自由波传播常数及自由波特征向量求解各周期性结构刚度矩阵, 在此基础上结合传统有限元、完美匹配层单元, 提出了多耦合周期性有限元法, 将该方法分别与2.5维有限元法及解析算法计算结果相对比, 验证了该算法的可靠性.
(3)多耦合周期性有限元法具有可考虑结构沿纵向变化特性的优点. 其与解析方法相比, 计算效率略低, 但与有限元-边界元方法相比, 该方法在计算效率上具有一定优势. 地基中结构的纵向变化特性会影响系统的动力响应, 因此, 进行动力响应预测时, 应考虑地基中结构沿纵向的变化特性.
附录
附录A
$ \;$ $$ \begin{split} & {B}_{11}=\frac{{a}_{1}\gamma -n{S}_{r}}{{K}_{s}} + \frac{n}{{K}_{l}}-{A}_{s}\frac{n}{{S}_{r}}\;\;\\ &{B}_{12}=\frac{{a}_{1}(1-\gamma )-n(1-{S}_{r})}{{K}_{s}} + {A}_{s}\frac{n}{{S}_{r}}\\ & {B}_{13}\text=1-n-\frac{{K}_{b}}{{K}_{s}}\;\;\\ &{B}_{14}\text=n,\;\;{B}_{21}={A}_{s} + \frac{({a}_{1}\gamma -n{S}_{r})(1-{S}_{r})}{n{K}_{s}}\\ & {B}_{22}\text=-{A}_{s} + \frac{[{a}_{1}(1-\gamma )-n(1-{S}_{r})](1-{S}_{r})}{n{K}_{s}} + \frac{(1-{S}_{r})}{{K}_{l}}\\ & {B}_{23}\text=\left(1-n-\frac{{K}_{b}}{{K}_{s}}\right)\frac{(1-{S}_{r})}{n}\;\;\\ &{B}_{24}\text=1-{S}_{r} \end{split} $$ 附录B
$$ \begin{split} & {c}_{11}=(1-n){\rho }_{\text{s}} + n{S}_{r}{\rho }_{l}\frac{{M}_{l}}{{M}_{l}-{\omega }^{2}{\rho }_{l}} + n(1-{S}_{r}){\rho }_{g}\frac{{M}_{g}}{{M}_{g}-{\omega }^{2}{\rho }_{g}}\\ & {c}_{12}=n{S}_{r}{\rho }_{l}\frac{1}{({M}_{l}-{\omega }^{2}{\rho }_{l})}\;\; \\ & {c}_{13}=n(1-{S}_{r}){\rho }_{g}\frac{1}{({M}_{g}-{\omega }^{2}{\rho }_{g})} \\ & {c}_{21}={S}_{r}\left(1-n-\frac{{K}_{b}}{{K}_{s}}\right)\;\;\\ &{c}_{22}=n{S}_{r}{M}_{l}/({M}_{l}-{\omega }^{2}{\rho }_{l}) \\ & {c}_{23}={S}_{r}\left(\frac{a\gamma -n{S}_{r}}{{K}_{s}} + \frac{n}{{K}_{l}}-{A}_{s}\frac{n}{{S}_{r}}\right)\;\;\\ &{c}_{24}=n{S}_{r}/({M}_{l}-{\omega }^{2}{\rho }_{l})\\ & {c}_{25}={S}_{r}\left[\frac{a(1-\gamma )-n(1-{S}_{r})}{{K}_{s}} + {A}_{s}\frac{n}{{S}_{r}}\right] \;\;\\ &{c}_{26}=0\\ & {c}_{31}=\left(1-n-\frac{{K}_{b}}{{K}_{s}}\right)(1-{S}_{r}) \;\; \\ &{c}_{32}=n(1-{S}_{r}){M}_{g}/({M}_{g}-{\omega }^{2}{\rho }_{g})\\ & {c}_{33}=n\left[{A}_{s} + \frac{(a\gamma -n{S}_{r})(1-{S}_{r})}{n{K}_{s}}\right] \end{split} $$ $$ \begin{split} &{c}_{34}=0 \;\;\\ & {c}_{35}=n\Bigg[-{A}_{s} + \frac{(a(1-\gamma )-n(1-{S}_{r}))(1-{S}_{r})}{n{K}_{s}} + \frac{(1-{S}_{r})}{{K}_{l}}\Bigg] \\ & {c}_{36}=n(1-{S}_{r})/({M}_{g}-{\omega }^{2}{\rho }_{g}) \end{split}$$ 附录C
$$ \begin{array}{l} {K_{11}} = - {({\boldsymbol{BN}})^{\rm{T}}}{\boldsymbol{DBN}}\left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\ {K_{21}} = {c_{21}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{{\boldsymbol{m}}^{\rm{T}}}{\boldsymbol{BN}}\left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {K_{31}} = {c_{31}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{{\boldsymbol{m}}^{\rm{T}}}{\boldsymbol{BN}}\left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\{K_{12}} = {\omega ^2}{c_{11}}{{\boldsymbol{N}}^{\rm{T}}}{\boldsymbol{N}}\left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {K_{22}} = - {c_{22}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{{\overset{\frown} {\boldsymbol{B}}}^{\rm{T}}}{\boldsymbol{N}}\left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\{K_{32}} = - {c_{32}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{{\overset{\frown} {\boldsymbol{B}}}^{\rm{T}}}{\boldsymbol{N}}\left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {L_{11}} = {a_1}\gamma {({\boldsymbol{BN}})^{\rm{T}}}{\boldsymbol{m}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\{L_{21}} = {c_{23}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {L_{31}} = {c_{33}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\{{\boldsymbol{L}}_{12}} = - {\omega ^2}{c_{12}}{{\boldsymbol{N}}^{\rm{T}}}{{\overset{\frown} {\boldsymbol{B}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {L_{22}}{\rm{ = }}{c_{24}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{{\overset{\frown} {\boldsymbol{B}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{B}}} {\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\ {L_{32}} = 0\\ {G_{11}} = {a_1}(1 - \gamma ){({\boldsymbol{BN}})^{\rm{T}}}m{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\{G_{21}} = {c_{25}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {G_{31}} = {c_{35}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \;\;\\{G_{12}} = - {\omega ^2}{c_{13}}{{\boldsymbol{N}}^{\rm{T}}}{{\overset{\frown} {\boldsymbol{B}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \\ {G_{22}} = 0\;\;\\{G_{32}} = {c_{36}}{{\overset{\frown} {\boldsymbol{N}}}^{\rm{T}}}{{\overset{\frown} {\boldsymbol{B}}}^{\rm{T}}}{\overset{\frown} {\boldsymbol{B}}} {\overset{\frown} {\boldsymbol{N}}} \left| J \right|{\rm{d}}\eta {\rm{d}}\zeta {\rm{d}}\xi \end{array} $$ -
![]()
图 1 一维压缩波位移衰减对比图
Figure 1. The comparison diagram of displacement attenuation of one-dimensional compression wave
![]()
图 4 周期性结构自由波传播示意图
Figure 4. Schematic diagram of free wave propagation in periodic structure
![]()
图 6 多耦合周期性有限元方法与2.5 D FE-PML模型计算结果竖向动位移ux对比图
Figure 6. Comparison of vertical dynamic displacement ux of multi coupling periodic finite element model and 2.5 D FE-PML model
![]()
图 7 多耦合周期性有限元方法与解析算法计算结果对比图
Figure 7. Comparison of calculation results of multi coupling periodic finite element model and analytical algorithm
![]()
图 9 不同工况下地表振动响应分布
Figure 9. The distribution characteristics of surface vibration response of double row pile wall
表 1 验证模型计算参数
Table 1 Calculation parameters of validation model
Medium Parameter Value tunnel E/Pa 5.0 × 1010 υ 0.3 ρt/(kg·m−3) 2500 r1/m 3 h/m 0.25 d/m 5 βt 0.03 air Kg/kPa 145 ρg/(kg·m−3) 1.29 ηg/(Pa∙s) 1.5 × 10−5 water Kl/MPa 2 ρl/(kg·m−3) 1000 ηl/(Pa∙s) 1.0 × 10−3 soil grain Ks/GPa 36 ρs/(kg·m−3) 2650 soil framework Kb/MPa 43.3 μg/GPa 20 n0 0.4 Sr 0.9 Sw0 0.05 γ 0.4 βs 0.04 κ 1.0 × 10−4 φ/(°) 27 [α1, α2, α3] [1 × 10−4, 0.5, 2] 
下载: 导出CSV
表 2 隔离桩计算参数
Table 2 The calculation parameters of isolation pile
elastic modulus
E/PaPoisson's ratio
νdensity
ρ/(kg·m−3)pile diameter
Rp/mpile length
Le/m3.0 × 1010 0.3 2200 1 14
下载: 导出CSV
-
[1] 翟婉明, 赵春发. 现代轨道交通工程科技前沿与挑战. 西南交通大学学报, 2016, 51(2): 209-226 (Zhai Wanming, Zhao Chunfa. Frontiers and challenges of sciences and technologies in modern railway engineering. Journal of Southwest Jiaotong University, 2016, 51(2): 209-226 (in Chinese) doi: 10.3969/j.issn.0258-2724.2016.02.001 [2] 雷晓燕, 王全金, 圣小珍. 城市轨道交通环境振动与振动噪声研究. 铁道学报, 2003, 25(5): 109-113 (Lei Xiaoyan, Wang Quanjin, Sheng Xiaozhen. Study on environment vibration and vibration noises induced by the urban rail transit system. Journal of The China Railway Society, 2003, 25(5): 109-113 (in Chinese) doi: 10.3321/j.issn:1001-8360.2003.05.021 [3] Di HG, Zhou SH, Yao XP, et al. In situ grouting tests for differential settlement treatment of a cut-and-cover metro tunnel in soft soils. Bulletin of Engineering Geology and the Environment, 2021, 80: 6415-6427
[4] Sheng X. A review on modelling ground vibrations generated by underground trains. International Journal of Rail Transportation, 2019, 7(4): 241-261
[5] Metrikine AV, Vrouwenvelder ACWM. Surface ground vibration due to a moving train in a tunnel: two-dimensional model. Journal of Sound and Vibration, 2000, 234(1): 43-66
[6] Koziol P, Mares C, Esat I. Wavelet approach to vibratory analysis of surface due to a load moving in the layer. International Journal of Solids and Structures, 2008, 45(7-8): 2140-2159 doi: 10.1016/j.ijsolstr.2007.11.008
[7] Hu A, Li Y, Deng Y, et al. Vibration of layered saturated ground with a tunnel subjected to an underground moving load. Computers and Geotechnics, 2020, 119: 103342 doi: 10.1016/j.compgeo.2019.103342
[8] Forrest JA, Hunt HEM. A three-dimensional model for calculation of train-induced ground vibration. Journal of Sound and Vibration, 2006, 294(4/5): 678-705
[9] Hussein M, Franois S, Schevenels M, et al. The fictitious force method for efficient calculation of vibration from a tunnel embedded in a multi-layered half-space- ScienceDirect. Journal of Sound & Vibration, 2014, 333(25): 6996-7018
[10] Di H, Zhou S, He C, et al. Three-dimensional multilayer cylindrical tunnel model for calculating train-induced dynamic stress in saturated soils. Computers & Geotechnics, 2016, 80(dec.): 333-345
[11] He C, Zhou S, Guo P, et al. Analytical model for vibration prediction of two parallel tunnels in a full-space. Journal of Sound and Vibration, 2018, 423: 306-321 doi: 10.1016/j.jsv.2018.02.050
[12] 郭慧吉, 狄宏规, 周顺华等. 非饱和土−隧道系统动力响应计算的波函数法. 力学学报, 2020, 52(2): 591-602 (Guo Huiji, Di Honggui, Zhou Shunhua, et al. Wave functions method for calculating the dynamic response of a tunnel in unsaturated soil. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(2): 591-602 (in Chinese) doi: 10.6052/0459-1879-19-280 [13] Balendra T, Chua KH, Lo KW, et al. Steady-state vibration of subway-soil-building system. Journal of Engineering Mechanics, 1989, 115(1): 145-162 doi: 10.1061/(ASCE)0733-9399(1989)115:1(145)
[14] 宫全美, 徐勇, 周顺华. 地铁运行荷载引起的隧道地基土动力响应分析. 中国铁道科学, 2005, 26(5): 47-50 (Gong Quanmei, Xu Yong, Zhou Shunhua. Dynamic response analysis of tunnel foundation by vehicle vibration in metro. China Railway Science, 2005, 26(5): 47-50 (in Chinese) doi: 10.3321/j.issn:1001-4632.2005.05.010 [15] 刘维宁, 夏禾, 郭文军. 地铁列车振动的环境响应. 岩石力学与工程学报, 1996, 15(s1): 586-593 (Liu Weining, Xia He, Guo Wenjun. Study of vibration effects of underground trains on surrounding environments. Chinese Journal of Rock Mechanics and Engineering, 1996, 15(s1): 586-593 (in Chinese) [16] Gardien W, Stuit HG. Modelling of soil vibrations from railway tunnels. Journal of Sound and Vibration, 2003, 267(3): 605-619 doi: 10.1016/S0022-460X(03)00727-2
[17] 和振兴, 翟婉明, 罗震. 地铁列车引起的地面振动. 西南交通大学学报, 2008, 43(2): 218-221, 247 (He Zhenxing, Zhai Wanming, Luo Zhen. Ground vibration caused by moving metro trains. Journal of Southwest Jiaotong University, 2008, 43(2): 218-221, 247 (in Chinese) doi: 10.3969/j.issn.0258-2724.2008.02.013 [18] Hwang RN, Lysmer J. Response of buried structures to traveling wave. Journal of the Geotechnical Engineering Division, 1981, 107(2): 183-200 doi: 10.1061/AJGEB6.0001096
[19] 谢伟平, 孙洪刚. 地铁运行时引起的土的波动分析. 岩石力学与工程学报, 2003, 22(7): 1180-1184 (Xie Weiping, Sun Honggang. FEM analysis on wave propagation in soils induced by high speed train loads. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(7): 1180-1184 (in Chinese) doi: 10.3321/j.issn:1000-6915.2003.07.025 [20] Degrande G, Clouteau D, Othman R, et al. A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite element-boundary element formulation. Journal of Sound and Vibration, 2006, 293(3-5): 645-666 doi: 10.1016/j.jsv.2005.12.023
[21] 刘卫丰, 刘维宁. 地铁振动预测的周期性有限元-边界元耦合模型. 振动工程学报, 2009, 22(5): 480-485 (Liu Weifeng, Liu Weining. A coupled periodic finite element-boundary element model for prediction of vibrations induced by metro traffic. Journal of Vibration Engineering, 2009, 22(5): 480-485 (in Chinese) doi: 10.3969/j.issn.1004-4523.2009.05.007 [22] Sheng X, Jones CJC, Thompson DJ. Prediction of ground vibration from trains using the wavenumber finite and boundary element methods. Journal of Sound and Vibration, 2006, 293(3): 575-586
[23] He C, Zhou S, Guo P, et al. Modelling of ground vibration from tunnels in a poroelastic half-space using a 2.5-D FE-BE formulation. Tunnelling and Underground Space Technology, 2018, 82: 211-221
[24] 姜忻良, 徐余. 地下隧道—土体系地震反应分析的有限元与无限元耦合法. 地震工程与工程振动, 1999, 19(3): 22-26 (Jiang Xinliang, Xu Yu. Finite element and infinite element coupling method for seismic analysis of soil-underground tunnel system. Earthquake Engineering and Engineering Vibration, 1999, 19(3): 22-26 (in Chinese) doi: 10.3969/j.issn.1000-1301.1999.03.004 [25] Yang YB, Hung HH. A 2.5D finite/infinite element approach for modelling visco-elastic bodies subjected to moving loads. International Journal for Numerical Methods in Engineering, 2001, 51(11): 1317-1336
[26] 雷晓燕, 徐斌, 徐满清. 移动荷载作用下半无限弹性空间中地铁隧道动力响应的频域-波数域比例边界有限元法分析. 中国铁道科学, 2017, 38(1): 77-85 (Lei Xiaoyan, Xu Bin, Xu Manqing. Analysis on the dynamic response of metro tunnel under moving load in semi-infinite elastic space by scaled boundary finite element method in frequency-wave domain. China Railway Science, 2017, 38(1): 77-85 (in Chinese) doi: 10.3969/j.issn.1001-4632.2017.01.11 [27] Bian X, Hu J, Thompson D, et al. Pore pressure generation in a poro-elastic soil under moving train loads. Soil Dynamics and Earthquake Engineering, 2019, 125(Oct.): 105711.1-105711.15
[28] Di, HG, Zhou, SH, Guo, HG, et al. Three-dimensional analytical model for vibrations from a tunnel embedded in an unsaturated half-space. Acta Mechanica, 2021, 232: 1543-1562
[29] Lopes P, Costa PA, Calçada R, et al. Influence of soil stiffness on building vibrations due to railway traffic in tunnels: Numerical study. Computers and Geotechnics, 2014, 61(3): 277-291
[30] Germonpré M, Degrande G , Lombaert G. A track model for railway-induced ground vibration resulting from a transition zone. Journal of Rail & Rapid Transit, 2018, 232(6): 1703-1717
[31] Lopes P, Ruiz JF, Costa PA, et al. Vibrations inside buildings due to subway railway traffic. Experimental validation of a comprehensive prediction model. Science of the Total Environment, 2016, 568: 1333-1343
-
期刊类型引用(10)
1. 吴正人，王琦峰，陆世佳，董帅，刘梅. 倾斜壁面上液膜流中孤立波及其内部涡的演化特性研究. 力学学报. 2024(07): 1983-1991 . 
本站查看
2. 何子豪，孙宏月，丁伟业，赵西增，李怡良. 海啸波作用下浮式水平板防波堤消浪性能试验研究. 水运工程. 2024(12): 11-19 . 百度学术
3. 王千，陆昊成，牟泽宇，尤天庆，刘桦. 气液两相自由表面流动的三维形态重构研究. 力学与实践. 2024(06): 1132-1144 . 百度学术
4. 林金波，毛鸿飞，田正林，纪然. SPH方法的孤立波与部分淹没结构物相互作用数值计算研究. 海洋学报. 2022(06): 116-127 . 百度学术
5. 刘玉娇，余明辉，田浩永. 带滩槽地形的连续弯道中纵向流速横向分布解析. 力学学报. 2021(02): 580-588 . 本站查看
6. 赵希宁，杨晓东，张伟. 含电学边界的压电层合梁的非线性弯曲波. 力学学报. 2021(04): 1124-1137 . 本站查看
7. 杨凯恩，王千，刘桦. 孤立波对淹没平板作用力的经验预报模型. 水动力学研究与进展(A辑). 2021(02): 288-294 . 百度学术
8. 高俊亮，张一兆，何志伟，王岗. 孤立波作用下水平板的水动力特性数值研究. 哈尔滨工程大学学报. 2021(09): 1287-1294 . 百度学术
9. 邓斌，王孟飞，黄宗伟，伍志元，蒋昌波. 波浪作用下直立结构物附近强湍动掺气流体运动的数值模拟. 力学学报. 2020(02): 408-419 . 本站查看
10. 邵奇，王千，房詠柳，刘桦. 孤立波中淹没圆盘的波浪力实验研究. 水动力学研究与进展(A辑). 2020(06): 681-687 . 百度学术
其他类型引用(7)
计量
- 文章访问数: 832
- HTML全文浏览量: 376
- PDF下载量: 122
- 被引次数: 17
