基于邊界面模型的砂土中板錨循環承載力數值分析

王 強，劉海笑，李 洲

(天津大學 建筑工程學院，天津 300072)

隨著海洋資源開發走向深水化，浮式油氣開發平臺、浮式波浪發電站和海上風力發電站等各類浮式海洋工程構筑物大量興建。這些結構除了承受正常生產活動所產生的工作荷載外，還會受到風、浪、流等引起的循環荷載，且與地震荷載不同，這種循環荷載大多具有長期、持續和低頻的特征，因此對系泊系統的要求極為嚴苛。錨固結構作為系泊系統中提供抗拔承載力的重要組成部分，是限制浮體漂移的關鍵所在。板錨是一種重要的嵌入式海洋工程錨固結構形式，盡管其體積小、重量輕，但仍然能夠提供很高的承載能力，且對系纜角度的要求不高，同時具有經濟實用、施工簡單等優點[1]，因此板錨在系泊系統中被廣泛應用。

海洋土的非線性動力特性，決定了嵌入其中的結構在循環荷載作用下會產生復雜的荷載—位移響應，也對結構的承載力設計和維護帶來諸多困難。現行設計規范如挪威船級社制定的DNV規范[2]和美國船級社制定的ABS規范[3]，在針對海洋結構基礎的循環承載性能評估中，都建議通過繪制土體循環強度輪廓圖(the cyclic contour diagram)的方法來考慮循環荷載效應。該方法需要考慮結構不同位置處土單元的受力狀態，通過一組單調三軸壓縮，拉伸和直接剪切試驗，以及系列循環試驗來獲得完整的循環強度輪廓圖，進而利用插值法獲得各土單元的循環強度[4]。

目前針對嵌入式海洋結構循環承載性能的研究常局限于黏土[5-6]，由于砂土更為復雜的應力—應變響應，針對板錨循環動力特性的研究還處于初步階段。Bemben等[7-8]對飽和松砂中的板錨進行了系列循環試驗，試驗中觀察到了明顯的“棘輪效應”。Petereit[9]針對循環荷載幅值對板錨循環承載特性的影響進行了系列模型試驗，研究發現循環幅值越高，板錨位移累積速率越大。Chow等[1，10]針對板錨在密實砂土中的行為特性做了大量試驗研究，考察了不同性質循環荷載作用下板錨的荷載響應特性。由于板錨在長期循環荷載作用下承載性能的變化，用試驗模擬耗時耗力，并且隨著深水化發展，現場試驗的開展難度越來越大。因此數值分析模型成為研究循環荷載作用下板錨承載特性的重要選擇。

現有海洋工程構筑物循環承載數值分析模型主要分為擬靜力分析模型[11-12]、顯式經驗模型[13-14]以及動力本構數值分析模型[15-17]。其中，動力本構數值分析模型基于能夠合理反映砂土動力響應特性的本構模型理論，是評估嵌入式海洋工程結構循環承載特性的重要手段，也受到了越來越多研究者的關注。現有動力本構數值模型具有模擬板錨循環動力響應的能力，但所采用的本構理論自身具有局限性。如基于廣義應力、應變建立的宏元模型[1，18]無法反映土體結構自身的動力響應，也無法考慮結構特征改變對承載力的影響。以及基于具有高頻、瞬時特點的地震荷載建立的本構模型[15-16]無法模擬在低頻、長期循環荷載作用下結構在砂土中的動力特性。目前，可以合理反映此類循環承載響應特性的本構模型仍然處于理論研究階段，沒有成熟的數值模型可用于嵌入式海洋工程結構的分析模擬。因此，一種能揭示長期循環荷載作用下板錨—海床土相互作用機理，進而可用以評估板錨循環承載力的數值分析模型，亟待研究、建立。

鑒于以上原因，基于一種適用于飽和砂土循環動力分析的邊界面模型[19]，通過修正模型塑性模量和硬化參數，使其能夠更加準確地描述循環荷載作用下土體的塑性變形。利用帶誤差控制的顯式積分算法，以張量形式，將上述邊界面塑性模型編寫成可供有限元軟件調用的用戶自定義材料子程序。建立土體單元數值模型以模擬砂土的排水靜、動三軸試驗，考察模型合理描述砂土在不同荷載條件下力學響應的能力。最后建立飽和砂土中板錨循環承載分析的數值模型，針對板錨在砂土中的單調抗拔特性和循環承載特性進行數值分析，并考察循環荷載要素對板錨循環承載特性的影響。

