基于多介質、多尺度離散元方法的冰載荷數值冰水池1)

季順迎 田于逵 )

* (大連理工大學工業裝備結構分析國家重點實驗室,遼寧大連 116024)

? (中國船舶科學研究中心,江蘇無錫 214082)

引言

近年來,中國在極地科學考察、北極商業航運、極地觀光郵輪、冰區油氣開發和新能源(風電、核電)等方面的發展,對極地船舶與海洋工程技術的研究提出了迫切的工程需求,其中工程結構冰載荷的合理確定是重要的核心內容[1].為確定船舶結構的冰載荷,芬蘭于20 世紀60 年代最早開展了實船測量和反演,隨后美國、加拿大、挪威、俄羅斯、韓國和中國均開展了大量的相關研究;對海洋工程結構的冰載荷測量則主要集中在中國渤海的導管架平臺[2-3]、波羅的海的燈塔[4-5]、加拿大的聯邦大橋[6]、沉箱式平臺[7]等.而對船舶與海洋工程結構冰載荷的模型試驗研究,則從20 世紀70 年代于德國漢堡冰水池(Hamburg Ship Model Basin,HSVA)開始,隨后芬蘭、美國、加拿大、俄羅斯、挪威、韓國和中國等寒區國家均開展了系統的研究[8-10].毫無疑問,現場測量和冰水池模型試驗是確定冰載荷的重要手段.但是相對于現場測量和模型試驗,數值方法具有準確性、詳實性、高效性、經濟性、場景性和開放性等優勢,也是國內外開展極地船舶與海洋工程研究的重要途徑[11-13].當然,適用于工程應用的數值方法需要依托于正確的理論模型、可靠的數值方法和合理的計算參數,并通過全面系統的試驗驗證.

在極地船舶與海洋工程的冰載荷研究中,國內外學者從20 世紀80 年代就開始注重采用數值方法進行計算分析.在對船舶結構冰載荷的數值模擬中,不僅考慮平整冰、碎冰和冰脊等不同海冰類型[14-15],還研究了船舶的直行、回轉以及Z型航行等不同航行狀態[16],從而確定了船舶結構的局部冰壓力和不同工況下的冰阻力.此外,對冰山撞擊下的船舶結構冰載荷和結構破壞也可進行數值分析.對于導管架式、自升式、浮式和半潛式等不同類型的海洋平臺結構,采用數值方法也可對其在不同海冰條件下的冰載荷進行系統的計算分析,從而確定冰厚、冰速、冰強度等海冰條件以及平臺結構的尺寸、幾何形式等因素對冰載荷的影響[15,17-18].通過對船舶與海洋平臺結構冰載荷的數值分析,不僅可以為結構抗冰設計、疲勞分析提供可靠的參考依據[19-20],還可對其在冰區的安全作業提供有力的數據支持.

在對冰載荷的數值分析中,人們廣泛地采用了有限元[21-24]、離散元[25-27]、黏聚力模型[28-29]、近場動力學[30-34]、環向裂紋法[35-36]和光滑粒子動力學[37]等不同方法,從而對海冰與結構物耦合作用時的破壞過程、冰載荷動力特性進行精確的分析.特別是離散元方法,自20 世紀80 年代被應用于海洋工程結構冰載荷分析以來,由二維模型向三維模型、由圓盤單元向塊體單元、由剛體碰撞計算向粘接破壞分析不斷發展,從而取得了最為廣泛的應用[38-41].同時,由于海冰與極地海洋工程結構的相互作用不僅涉及到海冰不斷發生斷裂、破碎和運行的動力過程,同時還要與海水、工程結構發生耦合作用,從而使該動力過程呈現出多介質、多尺度的復雜動力學行為.由此,基于冰-水-結構耦合的DEM-CFD-FEM數值方法在近年來取得了很大的研究進展[42-43],并將成為更有效確定極地船舶與海洋工程結構冰載荷的主要數值方法.季順迎等[44]、韓端鋒等[45]、徐瑩等[46]和Xue 等[47]從不同角度對目前海冰數值方法進行了相對全面的綜述,討論分析了目前存在的問題和發展的方向.

冰載荷數值方法的發展一直伴隨著試驗驗證工作的開展.通過現場測量和模型試驗,不僅可以明確海冰與不同類型海洋工程結構作用過程中的破壞模式[48-49]、冰載荷的局部分布和總冰力[50],同時還可獲得船舶與海洋工程結構在海冰作用下的力學響應[41,51].通過測量的海冰破壞模式、冰載荷和結構響應,可對數值方法及計算參數的可靠性進行驗證和改進[52].此外,海冰的物理力學性質是影響冰載荷特性的決定性因素,其包括冰晶結構、海冰溫度和鹽度、加載速率和方向等諸多因素[53-56].只有將海冰的物理力學性質、本構模型和破壞準則合理地引入到數值模擬中,才能準確地計算極地船舶和海洋工程結構的冰載荷.這也需要通過大量的試驗結果進行驗證.

極地船舶與海洋工程結構冰載荷的數值方法一直伴隨著相關學科的發展和新技術的應用.計算流體動力學的發展、近場動力學等數值方法的建立、高性能并行計算技術的應用均有力地促進了冰載荷數值研究的深入開展,并不斷地走向實際工程應用[27,47,57].此外,冰載荷數值分析的前后處理可有助于計算工作的簡便性、直觀性和實用性,其在很大程度上也取決于計算機技術的發展.特別是近幾年,機器學習和數據挖掘技術在力學計算中的廣泛應用[58-59],也將對冰載荷數值模擬研究產生深遠的影響.

在以上冰載荷的數值方法和試驗驗證中,不同數值方法對船舶與海洋工程結構冰載荷的研究相對獨立、分散性較強,從而導致目前數值研究工作的不系統和不規范,并產生較多的重復性研究.如何系統地考慮不同數值方法的優缺點、發展統一的前后處理系統,進而建立一個具有很強適用性的數值模擬系統平臺,則是目前冰載荷數值模擬工作所必須面對的問題.該工作可以借鑒近年來發展起來的數值水池研究[60-62].數值水池的概念雖然于20 世紀70 年代提出,但其快速發展則得益于在21 世紀初歐盟實施的“虛擬試驗水池計劃(VIRTUE)”,同時中國在近10 年也取得了顯著的研究成果[63-65].參考數值水池的研究模式,可建立數值冰水池或冰載荷數值試驗系統的研究框架.但是,由于數值冰水池需要考慮海冰在與船舶及海洋工程結構耦合作用中的斷裂、堆積等動力過程,其在數值實現上具有更大的復雜性和挑戰性.雖然,目前國際上也提出了數值冰水池的概念[42],但研究對象、工作場景、獲取信息等均過于單一,還遠遠未達到數值冰水池的范疇.同時,對數值冰水池的開發研究不僅限于對物理冰水池的數值再現,而是在此基礎上對其計算尺度、規模和場景進行擴展和延伸,從而發揮數值冰水池的優勢.

為此,本文將針對目前極地船舶與海洋工程結構對冰載荷及力學響應研究的工程需求,側重采用離散元方法探討數值冰水池的基本框架,介紹描述海冰物理力學性質和冰-水-結構耦合作用的數值方法,分析數值冰水池的軟件實現途徑和試驗驗證,由此提出數值冰水池的概念并探討其開發模式,并對其工程應用前景進行討論.

1 面向冰載荷的數值冰水池研究框架

下面提出面向極地船舶與海洋工程結構冰載荷預報的數值冰水池基本框架,簡要論述數值冰水池的主要研究目標、大體研究思路、主要研究內容和關鍵技術等,并將在后續研究中不斷改進和完善.

1.1 數值冰水池的基本概念和研究目標

