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

詹開宇，曹留帥，萬德成

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

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

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

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

在海冰覆蓋率較高的寒區(qū)，通常在海洋平臺立柱的水線面附近設(shè)置抗冰結(jié)構(gòu)，以增加平臺的抗冰性能，常見的抗冰結(jié)構(gòu)為抗冰錐體。由于海冰的抗剪強度遠小于其抗壓強度，因此在海冰與海洋平臺相互作用時，帶有特定錐角的抗冰錐體可以將海冰的主要破壞形式從擠壓破壞轉(zhuǎn)變?yōu)閺澢图羟衅茐?，從而極大的減少平臺所受載荷作用，過程如圖1所示?？贡F體是一種非常有效的抗冰方式，對于抗冰錐體的研究是十分必要的，許多學(xué)者都對冰區(qū)的錐體結(jié)構(gòu)進行了探索和研究[8-10]。

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

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

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

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

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

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

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

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

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

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

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

2 數(shù)值模擬參數(shù)設(shè)置

2.1 結(jié)構(gòu)物模型參數(shù)

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

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

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

表1 錐形立柱材料屬性

2.2 海冰模型參數(shù)

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

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

表2 海冰單元與黏聚單元模型參數(shù)

2.3 網(wǎng)格收斂性分析

網(wǎng)格尺度和時間步長對計算時間和數(shù)值結(jié)果有重要影響。在顯示動力學(xué)分析中，時間步長通過最小單元尺寸確定，因此本研究中僅對網(wǎng)格收斂性進行分析。分別選取尺寸為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的網(wǎng)格，計算層冰以0.4 m/s的速度與錐形立柱的碰撞過程。模擬時長20 s并統(tǒng)計水平方向平均冰力大小、相對誤差和計算所需時間，結(jié)果如表3所示?？梢钥吹剑S著網(wǎng)格尺寸的增加，計算結(jié)果逐漸收斂，同時計算時間迅速增加。權(quán)衡計算精度和計算時長，本研究選擇尺寸為0.05 m×0.05 m×0.01 m的網(wǎng)格進行計算。

表3 網(wǎng)格尺寸對計算結(jié)果和計算時間的影響

2.4 接觸模式

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

2.5 計算域與邊界條件設(shè)置

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

圖6 計算域設(shè)置Fig. 6 Computational domain settings

3 計算結(jié)果與分析

3.1 層冰破壞與碎冰堆積

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

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

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

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

3.2 立柱所受冰載荷

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

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

3.3.1 ISO冰力標(biāo)準(zhǔn)

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

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

(1)

VB=HB·hv

(2)

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

(3)

(4)

其中，f，gr和hv等系數(shù)由冰—結(jié)構(gòu)物間摩擦系數(shù)、錐徑、錐角大小和碎冰爬升高度等參數(shù)定義。

則在錐形結(jié)構(gòu)上的總冰力可以表示為：

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時間內(nèi)立柱所受水平方向和豎直方向的最大冰力與平均冰力大小，并將結(jié)果與ISO冰力標(biāo)準(zhǔn)計算得到的最大冰力以及狄少丞等[15]使用離散單元法得到的最大冰力結(jié)果進行對比。如圖10和圖11所示。

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

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

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

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

3.4 冰厚對冰載荷的影響

還研究了冰厚對冰載荷大小的影響。冰厚分別取0.2 m，0.25 m，0.3 m，0.35 m和0.4 m進行模擬。冰速取0.3 m/s，其余參數(shù)不變。圖12給出了t=60 s時間內(nèi)x方向最大冰力和平均冰力與冰厚的關(guān)系?？梢悦黠@看出，冰力大小與冰厚呈正相關(guān)，且增加速度接近線性。將計算結(jié)果與使用ISO冰力標(biāo)準(zhǔn)計算得到的結(jié)果相對比，具有較好的一致性。

圖12 冰厚與冰載荷的關(guān)系Fig. 12 The relation between ice thickness and ice force

4 結(jié) 語

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

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

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

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

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