寒區流冰撞擊輸水隧洞的數值模擬分析

王 鴻，貢 力，康春濤，楊軼群，王忠慧

(1.蘭州交通大學 土木工程學院，甘肅 蘭州730070；2.蘭州交通大學 調水工程及輸水安全研究所，甘肅 蘭州 730070)

1 研究背景

我國西北寒旱地區冬季氣溫低、冰期長且干旱缺水，為了緩解其干旱缺水的狀態，我國先后在該地區投資建設了諸多長距離輸調水工程，用于農業灌溉以及人畜飲水等[1-2]。然而由于該地區冬季氣候問題，輸水隧洞經常受到流冰的撞擊作用，長期的碰撞則會造成隧洞襯砌混凝土剝落或其他形式破壞，影響輸水隧洞的使用壽命。因此，開展流冰對輸水隧洞襯砌碰撞的研究迫在眉睫。

國內外學者在流冰與水工建筑物碰撞仿真模擬方面已取得了許多研究成果，Timco等[3]對冰的物理性質及力學性質進行了研究；單思鏑[4]在河冰對橋墩撞擊作用的數值模擬方面進行了分析；Rüdiger等[5]研究了在試驗過程中所使用冰的材料特性，并提出了考慮模型尺度的冰微觀結構和荷載作用下的各種物理效應的數值模型；張宿峰[6]對流冰的強度與厚度等對橋墩產生的等效應力方面進行了研究；Zhang等[7]利用LS-DYNA仿真模擬了流冰與T-rigid框架橋墩碰撞過程；貢力等[8]等在流冰對輸水隧洞撞擊作用力學特性方面進行了相應的研究與分析；Miryaha等[9]提出一種模擬浮冰對海洋結構影響的方法，討論了典型的浮冰作用模式和海上結構的壓力分布。綜上所述，國內外的學者對流冰與橋墩和海上結構撞擊研究較多，但針對流冰與寒旱地區引輸水隧洞的撞擊作用研究較少。

為了研究該類地區流冰對輸水隧洞撞擊作用的影響，本文采用理論分析結合數值模擬的方法，運用ANSYS/LS-DYNA有限元軟件及LS-PrePost后處理軟件完成流冰與隧洞撞擊的建模與計算分析，得到流冰速度及不同碰撞角度組合工況下對隧洞襯砌作用影響的規律，為西北寒旱地區冰期輸水安全提供支持。

2 流冰碰撞隧洞邊壁的數值模擬理論

碰撞力學按其一般分類可分為內部動力學和外部動力學。內部動力學主要研究內容為結構的變形、能量分配以及碰撞力等；外部動力學則研究兩個結構物碰撞后動能的變化、運動特性以及總能量耗散情況[10]。本文主要研究在碰撞過程中碰撞力及碰撞區變形量，因此采用內部動力學。ANSYS/LS-DYNA 是計算材料非線性以及結構非線性等問題比較精準的求解器，在爆炸和碰撞等領域內應用較為廣泛。Lagrange算法為該軟件的主要算法，且兼有Eluer和ALE兩種算法，其主要求解方式為顯式求解。對于固體力學中接觸-碰撞問題主要采用Lagrange描述方法建立有限元基本控制方程[11-12]，在Lagrange算法中，質點在隨著物體運動的整個過程中，其質量是始終保持恒定的。由連續介質力學理論可知，在流冰碰撞隧洞襯砌的過程中，需要滿足質量、動量以及能量三大守恒方程，本文采用Lagrange算法進行運動微分方程的計算。

2.1 質量守恒方程

ρ(X,t)J(X,t)=ρ0(X,0)

(1)

式中：ρ0(X,0)為t=0時刻初始碰撞模型中的介質密度，kg/m3；ρ(X,t)為t時刻碰撞模型中的介質密度，kg/m3;J(X,t)為碰撞模型的雅克比行列式。

2.2 動量守恒方程

(2)

式中：σij為現時構型的應力張量，Pa；bi為單位質量的體積力,N；üi為流冰和輸水隧洞模型質點的加速度，m/s2。

