一种冰区航行船舶流固耦合的CFD

一种冰区航行船舶流固耦合的cfd
–
fem
–
dem模拟方法
技术领域
1.本发明属于冰区船舶安全技术领域，具体涉及一种冰区航行船舶流固耦合的cfd
–
fem
–
dem模拟方法。


背景技术：

2.极地是世界各国未来发展必争的战略空间。破冰船是用于破碎水面冰层、开辟航道、保障舰船进出冰封港口锚地或引导舰船在冰区航行的勤务船。利用极地科学考察船对两极地区进行探索开发已成为世界各国极地战略中的重要发展方向。与常规的水面舰船不同，海冰对船体结构的冲击作用是影响冰区船舶航行安全性的重要因素。冰区船舶的设计建造与航行安全性问题是船舶工程领域当前的热点研究内容，相关技术亟需探索与发展。
3.常规水面舰船的水动力学问题当中，一般只涉及到船体结构与外部流体的相互作用，相关的理论研究方法起步较早并已取得了显著成绩。例如，通过采用势流理论或计算流体动力学(cfd)技术能够准确预报航行于静水或波浪中船舶的水动力性能。然而，冰区航行船舶在外部流体载荷和固态冰冲击载荷的共同作用下发生运动和结构变形，涉及到船体结构—冰结构—海洋流体三者之间的相互作用问题，且海冰和船舶结构相互作用现象非常复杂，常规的流固耦合方法难以模拟冰区航行船舶的水动力学与结构载荷等流固耦合问题。


技术实现要素：

