基于黏聚單元法計算分析海洋平臺錐形立柱冰載荷

詹開宇，曹留帥，萬德成

(上海交通大學 船海計算水動力學研究中心(CMHL) 船舶海洋與建筑工程學院，上海 200240)

在尋求新資源和新能源的大環境下，越來越多的國家將目光聚集到了極地。南北極地區雖然常年被冰雪覆蓋，但其遼闊冰原下覆蓋的豐富資源是人類的一筆寶貴財富。近年來，人類對于海洋資源的開發逐步向寒區延伸和發展。在極地探索研究中，需要考慮冰載荷對海洋結構物的影響。冰載荷是寒區的主要環境載荷。極地的海冰以多種形式存在，如層冰、碎冰、冰脊和冰山等，其中層冰是主要存在形式之一，也是研究海冰與結構物相互作用時最常考慮的載荷形式。

關于結構物層冰載荷的測量與預報，國內外學者最早使用的方法是直接測量和試驗模擬方法，并總結出了一系列用以估算層冰載荷的經驗公式[1-3]。但是直接測量的難度較大，且設備的安裝和維護成本較高；試驗方法對試驗環境的要求較高，且對于不同類型的海洋結構物和不同形式的海冰普適性不佳。近年來，學者們開始采用多種數值方法進行冰載荷模擬嘗試，其中得到廣泛應用的是離散單元法(discrete element method，簡稱DEM)和有限元法(finite element method，簡稱FEM)。離散單元法將海冰模擬成黏結在一起的顆粒，通過改變顆粒間作用形式和作用力的參數，能夠較準確模擬處平整冰的破壞和堆積效果，在模擬海冰離散特性方面有較明顯的優越性；有限元法的重點和優點在于可以較好地研究海冰結構的變形。季順迎等[4]利用離散元法對海冰與直立結構相互作用進行了數值模擬，獲得了不同柱徑下的冰載荷和結構冰振響應。Kim等[5]采用新開發的有限元模型和模型試驗方法，對一艘貨船在浮冰條件下的破冰阻力進行了數值和試驗研究，對比結果得出了一系列重要結論。

黏聚單元法是在有限元方法的基礎上逐漸發展起來的，該方法在有限單元間插入了黏聚單元，因此在研究中可以同時考慮材料的變形和破壞。黏聚單元法可以較準確模擬出海冰的變形、破碎和堆積現象，如今已經在海冰材料模擬研究中得到了廣泛應用。劉路平[6]使用黏聚單元法進行了平整冰和抗冰海洋平臺相互作用的數值模擬研究，根據以往的試驗結果對模擬結果進行了驗證。王峰等[7]基于非線性有限元數值方法，引入黏聚單元模型并結合線性軟化彈塑性本構模型，對平整冰與豎直固定圓柱體碰撞進行了數值模擬。

在海冰覆蓋率較高的寒區，通常在海洋平臺立柱的水線面附近設置抗冰結構，以增加平臺的抗冰性能，常見的抗冰結構為抗冰錐體。由于海冰的抗剪強度遠小于其抗壓強度，因此在海冰與海洋平臺相互作用時，帶有特定錐角的抗冰錐體可以將海冰的主要破壞形式從擠壓破壞轉變為彎曲和剪切破壞，從而極大的減少平臺所受載荷作用，過程如圖1所示。抗冰錐體是一種非常有效的抗冰方式，對于抗冰錐體的研究是十分必要的，許多學者都對冰區的錐體結構進行了探索和研究[8-10]。

圖1 層冰與錐形結構相互作用后的破壞與堆積Fig. 1 Failure and accumulation of ice after interaction with conical structure

基于黏聚單元法建立海冰模型，并進行層冰和帶有抗冰錐角的固定平臺立柱相互作用的數值模擬。該數值模擬將探究海冰的厚度等條件對結構物受力的影響。同時，對于帶有抗冰錐體的立柱，錐體的角度對立柱的受力也存在一定的影響，也會改變海冰的破碎和堆積現象，這也是本文研究的另一重點。將數值模擬的計算結果與近似環境條件下的試驗和觀測結果、使用離散單元法模擬得到的結果以及基于ISO(international organization of standardization)冰力標準計算得到的結果進行對比驗證，在一定范圍內具有較好的一致性。研究表明，基于黏聚單元法建立的海冰模型在研究層冰與帶有抗冰錐角的固定平臺立柱相互作用時具有較好的適用性，對于后續有關寒區結構物的冰載荷模擬以及抗冰結構的設計具有重要的參考意義。