2.3 能量守恒方程

(3)

2.4 邊界條件

(4)

利用加權余量法，將動量方程表示為方程式(5)，其中加權系數取虛速度值。

(5)

利用方程式(5)，可以得到碰撞系統的控制方程：

(6)

式中：At為模型的邊界面。

單元內任意一點的X方向位移可以近似地用單元結點位移來表示：

ui(X,t)=xi(X,t)-Xi=Na(X)uia(t)

(7)

式中：uia(t)為結點a的位移，uia=xia(t)-Xia。

同理可得，單元內在任意一點X方向的加速度、速度、變形率以及虛速度分別可表示為：

(8)

式中：δvia為a結點的虛速度。

將所得到的公式(7)和(8)轉換為矩陣的形式，并且將二式引入公式(6)中，經過整理可得：

Mü=fext-fint

(9)

其中：M=M(u1(t))

(10)

f=(u1(t))

(11)

式中：M為系統的質量陣，與時間無關;U為系統的位移列陣，即ü為系統的加速度列陣;fint為系統的內力矩陣;fext為系統的外力矩陣。

對于接觸類型而言，LS-DYNA 豐富的接觸類型可以針對不同碰撞類型的問題進行設置，有動態約束、罰函數法以及分布參數法3種，且其主要算法為罰函數法。

3 流冰碰撞隧洞邊壁有限元仿真模擬

3.1 流冰與輸水隧洞模型建立

本文選取引大入秦工程盤道嶺3號隧洞作為有限元仿真模型。引大入秦總干渠全長86.79 km，總干渠設計以隧洞為主，總長75.376 km，占渠線總長的86.85%，其他建筑物占總長的13.15%。盤道嶺 3 號隧洞位于中國甘肅省永登縣，全長 15.723 km，其加大流量為36 m3/s，設計流量為32 m3/s，設計水深為2.92 m，加大水深為3.37 m，縱坡 1/1000。該隧洞的結構為反拱底板式斷面以及圓拱直墻形結構，隧洞的凈高為4.40 m，凈寬為4.20 m，半圓形拱頂，其半徑為2.10 m。

水內冰形成的基本條件為紊動和過冷卻，除此之外還需要水中存在著冰晶核，當冰晶核與周圍水體交換熱量時，冰晶粒則會逐漸增大。當周圍氣溫使水體溫度下降到0℃以下時, 冰花和岸冰開始初步形成 , 進入結冰期[13-14]。大通河流域處在高寒山區，冰期較長，封凍最長天數達100 d(1984年)，流冰時間較長，一般在11月下旬開始至次年3月中、下旬，流冰時間120～130 d。冰科學屬于較為復雜的交叉科學，涉及固體力學、熱力學、流體力學以及自然地理等學科, 其最基礎的參數為厚度和強度[15-16]，綜合考慮到大通河河心最大冰厚度和盤道嶺3號隧洞所處地理位置，本次模擬選取厚度為0.2 m長寬均為0.5 m的中小體積流冰進行碰撞仿真模擬。隧洞及流冰模型網格劃分如圖1所示。

圖1 流冰與隧洞模型網格劃分

將隧洞邊壁碰撞區網格加密為2.5 cm六面體網格[17]。流冰網格劃分采用全局統一劃分原則，將流冰網格劃分為2.5 cm六面體單元。

3.2 流冰與隧洞材料參數及碰撞場景參數

漂浮在水面的流冰處于豎向平衡狀態，所以在仿真模擬時忽略流冰所受的豎向荷載，只考慮水流拖動流冰運動的水平荷載[18]，導致其水平運動的荷載主要是以流冰的初速度來體現，不同初速度的流冰以不同的碰撞角度與隧洞襯砌進行碰撞時，碰撞力以及隧洞碰撞區的位移也是不同的。

本文在模擬仿真時流冰和隧洞材料均為線彈性材料，冰的材料參數參考自于天來等[19]對河冰力學性能的研究成果，具體參數如表1所示。流冰以及隧洞所選取的單元均為Solid164 3D實體單元；采用映射網格劃分方式；接觸類型選用計算較為簡單的單面自動接觸(ASSC)；通過設定流冰初速度、邊界條件以及碰撞求解參數后，由ANSYS輸出K文件導入LS-DYNA進行求解計算。

