劉 順，唐小微，欒一曉

(大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗室，大連 116023)

重力式沉箱碼頭因其耐久性、方便建造而成為港口常見的濱水構(gòu)筑物之一，在神戶港(Japan)、阿爾戈斯托港(Greece)、斯里蘭卡港 (SriLanka)等重要國際港口都有應(yīng)用，我國青島港、廈門港、花蓮港(臺灣)等也有使用重力式碼頭。港口往往濱海或濱河建造，地基軟弱，為滿足承載力和穩(wěn)定性要求，需采取地基處理措施以支撐上部結(jié)構(gòu)物和工作載荷。圖1 為重力式碼頭的典型斷面形式，重力式碼頭的基礎(chǔ)和墻后通常使用砂土、碎石土等進(jìn)行處理。

圖1 重力式碼頭墻的典型斷面Fig. 1 Typical section of gravity wharf

飽和砂土地基在地震發(fā)生時可能發(fā)生液化，進(jìn)而引發(fā)結(jié)構(gòu)破壞，因此飽和砂土液化及減防災(zāi)措施研究始終是巖土工程領(lǐng)域研究的一個重點(diǎn)問題[1-6]。在過去幾十年，也有很多重力式碼頭墻遭遇地震液化破壞的記錄。1995 年，日本兵庫縣地震，神戶港碼頭的大量重力式沉箱結(jié)構(gòu)遭到破壞，整個港口陷入癱瘓。此外，中國臺灣Chi-Chi 地震、土耳其Kocaeli 地震都觀測到重力式碼頭液化破壞的現(xiàn)象。破壞原因往往是基底置換砂和墻后回填土液化引起的碼頭墻沉降、向海側(cè)位移、傾斜等。

Abu 等[7]研究了地震作用下重力式碼頭土和結(jié)構(gòu)相互作用，并分析了置換區(qū)和回填區(qū)密砂土實(shí)度對重力式岸壁結(jié)構(gòu)的影響，Alyami 等[8]探究了滲透系數(shù)對重力式碼頭結(jié)構(gòu)孔隙水壓力的影響以及回填土和地基土的相對密實(shí)度對岸壁殘余變形的影響。王麗艷等[9]詳細(xì)分析了液化地基中沉箱碼頭殘余變形隨各影響因素的變化規(guī)律，Tong 和Schaefer[10]研究了振沖加密優(yōu)化設(shè)計對重力式碼頭墻液化災(zāi)害的減弱作用。

上述研究均是基于有限元方法進(jìn)行的數(shù)值模擬，作為一種有效的工具，有限元方法在工程領(lǐng)域有重要的應(yīng)用。但在一些變形比較大的情況下，可能由于網(wǎng)格畸變導(dǎo)致計算中斷、精度喪失等。地震調(diào)查報告顯示，兵庫縣地震中神戶港區(qū)地面最大加速度達(dá)到0.6g，部分沉箱頂部水平殘余位移達(dá)5 m，殘余沉降達(dá)1 m～2 m，向海側(cè)傾角3°～5°。在強(qiáng)震作用下，基于傳統(tǒng)有限元方法的沉箱碼頭動力分析可能由于較大的變形引起分析中斷或者結(jié)果失真，不能合理反映強(qiáng)震作用下沉箱碼頭的最終變形，難以對沉箱碼頭的抗震優(yōu)化設(shè)計起到有效的參考作用?；谒阕臃至鸭夹g(shù)和自主開發(fā)的水土耦合動力液化數(shù)值分析平臺，筆者發(fā)展了一種ALE 有限元方法，并通過飽和土柱一維固結(jié)和離心機(jī)模型試驗?zāi)M進(jìn)行了驗證[11]，詳細(xì)實(shí)現(xiàn)過程可參考文獻(xiàn)[11]，這里將其用于強(qiáng)震作用下的沉箱碼頭地震液化分析。

1 飽和砂土控制方程

基于水土二相混合理論，以土骨架位移u和孔隙水壓力p作為變量，Uzuoka[12]推導(dǎo)了u-p形式的飽和砂土的場方程式：

