波浪作用下懸浮隧道動力響應(yīng)的有限元分析模型

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

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

引 言

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

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

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

1 懸浮隧道的有限元模型

1.1 隧道模型和計(jì)算模型的基本假設(shè)

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

圖1 懸浮隧道示意

建立模型前的基本假設(shè)如下：

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

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

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

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

1.2 控制方程

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

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

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

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

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

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

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

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

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

2）建立結(jié)構(gòu)動力平衡方程

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

即：

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

其中：M為整條隧道的質(zhì)量矩陣，為 6n×6n矩陣；K為整條隧道的剛度矩陣，為6n×6n矩陣；C為整條隧道的阻尼矩陣，為 6n×6n矩陣；，，X分別為隧道加速度、速度和位移矩陣，為6n×1矩陣；F為隧道受力矩陣，為6n×1矩陣。

a.剛度矩陣K

單元ij作為三維結(jié)構(gòu)，在局部坐標(biāo)下（單元軸向?yàn)榫植孔鴺?biāo)系x軸方向）有 12×12的單元剛度矩陣eij：

設(shè)Ax,Ay,Az分別是將整體坐標(biāo)系旋轉(zhuǎn)為局部坐標(biāo)系時，繞x、y、z軸旋轉(zhuǎn)的角度，則線性變換矩陣為：

局部坐標(biāo)系與整體坐標(biāo)系的轉(zhuǎn)換矩陣為：

整體坐標(biāo)下的單元剛度矩陣Keij為 12×12矩陣，可分成4個6×6的子矩陣Kii、Kij、Kjj、Kji，其中Kij表示j節(jié)點(diǎn)產(chǎn)生單位位移時i節(jié)點(diǎn)的受力矩陣。前文提到將結(jié)構(gòu)切分為m個單元，n個節(jié)點(diǎn)，矩陣K（6n×6n）為整個結(jié)構(gòu)的剛度矩陣，矩陣K由n×n個子矩陣Kij組裝而成。即：

m個單元共有4m個子矩陣，對應(yīng)填入上式即可得出不計(jì)錨索情況下整個結(jié)構(gòu)的剛度矩陣K。

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

b.質(zhì)量矩陣M

單元ij作為三維結(jié)構(gòu)，在局部坐標(biāo)下（單元軸向?yàn)榫植孔鴺?biāo)系x軸方向）有12×12的單元質(zhì)量矩陣eij：

c.阻尼矩陣C

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

d.受力矩陣F

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

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

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

2 波浪作用力

結(jié)構(gòu)單元上的波浪荷載可視為均布荷載，用Morison方程求出單元中點(diǎn)處的波浪力，將其作為整個單元上的均布線荷載。

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

Morison方程：

橢圓余弦波理論：

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

即可求出橢圓余弦波理論中水質(zhì)點(diǎn)的速度與加速度。

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

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

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

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

3 時域模擬

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

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

4）計(jì)算等效剛度，并作Cholesky分解：

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

6）tk+1時刻位移：

7）tk+1時刻響應(yīng)向量：

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

4 模型驗(yàn)證

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

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

波浪沿Y軸正向入射時，兩種方法管體中點(diǎn)在Y方向的位移響應(yīng)對比如圖2所示。

由對比圖知，本文使用模型與龍旭等采用的ANSYS軟件方法擬合程度良好，模型具有可信度。

圖2 兩種方法管體中點(diǎn)在Y方向的位移對比

5 算例分析

水下懸浮隧道全長 20 000 m，過渡隧道長40 m，每500 m設(shè)置一對錨索，每對包括兩根錨索對稱約束，錨索預(yù)拉力200 kN。計(jì)算模型中有限單元長度為 20 m。考慮極端天氣和海嘯等災(zāi)害的作用，本算例取波浪周期24.85 s（對應(yīng)波長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方程計(jì)算，質(zhì)量系數(shù)CM=2.0，拖曳力系數(shù)CD=1.2。阻尼矩陣C采用Rayleigh阻尼，α0=1.1042，α1=0.00165。

模擬80 s內(nèi)隧道的動力響應(yīng)，加入受力緩沖函數(shù)，前15 s屬于緩沖期。

入射波高改變，波浪沿Y軸方向入射，隧道Y方向與Z方向位移幅值見圖3。

圖3 隨入射波高改變，隧道Y方向與Z方向位移幅值

其他波浪要素保持不變，波浪沿Y軸向入射時，隨著波浪波高由6 m增加至20 m，懸浮隧道Y方向與Z方向的位移幅值呈線性增加，這一結(jié)論與直觀認(rèn)識相符合。結(jié)果表明隧道Y方向響應(yīng)要略小于Z方向響應(yīng)。

波高取20 m，波浪沿與X軸向不同角度入射，其他波浪要素保持不變，隧道Y方向與Z方向位移幅值見圖4。

圖4 隧道Y方向與Z方向位移幅值

其他波浪要素保持不變的情況下，波浪傳播方向與X軸向夾角越大，隧道動力響應(yīng)越明顯，當(dāng)夾角小于 50°時，隨入射角增加位移幅值增加程度不明顯，夾角在 50°～70°之間時，隨夾角增加位移幅值增加程度明顯。結(jié)果同時表明波浪不垂直于隧道軸向時，仍有Y方向響應(yīng)小于Z方向響應(yīng)的結(jié)論。

波高16 m沿Y方向入射，隧道位移幅值隨時間變化結(jié)果見圖5。

圖5 隧道Y方向與Z方向位移幅值隨時間變化

算例中，隧道最危險情況為波高20 m波浪沿Y軸正向入射，強(qiáng)度較小的 C25混凝土抗剪應(yīng)力為0.99 MPa，抗拉應(yīng)力為9.5 MPa，結(jié)構(gòu)安全。

6 結(jié) 論

1）其他波浪要素保持不變，波浪沿Y軸向入射時，隨著波浪波高由6 m增加至20 m，懸浮隧道Y方向與Z方向的位移幅值呈線性增加，這一結(jié)論與直觀認(rèn)識相符合。結(jié)果表明隧道Y方向響應(yīng)要略小于Z方向響應(yīng)。

2）其他波浪要素保持不變，波浪傳播方向與X軸向夾角越大，隧道動力響應(yīng)越明顯，當(dāng)夾角小于50°時，隨入射角增加位移幅值增加程度不明顯，夾角在 50°～70°之間時，隨夾角增加位移幅值增加程度明顯。結(jié)果同時表明波浪不垂直于隧道軸向時，仍有Y方向響應(yīng)小于Z方向響應(yīng)的結(jié)論。

3）前15 s過渡階段后，隧道響應(yīng)幅值隨時間變化較為規(guī)則。

4）研究結(jié)果表明，該方法能較準(zhǔn)確地計(jì)算出水下浮動隧道的動力響應(yīng)。