鹽水柱狀模型冰彎曲強度數值計算分析

剛旭皓，田于逵，余朝歌，季少鵬，寇瑩

中國船舶科學研究中心，江蘇 無錫 214082

0 引 言

船舶在冰區中航行時會受到諸多因素的影響，如海冰物理力學特性、海冰密集度、船體形狀、航行速度等，其中冰的物理力學特性是通過直接影響冰層的承載能力來對船舶冰阻力產生最本質的影響。在冰區船舶破冰過程中，彎曲破壞為冰層主要的破壞形式之一[1]，因此，研究冰的彎曲強度對于冰區航行船舶的破冰阻力具有重要意義。

自20 世紀70 年代以來，研究人員針對冰的彎曲強度陸續開展了研究，當前針對冰彎曲強度的研究方式主要有現場試驗、模型試驗和數值計算等方法。Maattanen[2]在波羅的海的現場試驗中就冰樣尺寸和加載速率對海冰彎曲強度的影響進行了研究，得到冰的彎曲強度會隨冰樣寬度的增加而增加，但幾乎不受冰樣長度影響的結論。Frederking和Timco[3-4]分別在不同的時間就模型冰的彎曲強度進行了試驗，測得冰的彎曲強度與加載速率無明顯相關性。李志軍等[5]對細粒酒精冰彎曲強度的影響因素進行研究，確定了加載速率對模型冰彎曲強度的影響與渤海海冰類似，均呈韌脆轉變的特性。在數值計算方面，季順迎[6]利用離散元方法（discrete element method，DEM）對冰的壓縮強度和彎曲強度進行了數值計算，并對計算模型的細觀參數予以了校準；韓端峰等[7]利用SPH（smoothed particle hydrodynamics）方法對冰的三點彎曲進行模擬，定義了粒子間的作用類型。DEM 法與SPH法均是采用粒子法來對冰層進行模擬，粒子法的關鍵是確定粒子之間接觸的搜索方法。魏龍海等[8]采用鏈式結構搜索法，利用全局搜索的方法完整記錄了鏈式結構體內的數據。劉璐[9]開發了鄰居列表法，在搜索中只需對相鄰顆粒進行搜索判斷即可，節省了計算時間。

本文將在前期中國船舶科學研究中心小型冰水池研究的基礎上，利用原位懸臂梁方法對鹽水柱狀模型冰的彎曲強度進行一系列試驗研究，分析彎曲強度的影響因素，然后在此基礎上利用冰力學特性作為參數輸入，基于DEM 法對彎曲強度進行數值計算，再利用試驗數據對其進行校核，驗證數值計算方法的準確性，用以為模型冰力學特性預報提供參考。

1 基 于DEM 法 的 模 型 冰 彎 曲 強 度數值計算

1.1 DEM 法基本原理

1） 單元失效判斷。

在離散元模型中，當內力平衡時，單元之間的相互作用可以看作是一個動態的過程，通過追蹤單元的運動，可以得到每個單元的接觸力和位移。考慮海冰的材料性質和離散特性，采用具有黏結?破碎特性的平行黏結單元對模型冰顆粒進行構造[7]。平行黏結模型是通過在兩相鄰顆粒接觸區域添加附加材料，并設想有一組具有恒定法向和切向剛度的線性彈簧，以兩個單元的接觸點為中心，均勻分布在黏結圓盤上，從而完成單元間作用力和力矩的傳遞，如圖1 所示。圖中：Fn，Fs分別為接觸力的法向力和切向力；Mn，Ms分別為法向力矩和切向力矩；xA，xB為單元中心點；R為黏結圓盤半徑；n為法向方向[10]。

圖1 單元之間的平行黏結模型Fig. 1 Parallel bonding model between elements

單元之間的黏結失效通過梁理論進行求解：

2） 冰層離散化。

在計算兩個單元的接觸力之前，需要確定單元之間的接觸方式，判定其接觸條件。本文首先對單元進行編號處理，采用元胞數組儲存顆粒的位置和運動信息，并隨著時間步的更新，利用元胞數組的特性對粒子的位置和運動信息進行實時更新，簡化顆粒搜索過程，同時可實時計算粒子間的相互作用力。矩形單元采用最密排列的方式對冰層進行構造，冰層離散單元通過元胞數組進行編號（k1，k2，…，kn），然后利用冰單元kn計算相互間的作用力以及與其他冰單元之間的相互作用，如圖2 所示。

圖2 冰單元的編號方式Fig. 2 Numbering of ice units

