基于多介质、多尺度离散元方法的冰载荷数值冰水池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 大学、美国船级社、新加坡国立大学、韩国船舶与海洋工程研究所等国内外单位诸多专家学者的支持和帮助,在此一并致谢;感谢大连理工大学刘璐、杨冬宝、吴捷等人对极地船舶与海洋工程结构冰载荷的离散元数值计算.