數值冰水池的建立源于當前國內外對極地船舶與海洋工程結構冰載荷、結構力學響應、抗冰設計和安全保障研究的工程需求,更得益于近幾十年來海冰物理力學性質、冰載荷特性、工程結構強度和疲勞分析等相關領域在現場測量、模型試驗中的不斷知識積累,并與當前計算機技術和高性能數值算法的發展密切相關.數值冰水池的建立和發展不僅與極地船舶與海洋工程的研究密不可分,更受到計算流體力學、計算固體力學、并行計算技術、軟件工程等基礎學科發展水平的影響.此外,近年來人工智能、深度學習、大數據和數字孿生等新興學科在工程應用中的發展和實踐,也必將對數值冰水池的研發和應用起到很大的促進作用.

數值冰水池的建設目標是面向極地船舶與海洋工程技術研發,采用先進、高效和準確的數值計算技術,對海冰力學行為、結構冰載荷和力學響應進行全面的數值分析和試驗驗證,形成具有獨立知識產權的數值計算分析軟件集成系統平臺,為極地船舶與海洋工程的冰載荷特性確定、結構抗冰設計、安全保障服務提供可靠的技術手段.作為一種數字化的虛擬試驗系統,數值冰水池將與冰水池模型試驗和冰區現場測量結合、互補,并在一定程度上代替物理模型試驗和現場測量,從而有效地降低研究成本、提高研發效率、擴展對象范圍并增強精細化,從而更好地服務于極地船舶與海洋工程.

類似于數值水池的內涵和技術特征[66],數值冰水池也是一種依托于基礎研究的工程應用型技術.作為一種綜合性的虛擬試驗系統,最重要的是數值模擬結果的可靠性.同時,它還具有計算模塊封裝化、可靠性定量化、操作遠程化、試驗過程精細化和情景化等基本特征.數值冰水池集成軟件系統平臺的研發首先需要進行頂層規劃和總體設計,并根據所涉及的研究內容進行模塊化分解.每個模塊既可獨立運行,又可協同工作,具有很強的靈活性.在整體設計時,需要考慮新技術的采用和新功能的擴展,使其具有很好的開放性和共融性.由于數值冰水池的研發涉及多個學科的交叉與融合,這就需要建立一個聯合的研究團隊,分工明確又密切合作.相關研究思路在中國數值水池建設中得到了很好的實踐[64].數值冰水池的研究和建設將借鑒中國數值水池的成功經驗,并結合其特點和難點,重點開展海冰物理力學性質、冰-水-結構耦合作用的數值分析和試驗驗證、高性能計算分析軟件及工程應用研究.

1.2 面向冰載荷預報的數值冰水池基本框架

為實現數值冰水池在極地船舶與海洋工程結構技術開發中的研究目標,需要從海冰材料的數學模型、海冰破壞模式、船舶與海洋工程結構數據庫、冰-水-結構耦合模型、高性能計算分析軟件、前后處理系統、試驗驗證和典型工程應用等幾個方面開展系統的研究(如圖1 所示).數值冰水池的研發涉及船舶與海洋工程、極地科學、固體力學、流體力學、計算力學、試驗力學、計算機科學、軟件工程等諸多學科,是一個典型的學科交叉型、綜合性研究系統.

圖1 數值冰水池的基本框架Fig.1 Basic framework of the numerical ice tank

下面對數值冰水池的主要研究內容進行簡要介紹,并分別討論其所面臨的主要關鍵技術.

(1) 海冰物理力學性質的數值模型.創建可準確描述極地海冰物理力學特征、本構模型和強度準則的不同數值方法和關鍵計算參數的選取方案,并形成平整冰、碎冰區和冰脊等不同冰類型的構造算法.

(2) 船舶與海洋工程結構數據庫.建立破冰船、極地運輸船和極地科考船等典型極地船舶以及自升式、導管架式、半潛式和浮式等典型海洋平臺的結構數據庫.

(3) 冰-水-固耦合作用的數值模擬.建立海冰與海水、船舶或海洋工程結構耦合作用的多介質、多尺度數值算法并發揮不同計算方法的優勢和特點,構造冰水池模型試驗或現場測量場景下的邊界約束條件.

(4) 高性能數值計算分析軟件及前后處理系統.建立基于多核GPU 或CPU 并行的數值冰水池高性能數值算法以實現工程尺度下的數值計算;研發數值冰水池模擬結果的三維再現技術,對數值模擬結果進行精細化立體呈現,發揮數值冰水池的情景化優勢.

(5) 數值冰水池模擬的試驗驗證系統及可靠性檢驗.基于典型冰水池物理模型試驗和冰區現場測量數據,建立數值冰水池的不同尺度試驗驗證系統,提高數值模擬結果的可靠性和工程應用性.系統地收集整理國內外相關測量結果,建立用于數值冰水池可靠性檢驗的數據庫,并考慮模型試驗和現場測量中測量結果的不確定性,發展相應的可靠性驗證技術;重點構建用于數值冰水池驗證的基準試驗,實現對其全面系統的可靠性驗證.計算結果的可靠性是數值冰水池的核心內容,也是其實現工程應用服務的基礎和前提條件.

(6) 典型極地船舶與海洋工程技術開發的工程應用.實現典型極地船舶與海洋工程結構的工程應用,通過極地裝備結構冰載荷和冰激結構振動的數值模擬結果,為結構抗冰設計提供可行的優化方案,為結構疲勞分析和結構完整性管理提供有力的數據支持和工程服務.

數值冰水池在極地船舶與海洋工程結構冰載荷研究中具有顯著的特色和優勢,主要表現為可靠性、經濟性、快速性、擴展性和情景化.為保證數值計算的可靠性,需要針對研究對象、研究方法進行知識提煉和封裝,從而避免人為因素對計算結果準確性的影響;相對于現場測量和模型試驗,數值冰水池的經濟性和快速性顯而易見;在數值冰水池計算中,由于不受試驗場地、試驗設備等因素的影響,其可對材料尺度、模型尺度、工程尺度,甚至是地球物理尺度下的海冰力學行為進行計算,從而具有出色的擴展性.此外,通過開發具有出色三維再現功能的后處理系統,可對數值模擬結果從任意視角,乃至結構內部的力學信息進行顯示,從而使其具有出色的情景化功能,以促進對相關內在機理的理解.

2 海冰物理力學性質的數值方法

由于海冰在不同環境條件下生成的內部冰晶結構尺寸和排列形式均有很大的差異,且受溫度、鹽度、加載速率和加載方向等因素的影響,其物理力學性質極為復雜[53,67-68].因此,對海冰的物理力學性質進行合理的數學建模是數值冰水池的核心內容.

2.1 海冰的典型物理力學性質

在自然條件下,海冰在細觀上呈現出形態各異的冰晶結構并導致其復雜的力學行為[67,69].海冰的冰晶細觀結構如圖2 所示[70].在冰水池模型試驗中,為使冰晶結構具有天然海冰的細觀結構形式,逐漸研發了溫播種法和噴灑法等模型冰模擬技術.然而,如何對海冰的細觀結構進行數值構造則是目前國內外海冰數值模擬的難點.采用相場模型可數值模擬冰晶的生長過程[71],通過海冰熱力學模型也可以計算海冰的凍結過程.然而,對海冰進行冰晶尺度的細觀構造以研究極地船舶與海洋工程問題,目前國內外還尚未開展.由于海冰內冰晶結構的不同,特別是對于柱狀冰,沿冰晶生長方向的強度明顯高于垂直于冰晶生長方向[53],海冰的力學行為表現為各向異性.海冰數值模型一般會考慮不同加載方向下海冰強度的差異并建立相應的參數化方案.

圖2 北極海冰的冰晶細觀結構[70]Fig.2 Micro-crystal structure of sea ice in the Arctic[70]