2 ALE 方法簡介

基于拉格朗日描述的有限元方法可能由于網(wǎng)格畸變導(dǎo)致計算中斷、精度喪失，而基于歐拉描述的有限元方法處理運(yùn)動的物質(zhì)邊界和相互作用問題是十分困難的。結(jié)合兩種方法的優(yōu)點(diǎn)，Donea 等[15]發(fā)展了一種任意的拉格朗日歐拉方法(Arbitrary Lagrangian Eulerian Method) 用于解決流體力學(xué)問題，隨后該方法被眾多學(xué)者用于解決固體力學(xué)中的非線性問題。

圖2 物質(zhì)域、網(wǎng)格域、空間域之間的映射關(guān)系Fig. 2 Mapping relationship among material domain, grid domain and spatial domain

根據(jù)ALE 方法的實(shí)現(xiàn)途徑，可以將其分為兩類：1) 完全耦合的(Fully Coupled)；2) 解耦合的(Decoupled)。完全耦合的ALE 方法將式(8)代入基于計算域的控制方程進(jìn)行線性化，將物質(zhì)運(yùn)動和網(wǎng)格運(yùn)動整合于離散化的有限元方程中，其基本形式如式(9)[16]所示：

式(10)中，耦合的ALE 方法需事先指定物質(zhì)點(diǎn)和網(wǎng)格點(diǎn)的運(yùn)動關(guān)系，這種關(guān)系通常比較簡單，如超限映射法[17]等，對于結(jié)構(gòu)形式復(fù)雜或變形較大的情況，往往不能作為有效的解決方案。在實(shí)踐中，解耦合的ALE 方法更為普遍，該方法基于算子分離技術(shù)將一個完整的ALE 分析步分解為[18]：1) 忽略網(wǎng)格和物質(zhì)相對運(yùn)動的拉格朗日分析步，即更新拉格朗日方法；2) 將場變量在物質(zhì)域和網(wǎng)格域之間傳遞的歐拉分析步。本文發(fā)展的ALE 有限元方法是在課題組開發(fā)的更新拉格朗日格式的水土耦合地震液化分析平臺基礎(chǔ)上，基于算子分離技術(shù)的解耦合ALE 有限元程序[11]。

3 有限元模型

3.1 模型簡介

本文考慮了邵琪[19]在飽和砂土地震液化網(wǎng)格自適應(yīng)的研究中采用的兩種重力式碼頭地基處理形式。如圖3 所示，沉箱位于海底黏土之上，沉箱底部用砂土置換，沉箱后部以砂土回填，兩種方案不同之處在于沉箱地基的置換范圍和形式。兩種方案中，模型兩側(cè)邊界設(shè)置寬高比較大(＞100)的單元和等位移邊界(EDB)模擬自由場邊界[11]，兩側(cè)邊界和底邊界為不透水邊界，模型頂部(粗虛線)可自由排水。各區(qū)域的尺寸見文獻(xiàn)[19]。

圖3 沉箱模型及有限元網(wǎng)格Fig. 3 Caisson model and Finite element meshes

3.2 輸入地震動

文中選擇Kobe (Japan，1995，PA=0.834g)，Kocaeli (Turkey，1999，PA=0.349g)，Northridge(USA，1994，PA=0.568g) 三種地震波作為輸入荷載，如圖4 所示。圖中，左側(cè)為地震波，右側(cè)為對應(yīng)傅里葉譜，其中Kocaeli 波的主要頻率范圍低于1 Hz，Kobe 和Northridge 波的主要頻率范圍在1 Hz～10 Hz?？紤]豎向地震動的影響，取水平向輸入荷載的1/2 作為豎向地震動輸入。

圖4 地震波加速度及傅里葉譜Fig. 4 Acceleration and Fourier spectra of earthquake waves

圖5 為近幾十年來，記錄的一些地震的峰值加速度分布，可以看出，強(qiáng)震動峰值加速度主要集中在0.5g左右。因此，本文對于強(qiáng)震動作用下沉箱碼頭的地震液化分析，所采用的三種地震波時程將按照其峰值強(qiáng)度與0.5g的比例進(jìn)行調(diào)整到0.5g。