然后，設置冰單元之間力的傳遞方式。根據冰層受力位置的不同，在冰單元之間力的傳遞也會有所不同。當左側冰單元（k1，k2）受力時，冰單元之間的接觸力會自左側向右和向后傳遞；當右側冰單元（k4，k5）受力時，接觸力將自右側向左和向后傳遞；當中部冰單元（k3）受力時，冰單元之間的接觸力自中間向兩側和向后傳遞。完成第1 步的傳遞之后，接觸力將逐漸拓展，直至覆蓋到整個冰層。最后，通過對接觸力的不同傳遞方式進行模擬，計算冰單元之間力的作用，以及冰層的裂紋生成與斷裂。

3） 單元間受力分析。

在單元黏結處，單元受力和位移的關系可以通過以下參數來表示：法向剛度Kn和切向剛度Ks、摩擦系數μ、法向阻尼和切向阻尼，以及粒子重疊時在重疊區域中心形成的接觸區。對接觸處進行力學分析，接觸力Fi可分解為法向力Fn和切向力Fs：

式中，ni和ti為單位向量。根據Hertz 接觸理論，法向力Fn為彈簧和法向阻尼器作用在單元上的彈性力Fe=Knxn與阻尼力Fν=Cnvn的合力，可以表達為

式中：xn為單元之間的重疊量；Cn為單元間的法向阻尼系數；vn為法向相對速度。單元間的切向力Fs一般以增量的形式進行計算。剪應力增量?Fs為每個時間步內切向方向的變化量與切向剛度的乘積：

式中：Δxs，vs分別為兩個單元在接觸點處的切向位移增量和相對切向速度；Cs為單元間的切向阻尼系數。假設單元間的接觸力滿足摩爾?庫倫定律，則有

1.2 模型冰彎曲強度數值計算

在數值計算中，冰樣的加載方式如圖3 所示。圖中，L為冰樣長度，L0為加載點至冰樣根部的距離，h為冰厚，b為冰樣寬度，P為冰樣上的作用力。在測試過程中，設置冰樣根部為約束邊界，其另一端為自由邊界，以平行黏結圓盤黏結冰單元，以不同的加載速率在冰樣端部向下施加載荷。冰單元的受力逐漸拓展至整個冰樣，直至超過單元間的黏結強度，冰樣發生彎曲破壞。根據懸臂梁的受力特點，冰樣的彎曲強度 σf為

圖3 用彎曲強度測試懸臂梁的方法示意圖Fig. 3 Schematic diagram of a test method for cantilever beams using flexure strength

當冰樣發生彎曲破壞時，彎曲測試力達到最大值Pmax，此時，計算得到的彎曲強度即為模型冰的彎曲強度 σf。

根據模型冰的力學參數試驗數據，包括彎曲強度、壓縮強度、彈性模量等，獲得數值計算的細觀參數如表1 所示。

表1 數值計算參數Table 1 Numerical calculation parameters

通過冰彎曲強度試驗研究，確定數值計算工況，變換冰樣的加載速率和冰樣尺寸。其中，冰樣加載速率的取值范圍為1～5 mm/s，冰樣寬度為2h～6h，冰樣長度為4h～10h。此試驗與冰的力學試驗工況相同，旨在研究彎曲強度的加載速率相關性以及尺寸效應。確定了計算參數和工況之后，開展冰樣彎曲強度的數值計算，圖4 為加載速率為5 mm/s、冰樣尺寸為240 mm×80 mm×40 mm時的計算時歷曲線。

圖4 彎曲強度計算時歷曲線Fig. 4 Time history curve of flexure strength by calculation

2 模型冰彎曲強度試驗測量

依托中國船舶科學研究中心小型冰水池（圖5）的試驗條件，針對鹽水柱狀模型冰開展彎曲強度試驗[13]。該模型冰是在0.5% 的氯化鈉溶液中通過自然冷凍形成。首先，將鹽水混合物降溫冷凍至冰點附近進行引晶，在水體表面形成一層極薄的細小顆粒冰晶層，進而在其下方引起以柱狀結構為主的冰顆粒的生長；然后，在?18 ℃的條件下逐漸冷凍形成包含鹵水泡的冰層，在達到目標厚度后，開始回溫。在此過程中，冰層中的鹵水慢慢流失，逐漸削弱冰層。達到目標強度后，維持溫度不變，開始試驗測試。

試驗中，冰樣厚度h的均值為40 mm。為研究冰樣尺寸與彎曲強度之間的關系，冰樣寬度取2h～6h，冰樣長度取4h～10h。在測試過程中，以1～5 mm/s 的加載速率在冰樣端部向下施加載荷直至冰樣發生彎曲破壞，然后記錄不同時刻冰樣受到的作用力。試驗時，每2～4 個冰樣為一組，以保證結果的準確性。圖6 所示為加載速率為5 mm/s、冰樣尺寸為240 mm×80 mm×40 mm 時的加載過程和時歷曲線。可見，冰樣受力呈線性增加，冰樣的破壞方式為脆性破壞。