由于海冰材料是由冰晶、鹵水和氣泡混合組成,孔隙率、溫度和鹽度等參數均對海冰強度有顯著的影響[57,72-73].海冰的單軸壓縮強度、彎曲強度、剪切強度、斷裂韌性和彈性模量等力學性能均隨鹽度和溫度的增加而明顯降低[55,57,68].考慮鹵水體積是海冰溫度和鹽度的函數,因此在對海冰的材料性能進行數值模擬時主要考慮鹵水體積的影響,并依據試驗結果將海冰強度設為鹵水體積平方根的負指數函數[56,74].

海冰在不同加載速率下具有不同的單軸壓縮強度,并隨加載速率的變化呈現出典型的韌-脆轉變特性,其轉變過程如圖3 所示[29,67,75-76].海冰的韌-脆轉變特性已被大量的現場試驗所證實,并被用于分析柔性直立海洋工程結構的冰激穩態振動現象[77].目前,為確定海冰在不同加載速率下的韌-脆力學行為和強度特性,直接定義其在不同加載速率下的強度是一種簡便可行的辦法.為將海冰韌-脆轉變特性的試驗結果進行全面的分析,最近基于機器學習的方法得以應用[58].然而,如何在數值模擬中描述海冰的韌-脆轉變特性還需要進一步對其發生的內在原因進行分析,并由此建立合理的數學模型.

圖3 海冰在單軸壓縮試驗中的韌-脆轉變過程[76]Fig.3 The ductile-brittle transition of sea ice in uniaxial compressive tests[76]

海冰與直立結構作用時會發生非同時破壞現象而產生高壓區(high pressure zone,HPZ),如圖4 所示[1,78].此外,海冰對直立結構的壓力隨接觸面積的增加而顯著降低,呈現出明顯的尺寸效應,如圖5 所示[79].在此基礎上依據現場實測結果進行改進,并分別確定了總體和局部冰壓力的變化趨勢[80-81].以上海冰的高壓區分布規律和尺寸效應已被國際規范ISO19906 所采用.海冰對結構的力學行為直接導致冰載荷的變化規律和極地海洋工程結構的力學響應,是海冰數值模擬中需要特別關注的問題[28,81-82].

圖4 海冰對直立結構擠壓破碎的高壓區[78]Fig.4 Sketch of high pressure zones on vertical structure[78]

圖5 冰壓力與接觸面積之間的關系[79]Fig.5 Relationship between ice pressure and contact area[79]

2.2 海冰力學行為的數值模擬方法

針對以上海冰的復雜力學行為,建立合理的海冰數學模型一直是國內外極地海洋工程領域關注的重要問題.針對海冰在自然條件下的離散分布特點及其與船舶和海洋工程結構相互作用中的破壞特性,離散元方法在海冰力學特性及冰載荷數值模擬中得到了廣泛應用[25,44,83].目前,采用離散元方法對海冰對斜面結構、直立和錐體海洋平臺、船舶結構的冰載荷進行了系統的數值分析[13,84-85].此外,采用離散元方法可以對平整冰、碎冰、冰脊以及航道內冰屑等不同冰類型進行數值分析(如圖6 所示)[86-88],也可以計算總冰力和局部冰壓力,具有很強的適用性[47,89-90].采用不同形態單元對冰-船/海洋平臺相互作用的計算結果如圖7 所示[29,84,91-92].為進一步提高離散元方法的計算效率,采用基于GPU 并行的數值算法可開展工程尺度下對船舶和海洋平臺結構冰載荷的數值計算[18,40,42,93].

圖6 采用球體離散元方法構造的不同海冰類型Fig.6 Different ice types constructed with the spherical discrete element method (SDEM)

圖7 海冰離散元方法中的不同單元形態Fig.7 Various element shapes in sea ice DEM simulations

在海冰離散元方法發展的同時,人們也采用有限元方法、黏聚單元法、環向裂紋法、近場動力學和光滑粒子動力學等數值方法對海冰的單軸壓縮、三點彎曲和巴西盤等物理試驗過程進行了系統的計算分析.以上方法同樣可以用于對海冰與船舶、海洋工程結構的相互作用進行數值計算,部分結果如圖8 所示.盡管以上數值方法的理論模型有很大的差異,但其中的一個共同點就是可以對海冰的破壞過程進行模擬,從而獲得海冰的動冰力過程.因此,如何建立合理的海冰材料失效準則并在數值方法中實現相關算法是決定計算可靠性的關鍵.當然,為更準確地模擬海冰的物理力學性質以及與工程結構作用時的破壞現象和冰載荷特性,還需要對相應的計算參數進行試驗驗證,同時也需要對其單元尺寸的依賴性進行系統的分析,以增強在實際工程應用中的可靠性[28,56].

圖8 不同數值方法模擬的冰-結構相互作用Fig.8 Interaction between sea ice and structure simulated with various numerical methods

3 冰-水-工程結構耦合作用的多介質、多尺度數值方法

在海冰與船舶及海洋工程結構的相互作用過程中,海冰的破壞模式、運動規律及由此導致的冰載荷特性均與流體、結構力學響應密切相關,因此在數值冰水池研究中需要考慮海冰與流體介質、工程結構的耦合效應[96].這就需要在海冰數值方法的基礎上,進一步發展冰-水-結構間的多介質、多尺度耦合數值模型,以提高冰載荷及結構力學響應的計算精度.以上研究已得到國內外相關學者高度重視.

3.1 冰載荷的多尺度計算模型

無論是對極地船舶,還是海洋工程結構,在具有顯著動力特性的冰載荷作用下均會引起結構的強烈振動,并由此產生耦合作用.例如,中國渤海的導管架式海洋平臺在海冰作用下發生強烈的冰激振動并導致疲勞損傷,引起上部設備和管線的振動破壞以及人員恐慌[97-98];加拿大的Molikpaq 沉箱平臺在持續冰激振動下導致砂心液化并引起結構下沉[99-100].船舶在冰區航行中,航速、航向和航態等因素對冰載荷有重要影響,同時船體會產生局部結構振動[84],且其螺旋槳也易受到海冰的影響發生葉片變形并導致推進系統的損傷[101-102].因此,冰激船舶與海洋工程結構的振動問題一直在通過現場觀測、模型試驗和數值模擬等不同途徑進行系統的研究.這些工作也將在數值冰水池研發中得以體現.

目前冰區海洋工程結構主要包括直立式和錐體導管架式海洋平臺、混凝土重力式海洋平臺、沉箱式海洋平臺等不同的結構類型,同時也在研發適用于極地冰區的自升式和浮式海洋平臺結構[103-105].以上研究中關注最多的是導管架式平臺結構,其中以慢冰速下直立導管架式平臺結構持續、劇烈且具有破壞力的自激振動或穩態振動為研究重點.對直立導管架平臺結構的冰激振動分析以往主要側重于強迫振動理論[81,106-107],而近年來則傾向于自激振動理論.采用有限元方法可對直立結構的冰載荷和冰激振動進行有效的數值分析[21].直立結構的穩態振動現象涉及到海冰的韌-脆轉變機理、微裂紋萌生與擴展特性、累積損傷模式等諸多細觀因素,是海冰與結構振動強耦合的過程,在數值模擬上具有很大的挑戰[108-109].目前,中國冰區的導管架式海洋平臺和風機基礎結構均在考慮如何避免自激振動現象[110-111].對錐體平臺結構的冰載荷時程和隨機冰激振動現象,則可采用DEM-FEM 耦合方法進行數值模擬,計算結果與現場測量結果相一致[43,112].由于DEM 的時間步長要遠遠小于FEM,因此在兩者耦合計算時需要采用相同的時間步長或進行時間上的插值,以實現兩種方法在計算信息傳遞的同步性.為提高DEM-FEM 耦合模擬的計算效率,可采用區域分解法并通過GPU 并行進行計算加速來擴大計算規模.以渤海遼東灣JZ20-2 油田單樁錐體NW 導管架海洋平臺為例(如圖9 所示),計算得到的冰激結構振動加速度如圖10 所示[51].