1 基于黏聚單元法的層冰模型建立

海冰是一種性質復雜的溫度敏感性復合材料，其物理和力學性質與其本身的鹽度、孔隙率，環境溫度加載速率等都有密切關系，在不同條件下會表現出脆性和塑性兩種特性。此外，在海冰與不同形式的結構物相互作用時，其破壞形式也是多種多樣的，包括彎曲破壞，擠壓破壞，剪切破壞等。因此在進行冰載荷的研究和數值模擬時，確定海冰的本構模型、材料參數和破壞形式，對得到準確的模擬結果十分關鍵。

有限元法是一種研究冰載荷的有效方法，可以模擬出海冰的變形和破壞，也能快速得到結構物的受力和響應情況，目前得到較廣的應用。但該方法存在局限性，即經過碰撞破壞后的冰單元會失效并被刪除，因此有限元法并不能很好地模擬出層冰破碎后的堆積現象以及碎冰對于結構物的后續作用，同時造成整個模擬過程的質量不守恒，使模擬結果與實際情況有偏差。為了解決這一問題，學者們提出了很多的解決方案，如FEM和SPH(Smoothed Particle Hydrodynamics)方法耦合，黏聚單元法等[11-13]。

黏聚單元法是在有限元法基礎上優化和發展得到的一種方法。如圖2所示，該方法將使用有限元方法建立的海冰單元離散化，并且在冰單元間插入黏聚單元。插入的黏聚單元與冰單元共用節點，且厚度為零。圖3是構造的海冰的8節點六面體單元，水平和垂直方向的黏聚單元將海冰單元分隔。

圖2 黏聚單元與冰實體單元Fig. 2 Cohesive element and bulk element

圖3 海冰六面體單元模型與黏聚單元模型Fig. 3 Ice hexahedral element model and cohesive element model

在層冰和結構物相互作用的過程中，原先的海冰有限單元發生變形，通過與黏聚單元共用的節點傳遞力和位移；而插入的黏聚單元會根據海冰材料的本構模型和材料屬性參數發生變形，當其位移達到設定的破壞臨界點時，黏聚單元破壞，而與之相連的海冰單元將會脫落。在整個層冰的黏聚單元模型中，破壞僅在黏聚單元處發生。

黏聚單元的破壞是以給定的牽引力—位移曲線作為依據進行判定的。該曲線規定了黏聚單元所受的牽引力與相鄰單元的分離位移之間的關系。牽引力—分離位移準則的重要參數有斷裂應力、斷裂韌性、最大分離位移、曲線形狀等。根據上述參數的不同可以將該模型分為不同種類，目前常用的黏聚單元破壞準則模型有三種：線性軟化模型；指數軟化模型以及延性軟化模型(又稱梯形軟化模型)[11]。圖4為三種模型的牽引力—分離位移曲線。其中，曲線與坐標軸所圍成的區域的面積即使黏聚單元破裂所需的能量。下文使用的是黏聚單元線性軟化模型。

圖4 黏聚單元破壞準則模型的牽引力—分離位移曲線Fig. 4 Common shapes of traction-separation law

2 數值模擬參數設置

2.1 結構物模型參數

結構物模型為帶有抗冰錐角的海洋平臺立柱，來源于渤海灣JZ20-2MUQ平臺[12]。取一根立柱與海冰發生作用，未加錐角之前立柱的直徑為1.7 m，錐體高度2.5 m，錐體最大直徑為4 m。定義錐體斜面與水線面夾角θ為錐角大小，并保證水線面處(即立柱與層冰作用位置)錐體直徑一致，為2.8 m。平臺和錐形結構物如圖5所示。

圖5 渤海灣JZ20-2 MUQ平臺與錐形結構物示意Fig. 5 JZ20-2 MUQ platform on Bohai Bay and the model of conical structure