3 數值計算與模型試驗比較分析

冰的彎曲強度對冰區船舶破冰航行時的冰阻力有著重要影響，充分研究冰的彎曲性能，可以更好地了解冰的破壞方式及其與船體作用的原理。本文將通過一系列數值計算，確定冰樣尺寸、加載速率和冰厚等因素的影響，并對數值計算模型進行校核，分析計算結果的誤差率，以確保計算結果的可靠性。

將數值計算時歷曲線與模型冰力學試驗時歷曲線進行對比，結果如圖7 所示。由圖可以發現，兩者在冰樣及彎曲過程中，時歷曲線均表現為線性變化，力的增加較為均勻，冰樣的破壞方式為明顯的脆性破壞。

圖7 數值計算與模型試驗時歷曲線對比Fig. 7 Comparison of time history curves between numerical calculation and model test

加載速率對冰彎曲強度的影響在不同的海域有著不同的表現。Timco 等[14]在統計分析大量的極地海冰研究數據后提出，加載速率對彎曲強度的影響并不明顯。隋吉學等[15]通過在渤海的實測發現，加載速率與彎曲強度之間實際有著明顯的相關性，且彎曲強度在不同的加載速率下有一定的韌脆性轉變關系。而在統計模型冰試驗數據后發現，彎曲強度與加載速率之間并無明顯的相關性，彎曲強度趨于某一恒定值，如圖8 所示。試驗和計算結果與Timco 等[14]和Frederking 等[3]的海冰彎曲強度試驗測量結果吻合，說明該模型冰與北極海冰的特性相似。對模型冰的破壞過程進行分析，顯示在不同加載速率下，冰樣晶體錯位程度不同，冰樣的破壞方式也不同。當加載速率較慢時，破壞裂紋逐漸生成，直至冰樣破壞；當加載速率變快時，則表現為明顯的脆性破壞。

圖8 不同加載速率下的計算結果比較Fig. 8 Comparison of calculation results at different loading rates

如圖9 所示，在試驗中改變冰樣長度后測量得到的模型冰彎曲強度其數據點十分離散，沒有明顯的規律性，冰樣長度的改變與彎曲強度之間沒有明顯的相關性。在數值計算中，改變冰樣的長厚比時，隨著冰樣長度的增加，發生彎曲破壞時冰樣根部的轉動幅度也會增加，導致冰晶顆粒錯位充分，變形量較大；數值計算結果呈現出與彎曲強度試驗結果相同的變化趨勢，其間的相對誤差均值約為12%。改變冰樣的寬度之后，如圖10所示，冰的彎曲強度隨著冰樣寬度的增加迅速增加，彎曲強度與冰樣的寬厚比呈明顯的線性相關，計算的相對誤差均值約為15%，能夠大致反映彎曲強度與冰樣寬厚比之間的相關關系。

圖9 冰樣長厚比變化時計算結果比較Fig. 9 Comparison of calculation results when the ice sample's length-to-thickness ratio changes

圖10 冰樣寬厚比變化時計算結果比較Fig. 10 Comparison of calculation results when the ice sample's width-to-thickness ratio changes

通過對不同工況下數值計算結果與試驗結果的比較發現，本文的數值計算方法能夠模擬出彎曲強度隨速度、冰樣尺寸等參數變化的趨勢，同時，也基本能夠計算得到彎曲強度值，驗證了計算模型的合理性。

4 結 論

冰的力學性質是研究冰區航行船舶冰阻力的重要參數。通過對模型冰力學特性的研究，可以為冰區航行船舶的設計、校核以及總冰阻力預報提供重要的參數依據。本文在利用原位懸臂梁方法對鹽水柱狀模型冰的彎曲強度進行測量的基礎上，采用離散元模型對其進行了數值模擬。通過對計算結果的處理和分析，并與試驗測量數據進行比較，得到以下主要結論：

1） 由于冰材料的特性較復雜，冰的彎曲強度受多種因素的影響，其中，冰樣的寬厚比對其影響較大，彎曲強度基本不受加載速率和長厚比變化的影響。

2） 基于DEM 法建立的計算模型能夠反映彎曲強度隨參數變化而變化的趨勢，得到了冰樣彎曲強度值的基本范圍，可為冰的力學特性預報提供參考。

在今后的研究中，將進一步考慮加載方式、加載方向以及浮力作用對彎曲強度的影響，從而系統研究模型冰的彎曲強度。