圖9 渤海JZ20-2 油田單樁錐體NW 導管架海洋平臺Fig.9 The NW jacket offshore platform of the JZ20-2 oil field in the Bohai Sea

圖10 采用DEM-FEM 方法模擬的錐體海洋平臺結構冰激振動[51]Fig.10 Ice-induced structure vibration of conical offshore platform simulated with DEM-FEM method[51]

船舶在冰區航行中,海冰與船舶結構碰撞的相對速度、接觸模式等因素對冰載荷影響顯著,并進一步反饋到航速和結構局部振動.無論是平整冰還是碎冰,船體冰阻力的經驗公式和模型試驗均表明,冰速是重要的影響因素[113-114].在數值模擬中,除考慮航速外,還要進一步考慮海冰的破壞、翻轉和滑動,及其對船體的作用位置和方向,從而合理地確定冰載荷的分布以及船體的力學響應.在采用環向裂紋法進行冰-船相互作用計算時,根據海冰破碎的物理特性定義其破碎長度,并由此確定具有隨機性的動冰載荷[14,36];類似于離散單元方法,Lubbad 和L?set[103]發展了具有粘接-破碎特性的塊體單元方法,對船舶結構的冰阻力進行了數值分析.采用近場動力學方法,對海冰與船體作用時的破壞現象及冰載荷變化時程進行了數值分析[33,40].如果考慮浮冰和冰山對船體的碰撞,則更需要考慮冰-船之間的耦合作用,進而從能量守恒的角度確定冰載荷的變化過程.目前對平整冰和碎冰區航行中的冰-船相互作用,更多地采用離散元方法進行數值計算,包括球體、圓盤和塊體離散元[13,39,84,115].在以上船體冰阻力的計算中,一般將船體視為剛體,但考慮其在海冰作用下的航速,并分析航速、冰厚、冰類型、冰塊尺寸等因素的影響.若進一步考慮船體在海冰作用下的應力分布、形變特性及局部振動,則需要建立船體結構的有限元模型,從而進行冰-船相互作用的全耦合分析.

3.2 冰載荷的多介質計算模型

作為在海洋中生成的海冰,其與海水之間無時無刻不在發生著熱力和動力耦合作用.特別是當海冰與船舶及海洋工程結構相互作用時,海冰的破碎、運動等均受海水影響,并導致冰載荷和結構響應的不斷變化.因此,在數值冰水池研發中必須要考慮海冰與海水間的耦合作用,以更加準確地進行數值模擬.以往,在對海冰與海洋結構物的相互作用數值模擬中,大多將海水設為定常流以簡化計算.然而對平整冰在斜面結構作用下的破壞模式研究表明,水動力學的影響對海冰的破壞模式有顯著影響[116].海冰也同時會導致船舶結構附近區域的流場更加復雜[117].在分析船舶和浮式平臺的動力定位時,也需充分考慮流體動力學的影響[104].因此,目前對海冰-海水間的多介質耦合計算得到越來越多的重視.

Hopkins[118-119]最早建立了三維圓盤和多面體離散元方法以模擬海冰的堆積、重疊及其與斜面結構的作用過程.該方法考慮了海水對海冰單元的浮力和拖曳力(矩),從而可以有效地計算冰塊在流體中的動力過程[120].在此基礎上,Huang 等[121]分別采用三維圓盤單元和有限體積法對海冰和波浪進行了耦合計算,分析了船舶在波浪作用下碎冰區航行中的冰阻力特性,如圖11 所示.研究平整冰與正/倒錐體海洋平臺結構相互作用的離散元模擬結果表明,受海水浮力影響,正錐體上的冰載荷要高于倒錐體[48].采用Star-CCM+軟件,Lou 等[88]對碎冰航道內的船舶航行阻力進行了數值分析,其中船舶結構采用剛體模型,碎冰屑采用組合球體單元,而流體采用有限體積法模擬.以上研究表明,在海冰與船舶及海洋工程結構相互作用時,流體動力學會有顯著的影響,需要在數值模擬時考慮海冰與流體的耦合作用,以更精確地模擬冰載荷和工程結構動力響應.

圖11 采用CFD-DEM 方法模擬的船舶在碎冰區航行[121]Fig.11 Ship navigation in broken ice field simulated with CFD-DEM coupled model[121]

海冰不僅與船舶及海洋工程結構存在多尺度耦合作用,同時還與海水存在多介質耦合,因此建立冰-水-結構間的多尺度、多介質全耦合模型是準確開展冰載荷及結構響應數值模擬的發展趨勢.這也是當前化工、巖土、交通和生物等領域所面臨的共同問題.為分析高壓水槍對巖石的沖擊破壞、泥石流對壩體沖擊時的運動特性及沖擊力、海底管道對顆粒-流體兩相介質輸運時的結構振動及損傷問題、均建立了CFD-DEM-FEM 耦合方法并取得了很好的數值結果[122-124].相對于其他領域,在極地船舶與海洋工程研究中,冰-水-結構間的全耦合模型研究剛剛開始,還需要借鑒其他相關領域的研究方法,發展海冰離散元、計算流體動力學與工程結構物有限元方法的DEM-CFD-FEM 全耦合模型.這里的海冰離散元方法可以采用球體或塊體單元構建合理的海冰類型,并考慮單元間的粘接-破碎功能以模擬海冰與結構物相互作用時的破壞現象.當然,也可采用近場動力學、黏聚單元方法、有限元方法等不同的數值方法模擬海冰材料;計算流體動力學則可以根據研究需要采用光滑粒子動力學方法、有限體積法等不同的數值方法;對于極地船舶與海洋平臺結構的有限元方法,則依據結構特點而采用梁單元、實體單元或殼單元.當然,更為重要的是建立不同數值方法間計算參數的精確、高效傳遞機制,并在此基礎上發展相應的高性能并行計算技術.

通過對海冰與海水、船舶及海洋工程結構間耦合作用的多介質、多尺度計算分析,可對冰激結構振動、形變等力學響應對海冰破壞模式的影響,進而更準確地計算冰載荷的動力特性.

4 數值冰水池的軟件實現

類似于數值水池的知識封裝、可靠性和情景化三大技術特征[66],數值冰水池也需要通過高性能計算機軟件實現相應的知識封裝和情景化.為此,需要發展相應的前后處理系統、工程結構數據庫和三維可視化顯示技術,并通過發展相應先進的并行算法實現高性能數值計算.

4.1 典型海冰與工程結構數據庫

考慮數值冰水池的使用對象主要為工程技術人員,其應為極地船舶和海洋工程的抗冰設計提供便捷的開發工具,因此它需要結合極地海冰的特點以及典型極地船舶與海洋工程結構,建立相對豐富的數據庫.該數據庫可從以下方面進行考慮:

(1) 針對工程海域或物理冰水池的實際情況,可自動定義相應的流體和海冰;根據計算區域的特點,可設計計算邊界處不同方向上海冰的剛度、速度等,也可設計周期性邊界條件.

(2) 針對極地海冰的類型,可自動定義平整冰、碎冰、冰脊等不同的海冰類型并設定相關的計算參數.對于碎冰,需自動設定其冰塊形狀、尺寸、密集度、厚度等參數;對于冰脊則需進一步設定其帆高及角度、龍骨深度及傾角、冰脊寬度和長度等幾何性能.因此,還需進一步自動設定海冰的單軸壓縮強度、彎曲強度、密集度、溫度和鹵水體積等物理力學性質.