4.本发明的主要目的在于克服现有技术的缺点与不足，提出一种冰区航行船舶流固耦合的cfd
–
fem
–
dem模拟方法。
5.为了达到上述目的，本发明采用以下技术方案：
6.一种冰区航行船舶流固耦合的cfd
–
fem
–
dem模拟方法，包括以下步骤：
7.s1、在fem求解器中建立船体结构模型；
8.s2、在cfd求解器中建立数值水池模型；
9.s3、在dem求解器中建立浮冰结构模型；
10.s4、设置cfd
–
fem
–
dem六向流固耦合的数据传递方法，使cfd、fem、dem三个求解器的任意两者之间都进行双向耦合并相互传递数据；
11.s5、程序运行求解，运行三种求解器进行数值模拟计算，输出结果数据。
12.进一步的，步骤s1包括：
13.建立船体结构的三维有限元模型；
14.船体结构有限元网格划分；
15.船体结构材料属性定义，定义各类结构单元的质量、密度、刚度、弹性模量、泊松比以及结构阻尼；
16.施加约束和载荷，结合船体运动自由度对网格节点施加位移约束，根据网格所受载荷类型施加相应载荷。
17.进一步的，步骤s2包括：
18.建立水池流体域，水池流体域为长方体空间域，由4个侧壁面、底面以及顶面，共6个面围成；
19.建立船体外壳的几何模型，设置船体参数及其在流体域中的初始位置；
20.建立浮冰的几何模型，设置浮冰参数及其在流体域中的初始位置；
21.采用流体体积法将水池流体域划分为水和空气两部分；
22.对流体域进行网格划分，并在自由面进行网格加密；
23.设置网格变形方式，分别采用重叠网格和变形网格技术模拟船体运动和结构变形，实现网格的混合变形；
24.设置边界条件，水池流体域的4个侧壁面和底面采用速度入口，顶面采用压力出口，船体物面采用无滑移壁面。
25.进一步的，步骤s3包括：
26.建立浮冰的几何模型，确定海冰类型和几何尺寸；
27.确定浮冰的材料模型，定义各项参数，包括密度、压缩强度、拉伸强度、弯曲强度、剪切强度、弹性模量以及泊松比；
28.将浮冰结构离散生成球形颗粒，定义球形颗粒单元的质量和大小并将球形颗粒单元进行粘结排列；
29.确定浮冰单元的接触模型和断裂模型。
30.进一步的，cfd求解器用于模拟数值水池波浪场以及海洋流体流动，fem求解器用于模拟船体结构在流体力与浮冰撞击力作用下的运动和结构变形，dem求解器用于模拟浮冰运动与受力情况。
31.进一步的，步骤s4具体包括：
32.设置cfd求解器与fem求解器之间的耦合方法；
33.设置cfd求解器与dem求解器之间的耦合方法；
34.设置fem求解器与dem求解器之间的耦合方法。
35.进一步的，设置cfd求解器与fem求解器之间的耦合方法具体为：
36.在cfd
–
fem双向耦合中，cfd求解器将作用于船体表面上的流体压力和剪切力传递给fem求解器；fem求解器根据外部流场力、刚体惯性力、结构弹性力对船体结构进行动力分析，并将所得到的船体运动和变形反馈给cfd求解器，用于更新流场信息；
37.fem模型中的各个网格上的流体力通过高斯面积分得到，每个高斯点处的流体压力定义为与其距离最近的cfd体网格节点处的压力值；cfd体网格节点的位移和变形则通过其周围fem网格节点的位移的形函数插值获取。
38.进一步的，设置cfd求解器与dem求解器之间的耦合方法具体为：
39.在cfd
–
dem双向耦合中，cfd求解器将作用于整个浮冰体积上的压力梯度力和作用于浮冰表面上的粘性应力传递给dem求解器；dem求解器根据流场的滞止阻力、附加质量力、船体撞击力、浮冰间的非线性接触力进行动力学分析，求解得到浮冰的运动和破碎情况，并将两相动量交换反馈到cfd求解器，用于更新流场信息；
40.dem模型内离散体系的浮冰粘结颗粒中心位置处的流体压力梯度、流速和空隙率根据其宿主cfd网格单元中心的相应物理量插值得到，cfd网格节点上的动量交换源项根据所有映射到该网格单元内浮冰粘结颗粒中心处的流固两相滞止阻力积分得到。
41.进一步的，设置fem求解器与dem求解器之间的耦合方法具体为：
42.在fem
–
dem双向耦合中，dem求解器将作用于船体表面上的浮冰撞击力和摩擦力传递给fem求解器；fem求解器根据外部浮冰撞击力和摩擦力、刚体惯性力、结构弹性力对船体结构进行动力分析，并将所得到的船体运动和变形反馈给dem求解器，用于下一时刻浮冰的运动和受力分析；在每个时间步内，对dem颗粒与fem网格之间的有效耦合界面进行接触判断，并将接触力等效为结构节点载荷，传递到fem网格单元中进行动力计算；与此同时，将fem求解器计算的结构变形量转化为耦合界面的位移，作为dem离散单元中的位移边界条件。
43.进一步的，步骤s5具体为：
44.设置联合模拟方法为隐式的分区交错耦合算法；
45.程序运行求解，设置计算时间步长、总模拟时长、每个时间步长内的数据交换次数，并同步运行三种求解器进行数值模拟计算；
46.输出结果数据，在cfd求解器中输出波面高程和船体表面压力，在dem求解器中输出浮冰运动和受力，在fem求解器中输出船体运动、剖面载荷、局部应力以及结构变形。
47.本发明与现有技术相比，具有如下优点和有益效果：
48.1、本发明提出了考虑船舶与流体、船舶与浮冰、浮冰与流体相互作用的一体化的时域模拟方法，将船舶结构、浮冰结构、海洋流体的相互作用综合考虑，充分考虑各类流固耦合问题，能够更精确地模拟极地船舶在冰区中航行时的运动和受力情况。
49.2、本发明综合三类数值方法的优势，采用cfd模拟外部海洋波浪场、采用fem模拟船舶结构、采用dem模拟浮冰结构，其中cfd
–
fem、cfd
–
dem和fem
–
dem均为双向耦合，共计六向耦合，完全考虑了各类流固耦合效应。
50.3、本方法除了能够计算海洋流体与浮冰共同作用下的船舶结构的运动和受力情况外，还能够考虑船体结构的变形效应及振动响应；因此，本方法能够分析冰块撞击频率与船体结构固有频率的相互关系对于船体运动、结构共振响应的影响，为冰区航行船舶的设计研发提供可靠方法。
附图说明
51.图1是本发明方法的流程图；
52.图2是cfd、fem、dem求解器在流固耦合问题中的作用区域示意图；
53.图3是数值水池计算域侧视图；
54.图4是数值水池计算域俯视图；
55.图5是cfd求解器中的船体外壳模型；
56.图6是实施例中浮冰颗粒单元的接触模型；
57.图7是六向耦合数据传递方案示意图；
58.图8是cfd和fem求解器之间的数据传递方法示意图；
59.图9是cfd和dem求解器之间的数据传递方法示意图；
60.图10是dem颗粒和三角形fem网格单元等效节点力计算模型示意图；
61.图11是dem颗粒和四边形fem网格单元等效节点力计算模型示意图；
62.附图标号说明：1-船体；2-水；3-空气；4-自由液面；5-浮冰；6-第一sph颗粒；7-第
二sph颗粒；8-法向阻尼；9-切向阻尼；10-法向刚度；11-切向刚度；12-摩擦系数；13-cfd网格单元；14-cfd网格节点；15-fem网格单元；16-高斯点；17-浮冰粘结颗粒；18-fem网格节点；19-等效节点力。
具体实施方式
63.下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。
64.术语说明：
65.cfd：计算流体动力学；
66.fem：有限元方法；
67.dem：离散单元法。
68.实施例
69.本实施例中，采用某集装箱船模为例，对其在冰缘区航行时的流固耦合问题进行模拟。计算模型缩尺比为1:40，模型水线长4.375m、排水量370.5kg。如图2所示，cfd求解器用于模拟数值水池波浪场以及海洋流体流动，fem求解器用于模拟船体1在流体力与浮冰撞击力作用下的运动和结构变形，dem求解器用于模拟浮冰运动与受力情况。
70.如图1所示，本发明，一种冰区航行船舶流固耦合的cfd
–
fem
–
dem模拟方法，包括以下步骤：
71.s1、在fem求解器中建立船体结构模型；包括以下步骤：
72.s11、建立船体1结构的三维有限元模型，包括船体外板、甲板、舱壁、骨材等结构。
73.s12、船体结构有限元网格划分，针对船体外板、甲板和舱壁进行壳单元网格划分，针对加强骨材进行梁单元网格划分。
74.s13、船体结构材料属性定义，定义各类结构单元的质量、密度、刚度、弹性模量、泊松比、结构阻尼。
75.本实施例中，根据实际的船体质量分布情况，设置各结构单元的质量；根据实际的船体刚度分布情况，设置结构单元的刚度。钢材密度为7850kg/m3、弹性模量为206gpa、泊松比0.3。采用瑞利阻尼模型，阻尼系数为0.05。
76.s14、施加约束和载荷，结合船体运动自由度对网格节点施加位移约束，根据网格所受载荷类型施加相应载荷。本实施例中，为了防止船体的纵向漂移，需要在船体重心位置处的节点上施加纵向位移约束。需要对所有网格节点施加重力载荷，对于与海洋流体或浮冰接触的船体表面的网格节点施加相应的cfd模拟的流体或dem模拟的固体载荷。
77.s2、在cfd求解器中建立数值水池模型；包括以下步骤：
78.s21、建立水池流体域，水池流体域为长方体空间域，由4个侧壁面、底面和顶面共6个面围成。数值水池流体域和船体模型相对位置的侧视图和俯视图分别如图3和图4所示。
79.s22、建立船体外壳的几何模型，船体外壳由船体外板和甲板组成。船体外壳的几何模型如图5所示。设置船体的重量、重心位置和转动惯量，设置船体在流体域中的位置。船模的重量为370.5kg，重心纵向位置距离艉柱2.125m、重心垂向位置距离基线0.213m，纵摇惯性半径为1.052m、横摇惯性半径为0.241m。船体中纵剖面、艉柱横剖面和静水面的交点处设置在坐标原点。
80.s23、建立浮冰的几何模型，设置浮冰参数及其在流体域中的初始位置。本实施例中，选取冰缘区碎冰航道，浮冰采用长方体单元，其边长分布符合正态分布规律，浮冰的厚度取20mm。浮冰随机分布在计算流体域中的静水面上。
81.s24、采用流体体积法将水池流体域划分为水2和空气3两部分，自由面的垂向位置位于坐标原点。当网格单元的体积分数γ＝0表示完全为空气，γ＝1表示完全为水，γ＝0.5表示由50％水和50％空气组成。
82.s25、对流体域进行网格划分，并在船体和自由面附近进行网格加密。采用六面体非结构化网格对计算域进行网格划分，采用重叠网格方案，计算域包括背景区域和重叠区域两部分。为精确捕捉自由液面和船体周围的湍流等物理量的剧烈变化，在自由液面4以及船体周围进行局部网格加密。波高范围内包含16层网格，波长范围内包含80层网格，船体外壳壁面处设置5层边界层网格。
83.s26、设置网格变形方式，分别采用重叠网格(overset grid)和变形网格(morphing grid)技术模拟船体运动和结构变形，实现网格的混合变形。重叠网格技术用于模拟船体在波浪中的大幅运动，变形网格技术用于模拟船体结构的弹性变形。
84.s27、设置边界条件，水池流体域的4个侧壁面和底面采用速度入口，顶面采用压力出口，船体物面采用无滑移壁面。
85.s3、在dem求解器中建立浮冰结构模型；包括以下步骤：
86.s31、建立浮冰5的几何模型，确定海冰类型、几何尺寸。本实施例中，选取冰缘区碎冰航道，浮冰采用长方体单元，其边长取50～200mm范围内并符合正态分布规律，浮冰的厚度取20mm。
87.s32、确定浮冰的材料模型，定义密度、压缩强度、拉伸强度、弯曲强度、剪切强度、弹性模量、泊松比等。本实施例中，浮冰密度取900kg/m3，弹性模量9gpa，泊松比0.3，阶段应力35mpa，最大主应力35mpa。
88.s33、将浮冰结构离散生成球形颗粒，定义球形颗粒单元的质量和大小，并将球形颗粒单元进行粘结排列。
89.s34、确定浮冰单元的接触模型和断裂模型。粘结颗粒单元之间的力与位移的关系通过接触模型来确定。粘弹性作用力通过计算颗粒之间的法向相对速度和法向位移来得到，剪切力则需要符合mohr-coulomb摩擦定律，颗粒单元的接触模型如图6所示，图6中，标号6为第一sph颗粒6，标号7为第二sph颗粒7，标号8为法向阻尼8，标号9为切向阻尼9，标号10为法向刚度10，标号11为切向刚度11，标号12为摩擦系数12。本实施例采用hertz-mindlin非线性粘弹性模型模拟颗粒碰撞过程。
90.s4、设置cfd
–
fem
–
dem六向流固耦合的数据传递方法，使cfd、fem、dem三个求解器的任意两者之间都进行双向耦合并相互传递数据，如图7所示；包括以下步骤：
91.s41、设置cfd求解器与fem求解器之间的耦合方法，具体为：
92.在cfd
–
fem双向耦合中，cfd求解器将作用于船体表面上的流体压力和剪切力传递给fem求解器；fem求解器根据外部流场力、刚体惯性力、结构弹性力对船体结构进行动力分析，并将所得到的船体运动和变形反馈给cfd求解器，用于更新流场信息。
93.由于cfd模型与fem模型中的船体表面网格是不匹配的(通常，fem模型中的船体表面网格要比cfd模型的粗糙)，数据传递不能直接进行。
94.fem模型中的各个网格上的流体力可通过高斯面积分得到，如图8所示，在四边形网格上引入9个高斯点16，每个高斯点处的流体压力定义为与其距离最近的cfd体网格节点处的压力值。
95.cfd体网格节点的位移和变形则通过其周围fem网格节点18的位移的形函数插值获取。
96.s42、设置cfd求解器与dem求解器之间的耦合方法，具体为：
97.在cfd
–
dem双向耦合中，存在两种运动体系，分别是基于欧拉观点的连续流体运动(用于cfd求解器)和基于拉格朗日观点的离散粘结颗粒运动(用于dem求解器)。通常，欧拉连续流体运动体系下的cfd网格单元尺度要比拉格朗日离散体系下dem模型中的浮冰粘结颗粒尺度大得多。因此这两种运动独自产生的物理量数据在体系间动量耦合过程中，无法直接传递，需要对压力梯度、流体速度和空隙率(或固含率)等物理量进行映射。
98.如图9所示，在基于拉格朗日离散体系的dem求解器中，先通过追踪浮冰粘结颗粒与cfd网格单元13的几何相对位置关系，精确计算cfd网格单元的空隙率(或固含率)；再通过最小二乘法或green-gauss方法分别重构cfd网格节点14上的空隙率梯度、流体压力梯度和流体速度梯度；进而利用广义taylor展开定理将欧拉网格节点处的流体压力梯度、空隙率(或固含率)和流体速度插值到浮冰粘结颗粒中心处，从而计算浮冰粘结颗粒所受流体的压力梯度力；再根据插值的空隙率和流体速度计算广义阻力系数、附加质量力及相间粘性滞止阻力。求得的流体压力梯度力、附加质量力、相间粘性滞止阻力将传递给dem求解器，与dem模型本身的非线性接触力合成作为载荷，施加在浮冰粘结颗粒之上，用于浮冰更新运动信息。
99.另一方面，在欧拉体系的cfd网格节点上，将网格单元内已算得的作用在每一个浮冰粘结颗粒17上的相间滞止阻力累加作为动量交换源项，反馈到cfd求解器，形成相间的动量耦合。
100.s43、设置fem求解器与dem求解器之间的耦合方法，具体为：
101.在fem
–
dem双向耦合中，dem求解器将作用于船体表面上的浮冰撞击力和摩擦力传递给fem求解器；fem求解器根据外部浮冰撞击力和摩擦力、刚体惯性力、结构弹性力对船体结构进行动力分析，并将所得到的船体运动和变形反馈给dem求解器，用于计算下一时刻的浮冰运动和受力情况。在每个时间步内，对dem颗粒与fem网格之间的有效耦合界面进行接触判断，并将接触力等效为结构节点载荷，传递到fem网格单元15中进行动力计算；与此同时，将fem求解器计算的结构变形量转化为耦合界面的位移量，作为dem离散单元中的位移边界条件。fem
–
dem耦合界面上的参数传递是耦合模型计算的关键问题。可将船体外壳的fem模型划分为三角形或四边形网格单元，进而计算dem颗粒与fem网格之间的接触力。采用inside-outside搜索算法判断dem颗粒与fem网格单元的接触类型，进而计算dem颗粒与fem网格单元的接触力。
102.对于dem计算部分，fem网格可看做三角形或四边形的几何边界，在其范围内进行搜索算法和接触力计算。对于fem计算部分，基于虚功原理，将dem颗粒单元与fem网格单元的接触力作为集中载荷等效施加到fem网格单元的节点上。
103.如图10和图11所示，分别为dem颗粒与三角形和四边形fem网格单元等效节点力19计算模型。由于fem网格单元的位移在节点计算，因此需要对接触点进行位移插值更新耦合
界面的位移信息。以dem颗粒单元与四边形fem网格单元的搜索为例，inside-outside算法在离散单元-结构搜索中的实现过程具体为：
104.首先，判断dem颗粒与fem网格单元的接触类型。dem颗粒与fem网格单元的接触类型主要有点面接触、点边接触和点点接触。通过inside-outside算法可判断dem颗粒单元在四边形fem网格单元边的内侧或外侧，通过对四条边分别进行判断即可得到dem颗粒与fem网格单元的接触类型。
105.s5、程序运行求解及数据导出；包括以下步骤：
106.s51、设置联合模拟方法为隐式的分区交错耦合算法；
107.s52、程序运行求解，设置计算时间步长、总模拟时长、每个时间步长内的数据交换次数，并同步运行三种求解器进行数值模拟计算。
108.s53、输出结果数据，在cfd求解器中输出波面高程、船体表面压力，在dem求解器中输出浮冰运动和受力，在fem求解器中输出船体运动、剖面载荷、局部应力、结构变形。
109.还需要说明的是，在本说明书中，诸如术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
110.对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其他实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。