碰撞過程為層冰以恒定速度撞擊固定立柱。碰撞不考慮平臺變形，因此選擇立柱的材料為剛性材料。具體的材料屬性如表1所示。

表1 錐形立柱材料屬性

2.2 海冰模型參數

海冰的有限元單元實體模型采用8節點六面體單元建模，取18 m×40 m的層冰進行計算。研究了層冰厚度對于冰力大小的影響，冰厚h分別取0.2 m、0.25 m、0.3 m、0.35 m和0.4 m。采用各向同性的線性彈塑性材料模擬海冰實體單元的變形。在冰實體單元間插入厚度為零的黏聚單元項，黏聚單元的牽引力—分離位移曲線選擇線性軟化曲線模型。

海冰材料的基本參數受到許多因素的影響，在不同工況下很難有確定的數值。國內外學者對海冰進行了大量的力學試驗，并考慮海冰的溫度、鹽度、孔隙率等對海冰性質的影響，給出了工程應用中所使用的海冰參數的推薦范圍[13-14]。這里也根據該推薦范圍選擇了海冰彈塑性模型的材料參數，以滿足本算例的工程需求。具體參數值見表2。

表2 海冰單元與黏聚單元模型參數

2.3 網格收斂性分析

網格尺度和時間步長對計算時間和數值結果有重要影響。在顯示動力學分析中，時間步長通過最小單元尺寸確定，因此本研究中僅對網格收斂性進行分析。分別選取尺寸為0.075 m×0.075 m×0.01 m、0.05 m×0.05 m×0.01 m和0.025 m×0.025 m×0.01 m的網格，計算層冰以0.4 m/s的速度與錐形立柱的碰撞過程。模擬時長20 s并統計水平方向平均冰力大小、相對誤差和計算所需時間，結果如表3所示。可以看到，隨著網格尺寸的增加，計算結果逐漸收斂，同時計算時間迅速增加。權衡計算精度和計算時長，本研究選擇尺寸為0.05 m×0.05 m×0.01 m的網格進行計算。

表3 網格尺寸對計算結果和計算時間的影響

2.4 接觸模式

計算前設定層冰和結構物碰撞的接觸模式。立柱和層冰之間的接觸采用面面接觸的侵蝕算法，在LS-DYNA中選擇算法“CONTACT_ ERODING_ NODES_ TO_ SURFACE”，該算法可以定義海冰單元的破壞失效標準。層冰與立柱之間的摩擦力系數設定為0.2，破碎堆積的碎冰之間的接觸采用單面接觸的算法，在LS-DYNA中選擇接觸算法“CONTACT_ ERODING_ SINGLE_ SURFACE”，該算法可以簡化碎冰之間的作用，節約計算時間。海冰間的摩擦系數設定為0.1。

2.5 計算域與邊界條件設置

如圖6所示的計算域中，海冰沿x軸負方向移動，z軸與海冰厚度方向平行，與水線面垂直，且定義向上為正方向。層冰與立柱碰撞的邊界設定為自由邊界，其余三邊約束三個方向自由度，并設置為無反射邊界條件，以模擬半無限大的層冰。本算例沒有考慮水對碰撞過程的影響，僅在海冰下方增加海平面，并在垂直方向定義海水對海冰的浮力作用。

圖6 計算域設置Fig. 6 Computational domain settings

3 計算結果與分析

3.1 層冰破壞與碎冰堆積

將使用黏聚單元法得到的模擬結果與國內外學者的實測、試驗以及數值模擬等不同方法得到的結果進行對比，從而驗證該方法的可行性。