(3) 對于簡單的錐體或直立結構,可建立相應的結構類型,并通過調整結構參數進行直接設定;對于相對復雜的工程結構,則可構建與相關通用商業軟件的數據接口,從而導入船舶及海洋工程結構數據.對于船舶及海洋工程結構的約束、運動方式及速度,也可通過交互式界面設定.對于典型的極地船舶及海洋工程結構,則建立相應的數據庫,以便于直接調用、計算.

以上相關功能已在我國自主研發的海冰離散元計算分析軟件IceDEM 中初步實現,部分模塊如圖12 所示[84].但還需要在后續數值冰水池研發中不斷豐富和完善.

圖12 離散元軟件IceDEM 中海冰及工程結構參數設定[84]Fig.12 Settings of sea ice and structure parameters in the software IceDEM[84]

4.2 面向工程尺度的高性能數值算法

在工程尺度下極地船舶與海洋工程結構冰載荷的數值模擬中,特別是采用離散元方法、近場動力學方法時,由于單元數量龐大且計算步長很小,對數值計算的效率提出了迫切的要求.此外在開展冰-水-工程結構的耦合模擬時,同時涉及三相介質的數據傳遞和同步計算,更需要強有力的高性能算法.目前提高數值計算性能的方式主要包括基于多核CPU的并行計算和基于GPU 的并行計算,其中后者在近年來得到更快的發展和更廣泛的應用.無論是采用球體離散元還是塊體離散元,目前均發展了基于CUDA-GPU 并行的數值算法,計算規模可達105～106個單元[27,93].在冰激海洋平臺結構振動的DEM-FEM 耦合,也采用了GPU 并行算法[51,125].在Jan?en 等[42]的數值冰水池研究中,為計算非規則碎冰與船體的耦合作用,也采用了GPU 并行算法,且效果良好.目前,基于GPU 的并行計算技術已成為多介質、多尺度數值計算的重要手段,并將在數值冰水池的冰-水-結構耦合計算中得到進一步的應用和發展.

近幾年,基于人工智能、機器學習、數據驅動的力學數值計算得到了迅速發展,并在航空航天、材料性能、巖土力學等領域取得了很好的工程應用[59,126].相對于傳統的數值方法,以上數值方法可以有效提高數值模擬的效率.基于數據驅動的計算力學可分為基于能量泛函和基于距離泛函的數據驅動算法[126].目前,基于能量泛函的數據驅動算法已應用于海冰物理力學性質的數值預測研究,其通過大量海冰物理力學性能的實測數據,建立海冰強度、破壞模式與冰晶結構、溫度、鹽度、加載速率、約束條件等因素的對應關系,從而將本構數據代替本構模型[58].隨著基于機器學習的力學數值計算技術不斷完善,可通過海冰參數、結構類型、作用模式等計算信息直接預測結構冰載荷和力學響應,或通過對典型工況下冰載荷的計算,采用機器學習的方法預測其他工況,從而有力地促進數值冰水池的發展和工程應用.

4.3 前后處理系統及三維可視化技術

作為數值冰水池的三大特征之一,情景化的實現需要基于出色的前后處理系統和可視化技術.數值冰水池的前處理系統需要同海冰與工程結構數據庫建設相結合,從而增強其靈活性.在前處理系統中,需要將數值計算參數與物理試驗相一致.特別是對于海冰類型和參數的輸入,一方面可以從統計意義上進行設定(如圖13 所示),另一方面也可通過對物理試驗參數識別的途徑實現真實數值再現.這也是當前數值冰水池亟需發展的內容.通過三維可視化技術可對數值模擬結果進行直觀顯示,進而從物理機制上對冰-水-結構的耦合作用機制、冰的破壞模式有一個深刻的認識,如圖14～ 圖16 所示.此外,還可對一些物理試驗中不易測量和觀測的現象進行補充顯示,從而獲得更全面的資料信息.

圖13 HSVA 冰水池中海冰密集度71%冰況及數值冰水池[127]Fig.13 The HSVA ice tank with 71% ice concentration and its numerical ice tank[127]

圖14 船舶在碎冰區及平整冰區中航行的離散元模擬Fig.14 DEM simulation of ship navigation in broken and level ice areas

圖15 平整冰與錐體海洋平臺上部、中部和下部作用時的離散元模擬三維再現[18]Fig.15 3D reconstruction of DEM simulation of level ice interacting with the upper,middle and lower parts of a conical offshore platform[18]

圖16 船舶在碎冰區航行的離散元模擬三維再現[12]Fig.16 3D reconstruction of DEM simulation of ship navigation in broken ice area[12]

在極地船舶與海洋工程結構冰載荷的高性能離散元計算分析軟件IceDEM 中,海冰對船舶結構的作用模式、船體冰載荷分布規律、船體冰阻力可以同步顯示(如圖17 所示),從而可直觀地分析其對應關系.該三維顯示模塊還具有縮放、旋轉、透視等顯示功能,從而可從不同視角更有側重地觀察冰-船作用模式.以上前期研究為數值冰水池的研發提供了很好的參考借鑒.

圖17 船舶在冰區航行的冰載荷分布、冰阻力及冰-船作用模式的三維顯示Fig.17 3D display of ice load distribution,ice resistance and ice-ship interaction model of ship navigation in level ice area

5 數值冰水池的試驗驗證

數值方法的工程實用性在很大程度上取決于計算結果的可靠性和準確性.這也是所有數值方法必須面對的問題.數值冰水池所涉及的海冰物理力學性質、海冰失效模式、冰-水-結構耦合作用、冰載荷分布特性等研究內容均基于試驗驗證的合理性.該試驗驗證不僅依據冰水池內開展的模型試驗,同時也基于船舶及海洋平臺結構的現場測量,并在必要時參考相應的理論模型或半經驗公式.

5.1 海冰物理力學性質的試驗驗證

在海冰的離散元、有限元、黏聚單元方法、近場動力學等不同海冰數值方法研究中,均注重對海冰物理力學性質的試驗驗證,以建立合理的海冰本構模型和強度準則、發展理想的數值方法并確定相應的計算參數[13,128].在海冰的物理力學性質數值模擬中開展最廣泛的是海冰的單軸壓縮和三點彎曲[52,56,92,129].這里采用離散元方法對顆粒排列、粒徑、摩擦系數等參數影響下的數值模擬結果進行分析.對于海冰的單軸壓縮試驗,海冰試件尺寸(a×a×h)設為15 cm × 15 cm × 37.5 cm,加載速率0.01 m/s,得到的海冰受壓破碎過程如圖18(a)所示;保持加載速度不變,海冰三點彎曲試驗中的試件尺寸(l×a×a)為140 cm × 15 cm × 15 cm,得到的海冰彎曲破壞情況如圖18(b)所示[52].

圖18 海冰單軸壓縮和彎曲強度的離散元模擬[52]Fig.18 DEM simulation of uniaxial compression and three-point bending tests of sea ice[52]

為確定單元粒徑對海冰力學性質離散元模擬結果的影響,分別設定不同的試樣尺寸和單元粒徑對海冰的單軸壓縮強度和彎曲強度進行了計算.這里引入無量綱尺寸L/D,其中L為海冰試件加載截面的尺寸,D為顆粒單元直徑.海冰單軸壓縮強度隨無量綱試樣尺寸的變化規律如圖19 所示.從圖中可以看出,海冰的單軸壓縮強度和彎曲強度均隨尺寸比L/D的增大而增大;但當L/D>20時,海冰的強度基本保持不變.這說明當離散單元的粒徑足夠小時,海冰力學性質的計算結果受單元的尺寸影響逐漸減小.

圖19 離散元模擬中 L/D 對海冰單軸壓縮強度的影響[52]Fig.19 Influence of L/D on the uniaxial compressive strength of sea ice in DEM simulations[52]

