沉箱碼頭地震液化任意拉格朗日歐拉方法大變形分析

劉 順，唐小微，趙 祥

(海岸和近海工程國家重點實驗室(大連理工大學)，遼寧 大連 116023)

重力式沉箱結構廣泛應用于港口碼頭和擋墻岸壁。港口碼頭常處于河流沖積區域，地基土軟弱，常需對沉箱基礎進行置換和回填處理，常用材料有砂土、碎石土等。圖1為典型的沉箱結構示意，沉箱底部和墻后一定范圍內的海底黏土需置換為強度較高的砂土或碎石土等以滿足承載力和正常工作要求。由于地下水位較高，在波浪或地震荷載作用下，基底和墻后置換砂土的動力液化可能導致沉箱碼頭出現傾斜、向海側移動、沉降等災害[1]。

圖1 沉箱碼頭示意

Taiyab等[2]研究了動力作用下沉箱碼頭結構與土的相互作用，并分析了置換區和回填區土體密實度對沉箱結構抗震能力的影響，發現沉箱趾部地基土體的致密化能顯著減輕沉箱碼頭的位移響應。Alyami等[3]探究了滲透系數對沉箱碼頭結構孔隙水壓力的影響以及回填土和地基土的相對密實度對沉箱碼頭抵抗地震殘余變形的影響。Dakoulas等[4]研究了沉箱碼頭在地震荷載作用下側向位移、塑性應變和孔隙水壓力的演化規律，發現了沉箱碼頭慣性力、墻后土壓力和孔隙水壓力之間復雜的相互作用。王麗艷等[5]詳細討論了地震動水平、結構寬高比等因素對液化地基中沉箱碼頭殘余變形的影響規律。佟彬等[6]分析了振沖加密優化設計對沉箱碼頭液化災害的減弱作用。

上述基于有限元方法的沉箱碼頭地震響應研究對工程領域沉箱結構減災設計具有重要指導作用。但對于存在較大變形的情況，傳統有限元方法可能由于網格畸變出現分析中斷、精度喪失等問題而不能正確反映結構的動力響應過程。對于重力式沉箱碼頭，在置換區和回填區土體密實度較小的情況下，嚴重的液化變形可能導致網格畸變從而引起分析中斷?；诖耍_發了一種解耦合任意拉格朗日歐拉方法(ALE方法)解決強震作用下沉箱碼頭大變形導致的網格畸變、分析中斷、精度喪失等問題，并將其用于沉箱碼頭的液化減災措施評價。

1 控制方程及ALE方法簡介

1.1 飽和砂土控制方程

基于水土二相混合理論，將土骨架位移u和孔隙水壓力p作為變量，Uzuoka[7]建立了u-p形式的飽和砂土場方程式：

(1)

采用虛位移原理對式(1)進行空間離散，其中超孔隙水壓力二階梯度項采用有限差分法進行空間離散[8]，得到FEM-FDM耦合的二相混合體空間離散方程，即

[KV]{t+ΔtpE}={F}

(2)

1.2 ALE方法簡介

Donea等[10]結合拉格朗日方法和歐拉方法的優點，發展了一種任意拉格朗日歐拉方法(Arbitrary Lagrangian Eulerian Method)用以解決流體力學問題，此后該方法被眾多研究者引入固體力學的非線性問題研究中。

ALE方法中，場方程式的建立和離散基于計算參考域，任意時刻t，式(1)的各變量為當前時刻網格域對應的值，而變量的更新基于材料的物質域。網格域和物質域的變量映射關系可由對流方程[11]表示

(3)

(4)

式中：C1、C2和K1均為完全耦合ALE方法中與對流項相關的矩陣，具體形式見文獻[11]，右上標m表示物質，g表示網格，其余各矩陣的含義均與式(2)相同。式(4)未考慮黏滯阻尼，如考慮可參考式(2)引入阻尼項。

