波浪作用下悬浮隧道动力响应的有限元分析模型

李若雨，滕 斌，丛龙飞，侯志莹，耿宝磊，金瑞佳

（1.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024；2.交通运输部天津水运工程科学研究院 工程泥沙交通行业重点实验室，天津 300456）

引 言

水下悬浮隧道（Submerged Floating Tunnel或阿基米德桥）是穿越海峡、内陆水域的一种新型结构方式，该结构利用自身封闭管道的浮力来维持使用载荷。正因为这种特性，SFT在跨越能力和建设成本方面具有巨大的潜在优势[1]。水下悬浮隧道由通车隧道和适当的锚固系统组成，根据隧道结构重力与浮力的相对大小，可选择锚索或浮筒等多种方式支撑[2-4]。通过适当调整设计参数，可使结构在多种环境条件下工作。此外，悬浮隧道位于水下20～50 m左右，不影响船舶航行，对生态环境影响不大。正是这些积极的因素，使得研究人员不遗余力地推动SFT技术的发展。由于水下悬浮隧道属于长细结构，有较大的变形性与可运动性，其动力响应在设计水下悬浮隧道时需要重点考虑。

世界范围内水下悬浮隧道技术兴起不久，现阶段其理论模型研究还不够深入，还在持续完善和试验中。作为一种海洋结构，SFT必须面对来自水下环境的各种威胁，如海浪、波浪、洋流、海啸和地震负荷。Kunisu[5]采用边界元法研究了圆形和椭圆形断面水下浮动隧道的波浪力特性。Remseth等[6]基于Navier-Stokes方程建立二维有限元模型，通过JONSWAP谱分析了随机风浪作用下SFT的动力响应。Seo等[7]建立了简化的动力学方程，并将数值计算结果与实验结果进行了比较。Mulas等[8]用ANSYS程序和二维刚体动力学方法求解了隧道-车辆的相互作用。秦延飞[9]用 ANSYS程序对不同水位的水下结构进行模态分析，考虑流固耦合效应。

现有的文章大多采用Airy线性波理论，悬浮隧道的设计需要考虑极端环境，极端环境下波浪很大，因此本文采用二阶Stokes波和椭圆余弦波理论进行计算。现有的文章大多采用商业软件或仿真程序对悬浮隧道进行有限元模拟计算，本文基于有限元理论，建立了波浪作用下悬浮隧道的时域耦合水弹性动力分析模型，将锚索作为无质量的静力结构进行分析，求解隧道和锚索的动力平衡方程，研究系统的动力响应。波浪引起的水动力载荷由运动物体的Morison方程计算。为悬浮隧道动力响应分析提供一种计算手段。

1 悬浮隧道的有限元模型

1.1 隧道模型和计算模型的基本假设

本文对双排隧道进行分析，隧道示意见图 1。双排隧道两端共有四个端点嵌固到岸边，双排隧道中间有两个过渡隧道将双排隧道联系在一起。取双排隧道的轴向为x方向，z轴垂直向上，y轴与双排隧道的轴向垂直。双排隧道沿程布置了多个锚索，将隧道与海底连接，一个节点处的两根锚索处在同一YOZ平面内。

图1 悬浮隧道示意

建立模型前的基本假设如下：

1）简化流体对隧道结构的作用，只考虑沿隧道轴线分布的均布荷载，不考虑水流带来的涡激振动作用，同时结构浮力简化为沿隧道轴线分布的均布荷载；

2）不考虑锚索的质量与动力响应，仅考虑锚索会对隧道施加力，将锚索简化为y、z方向的两根无质量弹簧，锚索仅在隧道振动时提供y、z方向的弹力。

本文隧道断面形式采用圆形，圆形是最稳定的流体静力学结构形式，用垂直于隧道轴向的平面将隧道切分为很短、等长的圆环体，每个圆环体作为一个单元，相邻单元间传递弯矩和转角，单元采用两端固定梁模型，单元受力为重力、浮力、缆索预拉力、Morison波浪力和振动状态下的锚索拉力，锚索预拉力和振动状态下的拉力处理方式相同，按节点处连接单元数均分到每个单元上，如锚索所在节点连接三个单元，则每个单元受拉力为总拉力的1/3，并转化为均布线荷载。重力、浮力均为线荷载。根据结构力学两端固定梁在均布荷载下的受力情况，可求出单元左右节点处六个自由度方向的受力。求解单元动力平衡方程，可得单元左右节点的受力及运动情况，所有节点的受力及运动情况即为结构的动力响应。