表1 流冰和輸水隧洞材料的物理力學參數

4 不同組合工況下碰撞仿真結果

流冰與隧洞襯砌的碰撞過程是一個較為復雜的冰與結構之間相互碰撞的問題，影響碰撞的因素主要包括流冰初速度、流冰體積、碰撞面積、幾何形狀、流冰厚度以及流冰與隧洞襯砌碰撞角度等[20-21]。本文主要研究流冰在不同流速、不同碰撞角度工況下對輸水隧洞襯砌的等效應力及隧洞碰撞區位移等。有限元仿真模擬中流冰與隧洞碰撞角度θ采用0°(正碰)、10°、20°、30°、45°共5種工況，碰撞角度示意圖如圖2所示(以θ=0°、θ=30°為例)。本次仿真模擬中根據引大入秦工程運營設計水深及最大水深選取隧洞控制水深為3 m，根據引大入秦設計流速及最大流速v選取隧洞內水流控制流速分別為1.5、2.2、2.6、3.0、3.5 m/s 5種工況，不考慮二次碰撞等問題[22]。

圖2 流冰與隧洞襯砌碰撞角度示意圖(單位：mm)

自然界中輸水工程所形成的流冰大多數屬于形狀不規則的板類結構，幾何尺寸大小與當地氣候條件密切相關，綜合考慮到大通河冰情和盤道嶺隧洞所處地理位置以及參照徐國賓等[23]和陳云飛[24]在流冰方面的研究，本文在模擬仿真時，流冰幾何形狀采用 0.5 m×0.5 m×0.2 m長方體板型結構模型冰。通過仿真模擬得到的流冰在碰撞角度為θ=0°、流速v=2.2 m/s時的內能和沙漏能時程圖、隧洞碰撞區X方向最大位移時程圖如圖3、4所示。

圖3 流冰碰撞的沙漏能和內能時程圖

冰作為一種非金屬的材料，在仿真模擬過程中非常容易出現沙漏情況，因此對網格劃分單元的規格有一定的要求。本文采用分區域劃分網格有效地降低了沙漏能，由圖3可以看出，在模擬過程中沙漏能約為內能峰值的5%，滿足在10%以內，所以模擬結果較為精確。由圖4可知，當t=0.004 59 s 時，碰撞角度為0°，流冰速度為2.2 m/s工況下的隧洞撞擊區X方向最大位移為9.2×10-5m，圖4中位移出現持續波動是因為在計算時考慮計算時間的問題，未設置阻尼。同理，可以得到碰撞角度為0°時，v=1.5、2.6、3.0、3.5 m/s 4種流速工況下對應的隧洞碰撞區X方向最大位移分別為2.23×10-5、1.45×10-4、1.89×10-4、2.35×10-4m。

圖4 隧洞碰撞區X方向最大位移時程圖

圖5為流冰在碰撞角度θ=0°、流速v=2.2 m/s時對隧洞撞擊的最大等效應力云圖。其碰撞過程屬于動態變化的過程，圖5表明，當t=0.004 59 s 時，可以得到等效應力瞬時最大值為8.995×106Pa。同理，可以得到碰撞角度為0°時，v=1.5、2.6、3.0、3.5 m/s 4種流速工況對應的最大等效應力分別為6.20×106、1.07×107、1.24×107、1.43×107Pa。其他20種不同組合工況下的相關數據如表2、3所示。

圖5 流冰對隧洞撞擊的最大等效應力云圖(θ=0°，v=2.2m/s)

表2 不同組合工況下隧道碰撞區X方向最大位移數據表 m

表3 不同組合工況下隧洞碰撞區最大等效應力數據表 Pa