海冰強度受海冰鹵水體積(溫度與鹽度的函數)、加載速率等因素的顯著影響.根據大量的海冰現場和室內測量數據,海冰宏觀強度與其鹵水體積呈負指數關系[53,55].由此,在離散元模擬時也將海冰單元間的粘接強度設為鹵水體積的函數.這里以海冰單軸壓縮強度的離散元模擬為例,其計算結果與試驗測量數據相一致,如圖20 所示.除考慮以上尺寸影響、鹵水體積外,在采用離散元方法模擬海冰的強度時,還需要進一步考慮單元在斷裂前的內摩擦系數、單元排列方式、加載速率等因素的影響,從而更好地數值再現海冰的物理力學性質.

圖20 鹵水體積對海冰單軸壓縮強度的影響Fig.20 Relationship between uniaxial compressive strength of sea ice and brine volume

5.2 船舶冰阻力及冰載荷的試驗驗證

船舶結構在平整冰區航行中的冰阻力已有多個經驗公式,其中應用最廣的是Lindqvist 公式、Riska 公式等,而在碎冰區的冰阻力估算,也有諸多學者提出了相應的計算公式.韓端峰等[45]對相關的經驗公式進行了詳細討論;Su 等[14,35]采用環向裂紋法對船舶冰區航行中的冰阻力進行了數值計算,并與Lindqvist 經驗公式進行了對比驗證.劉璐等[130]采用擴展多面體離散元方法模擬船體結構在冰區的航行過程,將計算的冰阻力與Lindqvist 經驗公式進行了對比分析.離散元模擬的平整冰區尺寸為600 m ×150 m,航速為1.0 m/s,冰厚為1.0 m,計算結果及冰阻力時程如圖21 和圖22 所示.將穩定階段的冰載荷均值作為冰阻力,并將其與Lindqvist 經驗公式計算結果對比,如圖23 所示.從圖中可以看出,離散元結果和Lindqvist 公式均隨冰厚增大而增大,且具有明顯的線性增加趨勢;同時兩者的結果在數值上十分接近.因此,通過Lindqvist 公式可以充分說明擴展多面體離散元方法在計算船舶結構冰載荷方面具有良好的準確性和可靠性.

圖21 雪龍號破冰船與平整冰相互作用的離散元模擬[130]Fig.21 DEM simulation of interaction between the Icebreaker Xuelong and level ice[130]

圖22 船體冰阻力時程曲線[130]Fig.22 Time history of ice resistance on ship hull[130]

圖23 船體冰載荷的離散元計算結果與Lindqvist 公式對比[130]Fig.23 Comparison between DEM results and Lindqvist empirical formula of ice loads on ship hull[130]

國際船級社組織(International Association of Classification Societies,IACS)對特定工況下船體結構的冰載荷給出了規范計算方法[131].為驗證離散元方法對船體結構與大塊浮冰碰撞過程中冰載荷計算的可靠性,劉璐等[132]將離散元計算的冰壓力與IACS 規范進行了對比.計算情景如圖24 所示,碰撞點到艏柱的距離分別為5 m,15 m,25 m 和35 m.當航速為2.25 m/s,海冰厚度為3.0 m,海冰彎曲強度為1.0 MPa 時,船體結構與大塊浮冰碰撞的離散元模擬過程如圖25 所示.圖26 為4 個算例中碰撞點處的冰壓力與規范對比情況,其中x為碰撞點到艏柱的距離.計算結果表明,離散元結果與規范值相比存在上下誤差,但壓力值保持在相同量級范圍內,誤差范圍為6.7%～ 18.1%之間.由此可見,離散元方法對船體結構的冰壓力可進行準確可靠的分析,具有良好的工程適用性.

圖24 冰壓力IACS 規范驗證的離散元模型[132]Fig.24 Sketch of DEM simulation for the validation with IACS standard[132]

圖25 船體結構與大塊浮冰碰撞的離散元模擬[132]Fig.25 DEM simulation of collision between ship hull and ice floe[132]

5.3 海洋平臺冰載荷的試驗驗證

海洋平臺結構的冰載荷研究也已開展了大量的現場試驗和模型試驗,并用于海冰數值方法的驗證.這里分別以球體離散元和擴展多面體離散元方法為例,介紹其對錐體結構冰載荷的模擬情況,采用渤海海洋平臺的實測數據、HSVA 模型試驗和ISO 標準對離散元結果進行對比驗證[26,18,48].

首先采用HSVA 內開展的錐體結構與平整冰的模型試驗對球體離散元方法進行驗證[3].離散元計算參數同模型試驗完全相同,其中海冰模型中彎曲強度為60.6 kPa,冰厚為32 mm,冰速為100 mm/s.HSVA 模型試驗與離散元模擬結果如圖27 所示,由此得到的冰載荷時程如圖28 所示.離散元模擬的冰載荷峰值的平均值與漢堡實測數據相一致,且均表現出明顯的周期性,并具有相近的冰載荷頻率[48].

圖27 海冰與錐體作用的模型試驗及離散元模擬[48]Fig.27 Model test and DEM simulation of interaction between sea ice and conical structure[48]

圖28 漢堡試驗與離散元模擬的冰載荷時程曲線對比Fig.28 Comparison of time history of ice loads in HSVA model test and DEM simulations

海冰的斷裂長度是表征其破壞模式的重要參數,還與錐體冰載荷的幅值和周期密切相關.為分析斷裂長度的影響因素,將離散元計算結果同時與HSVA 試驗和渤海現場實測結果進行對比分析.依據錐體結構的離散元模擬結果,對其冰力峰值進行了統計分析,并與相關的靜冰力理論結果、渤海現場實測結果進行了對比分析,如圖29 所示.該錐體結構靜冰力可表述為

圖29 錐體結構的靜冰力[48]Fig.29 Static ice loads on conical structure[48]

式中,σf為海冰彎曲強度;Dw為錐體直徑;hi為冰厚;Lb為斷裂長度.該錐體結構靜冰力公式不僅考慮了上述參數的影響,并與現場實測和模型試驗結果相一致.

采用塊體離散元方法模擬錐體結構的冰載荷,其計算結果如圖30(a)所示.作為對比,渤海JZ20-2 MUQ 海洋平臺現場情況如圖30(b)所示.從中可以看出,平整冰在與錐體結構接觸后發生彎曲破壞,破碎海冰會沿著錐體向上攀爬,從而在冰載荷時程曲線中形成規律性的峰值載荷.采用不同冰厚模擬平整冰與錐體結構的相互作用,通過ISO 標準和渤海現場實測數據分析計算結果的合理性.圖30(c)和圖30(d)分別是離散元模擬結果同ISO 標準、渤海的現場監測數據對比.從圖中可以看出,離散元計算結果在趨勢上與ISO 標準和現場監測結果保持一致,數值上基本處于ISO 標準和現場監測數據之間.由此驗證了擴展多面體離散元方法在海冰與海洋結構物相互作用模擬中的合理性.

6 數值冰水池的工程應用

基于以上分析,數值冰水池研究包括海冰數值模型、海冰與流體、工程結構相互作用的多介質多尺度數值算法、模型試驗和現場觀測的試驗驗證、數據庫及前后處理系統等多方面的研究內容.目前,我國在以上相關研究中取得了很好的初步成果,為后續建立完善的數值冰水池系統打下了很好的基礎.下面分別從船舶結構冰阻力及冰載荷、海洋工程結構冰載荷及力學響應等方面,對目前相關的數值冰水池研究基礎及工程應用進行簡要介紹.

6.1 船舶結構冰阻力及冰載荷的數值模擬