将结构切分为m个单元，n个节点，单元间约束条件为相邻单元接触面的位移、转角及受力情况相同，本文的锚索一律加在单元端部。

1.2 控制方程

隧道在波浪力作用下的动力响应属于多自由度系统的振动问题，单元与结构满足动力平衡方程。

1）建立单元动力平衡方程

结构离散为若干个单元之后，单元任意节点的位移是坐标与时间的函数，它对时间的一阶导数，即为节点速度，定义为：

设在单元节点上施加作用力F，利用虚功原理，则有：

式中：δ*为单元节点的虚应变；d*为单元节点的虚位移；σ为节点应力；v为阻尼系数；ρ为材料密度。

引入单元弹性矩阵D，单元应变矩阵B，单元形状函数矩阵N。节点划分足够小时，有：

d=Nδ,ε=Bδ,σ=DBδ

式（1）变为：

节点虚应变与应变以及速度与加速度均可视为常数，且对于任意值的虚应变，上式均成立。可得：

2）建立结构动力平衡方程

设式（4）中节点由局部坐标到单元整体坐标的变换阵为P，可得：

即：

将整个结构系统节点位移矢量记为X，并将各单元整体坐标系下变换矩阵，，和F进行组装，得整体坐标下相应矩阵M、K、C。根据达朗贝尔原理，整个系统的动力平衡方程为：

其中：M为整条隧道的质量矩阵，为 6n×6n矩阵；K为整条隧道的刚度矩阵，为6n×6n矩阵；C为整条隧道的阻尼矩阵，为 6n×6n矩阵；，，X分别为隧道加速度、速度和位移矩阵，为6n×1矩阵；F为隧道受力矩阵，为6n×1矩阵。

a.刚度矩阵K

单元ij作为三维结构，在局部坐标下（单元轴向为局部坐标系x轴方向）有 12×12的单元刚度矩阵eij：

设Ax,Ay,Az分别是将整体坐标系旋转为局部坐标系时，绕x、y、z轴旋转的角度，则线性变换矩阵为：

局部坐标系与整体坐标系的转换矩阵为：

整体坐标下的单元刚度矩阵Keij为 12×12矩阵，可分成4个6×6的子矩阵Kii、Kij、Kjj、Kji，其中Kij表示j节点产生单位位移时i节点的受力矩阵。前文提到将结构切分为m个单元，n个节点，矩阵K（6n×6n）为整个结构的刚度矩阵，矩阵K由n×n个子矩阵Kij组装而成。即：

m个单元共有4m个子矩阵，对应填入上式即可得出不计锚索情况下整个结构的刚度矩阵K。

构建锚索的刚度矩阵Km（6×6），设锚索y向与z向弹簧属性完全相同，其中E为锚索弹性模量，A为锚索截面面积，L为锚索长度，若i节点设置锚索，则用KII=Kii+Km替代Kii，可构造出考虑锚索作用的总体刚度矩阵K。

b.质量矩阵M

单元ij作为三维结构，在局部坐标下（单元轴向为局部坐标系x轴方向）有12×12的单元质量矩阵eij：

c.阻尼矩阵C

阻尼矩阵C采用Rayleigh阻尼，由刚度矩阵、质量矩阵和系数构成：C=a0M+a1K。

d.受力矩阵F

将结构切分为n个节点，第i个节点6个自由度方向的受力情况分别为，按节点编号顺序排列：……组成6n×1的列矩阵，即受力矩阵F。

元素排列顺序与F相同，为6n×1的矩阵。

隧道在波浪力作用下的动力响应属于多自由度系统的振动问题，隧道整体需满足动力平衡方程，解出隧道在波浪力的动力平衡方程，即可完成隧道的动力响应分析。

2 波浪作用力

结构单元上的波浪荷载可视为均布荷载，用Morison方程求出单元中点处的波浪力，将其作为整个单元上的均布线荷载。

Morison方程被广泛地用于估算作用在水下物体上的流体力。如果速度和流体的加速度是在物体的表示位置上给出的，可以很容易地求出包括阻力和力在内的波浪力。

Morison方程：

椭圆余弦波理论：

式中：u为波浪x方向速度；w为波浪z方向速度；η为椭圆余弦波波面升高；H为波高；c为波速；cn为椭圆余弦函数；m为椭圆积分参数；K和E分别为第一类椭圆积分和第二类椭圆积分；h为水深；ηx为波面对x求一阶导；ηxx和ηxxx类推。

即可求出椭圆余弦波理论中水质点的速度与加速度。

二阶 Stokes波理论：在有限振幅情况中，设φ(x,z,t)为速度势函数，有：

