李若雨，滕 斌，叢龍飛，侯志瑩，耿寶磊，金瑞佳

（1.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024；2.交通運輸部天津水運工程科學研究院 工程泥沙交通行業(yè)重點實驗室，天津 300456）

引 言

水下懸浮隧道（Submerged Floating Tunnel或阿基米德橋）是穿越海峽、內陸水域的一種新型結構方式，該結構利用自身封閉管道的浮力來維持使用載荷。正因為這種特性，SFT在跨越能力和建設成本方面具有巨大的潛在優(yōu)勢[1]。水下懸浮隧道由通車隧道和適當?shù)腻^固系統(tǒng)組成，根據(jù)隧道結構重力與浮力的相對大小，可選擇錨索或浮筒等多種方式支撐[2-4]。通過適當調整設計參數(shù)，可使結構在多種環(huán)境條件下工作。此外，懸浮隧道位于水下20～50 m左右，不影響船舶航行，對生態(tài)環(huán)境影響不大。正是這些積極的因素，使得研究人員不遺余力地推動SFT技術的發(fā)展。由于水下懸浮隧道屬于長細結構，有較大的變形性與可運動性，其動力響應在設計水下懸浮隧道時需要重點考慮。

世界范圍內水下懸浮隧道技術興起不久，現(xiàn)階段其理論模型研究還不夠深入，還在持續(xù)完善和試驗中。作為一種海洋結構，SFT必須面對來自水下環(huán)境的各種威脅，如海浪、波浪、洋流、海嘯和地震負荷。Kunisu[5]采用邊界元法研究了圓形和橢圓形斷面水下浮動隧道的波浪力特性。Remseth等[6]基于Navier-Stokes方程建立二維有限元模型，通過JONSWAP譜分析了隨機風浪作用下SFT的動力響應。Seo等[7]建立了簡化的動力學方程，并將數(shù)值計算結果與實驗結果進行了比較。Mulas等[8]用ANSYS程序和二維剛體動力學方法求解了隧道-車輛的相互作用。秦延飛[9]用 ANSYS程序對不同水位的水下結構進行模態(tài)分析，考慮流固耦合效應。

現(xiàn)有的文章大多采用Airy線性波理論，懸浮隧道的設計需要考慮極端環(huán)境，極端環(huán)境下波浪很大，因此本文采用二階Stokes波和橢圓余弦波理論進行計算?，F(xiàn)有的文章大多采用商業(yè)軟件或仿真程序對懸浮隧道進行有限元模擬計算，本文基于有限元理論，建立了波浪作用下懸浮隧道的時域耦合水彈性動力分析模型，將錨索作為無質量的靜力結構進行分析，求解隧道和錨索的動力平衡方程，研究系統(tǒng)的動力響應。波浪引起的水動力載荷由運動物體的Morison方程計算。為懸浮隧道動力響應分析提供一種計算手段。

1 懸浮隧道的有限元模型

1.1 隧道模型和計算模型的基本假設

本文對雙排隧道進行分析，隧道示意見圖 1。雙排隧道兩端共有四個端點嵌固到岸邊，雙排隧道中間有兩個過渡隧道將雙排隧道聯(lián)系在一起。取雙排隧道的軸向為x方向，z軸垂直向上，y軸與雙排隧道的軸向垂直。雙排隧道沿程布置了多個錨索，將隧道與海底連接，一個節(jié)點處的兩根錨索處在同一YOZ平面內。

圖1 懸浮隧道示意

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

1）簡化流體對隧道結構的作用，只考慮沿隧道軸線分布的均布荷載，不考慮水流帶來的渦激振動作用，同時結構浮力簡化為沿隧道軸線分布的均布荷載；

2）不考慮錨索的質量與動力響應，僅考慮錨索會對隧道施加力，將錨索簡化為y、z方向的兩根無質量彈簧，錨索僅在隧道振動時提供y、z方向的彈力。

本文隧道斷面形式采用圓形，圓形是最穩(wěn)定的流體靜力學結構形式，用垂直于隧道軸向的平面將隧道切分為很短、等長的圓環(huán)體，每個圓環(huán)體作為一個單元，相鄰單元間傳遞彎矩和轉角，單元采用兩端固定梁模型，單元受力為重力、浮力、纜索預拉力、Morison波浪力和振動狀態(tài)下的錨索拉力，錨索預拉力和振動狀態(tài)下的拉力處理方式相同，按節(jié)點處連接單元數(shù)均分到每個單元上，如錨索所在節(jié)點連接三個單元，則每個單元受拉力為總拉力的1/3，并轉化為均布線荷載。重力、浮力均為線荷載。根據(jù)結構力學兩端固定梁在均布荷載下的受力情況，可求出單元左右節(jié)點處六個自由度方向的受力。求解單元動力平衡方程，可得單元左右節(jié)點的受力及運動情況，所有節(jié)點的受力及運動情況即為結構的動力響應。

將結構切分為m個單元，n個節(jié)點，單元間約束條件為相鄰單元接觸面的位移、轉角及受力情況相同，本文的錨索一律加在單元端部。

1.2 控制方程

隧道在波浪力作用下的動力響應屬于多自由度系統(tǒng)的振動問題，單元與結構滿足動力平衡方程。

1）建立單元動力平衡方程