圖5 近幾十年地震最大加速度[22]Fig. 5 Peak acceleration of earthquakes in recent decades[22]

表1 模型參數(shù)Table 1 Model parameters

4 數(shù)值結(jié)果分析

4.1 網(wǎng)格形狀

圖6 和圖7 顯示峰值加速度為0.5g時，兩種方案在三種地震荷載作用下，最終的網(wǎng)格變形，左側(cè)為ALE 方法的網(wǎng)格，右側(cè)為UL 方法的網(wǎng)格。

圖6 方案1 網(wǎng)格變形圖(Peak=0.5 g)Fig. 6 Mesh deformation of Case 1 (Peak=0.5 g)

圖6 顯示，方案1 中使用UL 方法的沉箱趾部前端在Kobe 波和Kocaeli 波作用下發(fā)生大范圍的網(wǎng)格畸變，沉箱墻體后方的回填區(qū)在Northridge 波作用下存在局部的網(wǎng)格畸變；ALE 方法中網(wǎng)格在三種地震荷載作用下均處于健康的狀態(tài)。

圖7 顯示，方案2 中在三種地震波作用下，UL 方法在沉箱趾部前端局部(Kobe)和較大范圍(Northridge)以及墻體后部(Kocaeli)產(chǎn)生網(wǎng)格畸變，導(dǎo)致該方法失效；ALE 方法僅在沉箱趾部發(fā)生局部網(wǎng)格變形，整體網(wǎng)格處于健康狀態(tài)。需要說明的是，在Kocaeli 波作用下ALE 方法失效是由于方案2 中沉箱發(fā)生較大的傾斜，置換區(qū)產(chǎn)生較大隆起變形，沉箱趾部侵入置換區(qū)，導(dǎo)致局部網(wǎng)格由凸變凹，同時也反映出此種處理方案下，沉箱在強(qiáng)震作用下的變形特點(diǎn)。

圖7 方案2 網(wǎng)格變形圖(Peak=0.5 g)Fig. 7 Mesh deformation of Case 2 (Peak=0.5 g)

綜合圖6 和圖7，ALE 方法能夠在強(qiáng)震作用下保持沉箱結(jié)構(gòu)網(wǎng)格的整體質(zhì)量，從而保持計算的持續(xù)性和準(zhǔn)確性，反映沉箱結(jié)構(gòu)的變形特點(diǎn)；UL 方法在強(qiáng)震作用下可能由于結(jié)構(gòu)的大變形引起網(wǎng)格畸變導(dǎo)致數(shù)值方法失效。

4.2 殘余位移

圖8 為三種地震荷載作用下，方案1 沉箱頂點(diǎn)(圖3(a)A點(diǎn))在不同峰值加速度下的水平位移?？梢钥闯觯诜逯导铀俣萈A=0.1g時，加載初期變形較小，兩種方法的沉箱頂點(diǎn)位移基本一致，隨著加載進(jìn)行，UL 方法的頂點(diǎn)殘余位移稍大于ALE 方法。強(qiáng)震作用時(PA＞0.3g)，兩種方法對應(yīng)的沉箱頂點(diǎn)位移隨著加載進(jìn)行差異明顯，并且峰值加速度越大差異越大。圖8(a)、圖8(b)顯示，在強(qiáng)震作用下，由于沉箱變形太大引起網(wǎng)格畸變導(dǎo)致UL 方法失效而不能完整反映沉箱位移變化，而ALE 方法則完整反映了沉箱頂點(diǎn)位移變化，并提供了最終的沉箱殘余位移。圖8(c)中，在Northridge 波作用下，UL 方法和ALE 方法均能反映沉箱頂點(diǎn)位移變化過程，但二者的差異隨峰值加速度增大而變大。需要指出的是，強(qiáng)震時ALE 方法對應(yīng)的沉箱頂點(diǎn)最終位移高于UL 方法中斷時對應(yīng)的位移。如表2 所示，在Kobe 和Kocaeli 波作用下，二者的位移差值最大達(dá)16%。