式(4)中存在物質點和網格點兩組未知量，無法得到系統的唯一解，因此，需要引入反映物質點和網格點運動關系的補充方程，通常表達如下[12]：

ug=α+Bum

(5)

式中α和B為反映物質點和網格點相對運動的矢量和矩陣，在耦合的ALE方法中，可依據事先指定的網格運動方法和邊界限制條件確定其具體形式[12]。當α=0，B=I時，ALE方法退化為更新拉格朗日方法；當α=0，B=0時，ALE方法退化為歐拉方法。

耦合的ALE方法需事先指定物質點和網格點的運動關系，而預先指定的運動關系往往比較簡單，如超限映射法等。對于結構形式復雜或變形較大的情況，相對簡單的物質點和網格點運動關系可能限制ALE方法解決大變形問題的能力。因此在實踐中，解耦合的ALE方法更為普遍。該方法基于算子分離技術將一個完整的ALE分析步分解為：1) 一個忽略網格和物質相對運動的拉格朗日分析步，即傳統更新拉格朗日方法；2) 跟隨一個將場變量在物質域和網格域之間傳遞的歐拉分析步。在拉格朗日分析步，網格隨物質點變形；在歐拉分析步利用網格優化算法對計算域變形網格進行調整，并借助變量傳遞算法將場變量從物質域傳遞到調整后的網格域，詳細實現過程可參考文獻[13]。 本文基于有限變形理論和既有更新拉格朗日程序，利用算子分離技術，在每一個拉格朗日分析步后增加一個歐拉分析步，將其擴展為解耦的ALE分析步，發展了一種解耦合ALE有限元程序[13]并將其用于重力式沉箱碼頭的動力液化及減災措施分析。

2 數值模型

2.1 模型簡介

如圖3所示，沉箱模型位于深厚海底黏土層上，沉箱底部用砂土置換海底黏土以提高地基承載力，墻后回填砂土以滿足碼頭正常工作。模型底部完全固定，兩側邊界限制水平位移，頂部為可自由排水邊界(粗虛線)，底部和兩側邊界為不排水邊界。此外，在模型兩側設置寬高比較大(>100)的單元，并綁定兩側邊界單元上同一水平高度處節點的水平自由度模擬自由場邊界[14]，相關研究[15-16]也表明這種等位移邊界條件可以消除邊界效應的影響，在砂土動力液化分析中具有良好的效果。

圖3 沉箱結構及其有限元網格

采用Oka等[17]基于非線性移動硬化準則提出的一種循環彈塑性模型模擬飽和砂土。該模型通過引入超固結邊界面模擬超固結狀態下砂土的剪脹特性，同時考慮了初始各向異性的衰退特性從而能夠描述各項異性固結。模型中各部分土體的參數如表1所示，剪脹參數(D0，n)以調整抗液化強度線的斜率，參考應變(γrefp，γrefe)用來控制應力比超過相轉換線后的應變快速增長。由于黏土具有黏聚力且不會因液化迅速軟化，將黏土參數中的參考剪應變γrefp和γrefe均設定為一個較大的值，使得黏土剛度隨著剪應變的累積而逐漸衰減，但又不像砂土一樣在達到相變比后急劇減小[14]。此外，對于液化后的砂土，考慮本構模型中剪切模量的應力依賴，引入相應的液化機制，使液化土體表現出類似于牛頓黏性流體的特征。為反映地震作用下最不利的狀況，置換砂和回填土均為松砂，沉箱結構賦予較大的彈性模量和泊松比模擬剛體。

表1 土體參數

2.2 輸入地震動

Taiyab等[2]研究表明，相對于高頻荷載，重力式碼頭在低頻荷載作用下產生更大的位移。選擇1999年土耳其地震記錄的Kocaeli波作為輸入荷載，其加速度時程曲線和傅里葉譜如圖4所示。該波的峰值加速度為0.349g，作為強震輸入時將被調整至0.5g，持時25 s，記錄間隔為0.01 s。圖4(b)顯示，Kocaeli波的主頻小于1 Hz，可以使沉箱碼頭產生較大的動力響應。