1 適用于飽和砂土動力分析的邊界面塑性模型

董建勛等[19]和高源等[20]以相應于歷史所受最大應力的最大加載面作為邊界面，以此來記憶應力歷史。初始加載時，應力點始終位于邊界面上，邊界面隨著應力的增加而擴大；卸載時，應力點向內移動，最大應力保持不變，邊界面固定；再加載時，在應力超過歷史最大應力之前，邊界面一直固定，應力狀態達到最大應力狀態時，應力路徑與邊界面相交，在后續加載中，邊界面以初始加載的形式擴大，更新最大應力歷史。

采用由Bardet[21]提出的修正橢圓邊界面，邊界面圖形如圖1所示，邊界面方程如下：

圖1 邊界面塑性模型的應力更新流程Fig. 1 Stress update flow chart of boundary surface plasticity model

(1)

(2)

(3)

(4)

式中：J2、J3分別為偏應力第二、三不變量；c=Me/Mc，Me為三軸拉伸時的臨界狀態應力比，表達式為：

(5)

(6)

式中：i表示第i加載步；p0為初始圍壓；(opi，oSi)為第i步的映射中心坐標。

砂土的體積變化同時取決于所受荷載和材料屬性，為了合理反映砂土受力后的體積變化特性，采用非相關聯的流動法則，剪脹函數表達式如下：

(7)

式中：d0為模型常數；rij=sij/p為偏應力比張量；Md為剪脹應力比；nij為塑性偏應變方向。其中，

(8)

式中：m為模型參數。ψ為表征砂土屬性的狀態參數，表達式為：

ψ=e-ec

(9)

式中：e和ec分別為當前孔隙比和臨界狀態孔隙比，ec與當前平均有效應力p有關。ec采用Li等[23]提出的表達式：

(10)

式中：pat為標準大氣壓；eΓ、λ、ξ為土體固結壓縮曲線的擬合參數，描述土體的臨界狀態。

nij表征塑性偏應變方向，表達式為：

(11)

針對砂土的循環承載，土體結構的變化對土體剛度有顯著影響，該變化主要分為兩類，一類是循環荷載作用下土體逐漸密實，土體結構更加穩定，另一種是負塑性體積應變造成的破壞[24]。為了準確反映這兩種變化對土體剛度的影響，采用如下硬化公式：

(12)

式中：

(13)

(14)

式中：

(15)

式中：α、β、χ為控制土體結構變化速率的材料常數，L為加載指數。L表達式為：

(16)

彈性剪切模量G和彈性體積模量K采用如下表達式：

(17)

(18)

式中：G0為模型參數；υ為泊松比。

塑性模量表示式如下：

(19)

(20)

式中：n為模型參數。

應力—應變關系的增量方程可以表示為：

(21)

式中：

(22)

(23)

(24)

利用上述適用于飽和砂土動力分析的邊界面塑性模型得到土體在特定加載路徑下的應力—應變關系，并最終用于巖土工程邊值問題求解，需要采用有效的數值積分算法對模型進行數值實現。采用帶誤差控制的改進Euler積分算法，通過自動調整各子增量步的大小來控制積分誤差[25]。該方法是一種子增量步應力積分算法，通過監控局部截斷誤差來控制子增量步步長，將應變增量{Δε}離散成系列子增量步應變增量{Δεs}，即Δε=ΔT·Δεs。T為無量綱時間，初始時間T=0，ΔT為步長控制參數，值域為(0，1]，初始增量步中ΔT=1。通過監控局部截斷誤差來控制子增量步步長。截斷誤差通過每個子增量步應力增量的誤差值來計算，如果截斷誤差超過指定容許值SSTOL，則使用更小的子增量步進行積分[20]。圖1展示了積分算法的主要流程。

2 三軸試驗數值模型驗證

為了考察建立的本構模型描述砂土力學響應的能力，同時驗證帶誤差控制的顯式積分算法的有效性和準確性，首先將建立土體單元數值模型以模擬完全排水條件下的砂土靜、動三軸試驗。

2.1 靜三軸試驗模擬

Verdugo等[26]對Toyoura砂進行了系列三軸試驗，試驗所用Toyoura砂是一種均勻細砂，平均粒徑0.17 mm，不均勻系數1.7，最大和最小孔隙比分別為0.977和0.597，相對體積質量2.65。針對3組具有不同密度砂土的排水三軸試驗進行模擬，模型參數見表1，模型參數采用董建勛等[19，27]建議的方法進行確定。3組土樣固結壓力均為100 kPa，孔隙比分別為0.996、0917和0.831，試驗加載采用應變控制加載，加載過程中保持側向壓力恒定不變。