圖7給出了層冰以0.4 m/s的速度撞擊固定錐形立柱的模擬結果，圖中所示為模擬進行到t=10 s時的俯視圖和側視圖。觀察碰撞過程可知層冰與抗冰錐角斜面接觸和碰撞情況。海冰與斜面接觸時會在斜面上有一小段的爬升，同時海冰所受的主要載荷為彎曲載荷。當海冰的彎曲應力大于其抗彎強度時，海冰發生破壞，形成碎冰并在水線面附近堆積。隨著層冰的不斷推進，碎冰和層冰繼續與立柱相互作用。模擬中還觀察到碎冰在錐體斜面上爬升后滑落的現象。分析完整作用過程可知，層冰與固定錐形立柱相互作用時產生的冰載荷主要有兩個部分，一部分為層冰對立柱的作用，另一部分為破碎后堆積的碎冰對立柱的繼續作用。對比可知，本文模擬得到的一系列現象與渤海灣JZ20-2 MUQ平臺的實測結果以及HSVA(Hamburg Ship Model Basin，漢堡船模試驗水池)的模型試驗結果基本吻合[15]，如圖8所示。

圖7 t=10 s時刻俯視圖和側視圖Fig. 7 Top and side views at time t=10 s

圖8 模擬結果與實測和試驗結果對比Fig. 8 Simulation results compared with measured and experimental results

3.2 立柱所受冰載荷

圖9為層冰運動速度為0.4 m/s時錐形立柱所受x方向水平載荷的時歷曲線。從圖中可知，在與層冰碰撞的過程中，立柱所受冰載荷呈波動變化且幅值較大。這是因為層冰與立柱碰撞過程中，層冰一直在進行“接觸變形—破碎—堆積爬升—接觸變形”這一循環過程，海冰受到壓應力和彎曲應力作用，載荷上升，隨后黏聚單元破壞，應力卸載造成載荷下降。同時由圖9可知，碰撞前期冰載荷呈震蕩上升趨勢，隨后逐漸趨于穩定，在一定范圍內波動。這是因為在碰撞初期，隨著層冰不斷破壞，碎冰開始堆積，海冰與立柱的接觸面積逐漸增大，冰載荷開始上升。當碎冰堆積達到一定程度時，接觸面積不再有明顯變化，層冰碰撞與碎冰堆積的過程達到穩定，冰載荷開始在一定范圍內穩定波動。這一結果也與模型試驗得到的冰載荷結果相吻合。

3.3 抗冰錐角大小對冰載荷的影響

3.3.1 ISO冰力標準

ISO-19906標準[16]是廣泛應用于海洋結構物設計的行業標準，其中給出了不同類型的海冰與海洋結構物相互作用的冰載荷。使用冰力標準給出的層冰與錐形結構物相互作用過程中冰載荷計算公式計算得到了不同錐角條件下立柱所受最大冰載荷的大小，并將其與本文研究結果進行對比，以驗證該黏聚單元模型在冰載荷研究中的可行性。ISO標準中層冰與錐形結構物相互作用的冰載荷計算公式主要基于海冰的塑性破壞和脆性破壞。海冰的塑性破壞通常發生在載荷加載速率較小的情況下，海冰在外力作用下的變形量超過其最大彈性變形量，產生明顯的永久變形，進而破壞失效；當加載速率增加時，海冰破壞模式逐漸向脆性破壞轉變；在大速率外載荷沖擊下，海冰會發生無明顯變形的脆性破壞。選取海冰的塑性破壞模式進行計算和對比。

塑形破壞模式中，錐形結構物所受冰力由兩部分組成：層冰彎曲斷裂和碎冰爬升對結構物的作用力。層冰破壞時的水平冰力HB和垂直冰力VB表示為：

(1)

VB=HB·hv

(2)

碎冰在錐面爬升時的水平冰力HR與垂直冰力VR表示為：

(3)

(4)

其中，f，gr和hv等系數由冰—結構物間摩擦系數、錐徑、錐角大小和碎冰爬升高度等參數定義。

則在錐形結構上的總冰力可以表示為：

FH=HB+HR

(5)

FV=VB+VR

(6)

3.3.2 錐角影響

為了研究抗冰錐角θ的大小對冰載荷的影響，取錐角θ分別為40°，45°，50°，55°，60°和65°進行了一系列的模擬研究。為消除錐角改變時造成的錐徑變化帶來的影響，保持水線面處錐徑不變，為2.8 m。計算了層冰以0.4 m/s的速度碰撞錐形立柱90 s時間內立柱所受水平方向和豎直方向的最大冰力與平均冰力大小，并將結果與ISO冰力標準計算得到的最大冰力以及狄少丞等[15]使用離散單元法得到的最大冰力結果進行對比。如圖10和圖11所示。