对于流体中某一点（x，z），波动的水平方向速度u和竖直方向速度w分别为：

水平方向加速度ax和竖直方向加速度az分别为：

带入Morison方程中可求波浪力。

3 时域模拟

1）生成刚度阵K、质量阵M和阻尼阵C

3）计算积分常数：

4）计算等效刚度，并作Cholesky分解：

5）tk+1时刻的等效力：

6）tk+1时刻位移：

7）tk+1时刻响应向量：

重复上述步骤即可计算出每个时间步隧道各节点的位移、速度、加速度矩阵，即隧道的动力响应。本文中的大规模矩阵是稀疏矩阵（数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律），且稀疏度较高。因此采用有高效的稀疏矩阵存储格式：CSR格式（Compressed Sparse Row），可大幅减少运算量。

4 模型验证

为验证本模型可信度，对龙旭等[10]论文中的算例进行计算，隧道设计长度为100 m，轴线位于水下4.2 m，锚索连接于管体轴向坐标30 m，50 m，70 m处。隧道所处水域的水文数据和计算模型物理参数如表1所示。

表1 水文数据和物理参数

波浪沿Y轴正向入射时，两种方法管体中点在Y方向的位移响应对比如图2所示。

由对比图知，本文使用模型与龙旭等采用的ANSYS软件方法拟合程度良好，模型具有可信度。

图2 两种方法管体中点在Y方向的位移对比

5 算例分析

水下悬浮隧道全长 20 000 m，过渡隧道长40 m，每500 m设置一对锚索，每对包括两根锚索对称约束，锚索预拉力200 kN。计算模型中有限单元长度为 20 m。考虑极端天气和海啸等灾害的作用，本算例取波浪周期24.85 s（对应波长520 m），波高取为6 m，8 m，10 m，12 m，14 m，16 m，18 m和20 m，波浪入射方向与x轴向夹角取三组，分别为 30°，40°，50°，60°，70°，80°和 90°。作用于水下悬浮隧道的波浪力采用Morison方程计算，质量系数CM=2.0，拖曳力系数CD=1.2。阻尼矩阵C采用Rayleigh阻尼，α0=1.1042，α1=0.00165。

模拟80 s内隧道的动力响应，加入受力缓冲函数，前15 s属于缓冲期。

入射波高改变，波浪沿Y轴方向入射，隧道Y方向与Z方向位移幅值见图3。

图3 随入射波高改变，隧道Y方向与Z方向位移幅值

其他波浪要素保持不变，波浪沿Y轴向入射时，随着波浪波高由6 m增加至20 m，悬浮隧道Y方向与Z方向的位移幅值呈线性增加，这一结论与直观认识相符合。结果表明隧道Y方向响应要略小于Z方向响应。

波高取20 m，波浪沿与X轴向不同角度入射，其他波浪要素保持不变，隧道Y方向与Z方向位移幅值见图4。

图4 隧道Y方向与Z方向位移幅值

其他波浪要素保持不变的情况下，波浪传播方向与X轴向夹角越大，隧道动力响应越明显，当夹角小于 50°时，随入射角增加位移幅值增加程度不明显，夹角在 50°～70°之间时，随夹角增加位移幅值增加程度明显。结果同时表明波浪不垂直于隧道轴向时，仍有Y方向响应小于Z方向响应的结论。

波高16 m沿Y方向入射，隧道位移幅值随时间变化结果见图5。

图5 隧道Y方向与Z方向位移幅值随时间变化

算例中，隧道最危险情况为波高20 m波浪沿Y轴正向入射，强度较小的 C25混凝土抗剪应力为0.99 MPa，抗拉应力为9.5 MPa，结构安全。

6 结 论

1）其他波浪要素保持不变，波浪沿Y轴向入射时，随着波浪波高由6 m增加至20 m，悬浮隧道Y方向与Z方向的位移幅值呈线性增加，这一结论与直观认识相符合。结果表明隧道Y方向响应要略小于Z方向响应。

2）其他波浪要素保持不变，波浪传播方向与X轴向夹角越大，隧道动力响应越明显，当夹角小于50°时，随入射角增加位移幅值增加程度不明显，夹角在 50°～70°之间时，随夹角增加位移幅值增加程度明显。结果同时表明波浪不垂直于隧道轴向时，仍有Y方向响应小于Z方向响应的结论。

3）前15 s过渡阶段后，隧道响应幅值随时间变化较为规则。

4）研究结果表明，该方法能较准确地计算出水下浮动隧道的动力响应。