由動能定理和動量定理可知，當一定質量的物體以不同的速度運動時，其所具有的動能也是不同的，它們之間呈現正相關關系。但是，冰是一種較為特殊的材料，其在一定的條件下會存在韌性向脆性轉變的特性，當等效應力超過某個值時，流冰就會產生破碎現象，所以撞擊力與等效應力之間會存在細微的差異情況。當流冰在撞擊隧洞襯砌出現破碎情況時，破碎的流冰仍有可能會與隧洞襯砌產生一定的二次碰撞，由于二次碰撞過程涉及的因素較多且復雜，但是對隧洞襯砌的影響較小，因而在本文的仿真模擬過程中不計算二次碰撞的情況。

5 流冰碰撞隧洞邊壁結果分析

5.1 流冰不同碰撞角度對隧洞撞擊的影響

分析表2、3中碰撞角度與等效應力及位移的關系，繪制如圖6、7所示的關系圖。

圖6 不同流速下流冰碰撞角度與最大等效應力的關系 圖7 不同流速下流冰碰撞角度與X方向最大位移的關系

由圖6可以看出，在保持流冰速度不變的情況下，只改變流冰與隧洞襯砌碰撞的碰撞角度，碰撞角度的變化對最大等效應力有一定的影響，但較為不明顯。最大等效應力隨碰撞角度的增大呈現出近似增大的趨勢，其變化關系在流冰速度v=1.5、2.2、2.6、3.0、3.5 m/s時均呈現出非線性(多項式)關系，不同流速工況下碰撞角度與最大等效應力的非線性關系式如公式(12)所示，其多項式各項系數值見表4。

表4 不同流速工況對應的公式(12)中各項系數值

y=mx3-nx2+px+c

(12)

式中：y為最大等效應力，Pa；x為流冰與隧洞襯砌碰撞角度，(°)

由圖7可以看出，在保持流冰速度不變的情況下，只改變流冰與隧洞襯砌的碰撞角度，碰撞角度對隧洞襯砌撞擊區域X方向最大位移有一定的影響。X方向位移隨碰撞角度的增大呈現出增大的趨勢，其變化關系在流冰速度v=1.5、2.2、2.6、3.0、3.5 m/s時均呈現出非線性關系，不同流速工況下碰撞角度與X方向最大位移的非線性關系式如公式(13)所示，其多項式各項系數值見表5。

表5 不同流速工況對應的公式(13)中各系數值

y=m′x3-n′x2+p′x+c′

(13)

式中：y為隧洞碰撞區X方向最大位移，m；x為流冰與隧洞襯砌碰撞角度，(°)。

5.2 流冰不同流速對隧洞撞擊的影響

分析表2、3中流冰速度與等效應力及位移的關系，繪制如圖8、9所示的關系圖。

圖8 不同碰撞角度下流冰速度與最大等效應力的關系

由圖8可以看出，在不改變流冰與隧洞襯砌碰撞角度的情況下，只改變流冰速度，流冰速度對最大等效應力的影響顯著，其最大等效應力均隨著流冰撞擊隧洞速度的增大而增大，且呈現出近似的線性關系,不同碰撞角度工況下流冰速度與最大等效應力的線性關系如公式(14)所示，其對應的各系數值見表6。

表6 不同碰撞角度工況對應的公式(14)中各系數值

y=ax-b

(14)

式中：y為最大等效應力，Pa；x為流冰速度，m/s。

由圖9可以看出，在保持流冰與隧洞襯砌碰撞角度不變的情況下，只改變流冰速度，流冰速度對隧洞襯砌撞擊區域X方向最大位移的影響顯著，其X方向最大位移向隨著流冰撞擊隧洞速度的增大而呈現出增大的趨勢，與宋安等[25]對冰速與冰力研究結果基本一致。其對應關系在本文5種碰撞角度工況均呈現出近似的線性關系，不同碰撞角度工況下流冰速度與X方向最大位移的線性關系如公式(15)所示，其對應的各系數值見表7。

圖9 不同碰撞角度下流冰速度與X方向最大位移的關系

表7 不同碰撞角度工況對應的公式(15)中各系數值

y=a′x-b′

(15)

式中：y為隧洞碰撞區X方向最大位移，m；x為流冰速度，m/s。

5.3 流冰不同流速和碰撞角度對隧洞撞擊的影響