結構離散為若干個單元之后，單元任意節(jié)點的位移是坐標與時間的函數(shù)，它對時間的一階導數(shù)，即為節(jié)點速度，定義為：

設在單元節(jié)點上施加作用力F，利用虛功原理，則有：

式中：δ*為單元節(jié)點的虛應變；d*為單元節(jié)點的虛位移；σ為節(jié)點應力；v為阻尼系數(shù)；ρ為材料密度。

引入單元彈性矩陣D，單元應變矩陣B，單元形狀函數(shù)矩陣N。節(jié)點劃分足夠小時，有：

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

式（1）變?yōu)椋?/p>

節(jié)點虛應變與應變以及速度與加速度均可視為常數(shù)，且對于任意值的虛應變，上式均成立?？傻茫?/p>

2）建立結構動力平衡方程

設式（4）中節(jié)點由局部坐標到單元整體坐標的變換陣為P，可得：

即：

將整個結構系統(tǒng)節(jié)點位移矢量記為X，并將各單元整體坐標系下變換矩陣，，和F進行組裝，得整體坐標下相應矩陣M、K、C。根據(jù)達朗貝爾原理，整個系統(tǒng)的動力平衡方程為：

其中：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節(jié)點產(chǎn)生單位位移時i節(jié)點的受力矩陣。前文提到將結構切分為m個單元，n個節(jié)點，矩陣K（6n×6n）為整個結構的剛度矩陣，矩陣K由n×n個子矩陣Kij組裝而成。即：

m個單元共有4m個子矩陣，對應填入上式即可得出不計錨索情況下整個結構的剛度矩陣K。

構建錨索的剛度矩陣Km（6×6），設錨索y向與z向彈簧屬性完全相同，其中E為錨索彈性模量，A為錨索截面面積，L為錨索長度，若i節(jié)點設置錨索，則用KII=Kii+Km替代Kii，可構造出考慮錨索作用的總體剛度矩陣K。

b.質量矩陣M

單元ij作為三維結構，在局部坐標下（單元軸向為局部坐標系x軸方向）有12×12的單元質量矩陣eij：

c.阻尼矩陣C

阻尼矩陣C采用Rayleigh阻尼，由剛度矩陣、質量矩陣和系數(shù)構成：C=a0M+a1K。

d.受力矩陣F

將結構切分為n個節(jié)點，第i個節(jié)點6個自由度方向的受力情況分別為，按節(jié)點編號順序排列：……組成6n×1的列矩陣，即受力矩陣F。

元素排列順序與F相同，為6n×1的矩陣。

隧道在波浪力作用下的動力響應屬于多自由度系統(tǒng)的振動問題，隧道整體需滿足動力平衡方程，解出隧道在波浪力的動力平衡方程，即可完成隧道的動力響應分析。

2 波浪作用力

結構單元上的波浪荷載可視為均布荷載，用Morison方程求出單元中點處的波浪力，將其作為整個單元上的均布線荷載。

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

Morison方程：

橢圓余弦波理論：

式中：u為波浪x方向速度；w為波浪z方向速度；η為橢圓余弦波波面升高；H為波高；c為波速；cn為橢圓余弦函數(shù)；m為橢圓積分參數(shù)；K和E分別為第一類橢圓積分和第二類橢圓積分；h為水深；ηx為波面對x求一階導；ηxx和ηxxx類推。

即可求出橢圓余弦波理論中水質點的速度與加速度。

二階 Stokes波理論：在有限振幅情況中，設φ(x,z,t)為速度勢函數(shù)，有：

對于流體中某一點（x，z），波動的水平方向速度u和豎直方向速度w分別為：

水平方向加速度ax和豎直方向加速度az分別為：

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

3 時域模擬

1）生成剛度陣K、質量陣M和阻尼陣C

3）計算積分常數(shù)：

4）計算等效剛度，并作Cholesky分解：

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

6）tk+1時刻位移：

7）tk+1時刻響應向量：

重復上述步驟即可計算出每個時間步隧道各節(jié)點的位移、速度、加速度矩陣，即隧道的動力響應。本文中的大規(guī)模矩陣是稀疏矩陣（數(shù)值為0的元素數(shù)目遠遠多于非0元素的數(shù)目，并且非0元素分布沒有規(guī)律），且稀疏度較高。因此采用有高效的稀疏矩陣存儲格式：CSR格式（Compressed Sparse Row），可大幅減少運算量。

4 模型驗證

為驗證本模型可信度，對龍旭等[10]論文中的算例進行計算，隧道設計長度為100 m，軸線位于水下4.2 m，錨索連接于管體軸向坐標30 m，50 m，70 m處。隧道所處水域的水文數(shù)據(jù)和計算模型物理參數(shù)如表1所示。

表1 水文數(shù)據(jù)和物理參數(shù)

波浪沿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方程計算，質量系數(shù)CM=2.0，拖曳力系數(shù)CD=1.2。阻尼矩陣C采用Rayleigh阻尼，α0=1.1042，α1=0.00165。

模擬80 s內隧道的動力響應，加入受力緩沖函數(shù)，前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過渡階段后，隧道響應幅值隨時間變化較為規(guī)則。

4）研究結果表明，該方法能較準確地計算出水下浮動隧道的動力響應。