圖4 Kocaeli波及其傅里葉譜

3 結果分析

3.1 UL方法和ALE方法結果對比

ALE方法能完整反映結構變形并保持網格的整體質量。圖5為使用UL方法和ALE方法計算時，沉箱結構的網格變形，其中，使用UL方法計算時，墻后土體變形較大引起墻后局部網格畸變，導致計算在t=17.9 s時中斷。而ALE方法顯示，計算完成時整體網格仍處于健康狀態。圖5中兩種方法的變形網格在沉箱趾部前端區域有明顯的隆起，沉箱墻體后背有明顯的沉降。ALE方法計算結束時的位移矢量圖顯示了趾部前端隆起和墻后沉降的原因。如圖6所示，由于置換區和回填區砂土液化，沉箱結構沿著置換區和海底黏土交界面形成弧形的滑動帶，最終土體變形在沉箱趾部前端堆積，導致趾部前端隆起，墻體后方沉降。

圖5 兩種方法變形網格圖

圖6 計算結束時位移矢量圖(ALE)

ALE方法能反映最終的結構位移。圖7為沉箱頂點(圖3A點)使用UL方法和ALE方法的水平和豎向位移時程曲線。可以看出，在地震波加載的前10 s, 兩種方法計算的沉箱頂點位移基本保持一致。此后，兩種方法的計算結果出現偏差，水平位移比豎向位移偏差較大，直到UL方法由于網格畸變中斷。UL方法計算中斷時，由于地震波振動仍在持續，沉箱頂點位移無收斂趨勢，無法正確反映最終的位移響應，而ALE方法計算結束時的沉箱頂點水平位移達4.1 m，比UL方法中斷時大11%。

圖7 兩種方法對應沉箱頂點位移

ALE方法能提供相對誤差較小的數值解。為表征兩種方法的精度，引入基于單元應變的平均相對誤差[18]評估計算過程中的整體誤差。圖8為兩種方法的平均相對誤差曲線。兩種方法計算的平均相對誤差的發展趨勢與地震波的波動情況基本一致。在加載初期，荷載和變形很小，相對誤差的波動較小。在3～6 s時經歷地震荷載的峰值，總體的相對誤差不大，導致平均相對誤差經歷較大的波動，隨后出現下降。在8～10 s時，地震波出現一次較大波動，但由于此時總體誤差較大，平均相對誤差出現小幅增大，但未出現加載前期的劇烈波動。此外，ALE方法的平均相對誤差除在加載初期和UL方法相比差別不大之外，總體上遠低于UL方法的平均相對誤差。

圖8 兩種方法對應平均相對誤差

3.2 液化減災措施評價

在回填區和置換區為松砂的情況，強震作用下的沉箱趾部前端隆起嚴重，碼頭墻后方有明顯沉降，沉箱頂點位移達4 m，顯然不能保證經歷地震后的正常工作。沉箱抗液化減災措施通常有增加置換區和回填區土體的密實度、利用碎石樁加固回填區土體等，本節將利用ALE方法探究置換區和回填區密砂加固在液化減災中對沉箱碼頭位移、變形的影響，所用密砂材料參數見表1。

3.2.1 基底置換砂影響分析