表2 方案1 沉箱頂點(diǎn)位移Table 2 Displacement of caisson apex of Case 1

圖8 沉箱頂點(diǎn)水平位移時程(方案1)Fig. 8 Horizontal displacements of caisson apex of Case 1

圖9 為三種地震荷載下，方案1 沉箱頂點(diǎn)在不同加速度峰值下ALE 方法的殘余水平位移(圖9(a))和豎向位移(圖9(b))。從圖中可以看出，沉箱的殘余位移與地震荷載峰值加速度近似呈線性關(guān)系，與文獻(xiàn)[9]中殘余位移與輸入荷載相關(guān)性一致。Abu 等[7]指出低頻荷載對應(yīng)的沉箱碼頭向海側(cè)位移相對更大?？梢园l(fā)現(xiàn)，圖9 中Kobe 波和Northridge 波對應(yīng)的殘余位移增長率相似，且數(shù)值相近，而Kocaeli 波作用下，殘余位移隨峰值加速度的增長率高于其他兩種地震波，且強(qiáng)震時相同峰值加速度對應(yīng)的殘余位移遠(yuǎn)高于其他兩中地震波對應(yīng)的殘余位移?？紤]到輸入荷載的頻譜特性，沉箱碼頭的地震響應(yīng)受低頻地震波影響更大，這與文獻(xiàn)[7]的結(jié)論一致。

圖9 沉箱頂點(diǎn)殘余變形(方案1)Fig. 9 Residual displacements of caisson apex of Case 1

圖10 所示為三種地震荷載作用下，方案2 中沉箱頂點(diǎn)(圖3(b)B點(diǎn))在不同峰值加速度下的水平位移。在小震(PA=0.1g)和強(qiáng)震(PA＞0.3g)作用下，方案2 沉箱頂點(diǎn)位移均勻與方案1 的沉箱頂點(diǎn)位移發(fā)展趨勢相似。不同之處在于，方案2在強(qiáng)震(PA＞0.3g)作用下，UL 方法失效時間比方案1 有所提前，并且ALE 方法在一些情況下(Kobe波PA=0.5g和Kocaeli 波PA=0.4g、0.5g)也發(fā)生失效。不同于UL 方法導(dǎo)致的大范圍網(wǎng)格畸變，ALE 方法失效的共同特點(diǎn)是由于沉箱向海側(cè)傾斜，侵入置換砂地基導(dǎo)致沉箱趾部某些網(wǎng)格失效(圖7)，并未影響網(wǎng)格的整體質(zhì)量。此種情況下，除更換整體網(wǎng)格或計算方法外并無有效解決途徑。此外ALE 方法失效時，均已超過了輸入地震波的主震范圍，震動基本結(jié)束，且殘余位移水平已經(jīng)遠(yuǎn)超沉箱碼頭正常工作所允許的范圍，此時的殘余位移對于沉箱碼頭抗液化措施有效性評價仍具有參考意義。表3 為在峰值加速度為0.4g和0.5g的Kobe 和Northridge 波加載下，UL 方法和ALE 方法在計算中止時對應(yīng)的沉箱頂點(diǎn)位移，可以看出在方案2 情況下，ALE 方法的位移遠(yuǎn)大于UL 方法，最大差值達(dá)56%。

圖10 沉箱頂點(diǎn)水平位移(方案2)Fig. 10 Horizontal displacements of caisson apex of Case 2

表3 方案2 沉箱頂點(diǎn)位移Table 3 Displacement of caisson apex of Case 2