在北極航行中,采用破冰船為貨輪引航的方式可有效降低成本,并提高極地航行安全保障的專業性.中遠海運集團在北極貨運航行中的破冰船引航情況如圖31 所示.這里采用中國“雪龍”號科學考察船作為引航破冰船,以中遠某商船為計算模型.貨船的船長、船寬和吃水分別為破冰船的1.20 倍、1.23倍和1.37 倍,兩船的航速相同.圖32 為航速2.5 m/s、冰厚1.0 m 條件下破冰船引航的模擬結果.從圖中可以看出,破冰船的破冰過程與單船破冰沒有明顯區別.貨船在破冰船開辟的開闊水道中航行,與水道中的碎冰發生相互作用,不會發生明顯的海冰破碎現象.

圖32 破冰船引航下船舶在平整冰區航行的離散元模擬 [130]Fig.32 DEM simulation of ship navigation in level ice with icebreaker pilotage [130]

圖33(a)是引航條件下破冰船和貨船的冰載荷時程曲線,虛線是穩定階段冰載荷的均值.可以看出,引航的破冰船冰載荷與單船破冰條件下的冰載荷類似.船體完全進入冰區后,其冰載荷時程在穩定的水平附近上下波動.貨船則與破冰船差別較大,進入冰區后,其冰載荷沒有穩定的上升階段,而是會出現類似脈沖形態的波動,但是其整體水平比破冰船小.圖33(b)是無引航條件下貨船的冰載荷時程曲線,虛線是穩定階段的冰阻力.無引航條件下貨船冰阻力為5.9 MN,大于引航條件下的貨船冰阻力.但是由于貨輪不具備破冰能力,其艏柱傾角和進水角較小,水線處存在較為尖銳的艏部結構,容易造成海冰的擠壓破碎,這與破冰船艏部大傾角造成海冰彎曲破壞不同.因此,無引航條件下貨船冰載荷存在近100 MN 的峰值載荷,對貨船的結構安全構成嚴重威脅.由此可見,在破冰船引航條件下,貨船的冰阻力可明顯降低,且不會存在較大的峰值載荷.

圖33 有無破冰船引航條件下船舶結構冰阻力時程曲線Fig.33 Time history of ice resistances on ship hull with or without icebreaker pilotage

6.2 海洋工程冰載荷及結構力學響應數值模擬

為分析不同冰況下導管架海洋平臺結構的冰載荷和振動響應,采用具有粘接-破碎性能的等粒徑球體離散單元對海冰的破碎特性進行模擬,通過由梁單元構建的海洋平臺有限元模型獲得結構的振動響應.在離散元與有限元的接觸區域中實現了兩個模型間計算參數的傳遞.為提高該耦合模型的計算效率和規模,發展了基于動力子結構方法的DEM-FEM耦合模型[43].

為模擬平整冰的破碎特性,這里將具有粘接-破碎功能的等粒徑球體單元規則排列構成平整冰離散元模型.渤海四樁腿JZ20-2 MUQ 錐體導管架式海洋平臺主要由三部分組成,即上部模塊、導管架和樁基.在進行有限元動力分析時,為簡化計算,在保證主體結構幾何形狀以及結構振動頻率和振型真實性的基礎上,可對平臺結構的有限元模型進行一定的簡化,這里將上部結構簡化為梁單元同時樁基用等效6 倍的樁徑代替.為提高動力分析的計算效率,這里采用動力子結構法中的約束模態綜合法對結構的自由度進行縮減,并將該海洋平臺結構劃分為兩個子結構,如圖34 所示.采用DEM-FEM 耦合模型分析平臺冰激振動時,兩模型間的參數傳遞尤為重要.這里,將海冰離散單元對導管架式海洋平臺結構的冰載荷作為力邊界條件傳遞到有限元模型,由此計算海洋平臺的動力響應,再進一步將更新后的平臺位移作為離散元的位移邊界條件.在海洋平臺結構的有限元計算中,結構振動響應采用逐步積分法中的Newmark 方法計算.

圖34 渤海JZ20-2 MUQ 導管架式海洋平臺及其有限元模型Fig.34 The jacket offshore platform JZ20-2 MUQ and its FEM model in Bohai Sea

基于GPU 并行算法,對不同冰速和冰厚下的冰載荷和海洋平臺結構的振動響應進行了計算,由此得到的海冰與錐體海洋平臺的相互作用過程如圖35所示.該模擬結果與現場觀測到的現象有很好的相似性,并且可以看到平整冰在錐體作用下的破碎情況,即冰排呈現出初次斷裂、爬升、二次斷裂和清除的過程,并由此引起交變動冰力.當冰速0.2 m/s,冰厚0.2 m 時,海洋平臺樁腿水平方向的冰力時程如圖36 所示,其具有很強的隨機性.這與海冰現場監測和室內模型試驗結果相一致.為驗證該DEMFEM 耦合模型的可行性與適用性,這里將模擬得到的冰載荷峰值分別與ISO19906 標準以及我國《港口工程荷載規范》JTS 144?1?2010 標準進行了對比,如圖37 所示.可以發現,冰載荷與冰厚呈近似的二次非線性遞增關系.

圖35 DEM-FEM 模擬的海冰與JZ20-2 MUQ 平臺相互作用過程Fig.35 Interaction between sea ice and platform JZ20-2 MUQ simulated with DEM-FEM

圖36 離散元計算的樁腿冰載荷時程Fig.36 Ice forces on the platform simulated with DEM

圖37 DEM-FEM 耦合模型的冰載荷峰值與ISO19906,JTS 1441-1-2010 標準對比Fig.37 Comparison of peak ice forces in coupled DEM-FEM model with ISO19906 and JTS 144-1-2010 standards

現場實測和數值模擬得到的平臺冰激振動加速度時程如圖38 所示.從中可以發現兩者具有較好的相似性.根據渤海JZ20-2 MUQ 平臺現場冰激振動加速度測量數據的統計,發現冰激振動加速度與冰速和冰厚平方的乘積呈線性關系,對不同冰況下冰激結構響應計算結果的擬合如圖39 所示,其中給出了相應的現場測量結果.從中可以看到實測數據和數值模擬的結果較為接近.由此可見,DEM-FEM 耦合模型可以揭示渤海錐體導管架海洋平臺冰激振動加速度與冰速和冰厚的對應關系.

圖38 基于DEM-FEM 耦合方法的JZ20-2 MUQ 導管架式海洋平臺冰激振動結果分析[43]Fig.38 Analysis of ice-inducted vibration of jacket offshore platform JZ20-2 MUQ based on coupled DEM-FEM[43]

圖39 冰速、冰厚與JZ20-2 MUQ 平臺結構冰激振動加速度的關系[43]Fig.39 Relationship between ice velocity,ice thickness and iceinduced vibration of platform JZ20-2 MUQ[43]

此外,采用DEM-FEM 耦合模型,對渤海JZ20-2 NW 單錐海洋平臺結構的冰載荷和冰激振動特性進行了數值分析,并與渤海冰區現場實測結果進行了對比驗證[51].該海洋平臺的錐體部分采用平板型殼單元構造,其整體構架及錐體內部的加筋肋采用梁單元構造,如圖9 和圖10 所示.通過DEM-FEM耦合模型計算得到了海冰作用下錐體結構的Von-Mises 應力,如圖40 所示.圖40(a)為結構的Von-Mises 應力云圖,圖40(b)為相對應的法向壓力云圖.從中可以發現錐體結構應力分布與壓力分布相對應,壓力最大處也是應力最大的部分.此外,對不同冰厚和冰速下錐體結構應力的最大值進行統計分析,由此確定錐體結構發生破壞時的海冰條件.該DEM-FEM 耦合模型不僅可在細觀尺度上分析海冰與結構相互作用中的破壞狀態,還可揭示結構冰激振動的內在機理,可為冰區結構的冰載荷、結構動力響應分析提供有效的數值手段.

圖40 冰載荷作用下錐體結構的應力及壓力分布[51]Fig.40 Mises stress and pressure distributions of the conical structure under ice load[51]

6.3 數值冰水池與物理冰水池的比對研究