1)基底置換砂對沉箱向海側位移的影響?；字脫Q砂加密能有效減小沉箱的向海側位移，但存在有效加密深度。圖9(a)為沉箱結構墻后回填區不加密、置換區加密不同深度時沉箱頂點的位移曲線。可以看出，置換區加密能有效地減小和改變沉箱的位移響應。在置換區加密深度D=0時，沉箱的位移響應在5 s以后變化劇烈，且增長速度較快，加載結束時位移達4.1 m。當D=5 m時，沉箱的位移響應明顯發生改變，位移增長緩慢，加載結束時位移為1.8 m。當D增大至10 m時，沉箱的位移增長進一步放緩，加載結束時位移減小至1.26 m。D=15和 20 m時，位移增長與D=10 m時差別很小。以上分析表明，當D≥10 m時，置換區加密深度對于沉箱位移發展的影響很小。圖9(b)為沉箱墻后存在25 m寬的回填加密區時，置換區不同加密深度對應的沉箱頂點位移?？梢钥闯觯靥罴用軈^的存在不影響置換區加密深度增加對沉箱位移影響的規律，在D≥10 m時，置換加密區的增大對位移發展影響很小。圖9表明，置換區加密能顯著改變沉箱位移的發展趨勢，但存在有效加密深度，超過此深度，置換區加密深度的增加對沉箱位移的減弱作用影響很小。

圖9 沉箱頂點位移時程曲線

圖10為沉箱頂點的殘余位移?？梢钥闯觯贒≥10 m時，置換加密深度增加對殘余位移的影響很小，此時墻后回填加密區的存在對于殘余位移的影響很小。當D<10 m時，置換加密深度增加使沉箱頂點殘余位移迅速減小，此時墻后回填區加密能進一步減小殘余位移。以上表明，在置換加密深度達到臨界深度前，回填區加密能增強沉箱碼頭抗震液化變形能力，當置換加密區深度超過臨界值，墻后回填區加密對于沉箱碼頭的位移增長無明顯抑制作用。

圖10 沉箱頂點殘余位移

2)基底置換砂對基礎隆起的影響。基底置換區加密對沉箱趾部前端隆起有明顯抑制作用。圖11為基底置換砂加密條件下，加載結束時網格的最終變形。圖11(a)顯示加載結束，置換區不加密或加密深度不足時，沉箱趾部前端均有明顯的隆起現象，最大隆起量達1.7和1.6 m，但基底加密的情況下沉箱的位移顯著小于不加密情況。圖11(b)表明,置換區深度增加至10 m時，沉箱趾部前端的隆起現象得到了明顯的抑制，并且沉箱的位移進一步減小。圖11(c)和(d)顯示，基底置換區加密深度從10 m增加到20 m，趾部前端最大隆起量僅從0.49 m減小至0.4 m，表明基底置換區加密深度在超過臨界有效深度后(>10 m)，繼續增加對趾部前端隆起和沉箱位移的改善作用不明顯。綜上表明，基底置換區加密至臨界有效深度前，對沉箱趾部前端和沉箱位移有明顯的抑制作用，此后繼續增加置換深度對隆起和位移的改善作用有限。

圖11 基底置換區加密深度對基底隆起的影響

3.2.2 墻后回填土影響分析

1)回填區對沉箱向海側位移的影響。置換區不加密時，回填區加密能減小沉箱結構的位移，但存在有效加密范圍；置換區加密時，回填區加密對沉箱結構的位移影響較小。圖12(a)為基底置換砂土不加密時，回填區不同加密寬度對應的沉箱頂點位移。可以看出，在置換區不加密情況下，回填區加密寬度≤25 m時，回填區加密寬度的增加對沉箱頂點位移影響較大，沉箱位移時程曲線差異明顯，加密寬度越大，位移越小；回填區寬度≥40 m時，沉箱頂點位移時程曲線幾乎重合。這表明沉箱背后回填區也存在臨界加密寬度，超過此寬度，回填區的加密對沉箱位移無影響，而臨界加密寬度可能與荷載類型、沉箱結構高寬比、場地地基條件等因素有關，實際應用中應根據具體工況確定合理的臨界加密寬度。圖12(b)和(c)為基底置換區分別加密5和10 m時，回填區不同加密寬度對應的沉箱頂點位移時程?？梢钥闯?，在置換區加密的情況下，不同置換寬度對應的沉箱頂點位移呈現相似的發展趨勢，僅在數值上有較小的差別，回填區加密寬度的增加不能明顯減小沉箱頂點位移，說明回填區加密對沉箱頂點位移有抑制作用，但不是影響沉箱位移的主要因素。

