顯式算法中三維黏彈性人工邊界應用及輸水隧洞動力響應分析

陳旻銘，任旭華，張繼勛

(河海大學水利水電學院，江蘇南京210024)

0 引 言

隨著地下綜合管道的加密，既有水工隧洞結構的安全與穩定問題愈加復雜。采用有限元法模擬無限域地基中隧洞的動力響應問題時，人工截斷邊界單元參數加載的合理性與準確性難以控制。動力人工邊界中由于黏彈性人工邊界具有克服高頻失穩和低頻飄移、精度高、實用性強等特點，許多學者對此開展了研究。劉晶波等[1]分析了二維黏彈性人工邊界和半無限空間的靜力問題的基本解，建立了三維靜-動力統一人工邊界理論。谷音等[2]基于集中黏彈性邊界理論，提出了等效三維一致黏彈性人工邊界單元。章小龍等[3]通過Fortran 程序編寫命令流的方式，給計算區域批量定義黏彈性邊界單元，實現了黏彈性邊界理論在ABAQUS/Standard中的運用。在動力分析中，人工邊界區的穩定性比內部計算區域的穩定條件更為嚴格，目前引入黏彈性邊界時，多采用隱式的無條件穩定算法。顯式算法與之相比效率優勢明顯，但存在邊界參數設置較為復雜且容易發生失穩等問題。李述濤等[4]對此提出在二維黏彈性人工邊界近端增加集中質量的方法，提高了動力分析在顯式算法中的穩定性。劉晶波等[5]通過求解邊界區域子系統的穩定性系數，改善其影響條件，從而提高二維黏彈性邊界的穩定性寬容度及計算效率。

本文在總結前人研究的基礎上，提出一種適用于動力顯式計算的三維黏彈性人工邊界施加方法，并依托某輸水隧洞與臨近高鐵隧道的立交工程，在設有三維黏彈性邊界的基礎上，研究高速列車振動荷載對既有水工隧洞的振動影響特性與規律。

1 動力顯式中黏彈性人工邊界的實現

1.1 三維黏彈性人工邊界系數的建立

黏彈性人工邊界是一種設置在人工截斷面節點上的彈簧阻尼系統，用以吸收計算模型向邊界輻射的能量。為了改善黏彈性人工邊界在顯式動力分析中的逐步積分穩定性，本文參考文獻[6]穩定性改善方法，推導了三維黏彈性人工邊界的穩定性改善公式，建立整個模型的人工邊界子系統，包括：①頂角子系統、②棱邊子系統和③側面子系統。布置形式見圖1。圖1中，MB為邊界節點的集中質量系數。

圖1 人工邊界子系統

三維時域下，黏彈性人工邊界的穩定性改善以中心差分方法[7]為基礎，逐步積分遞推公式如下

(1)

U(t+Δt)=A·U(t)+P(t)

(2)

式中，U為復合矩陣；A為人工邊界子系統的傳遞矩陣；P為外力向量。U的具體表達式為

(3)

由文獻[3]可知，子系統的穩定性積分與向量P無關，實則是求解傳遞矩陣A的譜半徑。此時，邊界穩定性判別條件ρ(A)=max|λi|的取值應在[-1，1]之間(λi為傳遞矩陣的特征值)，即當λi的模不大于1時，積分格式滿足穩定性條件。臨界時間步長為

(4)

式中，γ為邊界子系統的穩定性系數；L為有限單元的最小邊長；cP為壓縮波波速。改善后人工邊界上的彈簧剛度與阻尼系數與改進前相同，人工邊界近端的集中質量系數MB取值為

(5)

式中，βB為黏彈性人工邊界質量參數。結合文獻[2]中三維黏彈性邊界系數的遞推思路可知，當βB≥3ρL/2時，三維邊界子系統的穩定性系數均大于內部區域，即模型的穩定性僅由內部區域決定，臨界時間步長可通過改善內部區域的單元條件來提高；ρ為質量密度。因此，三維黏彈性邊界系數的計算公式可表示為

(6)

