-
随着经济的发展，近海工程结构的发展也越来越受到关注，如海上风电、人工岛、跨海桥梁等[1-3]。地震作用下水-结构的相互作用会在结构的表面产生动水压力作用，该动水压力会影响结构的动力响应和动力特性[4-6]。因此研究水-柱体结构相互作用对近海工程结构的抗震设计具有重要的意义。
地震作用下水-结构的相互作用的研究可以追溯到1933年Westergaard对竖直坝面上的地震动水压力的求解。地震作用下水-柱体相互作用的研究则始于1965年椭圆截面刚性柱体上动水压力的推导[7]。之后，Chopra在讨论了水体压缩性和表面波对圆柱结构地震响应的影响[8]。水-柱体结构相互作用的分析方法主要有解析法、数值法、和物理实验法[9]等。解析法是通过首先求解水-柱体结构相互作用的定解问题的动水压力解析式，进而将动水压力公式代入结构的运动方程对结构进行动力响应分析[10-15]。解析法一般只能分析几何形状简单的等截面柱体结构，如圆柱[10-13]和椭圆柱[14-15]。数值法则是基于有限元[16-17]、边界元[18-19]、无限单元[20]等数值方法分析水-柱体结构相互作用问题的方法，数值方法对结构的几何形状没有那么高的要求，可以求解几何形状复杂的结构或者倾斜结构。
轴对称柱体结构是近海工程中常见的结构形式，如圆锥柱体、复合桶形基础等[17, 21-22]。数值方法可以用于求解地震作用下水-轴对称柱体的相互作用问题。近海工程中一般将水体介质视为无限域，数值模型中可通过引入人工边界条件将水体分为近场有限域和远场无限域。近场有限域可以用动力有限元法进行模拟，远场无限域则用人工边界条件模拟[23-26]。王丕光等[24]提出了一种三维圆柱形高精度时域人工边界条件求解三维复杂结构与水体的动力相互作用问题。本文则是针对地震作用下水-轴对称柱体的相互作用问题的特点，采用分离变量法将三维问题转变为一种环向解析、竖向和径向数值的二维分析模型；二维分析模型中，基于比例边界有限元法(SBFEM)[27-29]推导了一种模拟远场无限域的高精度人工边界。比例边界有限元法能够降低问题的维度，从而减小计算量、提高计算效率，并且能自动满足无限远处的辐射边界条件。
-
地震作用下水-轴对称柱体的动力相互作用问题如图1所示，图中h表示水深，H为柱体高度；在直角坐标系和柱坐标系中，z轴为沿着柱体轴线向上，坐标原点位于柱体底部；a(z)表示z处截面半径。假定地基为刚性，地面运动沿x方向运动，相应的位移时程为ug；结构-水体系统初始静止；水体假定为小扰动不可压缩流体，并忽略表面重力波的影响，水体的密度为
${\rho _{\rm{w}}}$ =1000 kg/m3。柱坐标系$(r,\;\theta ,\;{\textit{z}} )$ 下水体的控制方程为：$$\frac{{{\partial ^2}p}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial p}}{{\partial r}} + \frac{1}{{{r^2}}}\frac{{{\partial ^2}p}}{{\partial {\theta ^2}}} + \frac{{{\partial ^2}p}}{{\partial {{\textit{z}}^2}}} = 0$$ (1) 式中，p表示水体的动水压力。
图 1 轴对称结构-水相互作用示意图
Figure 1. Axisymmetric cylinder-water interaction
水底、水面和水与结构交界面的边界条件以及无穷远处的辐射条件分别为：
$$\frac{{\partial p}}{{\partial {\textit{z}}}} = 0,\;\;{\textit{z}} = 0\qquad \qquad $$ (2) $$p = 0,\;\;{\textit{z}} = h\qquad\qquad \; $$ (3) $${\left. {\frac{{\partial p}}{{\partial r}}} \right|_{r = a}} = - {\rho _{\rm{w}}}{\ddot u_{\rm{n}}}\cos \theta \;\;\;\;$$ (4) $$p = 0,\;\;r \to \infty\qquad\quad \; $$ (5) 式中，
${u_{\rm{n}}}$ 为结构水平位移在roz平面的法向分量，字母${u_{\rm{n}}}$ 上方的点表示对时间t的偏导。动力压力p的可以通过分离变量的方式并结合边界条件表示为
$p = {p_1}(r,\;{\textit{z}},\;t)\cos \theta$ ，则相应的轴对称结构-水相互作用问题的控制方程和边界条件可以用动水压力的环向分量${p_1}$ 表示为：$$\frac{{{\partial ^2}{p_1}}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {p_1}}}{{\partial r}} - \frac{1}{{{r^2}}}{p_1} + \frac{{{\partial ^2}{p_1}}}{{\partial {{\textit{z}}^2}}} = 0$$ (6) $$\frac{{\partial {p_1}}}{{\partial {\textit{z}}}} = 0,\;\;{\textit{z}} = 0 \qquad\qquad \qquad\quad\;\;$$ (7) $${p_1} = 0,\;\;{\textit{z}} = h \qquad \qquad\qquad\qquad$$ (8) $${\left. {\frac{{\partial p_1}}{{\partial r}}} \right|_{r = a}} = - {\rho _{\rm{w}}}{\ddot u_{\rm{n}}}\;\;\qquad\qquad \qquad\;\;$$ (9) 为方便叙述在下文的推导过程中
${p_1}$ 称为动水压力。在
$r = {R_0}$ 处引入人工边界条件，则无限域水体可以分为近场有限域${\Omega _{\rm{I}}}$ 和远场无限域${\Omega _{\rm{II}}}$ 。近场有限域可以用轴对称有限元离散。 -
对于域内任意单元的控制方程，利用加权余量法可得关于动水压力
${p_1}$ 的积分方程的弱形式为：$$\int {w\left( {\frac{{{\partial ^2}{p_1}}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial {p_1}}}{{\partial r}} - \frac{1}{{{r^2}}}{p_1} + \frac{{{\partial ^2}{p_1}}}{{\partial {{\textit{z}}^2}}}} \right){\rm{d}}{V^e}} = 0\;\;\;\;\;\;$$ (10) 式中，
$w$ 为任意权函数。对式(10)进行分部积分可得：
$$ \begin{split} & \int_{{r_1}}^{{r_2}} {\int_{{{\textit{z}}_1}}^{{{\textit{z}}_2}} {\left( {\frac{{\partial w}}{{\partial r}}\frac{{\partial {p_1}}}{{\partial r}} \!- \!\frac{1}{r}w\frac{{\partial {p_1}}}{{\partial r}}\! +\! \frac{1}{{{r^2}}}w{p_1}\! +\! \frac{{\partial w}}{{\partial {\textit{z}}}}\frac{{\partial {p_1}}}{{\partial {\textit{z}}}}} \right){\rm{d}}{\textit{z}}{\rm{d}}r} }\! = \\&\qquad \int {w\frac{{\partial {p_1}}}{{\partial n}}} {\rm{d}}{S^e}\\[-16pt] \end{split}$$ (11) 对动水压力和域内任意点的坐标选用相同的插值函数：
$${p_1} = {{N}}{{{p}}_1}\qquad\qquad$$ (12) $$r = {{NR}}\;\;\qquad\qquad$$ (13) $${\textit{z}} = {{NZ}}\;\;\qquad\qquad$$ (14) 式中：
${{N}}$ 为插值函数；${{{p}}_1}$ 、${{R}}$ 、${{Z}}$ 分别为单元节点处的动水压力和对应的r、z坐标的向量。对权函数也采用同样的插值：
$$w = {{Nw}}\qquad\qquad$$ (15) 将式(12)~式(15)代入式(11)，消去任意权函数
${{w}}$ ，整理得：$$( {{{K}} + {{{K}}^1}} ){{{p}}_1} = {{F}}\;\;\;\;$$ (16) 其中：
$$\tag{17a}{{K}} = \iint {\left( {\frac{{\partial {{{N}}^T}}}{{\partial r}}\frac{{\partial {{N}}}}{{\partial r}} + \frac{{\partial {{{N}}^T}}}{{\partial {\textit{z}}}}\frac{{\partial {{N}}}}{{\partial {\textit{z}}}}} \right){\rm{d}}{\textit{z}}{\rm{d}}r}\;\;$$ $$\tag{17b}{{{K}}^1} = \iint {\left( { - \frac{1}{r}{{{N}}^T}\frac{{\partial {{N}}}}{{\partial r}} + \frac{1}{{{r^2}}}{{{N}}^T}{{N}}} \right){\rm{d}}{\textit{z}}{\rm{d}}r}$$ $$\tag{17c}{{F}} = \int {{{{N}}^T}\frac{{\partial {p_1}}}{{\partial n}}} {\rm{d}}S\;\;\qquad \qquad\qquad\qquad $$ 采用与有限元相同的四节点四边形单元等参单元变换，再利用高斯积分即可得到刚度矩阵的值。
-
基于比例边界中心线的思想，本文将oz轴选取为比例边界中心线，如图2所示，图中比例边界坐标
$\xi $ 和$\eta $ ，为局部坐标。径向坐标$\xi $ 可以看作是比例系数，因此域内任一点的坐标(r, z)可以用局部比例边界坐标$\xi $ 和$\eta $ 表示为：$$\tag{18a}r = \xi {R_0}$$ $$\tag{18b}{\textit{z}} = {{N}}{{{z}}_{\rm{b}}}$$ 式中：
${{N}} = {{N}}\left( \eta \right) = [ { {{{\left( {1 - \eta } \right)}/2}}\;\;{{{\left( {1 + \eta } \right)}/2}} } ]$ 为形函数；${{{z}}_{\rm{b}}} = {[ {{{\textit{z}}_1}}\;\;{{{\textit{z}}_2}} ]^{\rm{T}}}$ 为边上的节点坐标。
图 2 SBFE坐标变换
Figure 2. The transformation of SBFE coordinates
两种坐标间导数之间的关系为：
$$\left\{ {\begin{array}{*{20}{c}} {\dfrac{\partial }{{\partial \xi }}} \\ {\dfrac{\partial }{{\partial \eta }}} \end{array}} \right\} = {{J}}\left\{ {\begin{array}{*{20}{c}} {\dfrac{\partial }{{\partial r}}} \\ {\dfrac{\partial }{{\partial {\textit{z}}}}} \end{array}} \right\}\;\qquad\qquad\qquad$$ (19) 式中，
${{J}}$ 为Jacobian矩阵，可以表示为：$${{J}} = \left[ {\begin{array}{*{20}{c}} {{R_0}}&0 \\ 0&{{{{N}}_{,\eta }}{{{z}}_{\rm{b}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{R_0}}&0 \\ 0&{{{{\Delta _{\textit{z}}}}/2}} \end{array}} \right]$$ (20) 式中，
${\Delta _{\textit{z}}} = {{\textit{z}}_2} - {{\textit{z}}_1}$ ，可以看出Jacobian矩阵只有依赖于边界形状。微分算子可以用局部比例边界坐标表示为：
$$\left\{ {\begin{array}{*{20}{c}} {\dfrac{\partial }{{\partial r}}} \\ {\dfrac{\partial }{{\partial {\textit{z}}}}} \end{array}} \right\} = {{{J}}^{ - 1}}\left\{ {\begin{array}{*{20}{c}} {\dfrac{\partial }{{\partial \xi }}} \\ {\dfrac{\partial }{{\partial \eta }}} \end{array}} \right\} = {{{b}}_1}\dfrac{\partial }{{\partial \xi }} + {{{b}}_2}\dfrac{\partial }{{\partial \eta }}$$ (21) 式中，
${{{b}}_1} = \dfrac{1}{{\left| {{J}} \right|}}\left\{ {\begin{array}{*{20}{c}}{{\Delta _{\textit{z}}}/2}\\0\end{array}} \right\}$ ，${{{b}}_2} = \dfrac{1}{{\left| {{J}} \right|}}\left\{ {\begin{array}{*{20}{c}} 0 \\ {{R_0}} \end{array}} \right\}$ ，因此关于径向坐标的微分可以表示为：$$\frac{\partial }{{\partial r}} = {b_1}\frac{\partial }{{\partial \xi }}\quad$$ (22) 式中，
${b_1} = {{{b}}_1}\left( {1,1} \right)$ 。对于域内一个微元体单元dV可以表示为：
$${\rm{d}}V = \left| {{J}} \right|{\rm{d}}\xi {\rm{d}}\eta $$ (23) 式中，
$\left| {{J}} \right| = {{{R_0}{\Delta _{\textit{z}}}}/2}$ 为Jacobian矩阵的行列式值。 -
对控制方程式(6)运用加权余量法可得：
$$\int {w\left( {{r^2}\frac{{{\partial ^2}{p_1}}}{{\partial {r^2}}} + r\frac{{\partial {p_1}}}{{\partial r}} - {p_1} + {r^2}\frac{{{\partial ^2}{p_1}}}{{\partial {{\textit{z}}^2}}}} \right){\rm{d}}V} = 0$$ (24) 将相应的坐标变换关系式(21)、式(22)和式(23)代入式(24)，并对其中的两项进行分部积分得：
$$ \begin{split} & \int_1^\infty {\int_{ - 1}^1 {{\xi ^2}R_0^2\left( {w{{b}}_1^{\rm{T}}{{{b}}_1}\frac{{{\partial ^2}{p_1}}}{{\partial {\xi ^2}}} + w{{b}}_1^{\rm{T}}{{{b}}_2}\frac{{{\partial ^2}{p_1}}}{{\partial \xi \partial \eta }}} \right)\left| {{J}} \right|{\rm{d}}\xi {\rm{d}}\eta } } - \\[-2pt]&\;\; \int_1^\infty \!\! {\int_{ - 1}^1 {{\xi ^2}R_0^2\left( {\frac{{\partial w}}{{\partial \eta }}{{b}}_2^{\rm{T}}{{{b}}_1}\frac{{\partial {p_1}}}{{\partial \xi }} + \frac{{\partial w}}{{\partial \eta }}{{b}}_2^{\rm{T}}{{{b}}_2}\frac{{\partial {p_1}}}{{\partial \eta }}} \right)\left| {{J}} \right|{\rm{d}}\xi {\rm{d}}\eta } } + \\[-2pt]&\;\; \int_1^\infty {\int_{ - 1}^1 {w\left( {\xi {R_0}{b_1}\frac{{\partial {p_1}}}{{\partial \xi }} - {p_1}} \right)\left| {{J}} \right|{\rm{d}}\xi {\rm{d}}\eta } } + \\[-2pt]&\;\; \int_1^\infty {w{{b}}_2^{\rm{T}}\left( {{{{b}}_1}\frac{{\partial {p_1}}}{{\partial \xi }} + {{{b}}_2}\frac{{\partial {p_1}}}{{\partial \eta }}} \right)\left| {{J}} \right|{\rm{d}}\xi } = 0 \qquad\qquad\quad(25) \end{split} $$ 对动水压力和权函数采用相同的形函数插值：
$${p_1} = {{N}}{{{p}}_1}$$ (26) $$w = {{Nw}}\;\;$$ (27) 式中：
${{{p}}_1} = {{{p}}_1}\left( \xi \right)$ 为单元节点动水压力向量；${{w}} = {{w}}\left( \xi \right)$ 为权函数向量。定义${{{B}}_1} = {{{b}}_1}{{N}}$ ，${{{B}}_2} = {{{b}}_2}{{{N}}_{,\eta }}$ ，下标撇号表示导数，如${{{N}}_{,\eta }}$ 即表示形函数N对坐标$\eta $ 的导数。将式(26)、式(27)代入式(25)，并消去任意权函数
${{w}}$ 和对$\xi $ 的积分可得：$$ \begin{split} & {\xi ^2}R_0^2( {{{E}}_0^{\rm{e}}{{{p}}_{1,\xi \xi }} + {{E}}_1^{\rm{e}}{{{p}}_{1,\xi }} - {{E}}_1^{{\rm{e}}T}{{{p}}_{1,\xi }} - {{E}}_2^{\rm{e}}{{{p}}_1}} ) + \\&\qquad {{E}}_3^{\rm{e}}\xi {R_0}{{{p}}_{1,\xi }} - {{E}}_4^{\rm{e}}{{{p}}_1} + {{{F}}^{\rm{t}}} = 0 \end{split} $$ (28) 其中：
$$\tag{29a}{{E}}_0^{\rm{e}} = \int_{ - 1}^1 {R_0^2{{B}}_1^T{{{B}}_1}\left| {{J}} \right|d\eta } = \frac{{{R_0}{\Delta _{\textit{z}}}}}{6}\left[ {\begin{array}{*{20}{c}} 2&1 \\ 1&2 \end{array}} \right]\;\;$$ $$\tag{29b}{{E}}_1^{\rm{e}} = \int_{ - 1}^1 {R_0^2{{B}}_1^T{{{B}}_2}\left| {{J}} \right|{\rm{d}}\eta } = {{0}}\qquad\qquad\qquad\;\; $$ $$\tag{29c}{{E}}_2^{\rm{e}} = \int_{ - 1}^1 {R_0^2{{B}}_2^T{{{B}}_2}\left| {{J}} \right|{\rm{d}}\eta } = \frac{{2R_0^3}}{{{\Delta _{\textit{z}}}}}\left[ {\begin{array}{*{20}{c}} 1&{ - 1} \\ { - 1}&1 \end{array}} \right]$$ $$\tag{29d}{{E}}_3^{\rm{e}} = \int_{ - 1}^1 {{R_0}{b_1}{{{N}}^T}{{N}}\left| {{J}} \right|{\rm{d}}\eta } = \frac{{{\Delta _{\textit{z}}}{R_0}}}{6}\left[ {\begin{array}{*{20}{c}} 2&1 \\ 1&2 \end{array}} \right]$$ $$\tag{29e}{{E}}_4^{\rm{e}} = \int_{ - 1}^1 {{{{N}}^T}{{N}}\left| {{J}} \right|{\rm{d}}\eta } = \frac{{{R_0}{\Delta _{\textit{z}}}}}{6}\left[ {\begin{array}{*{20}{c}} 2&1 \\ 1&2 \end{array}} \right]\qquad$$ 由式(29a)、式(29d)和式(29e)可以看出
${{E}}_0^{\rm{e}} = {{E}}_3^{\rm{e}} = {{E}}_4^{\rm{e}}$ 。将式(28)沿人工边界处进行装配，消去系数为零的项，并利用上下表面的边界条件可得：
$${\xi ^2}( {{{{E}}_0}{{{P}}_{1,\xi \xi }} - {{{E}}_2}{{{P}}_1}} ) + {{{E}}_0}\xi {{{P}}_{1,\xi }} - {{{E}}_0}{{{P}}_1} = 0$$ (30) 式中：
${{{E}}_0}$ 和${{{E}}_2}$ 为式(29)的单元系数矩阵装配的矩阵；${{{P}}_1}$ 为各节点的动水压力组成的向量。人工边界处的节点力可以表示为：
$${{{Q}}^{\rm{e}}} \!=\! \int {{{{N}}^T}r\frac{{\partial {p_1}}}{{\partial r}}{\rm{d}}{S^\xi }}\! = \!{{E}}_0^{\rm{e}}{{{p}}_{1,\xi }}\! + \!{{E}}_0^{\rm{e}}{{{p}}_1} \!=\! \frac{\xi }{{{R_0}}}{{E}}_0^{\rm{e}}{{{p}}_{1,\xi }}$$ (31) 式中，
${{{Q}}^{\rm{e}}}$ 为单元节点力向量。定义无限域的动力刚度为：
$${{R}} = - \frac{1}{r}{{Q}} = - {{S}}_{\rm{B}}^\infty {{{P}}_1}$$ (32) 式中，
${{Q}}$ 为离散的边界节点上的节点力向量。联立式(31)、式(32)和式(30)可得在人工边界(ζ=1)处动力刚度的比例边界有限元方程为：
$${{\tilde S}}{{{E}}_0}^{ - 1}{{\tilde S}} - {{{E}}_0} - {{{E}}_2} = {{0}}$$ (33) 其中，
${{\tilde S}} = R_0^2{{S}}_{\rm{B}}^\infty$ 。式(33)为关于${{\tilde S}}$ 的Riccati方程，可以直接用MATLAB中的‘care’函数求解。 -
将整个水域节点分为三块区域，即水-轴对称柱体交界面节点(区域
${S_1}$ )、人工边界节点(区域${S_3}$ )和其余节点(区域${S_2}$ )，则地震作用下内域和外域耦合的有限元方程可以用分块的形式表示为：$$\left[ {\begin{array}{*{20}{c}} {{{{S}}_{11}}}&{{{{S}}_{12}}}&{{{{S}}_{13}}}\\ {{{{S}}_{21}}}&{{{{S}}_{22}}}&{{{{S}}_{23}}}\\ {{{{S}}_{31}}}&{{{{S}}_{32}}}&{{{{S}}_{33}} + {{S}}_{\rm{B}}^\infty } \end{array}} \right]\left\{ \begin{array}{l} {{{P}}_{1,1}}\\ {{{P}}_{1,2}}\\ {{{P}}_{1,3}} \end{array} \right\}{{ = }}\left\{ \begin{array}{c} {{{D}}_1}\\ {{{D}}_2}\\ {{{D}}_3} \end{array} \right\}$$ (34) 式中：内域刚度矩阵
${{S}} = ( {{{K}} + {{{K}}^1}} )$ ，其中，${{S}}$ 为除了${{S}}_{\rm{B}}^\infty $ 以外的刚度矩阵；${{{D}}_2} = \left\{ {{0}} \right\}$ ；${{{D}}_3} = \left\{ {{0}} \right\}$ ；假定柱体各节点的水平位移向量为${{{U}}_1}\left( t \right)$ ，则${{{D}}_1}$ 的可以表示为：$$ \begin{split} {{{D}}_1}{\rm{ = }}&\int_{{S_1}} {{{N}}_{S1}^{\rm{T}}\frac{{\partial {P_1}}}{{\partial n}}{\rm{d}}S} = \\& - \int_{{S_1}} {{{N}}_{S1}^{\rm{T}}{{N}}_{S1}^{}\cos \left( {n,r} \right){\rm{d}}S} {{{{\ddot U}}}_1}\left( t \right) = \\& - {\rho _{\rm{w}}}{{{M}}^1}{{{{\ddot U}}}_1}\left( t \right) \end{split} $$ (35) 式中：
${{{M}}^1} = \displaystyle\int_{{S_1}} {{{N}}_{S1}^{\rm{T}}{{N}}_{S1}^{}\cos \left( {n,r} \right){\rm{d}}S}$ ，${{N}}_{S1}^{}$ 为交界面形函数向量；n表示在平面rOz内交界面的法向。整个水域的动水压力可以表示为：
$$ \begin{split} \left\{\! \begin{array}{l} {{{P}}_{1,1}}\\ {{{P}}_{1,2}}\\ {{{P}}_{1,3}} \end{array} \right\}\!=\!&{\left[\! {\begin{array}{*{20}{c}} {{{{S}}_{11}}}&{{{{S}}_{12}}}&{{{{S}}_{13}}}\\ {{{{S}}_{21}}}&{{{{S}}_{22}}}&{{{{S}}_{23}}}\\ {{{{S}}_{31}}}&{{{{S}}_{32}}}&{{{{S}}_{33}} + {{S}}_{\rm{B}}^\infty } \end{array}}\! \right]^{ - 1}}\left\{\! \begin{array}{c} {{{D}}_1}\\ {{{D}}_2}\\ {{{D}}_3} \end{array} \!\right\}\! = \\& \left[ {\begin{array}{*{20}{c}} {{{{G}}_{11}}}&{{{{G}}_{12}}}&{{{{G}}_{13}}}\\ {{{{G}}_{21}}}&{{{{G}}_{22}}}&{{{{G}}_{23}}}\\ {{{{G}}_{31}}}&{{{{G}}_{32}}}&{{{{G}}_{33}}} \end{array}} \right]\left\{ \begin{array}{c} {{{D}}_1}\\ {{{D}}_2}\\ {{{D}}_3} \end{array} \right\} \end{split} $$ (36) 结构表面的动水压力进一步可以求解为：
$${{{P}}_{1,1}}{\rm{ = }} - {\rho _{\rm{w}}}{{{G}}_{11}}{{{M}}^1}{{{\ddot U}}_1}\left( t \right)$$ (37) -
柱体结构可以用梁单元进行模拟，单元刚度和质量矩阵可以用一下公式表示[30]：
$$E{I_i} = \frac{1}{{{l_i}}}\int_{{{\textit{z}}_{i + 1}}}^{{{\textit{z}}_i}} {EI\left( {\textit{z}} \right)} {\rm{d}}{\textit{z}}\quad$$ (38) $${\rho _{\rm{s}}}{A_i} = \frac{1}{{{l_i}}}\int_{{{\textit{z}}_{i + 1}}}^{{{\textit{z}}_i}} {{\rho _{\rm{s}}}A\left( {\textit{z}} \right)} {\rm{d}}{\textit{z}}$$ (39) 进一步，轴对称柱体的有限元方程可写成：
$${{M\ddot U}}\left( t \right) + {{C\dot U}}\left( t \right) + {{KU}}\left( t \right){\rm{ = }} - {{M}}{{{\ddot u}}_{\rm{g}}}\left( t \right) + {{F}}\left( t \right)$$ (40) 式中：
${{M}}$ 、${{K}}$ 分别由单元质量矩阵和刚度矩阵装配得到；${{C}}$ 为结构的阻尼矩阵，本文用瑞利阻尼；${{F}}\left( t \right)$ 是由动水压力引起作用在柱体上动水力向量。该动水力向量可以表示为：$$ \begin{split} {{F}}\left( t \right){\rm{ = }}& - \int_{{S_1}} {{{N}}_{S1}^{\rm{T}}{p_{S1}}\cos \left( {n,r} \right){\rm{d}}S} = \\& - \int_{{S_1}} {{{N}}_{S1}^{\rm{T}}{{N}}_{S1}^{}\cos \left( {n,r} \right){\rm{d}}S} {{{P}}_{1,1}} \end{split}\quad\;\;\; $$ (41) 式中，
${p_{S1}}$ 为水-柱体交界面上的动水压力。将式(37)代入式(41)可得：
$${{F}}\left( t \right){\rm{ = }} - \rho {{{M}}^1}{{{G}}_{11}}{{{M}}^1}{{{\ddot U}}_1}\left( t \right) = - {{{M}}_{\rm{p}}}{{{\ddot U}}_1}\left( t \right)$$ (42) 式中，
${{{\ddot U}}_1}\left( t \right) = {{\ddot U}}\left( t \right) + {{{\ddot u}}_{\rm{g}}}\left( t \right)$ 。将式(42)代入式(40)可得水-轴对称柱体结构相互作用系统的方程为：
$${{{M}}_{\rm{s}}}{{\ddot U}}\left( t \right) + {{C\dot U}}\left( t \right) + {{KU}}\left( t \right){\rm{ = }} - {{{M}}_{\rm{s}}}{{{\ddot u}}_{\rm{g}}}\left( t \right)\quad$$ (43) 式中，
${{{M}}_{\rm{s}}} = {{M}} + {{{M}}_{\rm{p}}}$ 。式(43)所示的有限元方程，可以用Newmark-β法等数值积分方法进行求解。 -
本文中算例的结构的弹性模量、密度、泊松比和阻尼比分别为：21 GPa，7800 kg/m3，0.3和0.02，结构的厚度为0.06 m。图3为本文中选取的作用在刚性地基上的地震动时程曲线，图4为对应的地震动的傅里叶幅值曲线。
图 3 地震动时程
Figure 3. Time history of earthquakes
图 4 地震动的傅里叶幅值曲线
Figure 4. Furrier amplitude spectrum of ground motions
-
首先，对轴对称结构的等效梁单元进行验证。表1为本文轴对称结构的等效梁单元和ABAQUS中的实体梁单元计算的自振频率的比较，可以看出相对误差在5%以内，因此本文选用的等效变换方法是符合精度要求的。
表 1 等效变换梁单元和实体单元的自振频率
Table 1. The natural frequency of vibration comparison of solid element and equivalent beam
阶数 实体单元/Hz 等效梁单元/Hz 相对误差/(%) 一阶 1.2871 1.306 632 1.52 二阶 6.6814 6.995 571 4.70 三阶 17.4300 17.881 06 2.59 其次，对二维分析模型中的人工边界的验证。分别选取等截面圆柱和圆锥柱体两个模型：模型1等截面圆柱的半径和高度分别为3 m、80 m，模型2圆锥形柱体的半径、倾角和高度分别为
$a = 3\;{\rm{m}}$ 、${\theta _0} = a\tan (80/1)$ 和80 m。模型1的参考解为解析解，模型2的参考解为$20\;000 \times 80$ 的扩展有限元解，网格大小为$5 \times 5$ 。人工边界的设置在$r = 50$ m处。图5为刚性轴对称柱体上的动水力分布本文解和参考解的对比。由图5可以看出，采用本文提出的人工边界条件得到的本文解与参考解极好地吻合。
图 5 刚性柱体上无量纲动水力沿高度分布
Figure 5. The distribution of dimensionless hydrodynamic pressure
最后，对本文提出的水-轴对称相互作用的子结构分析方法的验证。等截面圆柱的参考解为将柱体结构与动水压力的解析解耦合求得的结构动力反应；圆锥柱体模型的参考解为柱体结构有限元方程与动水力的扩展有限元方程与柱体结构有限元方程耦合解得的结构动力反应。在模型1、模型2施加图3所示的编号为No.1地震动荷载，将结构的顶部位移作为参考量。图6为结构顶部位移响应本文解和参考解的对比。由图6可以看出本文方解与参考解吻合较好。
图 6 柱体顶部位移响应
Figure 6. Top displacement of the axisymmetric structure
-
计算了不同水深时水-轴对称结构相互作用系统的地震响应，分析水深对水-轴对称结构相互作用。引入无量纲系数
${R_{\rm{u}}} = {{{u_{\rm{im}}}}/{{u_{0 {\rm{m}}}}}}$ ，其中${u_{\rm{im}}}$ 表示不同水深时水-轴对称结构相互作用时位移响应的最大值，${u_{0 {\rm{m}}}}$ 表示轴对称结构在空气中时位移响应的最大值。以模型2为例，分别计算了水深为20 m、40 m、60 m和80 m时，荷载分别为图3所示的7条地震动时的结构的动力响应；图7为
${R_{\rm{u}}}$ 随水深的变化。由图7可以看出，地震动水压力对结构地震反应的影响随水深相对结构高度(h/H)的增加呈增大趋势，这主要是由于作用在结构上的地震动水压力随水深的增加而增大。需要指出的是，由于不同地震动的频谱特性不同，因此不同地震动作用下动水压力对结构地震响应的影响也不相同。
图 7
${R_{\rm{u}}}$ 与水深的关系Figure 7.
${R_{\rm{u}}}$ vary with different earthquakes定义自振频率降低率为：
$$Err = \frac{{{\omega _{\rm{s}}} - {\omega _{\rm{ws}}}}}{{{\omega _{\rm{s}}}}} \times 100 \text{%}$$ (44) 式中，
${\omega _{\rm{ws}}}$ 和${\omega _{\rm{s}}}$ 分别是指水-结构相互作用系统的自振频率和结构无水时的自振频率。表2为不同水深时水-轴对称相互作用系统的自振频率。表3为不同水深时自振频率的相对值。由表2和3可以看出当水深相对结构高度较小时(h/H≤0.25)，系统的基频和结构无水时的基频相差不大，由图7也可以看出此时结构的地震响应与无水时没有明显差异；这表明水深相对结构高度较浅时，地震动水压力对柱体结构地震响应的影响可忽略。随着水深相对结构的高度的增大(0.25<h/H≤1)，系统的频率降低率逐渐减小，自振频率的减低可能会使得系统的二阶自振频率出现在地震动的卓越频率范围内，图7中也可以看出此时结构的地震响应与无水时相差较大；因此水深相对结构高度较深时，地震动水力对柱体结构的动力响应和自振特性的影响不可忽略。
表 2 水-轴对称结构相互作用系统的自振频率
Table 2. The natural frequency of vibration of water-axisymmetric interaction system
水深h/m 自振频率/(rad/s) 一阶 二阶 三阶 0 1.31 7.00 17.88 20 1.30 6.65 14.45 40 1.21 4.63 12.51 60 0.94 4.26 10.10 80 0.66 3.47 8.96 表 3 自振频率降低率
Table 3. The reduction ratio of natural vibration frequency
水深h/m 降低率/(%) 一阶 二阶 三阶 20 0.76 5.00 19.18 40 7.63 33.86 30.03 60 28.24 39.14 43.51 80 49.62 50.42 49.89 -
基于不可压缩水体，本文提出了一种地震作用下水-轴对称柱体相互作用的动力响应的子结构分析方法。首先将水体三维模型转化为环向解析、竖向和径向离散的二维模型，其次用高精度人工边界条件模拟无限域水体的作用，在保证了精度的同时大大提高了计算效率。
数值算例中分别验证了结构的等效梁单元、动水压力和水-结构动力相互作用耦合模型的精确性。以圆锥柱体为例分析了地震作用下水深对水-轴对称结构相互作用的地震响应的影响，数值算例表明随着水深的增加，动水力对结构的自振频率和动力响应的影响呈增大的趋势。
由于本文提出的子结构可以将水-轴对称结构的相互作用通过附加质量的模拟，易于在商业有限元软件中实现和方便于工程应用。
A SUBSTRUCTURE MODEL FOR WATER-AXISYMMETRIC CYLINDER INTERACTION DURING EARTHQUAKES
-
摘要: 针对水-轴对称柱体动力相互作用问题，提出了一种地震作用下水-结构相互作用的时域子结构分析方法。基于三维不可压缩水体的波动方程和边界条件，利用分离变量法将其转换为环向解析、竖向和径向数值的二维模型；基于比例边界有限元推导了截断边界处无限域水体的动力刚度方程，并将水体内域有限元方程和人工边界处的动水压力进行耦合，从而得到结构表面的动水压力方程；将轴对称柱体结构的有限元方程与动水压力方程耦合，从而得到水-轴对称柱体结构系统的时域有限元方程；数值算例验证该文提出的水-轴对称动力相互作用的子结构方法，结果表明：该文方法具有很高的精度和计算效率。通过对水中轴对称结构地震响应和自振频率的分析表明：地震动水压力对结构自振频率和动力响应的影响随水深的增加而增大。
-
关键词:
- 地震 /
- 轴对称柱体 /
- 水-轴对称结构相互作用 /
- 有限元 /
- 比例边界有限元
Abstract: A substructure analysis method for water-axisymmetric cylinder dynamic interaction problems is presented to simulate the seismic response of water-structure interactions. According to the wave equation and boundary conditions for the three-dimensional incompressible water and by applying the method of separation of variables, the model is transformed into a two-dimensional model, which is analytical in the circumferential direction and numerical in the vertical and radial directions. The dynamic stiffness equation of the infinity water at the truncated boundary is derived by using the scaled boundary finite element method, then the hydrodynamic pressure on the structure can be obtained by coupling the finite element equation of the near field water with the hydrodynamic on the truncated boundary. The finite element equation in time domain of a water-axisymmetric cylinder interaction system is developed by coupling the finite element equation of the structure with the hydrodynamic pressure. The numerical examples demonstrated that the proposed method is accurate and efficient. Numerical results show that the effect of the hydrodynamic pressure on the natural frequency and the dynamic response of the axisymmetric structure is increasing with the increase of the water depth in general. -
图 5 刚性柱体上无量纲动水力沿高度分布
Figure 5. The distribution of dimensionless hydrodynamic pressure
表 1 等效变换梁单元和实体单元的自振频率
Table 1. The natural frequency of vibration comparison of solid element and equivalent beam
阶数 实体单元/Hz 等效梁单元/Hz 相对误差/(%) 一阶 1.2871 1.306 632 1.52 二阶 6.6814 6.995 571 4.70 三阶 17.4300 17.881 06 2.59 
下载: 导出CSV
表 2 水-轴对称结构相互作用系统的自振频率
Table 2. The natural frequency of vibration of water-axisymmetric interaction system
水深h/m 自振频率/(rad/s) 一阶 二阶 三阶 0 1.31 7.00 17.88 20 1.30 6.65 14.45 40 1.21 4.63 12.51 60 0.94 4.26 10.10 80 0.66 3.47 8.96
下载: 导出CSV
表 3 自振频率降低率
Table 3. The reduction ratio of natural vibration frequency
水深h/m 降低率/(%) 一阶 二阶 三阶 20 0.76 5.00 19.18 40 7.63 33.86 30.03 60 28.24 39.14 43.51 80 49.62 50.42 49.89
下载: 导出CSV
-
[1] 项海帆. 21世纪世界桥梁工程的展望[J]. 土木工程学报, 2000, 33(3): 1 − 6. doi: 10.3321/j.issn:1000-131X.2000.03.001 Xiang Haifan. Prospect of world’s bridge projects in 21st century [J]. China Civil Engineering Journal, 2000, 33(3): 1 − 6. (in Chinese) doi: 10.3321/j.issn:1000-131X.2000.03.001 [2] Zhang D, Zhang X, He J, Chai Q. Offshore wind energy development in China: Current status and future perspective [J]. Renewable and Sustainable Energy Reviews, 2011, 15: 467 − 484. [3] 闫静茹, 路德春, 杜修力, 等. 港珠澳大桥工程人工岛三维非线性地震反应分析[J]. 世界地震工程, 2016, 32(1): 161 − 168. Yan Jingru, Lu Dechun, Du Xiuli, et al. Three-dimensional nonlinear seismic response analysis of artificial island of Hong Kong-zhuhai-Macao Bridge Project [J]. World Earthquake Engineering, 2016, 32(1): 161 − 168. (in Chinese) [4] 黄信, 李忠献. 动水压力作用对深水桥墩地震响应的影响[J]. 土木工程学报, 2011, 44(1): 65 − 73. Huang Xin, Li Zhongxian. Influence of hydrodynamic pressure on seismic responses of bridge piers in deep water [J]. China Civil Engineering Journal, 2011, 44(1): 65 − 73. (in Chinese) [5] 李忠献, 黄信. 行波效应对深水连续刚构桥地震响应的影响[J]. 工程力学, 2013, 30(3): 120 − 125. Li Zhongxian, Huang Xin. Influence of traveling wave effect on seismic responses of continuous rigid-framed bridge in deep water [J]. Engineering Mechanics, 2013, 30(3): 120 − 125. (in Chinese) [6] 张文学, 黄荐, 陈盈, 杜修力. 渡槽结构考虑流固耦合的横向地震响应简化计算公式[J]. 工程力学, 2017, 34(8): 69 − 75. Zhang Wenxue, Huang Jian, Chen Ying, Du Xiuli. A simplified formula for the calculation of the transverse seismic response of aqueducts considering fluid-structure interaction [J]. Engineering Mechanics, 2017, 34(8): 69 − 75. (in Chinese) [7] Kotsubo S. Seismic force effect on submerged bridge piers with elliptic cross sections [J]. Proceedings of the Third World Conference on Earthquake Engineering. 1965, 1965(2): 342 − 356. [8] Liaw C Y, Chopra A K. Dynamics of towers surrounded by water [J]. Earthquake Engineering & Structural Dynamics, 1974, 3(1): 33 − 49. [9] 李乔, 刘浪, 杨万理. 深水桥梁墩水耦合振动试验研究与数值计算[J]. 工程力学, 2016, 33(7): 197 − 203. Li Qiao, Liu Lang, Yang Wanli. Experimental and numerical investigation on pier-water Coupling vibration of bridges in deep water [J]. Engineering Mechanics, 2016, 33(7): 197 − 203. (in Chinese) [10] Jiang H, Wang B, Bai X, et al. Simplified expression of hydrodynamic pressure on deepwater cylindrical bridge Piers during earthquakes [J]. Journal of Bridge Engineering, 2017, 22(6): 4017014. doi: 10.1061/(ASCE)BE.1943-5592.0001032 [11] Du X, Wang P, Zhao M. Simplified formula of hydrodynamic pressure on circular bridge piers in the time domain [J]. Ocean Engineering, 2014, 85(2014): 44 − 53. [12] Li Q, Yang W. An improved method of hydrodynamic pressure calculation for circular hollow piers in deep water under earthquake [J]. Ocean Engineering, 2013, 72: 241 − 256. [13] Han R P S, Xu H. A simple and accurate added mass model for hydrodynamic fluid-structure interaction analysis [J]. Journal of the Franklin Institute, 1996, 333(6): 929 − 945. doi: 10.1016/0016-0032(96)00043-9 [14] Wang P, Zhao M, Du X. A simple added mass model for simulating elliptical cylinder vibrating in water under earthquake action [J]. Ocean Engineering, 2019, 179: 351 − 360. doi: 10.1016/j.oceaneng.2019.02.046 [15] Zhang G J, Li T Y, Zhu X, et al. Free and forced vibration characteristics of submerged finite elliptic cylindrical shell [J]. Ocean Engineering, 2017, 129: 92 − 106. doi: 10.1016/j.oceaneng.2016.11.014 [16] Liao W. Hydrodynamic interaction of flexible structures [J]. Journal of Waterway Port Coastal and Ocean Engineering, 1985, 111(4): 719 − 731. doi: 10.1061/(ASCE)0733-950X(1985)111:4(719) [17] Liaw C Y, Chopra A K. Earthquake analysis of axisymmetric towers partially submerged in water [J]. Earthquake Engineering & Structural Dynamics, 1975, 3(3): 233 − 248. [18] Williams A N. Earthquake response of submerged circular cylinder [J]. Ocean Engineering, 1986, 13(6): 569 − 585. doi: 10.1016/0029-8018(86)90040-5 [19] Aviles J, Li X. Hydrodynamic pressures on axisymmetric offshore structures considering seabed flexibility [J]. Computers & Structures, 2001, 79(2001): 2595 − 2606. [20] Park W, Yun C, Pyun C. Infinite elements for evaluation of hydrodynamic forces on offshore structures [J]. Computers & Structures, 1991, 40(4): 837 − 847. [21] Wang P, Zhao M, Du X. Short-crested, conical, and solitary wave forces on composite bucket foundation for an offshore wind turbine [J]. Journal of Renewable and Sustainable Energy, 2018, 10(2): 023305. doi: 10.1063/1.4995649 [22] Wang P, Zhao M, Du X. Simplified formula for earthquake-induced hydrodynamic pressure on round-ended and rectangular cylinders surrounded by water [J]. Journal of Engineering Mechanics, 2019, 145(2): 329 − 355. [23] Wang P, Zhao M, Du X, et al. Simplified evaluation of earthquake-induced hydrodynamic pressure on circular tapered cylinders surrounded by water [J]. Ocean Engineering, 2018, 164: 105 − 113. doi: 10.1016/j.oceaneng.2018.06.048 [24] Zhao M, Li H, Du X, et al. Time-domain stability of artificial boundary condition coupled with finite element for dynamic and wave problems in unbounded media [J]. International Journal of Computational Method, 2018, 15(3): 1850099. [25] 王丕光, 赵密, 李会芳, 等. 一种高精度圆柱形人工边界条件: 水-柱体相互作用问题[J]. 工程力学, 2019, 36(1): 88 − 95. Wang Piguang, Zhao Mi, Li Huifang, et al. A high-accuracy cylindrical artificial boundary condition: water-cylinder interaction problem[J]. Engineering Mechanics, 2019, 36(1): 88 − 95. (in Chinese) [26] Wang P, Wang X, Zhao M, et al. A numerical model for earthquake-induced hydrodynamic forces and wave forces on inclined circular cylinder[J]. Ocean Engineering, 2020, 207: 107382. [27] Song C, Wolf J P. The scaled boundary finite-element method-alias consistent infinitesimal finite-element cell method-for elastodynamics [J]. Computer Methods in Applied Mechanics and Engineering, 1997, 147(3−4): 329 − 355. doi: 10.1016/S0045-7825(97)00021-2 [28] Song C, Wolf J P. The scaled boundary finite-element method: analytical solution in frequency domain [J]. Computer Methods in Applied Mechanics and Engineering, 1998, 164(1): 249 − 264. [29] Song C, Wolf J P. Semi-analytical representation of stress singularities as occurring in cracks in anisotropic multi-materials with the scaled boundary finite-element method [J]. Computers & Structures, 2002, 80(2002): 183 − 197. [30] 崔灿, 蒋晗, 李映辉. 变截面梁横向振动特性半解析法[J]. 振动与冲击, 2012, 31(4): 85 − 88. Cui Can, Jiang Han, Li Yinghui. Semi-analytical method for calculating vibration characteristics of variable cross-section beam [J]. Journal of Vibration and Shock, 2012, 31(4): 85 − 88. (in Chinese) -
点击查看大图
图(7) / 表 (3)
计量
- 文章访问数: 139
- HTML全文浏览量: 35
- PDF下载量: 40
- 被引次数: 0
