• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2021-11-10 09:48:50季順迎田于逵
    力學(xué)學(xué)報 2021年9期
    關(guān)鍵詞:海冰海洋工程模型試驗

    季順迎 田于逵 )

    * (大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,遼寧大連 116024)

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

    引言

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

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

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

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

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

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

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

    1 面向冰載荷的數(shù)值冰水池研究框架

    下面提出面向極地船舶與海洋工程結(jié)構(gòu)冰載荷預(yù)報的數(shù)值冰水池基本框架,簡要論述數(shù)值冰水池的主要研究目標(biāo)、大體研究思路、主要研究內(nèi)容和關(guān)鍵技術(shù)等,并將在后續(xù)研究中不斷改進和完善.

    1.1 數(shù)值冰水池的基本概念和研究目標(biāo)

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

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

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

    1.2 面向冰載荷預(yù)報的數(shù)值冰水池基本框架

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

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

    下面對數(shù)值冰水池的主要研究內(nèi)容進行簡要介紹,并分別討論其所面臨的主要關(guān)鍵技術(shù).

    (1) 海冰物理力學(xué)性質(zhì)的數(shù)值模型.創(chuàng)建可準(zhǔn)確描述極地海冰物理力學(xué)特征、本構(gòu)模型和強度準(zhǔn)則的不同數(shù)值方法和關(guān)鍵計算參數(shù)的選取方案,并形成平整冰、碎冰區(qū)和冰脊等不同冰類型的構(gòu)造算法.

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

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

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

    (5) 數(shù)值冰水池模擬的試驗驗證系統(tǒng)及可靠性檢驗.基于典型冰水池物理模型試驗和冰區(qū)現(xiàn)場測量數(shù)據(jù),建立數(shù)值冰水池的不同尺度試驗驗證系統(tǒng),提高數(shù)值模擬結(jié)果的可靠性和工程應(yīng)用性.系統(tǒng)地收集整理國內(nèi)外相關(guān)測量結(jié)果,建立用于數(shù)值冰水池可靠性檢驗的數(shù)據(jù)庫,并考慮模型試驗和現(xiàn)場測量中測量結(jié)果的不確定性,發(fā)展相應(yīng)的可靠性驗證技術(shù);重點構(gòu)建用于數(shù)值冰水池驗證的基準(zhǔn)試驗,實現(xiàn)對其全面系統(tǒng)的可靠性驗證.計算結(jié)果的可靠性是數(shù)值冰水池的核心內(nèi)容,也是其實現(xiàn)工程應(yīng)用服務(wù)的基礎(chǔ)和前提條件.

    (6) 典型極地船舶與海洋工程技術(shù)開發(fā)的工程應(yīng)用.實現(xiàn)典型極地船舶與海洋工程結(jié)構(gòu)的工程應(yīng)用,通過極地裝備結(jié)構(gòu)冰載荷和冰激結(jié)構(gòu)振動的數(shù)值模擬結(jié)果,為結(jié)構(gòu)抗冰設(shè)計提供可行的優(yōu)化方案,為結(jié)構(gòu)疲勞分析和結(jié)構(gòu)完整性管理提供有力的數(shù)據(jù)支持和工程服務(wù).

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

    2 海冰物理力學(xué)性質(zhì)的數(shù)值方法

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

    2.1 海冰的典型物理力學(xué)性質(zhì)

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

    圖2 北極海冰的冰晶細觀結(jié)構(gòu)[70]Fig.2 Micro-crystal structure of sea ice in the Arctic[70]

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

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

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

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

    圖4 海冰對直立結(jié)構(gòu)擠壓破碎的高壓區(qū)[78]Fig.4 Sketch of high pressure zones on vertical structure[78]

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

    2.2 海冰力學(xué)行為的數(shù)值模擬方法

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

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

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

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

    圖8 不同數(shù)值方法模擬的冰-結(jié)構(gòu)相互作用Fig.8 Interaction between sea ice and structure simulated with various numerical methods

    3 冰-水-工程結(jié)構(gòu)耦合作用的多介質(zhì)、多尺度數(shù)值方法

    在海冰與船舶及海洋工程結(jié)構(gòu)的相互作用過程中,海冰的破壞模式、運動規(guī)律及由此導(dǎo)致的冰載荷特性均與流體、結(jié)構(gòu)力學(xué)響應(yīng)密切相關(guān),因此在數(shù)值冰水池研究中需要考慮海冰與流體介質(zhì)、工程結(jié)構(gòu)的耦合效應(yīng)[96].這就需要在海冰數(shù)值方法的基礎(chǔ)上,進一步發(fā)展冰-水-結(jié)構(gòu)間的多介質(zhì)、多尺度耦合數(shù)值模型,以提高冰載荷及結(jié)構(gòu)力學(xué)響應(yīng)的計算精度.以上研究已得到國內(nèi)外相關(guān)學(xué)者高度重視.

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

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

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

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

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

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

    3.2 冰載荷的多介質(zhì)計算模型

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

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

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

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

    通過對海冰與海水、船舶及海洋工程結(jié)構(gòu)間耦合作用的多介質(zhì)、多尺度計算分析,可對冰激結(jié)構(gòu)振動、形變等力學(xué)響應(yīng)對海冰破壞模式的影響,進而更準(zhǔn)確地計算冰載荷的動力特性.

    4 數(shù)值冰水池的軟件實現(xiàn)

    類似于數(shù)值水池的知識封裝、可靠性和情景化三大技術(shù)特征[66],數(shù)值冰水池也需要通過高性能計算機軟件實現(xiàn)相應(yīng)的知識封裝和情景化.為此,需要發(fā)展相應(yīng)的前后處理系統(tǒng)、工程結(jié)構(gòu)數(shù)據(jù)庫和三維可視化顯示技術(shù),并通過發(fā)展相應(yīng)先進的并行算法實現(xiàn)高性能數(shù)值計算.

    4.1 典型海冰與工程結(jié)構(gòu)數(shù)據(jù)庫

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

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

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

    (3) 對于簡單的錐體或直立結(jié)構(gòu),可建立相應(yīng)的結(jié)構(gòu)類型,并通過調(diào)整結(jié)構(gòu)參數(shù)進行直接設(shè)定;對于相對復(fù)雜的工程結(jié)構(gòu),則可構(gòu)建與相關(guān)通用商業(yè)軟件的數(shù)據(jù)接口,從而導(dǎo)入船舶及海洋工程結(jié)構(gòu)數(shù)據(jù).對于船舶及海洋工程結(jié)構(gòu)的約束、運動方式及速度,也可通過交互式界面設(shè)定.對于典型的極地船舶及海洋工程結(jié)構(gòu),則建立相應(yīng)的數(shù)據(jù)庫,以便于直接調(diào)用、計算.

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

    圖12 離散元軟件IceDEM 中海冰及工程結(jié)構(gòu)參數(shù)設(shè)定[84]Fig.12 Settings of sea ice and structure parameters in the software IceDEM[84]

    4.2 面向工程尺度的高性能數(shù)值算法

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

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

    4.3 前后處理系統(tǒng)及三維可視化技術(shù)

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

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

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

    圖15 平整冰與錐體海洋平臺上部、中部和下部作用時的離散元模擬三維再現(xiàn)[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 船舶在碎冰區(qū)航行的離散元模擬三維再現(xiàn)[12]Fig.16 3D reconstruction of DEM simulation of ship navigation in broken ice area[12]

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

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

    5 數(shù)值冰水池的試驗驗證

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

    5.1 海冰物理力學(xué)性質(zhì)的試驗驗證

    在海冰的離散元、有限元、黏聚單元方法、近場動力學(xué)等不同海冰數(shù)值方法研究中,均注重對海冰物理力學(xué)性質(zhì)的試驗驗證,以建立合理的海冰本構(gòu)模型和強度準(zhǔn)則、發(fā)展理想的數(shù)值方法并確定相應(yīng)的計算參數(shù)[13,128].在海冰的物理力學(xué)性質(zhì)數(shù)值模擬中開展最廣泛的是海冰的單軸壓縮和三點彎曲[52,56,92,129].這里采用離散元方法對顆粒排列、粒徑、摩擦系數(shù)等參數(shù)影響下的數(shù)值模擬結(jié)果進行分析.對于海冰的單軸壓縮試驗,海冰試件尺寸(a×a×h)設(shè)為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]

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

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

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

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

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

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

    圖21 雪龍?zhí)柶票c平整冰相互作用的離散元模擬[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 船體冰載荷的離散元計算結(jié)果與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)對特定工況下船體結(jié)構(gòu)的冰載荷給出了規(guī)范計算方法[131].為驗證離散元方法對船體結(jié)構(gòu)與大塊浮冰碰撞過程中冰載荷計算的可靠性,劉璐等[132]將離散元計算的冰壓力與IACS 規(guī)范進行了對比.計算情景如圖24 所示,碰撞點到艏柱的距離分別為5 m,15 m,25 m 和35 m.當(dāng)航速為2.25 m/s,海冰厚度為3.0 m,海冰彎曲強度為1.0 MPa 時,船體結(jié)構(gòu)與大塊浮冰碰撞的離散元模擬過程如圖25 所示.圖26 為4 個算例中碰撞點處的冰壓力與規(guī)范對比情況,其中x為碰撞點到艏柱的距離.計算結(jié)果表明,離散元結(jié)果與規(guī)范值相比存在上下誤差,但壓力值保持在相同量級范圍內(nèi),誤差范圍為6.7%~ 18.1%之間.由此可見,離散元方法對船體結(jié)構(gòu)的冰壓力可進行準(zhǔn)確可靠的分析,具有良好的工程適用性.

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

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

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

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

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

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

    圖29 錐體結(jié)構(gòu)的靜冰力[48]Fig.29 Static ice loads on conical structure[48]

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

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

    6 數(shù)值冰水池的工程應(yīng)用

    基于以上分析,數(shù)值冰水池研究包括海冰數(shù)值模型、海冰與流體、工程結(jié)構(gòu)相互作用的多介質(zhì)多尺度數(shù)值算法、模型試驗和現(xiàn)場觀測的試驗驗證、數(shù)據(jù)庫及前后處理系統(tǒng)等多方面的研究內(nèi)容.目前,我國在以上相關(guān)研究中取得了很好的初步成果,為后續(xù)建立完善的數(shù)值冰水池系統(tǒng)打下了很好的基礎(chǔ).下面分別從船舶結(jié)構(gòu)冰阻力及冰載荷、海洋工程結(jié)構(gòu)冰載荷及力學(xué)響應(yīng)等方面,對目前相關(guān)的數(shù)值冰水池研究基礎(chǔ)及工程應(yīng)用進行簡要介紹.

    6.1 船舶結(jié)構(gòu)冰阻力及冰載荷的數(shù)值模擬

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

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

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

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

    6.2 海洋工程冰載荷及結(jié)構(gòu)力學(xué)響應(yīng)數(shù)值模擬

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

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

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

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

    圖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 標(biāo)準(zhǔn)對比Fig.37 Comparison of peak ice forces in coupled DEM-FEM model with ISO19906 and JTS 144-1-2010 standards

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

    圖38 基于DEM-FEM 耦合方法的JZ20-2 MUQ 導(dǎo)管架式海洋平臺冰激振動結(jié)果分析[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 平臺結(jié)構(gòu)冰激振動加速度的關(guān)系[43]Fig.39 Relationship between ice velocity,ice thickness and iceinduced vibration of platform JZ20-2 MUQ[43]

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

    圖40 冰載荷作用下錐體結(jié)構(gòu)的應(yīng)力及壓力分布[51]Fig.40 Mises stress and pressure distributions of the conical structure under ice load[51]

    6.3 數(shù)值冰水池與物理冰水池的比對研究

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

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

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

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

    6.3.1 碎冰區(qū)航行模擬

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

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

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

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

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

    圖45 浮冰區(qū)船體結(jié)構(gòu)冰載荷時程Fig.45 Time history of ice loads on ship hull in broken ice area

    6.3.2 平整冰區(qū)航行模擬

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

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

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

    圖48 平整冰區(qū)船體結(jié)構(gòu)冰阻力時程Fig.48 Time history of ice resistance on ship hull in level ice

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

    7 結(jié)語

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

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

    8 致謝

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

    猜你喜歡
    海冰海洋工程模型試驗
    末次盛冰期以來巴倫支海-喀拉海古海洋環(huán)境及海冰研究進展
    海洋通報(2021年3期)2021-08-14 02:20:38
    反推力裝置模型試驗臺的研制及驗證
    基于SIFT-SVM的北冰洋海冰識別研究
    海洋工程專家 劉培林
    臺階式短加筋土擋墻行為特征的離心模型試驗
    巨厚堅硬巖漿巖不同配比的模型試驗研究
    《海洋工程》第二屆理事會
    海洋工程(2015年1期)2015-10-28 01:36:21
    海洋工程學(xué)會第四屆理事會
    海洋工程(2015年1期)2015-10-28 01:29:14
    電滲—堆載聯(lián)合氣壓劈烈的室內(nèi)模型試驗
    應(yīng)用MODIS數(shù)據(jù)監(jiān)測河北省近海海域海冰
    河北遙感(2014年4期)2014-07-10 13:54:59
    亚洲精华国产精华液的使用体验| 国产成人午夜福利电影在线观看| 国产午夜精品一二区理论片| 岛国毛片在线播放| 久久亚洲国产成人精品v| 男人爽女人下面视频在线观看| 卡戴珊不雅视频在线播放| 97在线视频观看| 久久韩国三级中文字幕| 日韩 亚洲 欧美在线| 三级国产精品片| 精品人妻熟女av久视频| 亚洲综合色惰| 午夜日本视频在线| 这个男人来自地球电影免费观看 | 日本爱情动作片www.在线观看| 精品久久久精品久久久| 亚洲国产最新在线播放| 中国三级夫妇交换| 国产成人精品久久久久久| 女人久久www免费人成看片| 国产永久视频网站| 亚洲怡红院男人天堂| 欧美少妇被猛烈插入视频| 热re99久久国产66热| 丝袜脚勾引网站| 91久久精品国产一区二区成人| 大码成人一级视频| 欧美日韩在线观看h| 国产欧美日韩综合在线一区二区| a级毛片免费高清观看在线播放| 男男h啪啪无遮挡| 老司机影院成人| 中文字幕免费在线视频6| 97在线人人人人妻| 大香蕉97超碰在线| 国产精品三级大全| 亚洲国产av影院在线观看| 午夜久久久在线观看| 一二三四中文在线观看免费高清| 狂野欧美激情性bbbbbb| 三上悠亚av全集在线观看| 日韩一区二区视频免费看| 少妇的逼水好多| 国产精品久久久久久精品古装| 亚洲高清免费不卡视频| 午夜91福利影院| 青春草国产在线视频| 国产日韩欧美在线精品| 国产白丝娇喘喷水9色精品| 国产免费福利视频在线观看| 国产精品秋霞免费鲁丝片| 啦啦啦在线观看免费高清www| 日韩精品有码人妻一区| 成人毛片a级毛片在线播放| 国产精品欧美亚洲77777| 人人澡人人妻人| 国产精品免费大片| 欧美亚洲 丝袜 人妻 在线| 99精国产麻豆久久婷婷| 成人亚洲精品一区在线观看| 只有这里有精品99| 国产国拍精品亚洲av在线观看| 高清不卡的av网站| 国产69精品久久久久777片| 尾随美女入室| 人成视频在线观看免费观看| 久久狼人影院| 夜夜骑夜夜射夜夜干| 午夜福利视频精品| 午夜91福利影院| 亚洲图色成人| 99视频精品全部免费 在线| 亚洲国产精品专区欧美| 国语对白做爰xxxⅹ性视频网站| 国产在线一区二区三区精| 国产精品一区二区三区四区免费观看| 少妇高潮的动态图| 人妻制服诱惑在线中文字幕| 女人久久www免费人成看片| 天堂中文最新版在线下载| xxx大片免费视频| 在现免费观看毛片| 亚洲精品久久成人aⅴ小说 | 成人影院久久| 婷婷色av中文字幕| 久久国产精品男人的天堂亚洲 | 中文字幕人妻熟人妻熟丝袜美| 伊人久久国产一区二区| 你懂的网址亚洲精品在线观看| 久久久久视频综合| 久久久久精品性色| 在线播放无遮挡| 欧美精品一区二区免费开放| 一区二区三区四区激情视频| 免费看av在线观看网站| 黄色欧美视频在线观看| 国产不卡av网站在线观看| 亚洲国产日韩一区二区| 国产精品久久久久久久电影| 自拍欧美九色日韩亚洲蝌蚪91| 熟女人妻精品中文字幕| a级片在线免费高清观看视频| 大香蕉久久成人网| 9色porny在线观看| 婷婷色综合www| 丰满少妇做爰视频| 三级国产精品片| 内地一区二区视频在线| 你懂的网址亚洲精品在线观看| 中文精品一卡2卡3卡4更新| 国产片特级美女逼逼视频| 午夜视频国产福利| 免费观看a级毛片全部| 在线观看免费高清a一片| 美女视频免费永久观看网站| 免费黄频网站在线观看国产| 男女无遮挡免费网站观看| 少妇被粗大的猛进出69影院 | 极品少妇高潮喷水抽搐| 日本-黄色视频高清免费观看| 成人手机av| 成人漫画全彩无遮挡| 欧美 日韩 精品 国产| 久久av网站| 欧美三级亚洲精品| 美女大奶头黄色视频| 观看av在线不卡| 黑人高潮一二区| 黄色视频在线播放观看不卡| 精品少妇黑人巨大在线播放| 七月丁香在线播放| 九色亚洲精品在线播放| 99国产精品免费福利视频| 成人国产麻豆网| 美女福利国产在线| 久久久久国产精品人妻一区二区| 亚洲av不卡在线观看| 一个人免费看片子| 大香蕉久久网| 18禁裸乳无遮挡动漫免费视频| 亚洲四区av| 国产成人午夜福利电影在线观看| 最新中文字幕久久久久| 在现免费观看毛片| 91精品一卡2卡3卡4卡| 一区二区三区精品91| 日日撸夜夜添| 狂野欧美激情性xxxx在线观看| 国产成人a∨麻豆精品| 精品人妻熟女av久视频| 国产探花极品一区二区| www.色视频.com| 欧美日韩av久久| 色5月婷婷丁香| 成人国产麻豆网| 免费观看的影片在线观看| 一区二区日韩欧美中文字幕 | 国产爽快片一区二区三区| 一级爰片在线观看| 国产免费一级a男人的天堂| 在线观看国产h片| 蜜桃久久精品国产亚洲av| 亚洲无线观看免费| 高清在线视频一区二区三区| 国产伦理片在线播放av一区| 22中文网久久字幕| 国产黄片视频在线免费观看| 又黄又爽又刺激的免费视频.| 男女啪啪激烈高潮av片| 91精品国产国语对白视频| 午夜福利视频在线观看免费| 国产精品国产三级专区第一集| 久久久精品94久久精品| 精品卡一卡二卡四卡免费| 春色校园在线视频观看| 热re99久久精品国产66热6| 国产免费又黄又爽又色| 超色免费av| 一边亲一边摸免费视频| 免费观看的影片在线观看| 日韩欧美一区视频在线观看| 69精品国产乱码久久久| 免费日韩欧美在线观看| 欧美亚洲日本最大视频资源| 大又大粗又爽又黄少妇毛片口| 国产日韩欧美亚洲二区| 秋霞伦理黄片| 国产av精品麻豆| 少妇人妻久久综合中文| 黄色视频在线播放观看不卡| videosex国产| 国产亚洲最大av| 亚洲人成网站在线播| 久久99一区二区三区| 色视频在线一区二区三区| 狠狠精品人妻久久久久久综合| 国产成人一区二区在线| 国产av精品麻豆| 蜜臀久久99精品久久宅男| 色哟哟·www| 天堂俺去俺来也www色官网| 老司机影院毛片| √禁漫天堂资源中文www| 特大巨黑吊av在线直播| 亚洲av欧美aⅴ国产| 国产av国产精品国产| 成人综合一区亚洲| 日本av免费视频播放| 成人18禁高潮啪啪吃奶动态图 | 九色亚洲精品在线播放| 国产日韩欧美亚洲二区| 边亲边吃奶的免费视频| 91精品国产国语对白视频| 在线观看免费视频网站a站| 精品国产乱码久久久久久小说| 精品一区二区三区视频在线| 少妇高潮的动态图| 国产片内射在线| 妹子高潮喷水视频| 亚洲精品一二三| 精品熟女少妇av免费看| 日本色播在线视频| 免费大片18禁| 欧美97在线视频| 美女主播在线视频| 国产精品99久久99久久久不卡 | 国产成人av激情在线播放 | 久久国产精品男人的天堂亚洲 | 亚洲成色77777| 一级a做视频免费观看| 国产爽快片一区二区三区| 久久热精品热| 亚洲综合精品二区| 欧美日本中文国产一区发布| 一边亲一边摸免费视频| 亚洲国产精品国产精品| 欧美xxⅹ黑人| 午夜激情久久久久久久| 亚洲欧美清纯卡通| 国产精品久久久久久久久免| 国产精品不卡视频一区二区| 亚洲国产av影院在线观看| 一级黄片播放器| 高清午夜精品一区二区三区| 亚洲色图综合在线观看| 国产成人午夜福利电影在线观看| 成人国语在线视频| 亚洲人成网站在线播| 亚洲精品乱码久久久久久按摩| 一区二区三区四区激情视频| 久久久久久久大尺度免费视频| 另类亚洲欧美激情| 91久久精品电影网| 2018国产大陆天天弄谢| 婷婷成人精品国产| 交换朋友夫妻互换小说| 美女内射精品一级片tv| 欧美+日韩+精品| 爱豆传媒免费全集在线观看| 熟妇人妻不卡中文字幕| 欧美日韩成人在线一区二区| 简卡轻食公司| 成人免费观看视频高清| 少妇丰满av| 天堂8中文在线网| 99热全是精品| 欧美+日韩+精品| 一区二区av电影网| 韩国高清视频一区二区三区| 22中文网久久字幕| 中文字幕精品免费在线观看视频 | 亚洲av中文av极速乱| 最后的刺客免费高清国语| 午夜免费男女啪啪视频观看| 超色免费av| 国产高清国产精品国产三级| 少妇 在线观看| 久久婷婷青草| freevideosex欧美| 最后的刺客免费高清国语| 亚洲在久久综合| 2018国产大陆天天弄谢| 午夜福利视频在线观看免费| 少妇人妻 视频| 韩国高清视频一区二区三区| av网站免费在线观看视频| 日韩欧美精品免费久久| 中国美白少妇内射xxxbb| 精品久久久久久久久亚洲| 亚洲精品成人av观看孕妇| 五月开心婷婷网| 一级二级三级毛片免费看| 美女中出高潮动态图| 国产成人免费观看mmmm| 亚洲精品456在线播放app| 中文字幕最新亚洲高清| 熟女人妻精品中文字幕| 国产成人精品无人区| 日日爽夜夜爽网站| 国产亚洲av片在线观看秒播厂| 色94色欧美一区二区| 极品人妻少妇av视频| 亚洲欧美清纯卡通| 成人午夜精彩视频在线观看| av又黄又爽大尺度在线免费看| 老熟女久久久| 亚洲av成人精品一二三区| 在线观看www视频免费| 大话2 男鬼变身卡| 欧美 亚洲 国产 日韩一| 国产在线免费精品| 亚洲色图综合在线观看| 国产精品不卡视频一区二区| 菩萨蛮人人尽说江南好唐韦庄| 国产在视频线精品| 丁香六月天网| 插逼视频在线观看| 少妇的逼好多水| 亚洲美女搞黄在线观看| 日本午夜av视频| a级毛色黄片| 91午夜精品亚洲一区二区三区| 五月天丁香电影| 日本免费在线观看一区| av电影中文网址| 成人亚洲精品一区在线观看| 国产成人a∨麻豆精品| 中国美白少妇内射xxxbb| 色婷婷av一区二区三区视频| 少妇猛男粗大的猛烈进出视频| videos熟女内射| 日本黄大片高清| 男女国产视频网站| 久久久精品免费免费高清| 国产欧美另类精品又又久久亚洲欧美| 久久久久久久亚洲中文字幕| 精品久久蜜臀av无| 国产高清三级在线| 日韩大片免费观看网站| 国产精品麻豆人妻色哟哟久久| 亚洲美女视频黄频| 国产在线一区二区三区精| 中文欧美无线码| av电影中文网址| 人妻 亚洲 视频| 成人手机av| 一区二区三区免费毛片| 美女中出高潮动态图| 一区二区av电影网| 久久久久久久久久久丰满| av福利片在线| 免费黄网站久久成人精品| 精品国产乱码久久久久久小说| 国产精品久久久久久久久免| 欧美日韩视频精品一区| 美女视频免费永久观看网站| 亚洲精品日韩av片在线观看| 国产日韩欧美视频二区| 九九在线视频观看精品| 午夜激情久久久久久久| 秋霞在线观看毛片| 麻豆成人av视频| 亚洲精品乱久久久久久| 乱码一卡2卡4卡精品| 少妇被粗大猛烈的视频| 亚洲国产精品成人久久小说| 在线播放无遮挡| 又大又黄又爽视频免费| 亚洲国产最新在线播放| 亚洲丝袜综合中文字幕| 久久人人爽人人爽人人片va| 亚洲丝袜综合中文字幕| 91久久精品国产一区二区三区| 午夜免费男女啪啪视频观看| 久久人人爽av亚洲精品天堂| 精品少妇黑人巨大在线播放| 久久国产亚洲av麻豆专区| 女性生殖器流出的白浆| 一级毛片我不卡| 99国产综合亚洲精品| 天美传媒精品一区二区| 国产高清三级在线| 国产黄色免费在线视频| 黄色毛片三级朝国网站| 午夜影院在线不卡| 亚洲精品日韩av片在线观看| 亚洲精品中文字幕在线视频| 久久精品国产自在天天线| 99热这里只有精品一区| 99热网站在线观看| 夜夜爽夜夜爽视频| 男女边吃奶边做爰视频| 熟女av电影| 久久人人爽av亚洲精品天堂| 大香蕉久久网| 插阴视频在线观看视频| 夫妻性生交免费视频一级片| 美女xxoo啪啪120秒动态图| 亚洲人成网站在线播| 美女大奶头黄色视频| 国产高清有码在线观看视频| 亚洲国产毛片av蜜桃av| 色吧在线观看| 久久婷婷青草| 99九九线精品视频在线观看视频| 人妻制服诱惑在线中文字幕| 午夜91福利影院| 欧美日韩一区二区视频在线观看视频在线| 亚洲国产精品成人久久小说| 免费高清在线观看视频在线观看| 欧美 日韩 精品 国产| 又大又黄又爽视频免费| 国产色爽女视频免费观看| 26uuu在线亚洲综合色| a级毛片黄视频| 性高湖久久久久久久久免费观看| 亚洲精品色激情综合| 18禁动态无遮挡网站| 永久免费av网站大全| 亚洲欧美中文字幕日韩二区| 免费人妻精品一区二区三区视频| 纯流量卡能插随身wifi吗| 搡女人真爽免费视频火全软件| 中文字幕人妻熟人妻熟丝袜美| 五月玫瑰六月丁香| 久久久欧美国产精品| 国产综合精华液| 精品酒店卫生间| 一级毛片我不卡| 国产有黄有色有爽视频| 日日爽夜夜爽网站| 夜夜爽夜夜爽视频| 一区在线观看完整版| 99久久精品国产国产毛片| 大香蕉久久网| 国产综合精华液| 黄色配什么色好看| 插阴视频在线观看视频| 日韩人妻高清精品专区| 国产一级毛片在线| 一级毛片 在线播放| 青春草视频在线免费观看| 九色成人免费人妻av| a级片在线免费高清观看视频| av国产精品久久久久影院| 免费高清在线观看日韩| 国产精品久久久久成人av| 国产精品一区二区三区四区免费观看| 国产精品一区二区在线不卡| 国产日韩欧美亚洲二区| 亚洲美女黄色视频免费看| 久久久久国产精品人妻一区二区| 母亲3免费完整高清在线观看 | 伦理电影大哥的女人| 3wmmmm亚洲av在线观看| 国产日韩欧美亚洲二区| 国产精品三级大全| 一个人看视频在线观看www免费| 黄色毛片三级朝国网站| 欧美xxⅹ黑人| 麻豆乱淫一区二区| 看非洲黑人一级黄片| av线在线观看网站| 我的女老师完整版在线观看| 亚洲av在线观看美女高潮| 日本色播在线视频| 亚洲欧美日韩另类电影网站| 高清午夜精品一区二区三区| 国产高清有码在线观看视频| 人妻一区二区av| 美女xxoo啪啪120秒动态图| 少妇高潮的动态图| 女性被躁到高潮视频| 欧美精品人与动牲交sv欧美| 国产精品久久久久久精品古装| 精品久久久久久电影网| 秋霞伦理黄片| 亚洲精品av麻豆狂野| 日本免费在线观看一区| 国产精品国产三级国产专区5o| 又大又黄又爽视频免费| 亚洲精品国产av蜜桃| 少妇丰满av| 熟妇人妻不卡中文字幕| 精品久久久久久久久亚洲| 狂野欧美白嫩少妇大欣赏| 男女啪啪激烈高潮av片| 欧美日韩av久久| 一级毛片黄色毛片免费观看视频| 午夜免费鲁丝| 久久毛片免费看一区二区三区| 国产成人精品婷婷| 国产精品麻豆人妻色哟哟久久| 少妇的逼好多水| 国产伦理片在线播放av一区| 最后的刺客免费高清国语| 天天躁夜夜躁狠狠久久av| 国产精品蜜桃在线观看| 亚洲欧洲国产日韩| 十八禁高潮呻吟视频| 日韩大片免费观看网站| 国产深夜福利视频在线观看| 国产极品天堂在线| 97精品久久久久久久久久精品| 国产女主播在线喷水免费视频网站| 婷婷色综合www| 亚洲欧美精品自产自拍| 日韩欧美精品免费久久| 日本免费在线观看一区| 曰老女人黄片| av免费观看日本| 熟妇人妻不卡中文字幕| 日韩制服骚丝袜av| 涩涩av久久男人的天堂| 亚洲欧美成人综合另类久久久| 成人二区视频| 看非洲黑人一级黄片| 内地一区二区视频在线| 国产亚洲欧美精品永久| 成年人免费黄色播放视频| 婷婷色综合www| 亚洲精品,欧美精品| av黄色大香蕉| 久久影院123| 久久精品人人爽人人爽视色| 久久精品夜色国产| 精品一区二区免费观看| 18+在线观看网站| 大片免费播放器 马上看| 亚洲精品乱码久久久久久按摩| 高清毛片免费看| 久久久精品94久久精品| 国产男女超爽视频在线观看| 一边摸一边做爽爽视频免费| 女人精品久久久久毛片| 久久av网站| 中文字幕精品免费在线观看视频 | 一级毛片 在线播放| 日日摸夜夜添夜夜添av毛片| av电影中文网址| 大香蕉久久成人网| 国产一区有黄有色的免费视频| 亚洲国产精品一区三区| 免费观看在线日韩| 国产无遮挡羞羞视频在线观看| 国产欧美另类精品又又久久亚洲欧美| 国产探花极品一区二区| 免费日韩欧美在线观看| 3wmmmm亚洲av在线观看| 99久久精品一区二区三区| 嫩草影院入口| 超色免费av| 亚洲国产毛片av蜜桃av| 亚洲精品久久午夜乱码| 男女免费视频国产| 狠狠精品人妻久久久久久综合| 国产精品欧美亚洲77777| 免费av不卡在线播放| 麻豆乱淫一区二区| 亚洲精品色激情综合| 日韩,欧美,国产一区二区三区| 丝袜脚勾引网站| 国产伦精品一区二区三区视频9| av视频免费观看在线观看| 亚洲高清免费不卡视频| 国产精品偷伦视频观看了| 国产日韩欧美在线精品| 久久国产亚洲av麻豆专区| 亚洲精品视频女| 新久久久久国产一级毛片| 欧美日韩视频高清一区二区三区二| 大话2 男鬼变身卡| 成年av动漫网址| av天堂久久9| 免费观看在线日韩| 在线亚洲精品国产二区图片欧美 | 看免费成人av毛片| 狂野欧美激情性bbbbbb| 一级毛片 在线播放| 国产探花极品一区二区| 国产 精品1| 极品人妻少妇av视频| 国产成人一区二区在线| 亚洲国产精品成人久久小说| 亚洲怡红院男人天堂| 午夜福利,免费看| 久久久久网色| 日韩制服骚丝袜av| 欧美+日韩+精品| 丝袜脚勾引网站| 亚洲综合色惰| 国产成人freesex在线| 国产精品久久久久久精品电影小说| 国产精品秋霞免费鲁丝片| 色5月婷婷丁香| a级毛片黄视频| 日本vs欧美在线观看视频| 亚洲国产精品999| 五月天丁香电影| 菩萨蛮人人尽说江南好唐韦庄| 欧美人与善性xxx| 水蜜桃什么品种好| 日韩大片免费观看网站| 免费大片黄手机在线观看| 亚洲精品成人av观看孕妇| av又黄又爽大尺度在线免费看| videosex国产| 欧美丝袜亚洲另类| 有码 亚洲区| 在线天堂最新版资源| 纯流量卡能插随身wifi吗| 在线观看国产h片| 永久免费av网站大全| 中文字幕人妻熟人妻熟丝袜美|