通過上述對流冰速度以及與隧洞碰撞角度對隧洞作用的單一影響因素分析發現，等效應力及隧洞碰撞區X方向位移均會隨著兩種因素的變化產生較為明顯的變化，為了綜合分析兩種因素對隧洞作用的共同影響，通過Matlab軟件建立流冰速度-碰撞角度-等效應力和位移關系圖如圖10、11所示。

圖10 流冰速度-碰撞角度-最大等效應力關系圖

由圖10可以發現，在本文所模擬仿真的兩種組合工況下，碰撞角度在0°～10°過程中等效應力有明顯的增大，且增大率為0°～45°碰撞角度變化過程中的最大值，而在10°～45°過程中等效應力雖然也呈現出增大狀態，但增幅相對于0°～10°過程為較平緩，且在20°～45°過程中曲面斜率逐漸趨平；而流冰速度與等效應力則呈現近似線性關系，流冰速度較大時，等效應力相應的較大，且其變化過程中等效應力增大率沒有明顯的變化。綜合上述分析可認為，撞擊角度對最大等效應力的影響較不明顯，但流冰速度對最大等效應力的影響則較為明顯，且流冰速度對最大等效應力的影響要強于碰撞角度對其的影響。由圖10還可以發現，對等效應力影響最大的組合工況為θ=45°、v=3.5 m/s的組合工況。

由圖11可以發現，在本文所模擬仿真的兩種組合工況下，流冰速度和碰撞角度較小時，X方向最大位移增幅也較小，例如當v=2.2 m/s、θ=0°時，X方向最大位移增幅明顯小于較大流速和較大碰撞角度的趨勢；碰撞角度在0°～10°變化過程中隧洞襯砌撞擊區域X方向最大位移有明顯的增長狀態，且其增長率為0°～45°碰撞角度變化過程中的最大值；而流冰速度與X方向最大位移則呈現近似線性關系，流冰速度較小時，X方向最大位移相應的較小，其變化過程中X方向最大位移增長率的變化較為均勻。綜合上述分析可認為，撞擊角度對隧洞碰撞區X方向最大位移的影響較為不明顯，但流冰速度對最大位移的影響則較為明顯，且流冰速度對最大位移的影響要強于碰撞角度對最大位移的影響。由圖11可以進一步發現，對隧洞襯砌X方向最大位移影響最大的組合工況為θ=45°、v=3.5 m/s的組合工況。

圖11 流冰速度-碰撞角度-X方向最大位移關系圖

6 結 論

本文忽略流冰碰撞破碎后形成二次碰撞的因素，應用有限元分析軟件ANSYS/LS-DYNA以線彈性材料主要研究了碰撞角度及流冰速度組合工況下流冰對隧洞襯砌的碰撞影響，根據流冰的速度與不同碰撞角度的組合工況變化對輸水隧洞襯砌撞擊的影響得到了以下結論：

(1)不改變流冰速度的情況下，只改變流冰與隧洞襯砌的碰撞角度時，隧洞撞擊區X方向最大位移和最大等效應力等均隨著碰撞角度的增大而增大，且均呈現出非線性(多項式)關系。

(2)在不改變碰撞角度的情況下，只改變流冰速度時，隧洞撞擊區X方向最大位移和最大等效應力等均隨著流冰速度的增大而增大，且均呈現出近似線性的關系。

(3)綜合分析本文所模擬的25種組合工況可以發現對隧洞襯砌作用影響最大的組合工況為θ=45°、v=3.5 m/s的組合工況，且流冰與隧洞襯砌碰撞角度對隧洞襯砌作用影響較為不明顯，而流冰速度對其的影響卻較為明顯。

在實際工程中,漂浮在水面的流冰的運動是受多因素影響的，例如水與流冰之間的黏滯力和水流拖動流冰運動的拖曳力等，而對于中小型流冰長期碰撞隧洞襯砌導致其作用失效等問題，還需要做進一步的研究。因此，本文研究內容可以為后續流冰長期碰撞導致隧洞襯砌失效的研究提供一定的支持與參考。