圖10 水平方向冰載荷與θ大小的關系Fig 10. The relation between θ and ice force of horizontal direction

圖11 豎直方向冰載荷與θ大小的關系Fig 11. The relation between θ and ice force of vertical direction

由計算結果可知，水平方向的最大冰力和平均冰力隨著抗冰錐角的增大而增大；豎直方向的最大冰力在一定錐角范圍內變化不大，而平均冰力隨著錐角的增大而減小。與其他結果對比可知，使用黏聚單元法得到的結果與其他方法得到的結果吻合較好，在變化趨勢上保持一致，但整體的冰載荷結果偏小。相較ISO冰力標準，造成這一結果的原因可能是黏聚單元和冰體單元參數的選擇和實際海冰材料不完全一致；另外，實際海冰材料內部結構復雜，破碎形狀隨機，而使用黏聚單元法建立的海冰模型是規則的各向同性結構，這也是造成誤差的主要原因。相較DEM，黏聚單元法在模擬海冰結構的變形和破壞時具有一定優勢，由于黏聚單元是插入在海冰單元之間的零厚度片狀單元，可以承受拉應力和剪應力但幾乎無法承受壓應力，而DEM中顆粒間的連接基于彈簧阻尼系統，可以更好地模擬海冰多個方向的受力，因此模擬結果更接近實際情況。這里使用的黏聚單元模型需要在這一方面進一步改進和優化。

考慮層冰與錐形結構物作用過程，錐角對冰載荷造成影響的主要原因可能是海冰與錐形結構物作用時，海冰的破壞形式發生了變化。海冰與錐體作用時的破壞形式主要由彎曲破壞和擠壓破壞，兩種破壞形式同時存在。當錐角較小時，海冰的破壞以彎曲破壞為主，隨著錐角增大，彎曲和擠壓同時發生且擠壓破壞逐漸稱為主要破壞形式，而海冰的抗壓強度要明顯大于抗彎強度。因此，隨著錐角增大，x方向上的冰力逐漸增大；同時，隨著錐角增大，碎冰在斜面上的爬升越來越困難，爬升高度降低，因此z方向的冰力呈減小趨勢。

3.4 冰厚對冰載荷的影響

還研究了冰厚對冰載荷大小的影響。冰厚分別取0.2 m，0.25 m，0.3 m，0.35 m和0.4 m進行模擬。冰速取0.3 m/s，其余參數不變。圖12給出了t=60 s時間內x方向最大冰力和平均冰力與冰厚的關系。可以明顯看出，冰力大小與冰厚呈正相關，且增加速度接近線性。將計算結果與使用ISO冰力標準計算得到的結果相對比，具有較好的一致性。

圖12 冰厚與冰載荷的關系Fig. 12 The relation between ice thickness and ice force

4 結 語

基于黏聚單元法建立了海冰模型，并模擬了帶有錐形結構的立柱與層冰相互作用的動力過程，分析不同參數對碰撞過程和冰載荷的影響，得到的主要結論如下：

1) 黏聚單元法可以較準確地模擬層冰與錐形立柱的碰撞過程，該方法在模擬海冰的破碎、堆積以及碎冰在立柱上的爬升等方面有較好的表現，模擬結果與實測結果以及試驗結果基本吻合。

2) 抗冰錐角的大小對冰載荷有明顯影響。研究結果表明，錐角角度在40°左右時有較好的抗冰效果。隨著錐角增大，水平方向的冰載荷有明顯增加，而豎直方向冰載荷逐漸減小。造成這一結果的主要原因為錐角改變了海冰的破壞形式。該模擬結果與ISO冰力標準以及其他學者得到的結果吻合。

3) 研究了層冰冰厚對冰載荷影響。結果表明，錐形立柱所受層冰載荷與冰厚呈正相關，且增長速度接近線性。該結果與ISO冰力標準吻合較好。

在下一步工作中，將增加計算工況，進一步研究影響冰載荷的因素。此外，將對海冰的黏聚單元模型進行改進，同時優化網格，提高計算結果的準確性和可信度。