式中，KBT1、KBT2和KBN為第一切向、第二切向和法向的彈簧的剛度系數；CBT1、CBT2和CBN為阻尼元件第一切向、第二切向和法向的阻尼系數；MBT1、MBT2和MBN為邊界節點第一切向、第二切向和法向的集中質量系數；αT1、αT2和αN為黏彈性邊界第一切向、第二切向和法向的修正系數，三維取值參考李金友[8]等研究，取法向αN=1.3，切向αT1=αT2=0.7；G為巖土體的剪切模量；R為黏彈性邊界節點到軌道中心線或固定波源的最短距離；cS為剪切波的波速。cP和cS可表示為

(7)

(8)

式中，E為彈性模量；μ為泊松比。

1.2 黏彈性人工邊界施加

在計算黏彈性人工邊界的彈簧剛度、阻尼及附加質量時，將式(6)中的各系數乘以節點控制面積[9]即可。在ABAQUS/Explicit模塊，動力隱式階段設置的接地彈簧將會失效，此時需要重新定義接地端節點。將人工邊界的所有節點分別沿X、Y、Z軸復制并移動一段距離，作為動力顯式計算的遠端節點，設置遠端節點的Point mass/inertia，并將節點完全固定，使其等價于接地端，最后在遠端節點與原人工邊界的近端節點之間通過ABAQUS設置Springs/Dashpots單元即可，批量處理可通過Matlab程序完成。遠端節點的設置與黏彈性邊界施加見圖2。ABAQUS顯式動力分析中黏彈性邊界加載的整體過程見圖3。

圖2 遠端節點的定義與黏彈性邊界的施加

圖3 ABAQUS顯式計算流程

1.3 黏彈性人工邊界驗證

建立尺寸為50 m×50 m×100 m三維模型，見圖4。圖4中，A、B、C點分別為模型的頂面中心、幾何中心和底面中心。采用各邊邊長為5 m的單元離散，頂面為自由邊界面，其他面均為黏彈性人工邊界面。土體質量密度為2 000 kg/m3，彈性模量為208 MPa，泊松比為0.3。從底部邊界垂直入射1個前進方向沿Z軸的P波和1個前進方向沿X軸的SV波，P波和SV波的振動方向均為Z軸正向，持時0.25 s，其波形見圖5。由式(7)和(8)計算cP和cS分別為374.17、200 m/s。

圖4 顯式黏彈性邊界驗算模型

圖5 振動波位移曲線

圖6給出了顯式計算的振動波位移數值解。從圖6可知，P波、SV波分別在0.27、0.25 s時到達A點，經過自由面反射后向下傳播，分別于0.53、0.50 s時回到C點，0.78、0.75 s后，P波、SV波完全穿過模型底部并向無限遠處區域傳播，不再向上反射。振動波在A點的最大位移是B、C點的2倍，與文獻[10]中隱式求解的脈沖波位移時程曲線的變化規律完全一致，說明人工邊界具有良好的吸能效果，滿足遠域地基彈性恢復的性能，同時也說明了上述的黏彈性人工邊界在顯式求解中具有的穩定性，有效地削弱了波在傳遞過程中的振動疊加效應[11]。

圖6 振動波位移時程數值解

2 輸水隧洞動力響應分析

2.1 動力計算模型及材料參數

以某輸水隧洞與高鐵隧道的立體交叉工程為例，研究高鐵振動荷載對既有輸水隧洞的動力響應。上行線為馬蹄形無壓隧洞，凈斷面尺寸7.4 m×7.4 m，襯砌厚0.4 m，拱頂埋深40.3 m；下穿線為雙洞單線的列車隧道，內側直徑5.5 m，襯砌厚0.35 m。2條隧道中心線的間距為12.5 m，交叉洞線以90°立交，豎向最小凈距5.4 m，交叉段為Ⅳ級圍巖。截取三維模型的尺寸為90 m×85 m×85 m，見圖7。

圖7 交叉隧洞整體結構

數值計算僅考慮最不利工況，即下穿隧道左右線并排通行同向列車對輸水隧洞運行期的動力響應。列車最高時速為350 km/h，無壓輸水隧洞運行期水深約5 m，振動會導致水體對襯砌產生較大的作用，動水壓力參考文獻[12]以附加質量的形式施加。輸水隧洞監測斷面見圖8。觀測點布置見圖9。模型材料參數見表1。