圖11 為三種地震荷載作用下，方案2 中沉箱頂點(diǎn)在不同峰值加速度下使用ALE 方法的殘余水平位移(圖11(a))和豎向位移(圖11(b))。與方案1中結(jié)果相似，沉箱殘余位移與地震波峰值加速度近似呈線性關(guān)系，并且Kocaeli 波的殘余位移大于Kobe 波和Northridge 波，反映出沉箱地震響應(yīng)的低頻敏感性。文獻(xiàn)[7]中顯示地基處理方案對將沉箱碼頭系統(tǒng)的模態(tài)頻率有較大影響，具體表現(xiàn)為隨著置換區(qū)地基區(qū)域的變大，模態(tài)頻率逐漸變大；沉箱碼頭墻底部置換區(qū)加固對模態(tài)頻率的影響高于沉箱碼頭墻后回填區(qū)。由圖9 和圖11 可以看出，在相同峰值加速度下，方案2 中Kobe 波和Northridge 波對應(yīng)的殘余位移大于方案1，這是由于地基處理方案的不同，改變了沉箱碼頭的自然頻率，使得沉箱碼頭對Kobe 波和Northridge 波的地震響應(yīng)變大，也減小了Kocaeli 波的地震響應(yīng)與Kobe 波Northridge 波地震響應(yīng)的差距。

圖11 沉箱頂點(diǎn)殘余變形(方案2)Fig. 11 Residual displacements of caisson apex of Case 2

4.3 相對誤差分析

在地震液化的數(shù)值分析中，土體和結(jié)構(gòu)由于液化引起的變形是關(guān)注的重點(diǎn)，因此，本文將單元應(yīng)變作為誤差評估對象來評估UL 方法和ALE方法在地震液化計算中對計算結(jié)果的影響?；趹?yīng)變的平均相對誤差定義形式可參考文獻(xiàn)[19]。

圖12 顯示，方案1 中，在小震(PA=0.1g)作用時，除初期兩種方法的相對誤差相差不大以外，UL 方法平均相對誤差隨著變形增大而持續(xù)增大，ALE 方法對應(yīng)的平均相對誤差處于很低的水平。強(qiáng)震(PA＞0.3g)作用下，計算初期兩種方法的平均相對誤差亦差別不大。隨著變形的增大，UL 方法對應(yīng)的平均相對誤差持續(xù)增大，在計算失效時，相對誤差出現(xiàn)陡升，對應(yīng)圖6 中大范圍的網(wǎng)格畸變，而ALE 方法則保持平均相對誤差始終處于較低的水平，直至計算結(jié)束。

圖13 顯示，方案2 中，小震作用下，UL 方法和ALE 方法的相對誤差發(fā)展趨勢與方案1 相似。強(qiáng)震(PA＞0.3g)作用下，UL 方法計算失效，對應(yīng)了圖7 中部分網(wǎng)格畸變，ALE 方法在Kobe 波和Kocaeli 波作用下盡管比UL 方法計算持時較長，也未能完成整個計算過程。如第4.2 節(jié)所述，ALE 方法失效但并未影響整體網(wǎng)格質(zhì)量，從圖13可以看出UL 方法失效時，由于較多網(wǎng)格畸變，平均相對誤差出現(xiàn)陡升，而ALE 方法雖然失效，平均相對誤差始終處于很低的水平，且未出現(xiàn)由于大范圍網(wǎng)格失效導(dǎo)致的誤差陡增。

綜合圖12 和圖13，ALE 方法用于沉箱碼頭的地震液化分析，在保證網(wǎng)格質(zhì)量的基礎(chǔ)上，提供了相比UL 方法精度更高的結(jié)果。

圖12 平均相對誤差曲線 (方案1)Fig. 12 Average relative error curves of Case 1

圖13 平均相對誤差曲線 (方案2)Fig. 13 Average relative error curves of Case 2

4.4 超孔隙水壓力

為表征飽和砂土在地震作用下的孔壓水平，這里定義超孔隙水壓力比(excess pore water pressure ratio)為EPWPR=Δu/σ′y0，其中， Δu為單元超孔隙水壓力， σ′y0為單元當(dāng)前豎向有效應(yīng)力。