圖12 沉箱頂點位移時程曲線

圖13為加載結束時，回填區不同加密寬度對應的沉箱頂點殘余位移。可以看出，在置換區加密深度D=0 m時，隨著回填區加密寬度增加，沉箱頂點殘余位移有較大幅度的減小。而當D=5和10 m時，回填區加密寬度增加對沉箱頂點殘余位移減小幅度很小，這也說明回填區的液化變形不是引起沉箱位移的主要原因。此外，可以發現，不同置換區加密深度下，回填區的臨界有效加密寬度基本保持不變，這也說明回填區不是影響沉箱位移的主要原因。

圖13 沉箱頂點殘余位移

2)回填區對基礎隆起的影響?；靥顓^加密對基礎的隆起無影響。圖14為置換區加密深度為5 m時，回填區不同加密寬度對應的網格變形?？梢钥闯觯煌用軐挾葘某料渲翰壳岸说穆∑鹦螒B基本不變，這表明沉箱趾部隆起主要受置換區液化變形影響，與回填區變形無關。此外，觀察沉箱變形可以看出，置換區加密情況下，回填區加密對沉箱結構的變形影響很小。

圖14 回填區加密寬度對基底隆起的影響

3.2.3 碎石樁加固影響分析

在沉箱結構的液化減災措施中，碎石樁加固回填區土體也是一種有效的措施[1，19]。前兩節的分析中沉箱結構的置換區和回填區加密都存在一個有效的臨界加密范圍，超過此范圍，增大加密范圍對沉箱結構位移響應減小作用基本不變。本節考慮在臨界加密范圍基礎上，加密置換區和回填區土體聯合碎石樁加固回填區對沉箱結構抗震性能的影響。

聯合碎石樁加固能進一步增強沉箱結構抗震性能。圖15為聯合碎石樁加固條件下，沉箱頂點位移時程曲線?？梢钥闯?，在加載前10 s左右，兩種措施對應的沉箱頂點位移基本一致；隨著加載進行和液化區的發展，聯合碎石樁加固的條件下，沉箱位移相比無碎石樁的進一步減小，這表明聯合碎石樁加固方案能在置換區和回填區臨界加密范圍的基礎上進一步提高沉箱結構的抗震性能。

圖15 沉箱頂點位移時程曲線

圖16 超孔隙水壓力比云圖

4 結 論

1)松砂地基沉箱碼頭在強震作用下由于液化大變形造成網格畸變而導致UL方法失效，不能完整反映沉箱碼頭的強震響應；本研究開發的ALE方法能夠在保證網格質量的情況下，反映沉箱碼頭的強震響應過程，并提供精度較高的數值結果。

2)在松砂地基沉箱碼頭的抗液化變形措施中，置換區和回填區的加密均存在一個臨界的有效加密范圍，超過此范圍，對于沉箱碼頭的抗液化變形能力影響不大，與文獻[2]的結論一致。而臨界有效范圍通常與荷載類型、結構高寬比、場地地層條件等因素有關，在實際應用中應根據具體工況加以確定。

3)沉箱趾部前方的隆起主要受置換區砂土密實度的影響，而與回填區無關，在臨界范圍內，增加置換區砂土的密實度能顯著改善趾部前方的隆起狀況。在置換區不加密的情況下，回填區加密能降低沉箱結構的位移；而在置換區加密時，回填區加密處理對沉箱結構的位移影響很小，這也表明置換區土體變形是影響沉箱結構位移的主要原因。這與傳統有限元方法結論一致，也表明ALE方法應用于沉箱碼頭動力液化分析的有效性，可以作為解決砂土動力液化大變形問題的一種有效途徑。

4)碎石樁聯合置換區和回填區加密措施能夠進一步減小沉箱結構的位移，提高沉箱結構的抗液化潛能。