表1 交叉隧洞材料參數

圖8 輸水隧洞計算斷面

圖9 輸水隧洞測點布置

2.2 振動荷載模型建立

軌道結構與列車既是組合體也是離散體，振動與輪軌不平順有關，可以歸納為：①軌道幾何不平順；②軌下基礎欠缺；③懸掛與車軸發生豎向相對運動。英國鐵路技術中心通過諸多理論與試驗分析了軌道隨機不平順與軌面局部波形磨損，對于時速200 km/h的列車建立軌道不平順管理標準，見表2[13]。

表2 英國軌道不平順管理值

梁波等[14]考慮軌道行駛的車輪間存在振動疊加現象與荷載在軌道上傳遞的分散性因素，提出了列車激振力函數F(t)表達式，即

F(t)=k1k2(f0+f1sin(ω1t)+f2sin(ω2t)+f3sin(ω3t))

(9)

以CRH3型車為例，列車單邊軸重8.5 t，簧下質量M0取750 kg。式(9)的函數應用在時速較高的列車上時，其振動荷載會出現負值現象，不符合列車-軌道耦合模型。于鶴然等[15]對此分析了軌面不平順波形，將表2中典型矢高進行了40%～60%的折減，得到3種控制條件下的典型波長及相應的正矢，分別為L1=10 m、α1=1.75 mm；L2=2 m、α2=0.2 mm；L3=0.5 m、α3=0.04 mm。當v=350 km/h時，激振力時程函數曲線見圖10。

圖10 激振力時程函數曲線

車輪荷載的位置和大小都是與時間相關的函數。為提高多軸移動荷載的施加效率，本文結合ABAQUS二次開發功能，采用Fortran程序編寫移動載荷。為使該程序便于三維模型的計算，作如下簡化：①各節車箱長度、輪對之間的間距均相等；②各節車的單邊軸重、簧下質量、車輪作用面積S均相等；③位于兩軌道上的同軸車輪的激振力時程完全相同。本文分析的列車由8節車箱編組，單邊車輪荷載分布見圖11。

圖11 高鐵列車單邊輪荷載分布(單位：m)

2.3 計算結果

振動荷載從駛入模型到完全駛出模型用時3.34 s，各斷面測點最大位移見圖12。從圖12可知，位于軌道左右線中部的S1斷面的豎向位移最大，距離右線豎直凈距最小的S2斷面其次，說明S1斷面同時受到左右線振動荷載的影響，其疊加響應效果大于主要受單線影響的S2斷面；隨著監測斷面與右線軌道距離的增加，S3、S4斷面的響應程度逐漸減弱，滿足振動波在土體介質中傳播時發生折減與發散的工程實況。各斷面監測點的沉降值由大到小依次為：底板中心>拱腳>拱肩>拱頂。

圖12 測點最大位移

分析各斷面底板中心加速度的響應情況。圖13為S1～S4斷面的底板中心加速度曲線。從圖13可知，S1～S4斷面加速度正向極值依次為0.079、0.044、0.024 m/s2和0.002 m/s2，負向極值依次為-0.076、-0.043、-0.024 m/s2和-0.002 m/s2。由此可見，下穿列車對上部輸水隧洞的動力響應主要以正向(豎直向上)作用為主，對于S4斷面的振動影響極其微弱，說明了黏彈性人工邊界充分吸收了振動波能，具有較強的實用性。

圖13 輸水隧洞底板中心加速度

3 結 語

本文基于二維黏彈性邊界理論及其在動力顯式分析中穩定性的改善方法，對ABAQUS顯式算法中三維黏彈性邊界的實現及既有輸水隧洞動力響應進行研究，得出以下結論：

(1)通過重新定義黏彈性人工邊界的遠端節點和近端節點，可實現邊界彈簧和阻尼在ABAQUS/Explicit中的應用。

(2)下行振動荷載對既有隧洞的動力響應隨著距離的增大而減小，同一斷面的響應程度由強到弱依次為底板中心、拱腳、拱肩、拱頂。

(3)三維模型顯式動力計算的穩定性是影響計算效率的主要因素，控制穩定性的條件仍需進一步研究。