圖14 為方案1 在不同峰值的Kocaeli 波(0.3g～0.5g)作用下，ALE 方法中E1(圖3(a))的超孔壓比曲線和云圖。從圖14(a) 可以看出在不同峰值作用下，E1 的超孔壓水平都較低，未達(dá)到液化水平。整體來看，孔壓發(fā)展在該地震波激勵下波動較為劇烈，0.5g對應(yīng)的孔壓曲線加載后期與0.3g和0.4g相差較大，達(dá)到較高的水平。從圖14 超孔壓云圖可以看出，隨著峰值加速度的增大，沉箱底部雖未達(dá)到液化，但高超孔壓區(qū)逐漸擴(kuò)展。回填區(qū)液化主要發(fā)生在距離墻后一定距離的區(qū)域，緊鄰碼頭墻后方區(qū)域未發(fā)生液化，這種現(xiàn)象與實(shí)際地震中墻后液化現(xiàn)象一致[23]。隨著峰值加速度增大，回填區(qū)和置換區(qū)的高孔壓區(qū)逐漸產(chǎn)生連通，這將加大沉箱的殘余變形，因此，設(shè)計中可考慮設(shè)置阻斷置換區(qū)和回填區(qū)孔壓轉(zhuǎn)移的措施，如振沖加密置換回填區(qū)土體[10]或設(shè)置水平止水層。

圖15 為方案2 在不同峰值的Northridge 波(0.3g～0.5g)作用下，ALE 方法中E1(圖3(b))的超孔壓比曲線和云圖。方案2 中E1 孔壓隨加載進(jìn)行逐漸累積，且峰值加速度增大，相應(yīng)超孔隙水壓力越大。圖15 云圖顯示在方案2 的沉箱下部，有較大范圍的液化區(qū)，隨著峰值加速度的增大，液化區(qū)范圍逐漸增大，并且有與圖14 類似置換區(qū)和回填區(qū)高孔壓區(qū)貫通擴(kuò)展的現(xiàn)象。方案1 沉箱底部未出現(xiàn)明顯液化區(qū)，而方案2 沉箱底部存在較大范圍的液化區(qū)，說明方案2 的形式在沉箱底部易產(chǎn)生孔壓積聚。

圖15 超孔隙水壓力曲線和云圖(方案2)Fig. 15 Curves and contours of EPWPR (Case 2)

5 結(jié)論

基于飽和砂土控制方程和算子分裂技術(shù)，發(fā)展了一種解耦和的ALE 有限元方法，并將其用于飽和砂土地基沉箱碼頭地震液化分析。使用UL 方法和ALE 方法對比了兩種沉箱地基處理方案在三種不同地震波激勵下的地震響應(yīng)，主要結(jié)論如下：

(1) 小震作用下，UL 方法和ALE 方法均能完成分析，并且UL 的計算結(jié)果略大于ALE 計算結(jié)果；強(qiáng)震作用下，UL 方法失效引起計算中止時的沉箱頂點(diǎn)位移與ALE 方法完成計算時相差較大。因此，文中發(fā)展的ALE 方法可用于對于UL方法不能夠提供完整的沉箱碼頭地震響應(yīng)的情況。此外，相對誤差分析結(jié)果顯示，無論UL 方法是否能完成計算，ALE 方法的平均相對誤差均低于UL 方法，表明ALE 方法能夠提供精度可靠的數(shù)值解。

(2) 沉箱碼頭的殘余位移與地震加速度峰值近似呈線性關(guān)系，且對低頻的荷載較為敏感，相應(yīng)的地震響應(yīng)也較大。地基處理方案不同將改變結(jié)構(gòu)的自然頻率，進(jìn)而影響沉箱碼頭對地震荷載的敏感頻率范圍。因此，在沉箱碼頭的設(shè)計中，應(yīng)考慮地基的處理形式對自然頻率的影響，進(jìn)而影響到沉箱碼頭的地震響應(yīng)。

(3) 置換區(qū)和回填區(qū)在地震荷載作用下，可能存在高孔壓連通區(qū)，并且地震峰值加速度越大，連通區(qū)范圍越大。連通區(qū)將增大沉箱碼頭地震響應(yīng)的殘余位移，因此，在設(shè)計中，可考慮設(shè)置阻斷置換區(qū)和回填區(qū)孔壓轉(zhuǎn)移的措施。