物理冰水池中開展的冰區船舶及海洋結構物的物理模型試驗為極區船舶的船型設計、性能預報和結構強度評估提供技術支撐[133].模型試驗設計與實際船舶航行工況通常采用佛汝德相似準則和柯西相似準則,其縮尺比 λ 考慮冰物理力學性質與船舶模型尺寸,表1 列出模型試驗中各個變量的縮尺關系.

表1 海冰模型試驗中的主要物理量比尺Table 1 Scale ratio of primary physical parameters in sea ice model tests

數值冰水池的一個重要功能是對物理冰水池的虛擬再現,其數值試驗設計繼承了模型試驗中物理量的縮尺,各個試驗設計變量包括了海冰類型(平整冰、碎冰和冰脊等)和航行狀態(直航、回轉).此外,數值試驗作為物理冰水池模型試驗的補充,低成本的大量重復性數值試驗可以避免物理試驗中不可控因素的影響.參考國內外冰水池,采用離散元方法建立數值冰水池模型,圖41 給出了HSVA 和對應建立的數值冰水池模型.

圖41 物理冰水池與數值冰水池Fig.41 Physical ice tank and numerical ice tank

6.3.1 碎冰區航行模擬

船舶碎冰區航行的模型試驗重點在于碎冰場的生成,經典的Voronoi 切割方法優勢在于可定量控制碎冰塊幾何形狀,保證分割的碎冰形狀符合隨機分布到規則分布的連續變換,并且分割后形狀與天然海冰存在高度的幾何形態相似性[86].采用Voronoi 切割方法可以定量(冰塊面積、密集度和幾何規則度)生成數值冰水池中的碎冰場并相繼開展碎冰場幾何特性對數值試驗結果的影響研究.但是Voronoi 切割方法所構造的碎冰場與物理冰水池碎冰場依然存在差距.相關研究表明,碎冰初始位置和大小差異性對于船舶的模型試驗具有一定影響[12].為了更加真實地反應真實情況下冰水池內海冰初始位置和大小,基于圖像處理方法實現數值冰水池的碎冰場的快速生成,圖42 分別給出了由Voronoi 方法和數值圖像方法所生成的碎冰場與物理冰水池中的碎冰場.

圖42 數值冰水池與物理冰水池中的碎冰場Fig.42 Broken ice field in numerical and physical ice tank

圖43 給出了在數值冰水池中模擬航速5 kn 的某極地船舶在冰厚1.2 m 碎冰區的航行狀態,其模型縮尺采用1∶25,即航速和冰厚分別為0.10 m/s和0.048 m.船舶在直航過程與大塊碎冰發生碰撞,海冰被船舶排開并發生彎曲/擠壓破壞,船艉后側出現明顯的冰間水道.圖44 為船舶結構與浮冰塊相互作用時數值模型試驗和物理模型試驗的對比.大塊碎冰與船艏作用而發生彎曲破壞,船肩處海冰發生局部擠壓.船舶碎冰區航行中在3 個方向上選取的50 s 穩定階段冰載荷如圖45 所示,由此可見冰載荷具有很強的脈沖性,這種脈沖性載荷是由大塊碎冰與船碰撞發生破壞而不能持續作用所引起.

圖43 數值冰水池中船與碎冰相互作用模擬Fig.43 Simulation of interaction between ship and float ice in numerical ice tank

圖44 數值與物理模型試驗中海冰局部破壞現象對比[12]Fig.44 Comparison of local damage of sea ice in numerical and physical model tests[12]

圖45 浮冰區船體結構冰載荷時程Fig.45 Time history of ice loads on ship hull in broken ice area

6.3.2 平整冰區航行模擬

對于船舶在平整冰區航行模型試驗的數值模擬,數值冰水池中海冰初始場采用球體單元規則排列快速生成,并參考天津大學冰水池的相關試驗參數[134].圖46 為模擬破冰船在平整冰區航行過程,其試驗條件為航速5 kn,冰厚1.2 m,數值冰水池試驗設計縮尺為1∶25.為了驗證數值模擬的合理性,圖47 給出了模型試驗中海冰破碎現象與數值模型的對比情況.在船舶破冰航行中,船-冰碰撞使海冰發生彎曲和擠壓破壞,其中在船肩兩側主要發生彎曲破壞,彎曲破壞形成的海冰尺寸較大;海冰的擠壓破壞主要發生在艏柱附近區域,海冰破碎尺寸較小.此外,船舶航行過程中,在船舯兩側會出現明顯的環形裂紋,船艉后側也出現明顯的冰間水道.圖48 為數值計算得到的冰阻力時程曲線.它與碎冰區航行冰阻力時程不同,平整冰區航行冰載荷具有很好的穩定階段,但是也會出現脈沖形態波動.這主要是由船-冰作用過程中海冰彎曲破壞所導致.

圖46 數值冰水池中船與平整冰相互作用模擬Fig.46 Simulation of interaction between ship and level ice in numerical ice tank

圖47 數值冰水池與物理冰水池試驗現象對比[134]Fig.47 Comparison of experimental phenomenon in numerical and physical ice tanks[134]

圖48 平整冰區船體結構冰阻力時程Fig.48 Time history of ice resistance on ship hull in level ice

通過對冰水池中船舶與海洋工程結構冰載荷模型試驗的數值分析,可對冰載荷離散元方法及計算分析軟件的可靠性進行驗證.在后續研究中將進一步考慮不同海冰類型和結構形式,對海冰破壞模式、冰載荷和結構響應開展更加系統的對比驗證,從而保證冰載荷數值模擬系統的可靠性.當然,對物理冰水池的數值模擬是數值冰水池研究的一部分,由此可對數值冰水池的可靠性進行驗證.在此基礎上,進一步發展不同尺度下海冰與船舶及海洋工程結構的耦合作用則可克服模型試驗尺度下的尺寸影響,進而服務于原型尺度下的結構冰載荷數值預測.

7 結語

數值冰水池的建設是當前國內外極地船舶與海洋工程發展的需求,也依托于海冰物理力學性質、冰載荷特性、高性能數值方法和計算技術等相關學科的發展,同時又與冰載荷現場測量和冰水池模型試驗驗證密切相關.本文針對數值冰水池的發展思路、框架、研究目標和內容進行了相應的介紹,并重點對海冰物理力學性質的數值建模、冰-水-工程結構的耦合數值算法、數值冰水池的軟件實現和試驗驗證進行了較為詳細的論述.最后,針對典型的工程應用實例,對采用離散元方法對極地船舶和導管架式海洋平臺結構的冰載荷和力學響應進行了數值分析和驗證,并與相應的冰水池模型試驗進行了初步的對比分析.作為一個具有獨立知識產權的軟件系統平臺,數值冰水池將與冰載荷的理論分析、現場測量和冰水池物理模型試驗研究相互補充、密切結合,共同解決極地船舶與海洋工程結構的抗冰設計、強度分析、疲勞評估和安全保障問題.

本文初步提出了基于離散元方法的數值冰水池基本概念、研究思路和初步的工程應用,在后續建設中還將密切結合具體的極地船舶與海洋工程問題進行完善,形成具有獨立知識產權的數值冰水池計算分析軟件系統.將數值冰水池研發與中國冰水池建設相結合,不斷發展和采用相關的先進技術,可有力地促進極地船舶與海洋工程的發展.

8 致謝

本文工作得到中國船級社、哈爾濱工程大學、天津大學、大連海事大學、國家海洋環境預報中心、中國極地研究中心、中船集團第708 研究所和第719 研究所、美國Clarkson 大學、美國船級社、新加坡國立大學、韓國船舶與海洋工程研究所等國內外單位諸多專家學者的支持和幫助,在此一并致謝;感謝大連理工大學劉璐、楊冬寶、吳捷等人對極地船舶與海洋工程結構冰載荷的離散元數值計算.