表1 模型參數Tab. 1 Model parameters

圖2中，ε1為軸應變，q為偏應力，e為孔隙比，由圖2可以看出，剪切荷載作用下，3組不同孔隙比土樣的應力—應變和應力—孔隙比計算結果與三軸試驗結果相吻合，說明邊界面塑性模型所采用的塑性模量和剪脹函數能有效反映在單調荷載作用下砂土的力學特性。

圖2 單調荷載下Toyoura砂三軸壓縮試驗數據與模型對比Fig. 2 Comparison of Toyoura sand under monotonic loading between data from triaxial compression tests and results from modelling

2.2 動三軸試驗模擬

Pradhan等[28]在循環荷載下針對砂土的應力—應變關系做了一系列排水三軸試驗，試驗所用砂為Toyoura砂。選取其中2組正常固結試驗進行模擬分析，2組試驗土樣初始孔隙比分別為0.832和0.845，平均主應力P′為100 kPa，且在試驗過程中保持恒定。

圖3和圖4分別對比了孔隙比為0.832的Toyoura砂在平均主應力和循環應力幅值都保持不變時的循環荷載作用下循環動力特性以及孔隙比為0.845的Toyoura砂在平均主應力恒定且循環應力幅值不斷增大的循環荷載作用下循環動力特性。從對比模型計算結果與試驗數據可以看出，邊界面模型可以合理反映循環荷載作用下砂土的排水動力響應特性。

圖3 循環荷載下Toyoura砂應力比—應變曲線的循環三軸試驗數據與數值模擬對比Fig. 3 Comparison of stress-strain curve of Toyoura sand under cyclic loading between data from triaxial compression tests and results from modelling

圖4 循環荷載下Toyoura砂循環三軸試驗數據與數值模擬對比Fig. 4 Comparison of Toyoura sand under cyclic loading between data from triaxial compression tests and results from modelling

3 板錨承載力分析

3.1 板錨極限承載力分析

Dickin等[29]針對板錨在砂土中的單調抗拔效應做了系列離心機試驗，考察了不同板錨尺寸和埋深下板錨的荷載—位移響應。試驗采用的是干燥、均勻的Erith石英砂，80%的粒徑在0.125 mm到0.25 mm之間，不均勻系數為1.5，密砂和松砂的容重分別為16 kN/m3和14.5 kN/m3，離心機加速度為40g，試驗所用板錨為25 mm寬、3 mm厚的不銹鋼模型板錨，模型參數見表1。

由于錨和荷載都是對稱的，可利用對稱性直接建立四分之一的三維軸對稱數值模型進行分析，垂直對稱軸通過板錨的中心，選取適當的邊界位置(3B)以消除邊界效應，有限元模型如圖5所示，其中B為板錨寬度，H為板錨埋深。為板錨施加剛體約束，將板錨的運動約束為參考點的運動，給土體底部施加固定約束，在對稱面上施加對稱邊界條件，板錨和砂土之間采用綁定接觸。土體單元以及板錨單元的單元類型為C3D8，網格尺寸對于模型的計算效率和精度有很大的影響，通常網格尺寸越小，單元個數越多，計算精度越高，但計算效率越低。為了提高模型的計算精度，同時又保證計算效率，需要對網格尺寸進行考察，以獲得最優化的網格尺寸。研究考察了3組網格密度，最小網格尺寸分別為0.05B、0.025B、0.0125B，網格數分別為66 781(Mesh1)、247 491(Mesh2)、963 200(Mesh3)，如圖6所示，W是板錨上拔位移，B為板錨寬度，F為四分之一板錨所受抗拔力。結果顯示3組網格密度所計算的結果都已經收斂，但是網格尺寸越小，曲線越光滑，綜合考慮計算精度和計算效率，研究選取最小網格尺寸為0.025B。

圖5 數值模型示意 Fig. 5 Schematic diagram of numerical model

圖6 網格考察結果Fig. 6 Results of grid survey

圖7、圖8展示了在不同埋深下，模型試驗與數值模擬獲得的不同長寬比板錨在松砂和密砂中的荷載—位移響應曲線。由圖可知，隨著板錨長寬比、埋深和土體密度的增大，最大拉拔承載力也隨之增大，數值計算得到的最大抗拔承載力、對應最大抗拔承載力的位移、失效后殘余抗拔承載力均與實測試驗結果一致。

圖7 板錨在砂土中的荷載—位移曲線(L/B=1)Fig. 7 Load-displacement curve of plate anchors in sand (L/B=1)

圖8 板錨在砂土中的荷載—位移曲線(L/B=8)Fig. 8 Load-displacement curve of plate anchors in sand (L/B=8)

針對淺埋工況(H/B=3)，在密砂中，初始階段承載力隨著板錨位移的增加而迅速增大，增長速率隨著接近最大值(峰值)而逐漸減小，出現明顯峰值后，承載力逐漸減小。在松砂中，不同長寬比的板錨所受荷載以不同的速率單調增加到一個最大值，之后板錨繼續上拔，但承載力幾乎不再發生變化。針對深埋工況(H/B=7)，在密砂和松砂中的板錨荷載—位移曲線都未觀察到明顯的峰值。

圖9展示了在不同埋深下板錨上方土體位移的差異，可以看出，淺埋時，土體位移會延伸到土體表面，而深埋時，只在錨上方的局部區域觀察到位移變化。上述現象Ilamparuthi等[30]在圓形板錨的單調抗拔承載力試驗研究中也有報道，通過對砂土中淺埋和深埋板錨的單調抗拔模型試驗，總結出兩類典型破壞模式：淺埋時，從板錨到砂土表面會形成一個向上隆起的土錐臺，即整體剪切破壞；深埋時，板錨上方會形成一個氣球狀的破裂帶，即局部剪切破壞，兩種破壞模式如圖10所示。隨著埋入深度的增加，板錨的破壞行為由整體剪切破壞轉變為局部剪切破壞。

圖9 板錨在松砂中的位移云圖Fig. 9 Displacement contours for a plate anchor in sand

圖10 松砂中板錨的失效機理Fig. 10 Failure mechanism in sand above plate anchor

綜上，建立的適用于飽和砂土中板錨排水循環動力分析的數值模型能夠很好地模擬板錨在砂土中的行為特性，并有效反映板錨尺寸、埋深和土體密度對板錨極限抗拔承載力的影響。

3.2 板錨循環承載力分析

Petereit[9]對圓形板錨在砂土中的循環承載特性做了大量模型試驗，試驗所用砂為Ottawa砂，其相對體積質量為2.65，不均勻系數為1.9，摩擦角為40°，最大密度1.7 g/cm3，最小密度1.47 g/cm3，相對密度80%，孔隙比為0.61。板錨直徑為10.16 cm，厚度為0.635 cm，板錨初始埋深為914 cm，試驗中此埋深下板錨單調抗拔承載力Pu=4 895 N，試驗用循環荷載為正弦荷載，如圖11，周期T范圍為0.5～10.0 s。模擬的試驗組為：循環荷載周期T為4 s，循環幅值Qc=25%Pu、50%Pu和75%Pu，試驗用砂密度1.65 g/cm3，試驗工況為板錨水平埋置，施加豎向荷載向上拔出板錨，模型參數見表1。

圖11 循環荷載參數Fig. 11 Parameters of cyclic load

3.2.1 板錨極限承載力分析

為了驗證Ottawa砂模型參數的準確性，首先針對板錨在Ottawa砂中的單調抗拔過程進行模擬。圖12給出了模型試驗和數值模擬中板錨在單調荷載作用下的荷載—位移曲線，該工況埋深比H/B=9，符合3.1節中介紹的深埋條件，其荷載—位移響應也與深埋板錨上拔行為類似，因此再次證實了所建立的適用于飽和砂土中板錨排水循環動力分析數值模型的有效性。

圖12 板錨在單調荷載作用下的荷載—位移曲線Fig. 12 Load-displacement curve of a plate anchor under monotonic load

3.2.2 板錨循環承載力分析

圖13中，U表示板錨位移，F為板錨所受荷載。可以看出，初始加載時，板錨在荷載作用下會產生一定位移，隨后卸載，板錨位移會恢復一部分，但是無法回到初始位置，即一次加—卸載之后土體產生了塑性變形。再加載時，板錨再次經歷移動和回彈，由圖13(a)可以看出，在后續加—卸載過程中，每一周次產生的塑性變形逐漸減小直到一個穩定值。Bemben等[8]和Hanna等[31]對長周期循環荷載下板錨在飽和砂土中的承載力做了模型試驗，觀察到了類似現象：當承受循環荷載時，板錨會產生不可恢復的位移，發生明顯的“棘輪效應”，即在應力控制循環加載下產生循環塑性應變累積[32]。隨著荷載周次的增加，每周次產生的位移會逐漸減小，累積位移速率雖然減小，但并沒有停止。這種永久位移累計速率隨著循環次數的增加而逐漸減小的特性，與荷載條件有關，但文中本構模型中的硬化公式C(ξ)利用應變累積來考慮應力歷史的影響，得到的土體剛度隨加載周次的變化速率偏小，所以針對大周次“棘輪效應”的模擬效果有限。

圖13 板錨在循環荷載作用下的荷載—位移曲線(Qc=75%Pu)Fig. 13 Load-displacement curve of plate anchor under cyclic load (Qc=75%Pu)

3.2.3 循環幅值影響

圖14給出了不同循環幅值的循環荷載作用下，板錨的位移—周次響應，圖中N為循環周次，U為板錨位移。由圖可以看出，循環幅值越大，初始位移和位移變化率越大，這是因為更高的荷載會給周圍的土體帶來更大的應變。從圖中可以看出，每組試驗中，數值計算的位移值一開始都先小于試驗測量值，但隨著循環周次的增加逐漸增大并超過測量值，且循環幅值越大，越早完成超越。這是由于數值計算中采用的幾乎恒定的塑性變形率和實際變化的塑性應變率的差別所導致的，在初始階段，數值計算的累積塑性變形略小于模型試驗所測得的塑性變形，但隨著模型試驗中累積塑性變形速率的減小，數值計算的累積變形會逐漸超過模型試驗所測得的結果。

圖14 板錨在不同幅值循環荷載作用下位移—周次響應曲線(Qm=40%Pu)Fig. 14 Displacement-cycle number response curve of plate anchors under cyclic loads of different amplitudes (Qm=40%Pu)

3.2.4 荷載均值影響

圖15給出了Qc相同、Qm不同的工況對應的位移—周次響應曲線，可以看出，荷載均值Qm越大，初始位移越大，但是位移變化率越小，上述結論與離心機試驗結論[33]相同。因此，對于長期承受循環荷載的嵌入式海洋工程結構，可以通過增大循環均值Qm的方法來減小循環位移累積，延長服役壽命。

圖15 板錨在不同荷載均值循環荷載作用下位移—周次響應曲線(Qc=10%Pu)Fig. 15 Displacement-cycle number response curve of plate anchors under cyclic loads with different mean values (Qc=10%Pu)

4 結 語

利用帶誤差控制的顯式積分算法，將一種適用于飽和砂土排水循環動力分析的邊界面塑性模型編寫成可供有限元軟件調用的用戶自定義材料子程序，從而建立起適用于板錨在飽和砂土中循環承載分析的數值模型。得到以下結論：

1) 針對Toyoura砂的3組排水靜三軸試驗和2組排水動三軸試驗進行了模擬，數值計算結果與模型試驗結果相吻合，驗證了模型具備合理描述砂土在不同荷載條件下的力學響應的能力；

2) 建立了飽和砂土中板錨承載分析的數值模型，并模擬了板錨單調上拔過程。數值計算得到的最大抗拔承載力、對應最大抗拔承載力的位移、殘余抗拔承載力與離心機試驗結果一致，表明建立的動力本構數值模型可以有效模擬板錨在砂土中的行為特性；

3) 利用以上建立的數值分析模型，模擬板錨在砂土中的循環承載過程。結果表明，隨著循環荷載的施加，板錨永久位移逐漸累積，循環荷載會導致板錨產生累積位移。通過考察循環荷載要素對板錨循環承載特性的影響，得到了與試驗結果相吻合的結論，即循環幅值越大，初始位移和位移變化率越大；循環均值越大，初始位移越大，但位移變化率越小。

4) 盡管利用建立的適用于板錨在飽和砂土中循環承載分析的數值模型能夠獲得上述結論，但仍然存在以下3個問題：① 迭代計算時間長；② 循環周次較大時計算累積誤差大；③ 確定模型參數需要較多室內試驗。動力本構模型是根據應變或應力增量進行計算，每一個循環周次需要大量的增量步，如果循環周次較大，計算時間就會很長，并且每一個增量步都會有一定的計算誤差，隨著循環周次的增加，累計誤差也會越來越大。而且動力本構模型參數較多，用于循環承載計算需要大量的室內土動力試驗來標定參數。這3方面原因制約了文中建立的適用于板錨在飽和砂土中循環承載分析的數值模型廣泛應用于工程計算。