基于双向流固耦合模型的溃坝水流数值模拟

2023-10-09王春正孙亮黄玲玲陈丽芬

科学技术与工程 2023年26期

王春正, 孙亮*, 黄玲玲, 陈丽芬

(1. 武汉理工大学船海与能源动力工程学院, 武汉 430063; 2. 大连理工大学港口海岸及近海工程国家重点实验室, 大连 116024; 3. 大连市海洋可再生能源重点实验室, 大连 116024)

流固耦合是一类多学科交叉融合的问题,广泛存在于工程、科学和自然界中。例如,液化天然气船在运输过程中遭遇的液舱晃荡问题;海洋立管在复杂海洋环境下发生的涡激振动问题;风电机组叶片在运转过程中面临的气动弹性问题;超大型浮式结构物在复杂波浪条件下产生的水弹性响应问题以及盾构隧道施工[1]等场景。这些问题的共性是流体作用于结构而引起弹性结构发生位移或者变形,同时结构的动力响应又对流体的运动产生影响。

流固耦合问题较为复杂且涉及结构变形破坏、流体剧烈变化等因素,往往难以通过严格的公式推导获得解析解,因此,模型试验和数值模拟成为研究流固耦合问题的两种常用方法。然而模型试验受试验成本、仪器精度、环境因素干扰的制约,常难以获得较为精确的结果。与模型试验相比,数值模拟具有成本低、周期短、物理过程清晰等优势,从而成为解决流固耦合问题的主流。

对于流固耦合问题的数值模拟,根据是否在同一组控制方程中求解流体域和固体域,可以将研究方法分为强耦合和弱耦合两种。其中强耦合方法是在同一组控制方程中求解,又被称为直接耦合法。Walhorn等[2]提出了时空有限元法来分析具有自由液面流动的流固耦合问题,并对溃坝水流冲击弹性结构问题进行了数值模拟;Idelsohn等[3]提出了一个以统一形式处理流体子域和固体子域的拉格朗日公式,并基于PFEM(particle finite element method,粒子有限元法)对若干流固耦合问题进行了数值模拟;Rafiee等[4]提出了一种基于光滑流体动力学(smoothed particle hydrodynamics,SPH)的数值方法,并对弹性板在水压下变形和溃坝水流冲击弹性结构等问题进行了模拟;Sun等[5]基于材料点法(material point method,MPM)同样提出了一种模拟流固耦合的方法,并对液仓晃荡、弹性水坝溢流等问题进行了数值模拟。

弱耦合方法是在不同组控制方程中分别求解流体和固体域,又被称为间接耦合法。Cerquaglia等[6]基于PFEM和有限元法(finite element method,FEM)相耦合,以分区的方法对鸟类撞击、溃坝等流固耦合问题进行了数值模拟;张友林等[7]基于移动粒子半隐式方法(moving particle semi-implicit method,MPS)和FEM,开发出了MPS-FEM耦合求解器,并使用该求解器数值研究了溃坝流动对弹性结构的冲击问题;张智琅等[8]基于解耦有限粒子法(decoupled finite particle method,DFPM)和光滑有限元(smoothed finite element method,SFEM)相耦合,以分区的方法对结构入水,溃坝水流冲击弹性结构等流固耦合问题进行了数值模拟;赵冉等[9]基于格子玻尔兹曼法(lattice Boltzmann method,LBM)和浸入边界法(immersed boundary method,IBM)相耦合的方法,对重力场内单根柔性纤维与均匀流场的耦合行为进行了数值模拟;姚学昊等[10]提出了一种基于SPH和虚粒子和排斥力的近场动力学(peridyna-mics,PD)相耦合的方法,利用流体粒子-虚粒子接触算法处理流-固界面,并使用该方法对弹性结构在水压下变形和溃坝冲击弹性结构等问题进行分区求解。

强耦合方法因流固控制方程组同时收敛较为困难,故难以针对复杂的物理问题建立统一的控制方程组且消耗计算资源巨大,而弱耦合方法使用不同的控制方程分离求解流体域和固体域,因此,可以组合使用各种成熟的流固求解器,灵活分析复杂的流固耦合问题,对计算机性能的需求也大幅降低。故而,本文研究采用弱耦合的方法对流固耦合问题进行分析。

长期以来,有限体积法(finite volume method,FVM)在解决带有自由表面的两相流问题中占据了主导地位,FEM可以精确、稳定、高效地处理结构的变形问题。现使用基于FVM的开源计算流体动力学(computational fluid dynamics,CFD)软件OpenFOAM和基于FEM的开源计算固体力学(computational solid mechanics,CSM)软件CalculiX,通过开源耦合库preCICE对双向流固耦合问题展开分析。首先对代数流体体积法(volume of fluid,VOF)和几何VOF两种自由液面捕捉的方法进行对比,随后使用几何VOF方法对溃坝水流冲击弹性结构的双向流固耦合问题进行模拟,通过与已有的数据进行对比,验证本文所使用分析方法的可靠性,为今后研究相关问题提供参考。

1 数学模型

1.1 FVM方法

在本文所采用的方法中,将气体和液体视为单一可变密度的连续流体,并考虑黏性,其控制方程和自由液面捕捉的VOF方法如下所述。

1.1.1 控制方程

开源CFD软件OpenFOAM采用有限体积法求解连续方程[式(1)]和N-S方程[式(2)][11]。

(1)

(2)

式中:U为速度矢量;ρ为流体密度;μ为动力黏性系数;g为重力加速度;p为流体压力;fσ为表面张力;t为时间。

1.1.2 VOF方法

在使用OpenFOAM中的两相流求解器分析溃坝水流冲击结构物时,需要追踪流体区域内液体和气体的交界面,即自由液面。Hirt等[12]提出的VOF方法根据每个时刻网格中流体体积的变化来追踪自由液面,是目前使用最广泛的自由液面追踪方法。在分析中引入流体体积分数α表示流体计算域内气体和液体的分布:当α=1时,表示流体区域内全部为液体;当α=0时,表示流体区域内全部为气体;当0<α<1时,表示流体区域内既有气体也有液体。根据求解方法的不同,VOF方法可以分为代数VOF方法和几何VOF方法。

(1)代数VOF方法。代数VOF方法为获取更小的界面厚度,引入了人工压缩项来建立体积分数运输方程[13],即

(3)

(2)几何VOF方法。几何VOF方法(isoAdvector方法)是Roenby等[14]提出的一种新型几何VOF方法。在每个分析步内,该方法首先会在气体和液体交界的网格单元中,找到Isoface(最适合划分气体和液体的平面)[14],并基于这些Isoface对自由液面进行几何重构。之后,会使用这些重构的自由液面对流体体积分数场进行更新,下一分析步的流体体积分数场由式(4)[14]给出。

(4)

(5)

式中:αi为第i网格处的流体体积分数;Vi为第i个网格处的体积;Bi为组成网格单元i的所有面集合;Sij为确定法向量方向的辅助因子;Sj为网格表面的面积法向量;Aj(τ)为τ时刻网格表面Fj处的瞬时淹没面积;φj(t)为通过Fj处的容积通量;x为笛卡尔坐标;u(x)为网格表面Fj中x处的流速矢量。

1.2 FEM方法

使用基于FEM的开源CSM软件CalculiX[15]对结构场分析,其控制结构单元变形的动力学方程由式(6)给出。

(6)

C=α1M+α2K

(7)

式中:M为结构的质量矩阵;C为瑞丽阻尼矩阵;K为结构刚度矩阵;F(t)为施加在结构上随时间t变化的外力;y为结构单元的位移矢量;系数α1和α2的取值与结构的阻尼比和固有频率有关。

(8)

(9)

(10)

式中:参数α、β、γ控制α-Method求解的精度和稳定性。在CalculiX中参数β、γ由参数α表示。

(11)

(12)

1.3 FVM-FEM耦合方法

图1 双向流固耦合模型

2 算例分析

2.1 溃坝水流冲击刚性结构

为了比较代数VOF方法和几何VOF方法捕捉自由液面在实际应用中的效果,首先分别使用基于代数VOF法的interFoam求解器和基于几何VOF法的interIsoFoam求解器对该溃坝水流冲击刚性结构的过程[18]进行研究。如图2所示,初始时刻左侧区域有宽0.146 m、高0.292 m的水柱,在下壁面的中间位置有一个宽0.024 m、高0.048 m的刚性结构,A点位于刚性结构的左侧中点(距离底边界0.024 m)。在该算例中,空气的密度为1 kg/m3,液体的密度为1 000 kg/m3。计算过程中经过收敛性分析最终使用的网格大小为0.002 m,时间的步长设置为自适应时间步长,该时间步长由最大库朗数控制,本文采用的最大库朗数为0.5。

图2 溃坝水流冲击刚性结构

图3给出了分别使用基于代数VOF法和几何VOF法的求解器进行数值模拟,得到的不同时刻(t=0.2、0.3、0.4、0.5 s)液体的界面形态与试验[18]的对比。从图3中可以看出,随着坍塌的水流抨击刚性结构,使用代数VOF方法捕捉的自由液面可以观察到非常明显的液体体积损失,并且由于碰撞本应该出现的飞溅液珠和复杂尖锐的自由液面没有被精确捕捉,反而因为数值模拟的误差部分水体消失,而使用几何VOF法可以将自由液面很好地捕捉并且没有观察到有明显的水体损失。因此,几何VOF方法对于剧烈变化的自由液面具有更好的捕捉能力。

图3 水流冲击刚性结构物的自由表面分布

图4为基于代数VOF方法和几何VOF方法得到的计算域内水体质量的变化情况。初始时,积分得到水体的质量为0.085 264 kg,与计算得到的理论值相同,即初始误差为0。在3.0 s时,使用几何VOF法得到的水体质量稳定在0.084 71 kg,误差为0.649 7%,而使用代数VOF法得到的水体质量稳定在0.083 48 kg,误差为2.092 3%。因此,基于几何VOF方法得到的结果具有更好的质量守恒性。

图4 计算域中液体质量的变化

图5给出了使用两种VOF方法模拟左边壁上的液面变化。从图5可以看出,与试验值[18]和Issakhov等[19]的数值模拟值对比,使用两种VOF方法对左边壁上液面的模拟结果没有较大差别,原因在于左边壁处液体的运动比较平缓,并没有出现碰撞、飞溅等强非线性现象,这也说明了对于比较缓慢变化的冲击过程几何VOF方法和代数VOF方法均可以给出高精度的预测值。

图5 左侧壁处的液面变化

图6给出了使用两种VOF方法模拟图2中A点处的压力随时间变化的对比。从图6可以看出当水柱突然坍塌冲击刚性结构的A点时,该处的压力突然增大并出现最大峰值,之后刚性障碍物左侧的水体会继续冲击A点,并在0.6 s内逐渐形成较为明显的三个峰值。使用几何VOF方法所得的结果与Issakhov等[19]的结果基本吻合,而使用代数VOF方法所得的结果与Issakhov等[19]的结果差别较大。这在一定程度上与图3中观察到的自由液面破碎现象以及图4中使用基于代数VOF法模拟会出现较大的液体质量损失有关。

图6 A点处压力变化

以上的分析表明:基于代数VOF方法的流体求解器在模拟溃坝水流冲击等强非线性问题时的精度要低于基于几何VOF方法所得到的相应结果(自由表面分布、压力和质量守恒)。在分析双向流固耦合问题时,对流体区域压力的较大计算误差会影响结构的响应,而结构响应的分析误差又会影响流体的运动,这些误差会在耦合过程中不断累积,从而导致最终的模拟结果产生较大的差异。为保证双向流固耦合的分析精度,进一步的研究中将使用具有更强液面捕捉能力、质量守恒性更好的几何VOF方法来模拟双向流固耦合问题。

2.2 溃坝水流冲击弹性结构

为验证本文基于FVM-FEM方法在模拟带有自由液面的双向流固耦合问题时的准确性,同时进一步揭示此类流固耦合问题的内在机理,对溃坝水流冲击弹性结构问题[2]进行研究。如图7所示,水箱宽0.584 m,水箱左下脚存有宽0.146 m,高0.292 m的水柱,在下壁面的中间位置有一个宽0.012 m、高0.08 m的弹性结构。在该算例中,空气的密度为1 kg/m3,液体的密度为1 000 kg/m3,弹性结构的密度2 500 kg/m3,弹性模量为1 MPa,泊松比为0。

图7 溃坝水流冲击弹性结构

图8给出了本文方法(FVM-FEM)与Walhorn等[2]、Idelsohn等[3]、Rafiee等[4]所计算得出的弹性结构左上角水平位移的对比。从图8可以看出,使用本文的方法FVM-FEM与其他3种方法所得的结果曲线,在总体趋势上具有较好的一致性,且在关键点处的数值基本吻合。在0.14 s左右,4种方法的结果中均出现向左的水平位移(出现负值),本文结果为0.068 cm;在0.24 s左右时,4种方法所模拟的弹性障碍物的水平位移均达到向右的最大值,本文结果为4.88 cm,与Idelsohn等[3]的结果几乎完全一致;在0.65 s左右,弹性障碍物的水平位移达到第二个峰值,本文结果为1.68 cm,与Walhorn等[2]和Idelsohn等[3]的结果吻合较好,优于Rafiee等[4]的结果,验证了本文方法的可靠性;在第二个峰值之后,4种数值模拟方法得到的水平位移曲线具有较为明显的差别,这是因为溃坝水流冲击弹性结构是一个流体和固体相互耦合的过程,同时液体区域具有更强的非线性,求解难度更大,对求解的精确度和稳定性要求高,但可以看到的是本文的方法与其他3种方法所得的结果在总体趋势上保持一致,再次验证了本文所采用的方法在模拟具有自由液面的双向流固耦合问题时具有较好的可靠性。

图8 弹性结构左上角的位移变化

图9给出了不同方法对于弹性结构变形和流体分布在不同时刻的比较。可以看出:在t=0.14 s时,弹性结构下端受到水流的冲击作用,结构下部向右弯曲,顶部出现向左的偏转,因此结构左上角的位移为负值(图8);在t=0.26 s时,水流完全覆盖弹性结构左侧,出现明显的“射流”现象,此时弹性结构具有最大向右的弯曲幅度,同时在计算域内可以观察到飞溅的液珠和尖锐的自由液面形状,这与参考文献[3- 4]中的结果非常接近,本文方法在捕捉剧烈变化的自由液面时具有较高的精度;t=0.34 s时,水体越过弹性结构的抨击右侧壁面,弹性结构开始出现向左的回弹;t=0.42 s时,水流越过弹性结构继续砰击右侧壁面,形成较为明显的“分流”现象,部分液体沿右侧壁面开始回流,弹性结构继续向左回弹。

图9 不同时刻的弹性板变形和水流形态

图10给出了不同时刻(对应图8中位移的峰值和谷值,即t=0.14、0.23、0.53、0.65、0.84、0.94 s)液体压力和弹性板应力的分布。从图10可以看出,t=0.14 s时,溃坝水流沿着弹性结构左侧向上流动,水流受到弹性结构的阻碍作用,汇聚在弹性结构的左下角,水体在该处产生压力的最大值(约为4.2 kPa),弹性结构在中部偏下向右侧弯曲,而弹性结构的顶端受到的压力很小,因而发生了向左侧的较小偏移(图8);t=0.23 s时,溃坝水流越过弹性结构顶端冲出去,形成“射流”现象,此时弹性结构的顶端出现向右水平位移的最大值,同时弹性结构左下角的压力因为液体的流动而变小,弹性结构左右两侧的应力均出现了整个过程中的最大值,约为221 kPa和-221 kPa,此时也是弹性结构最可能发生破坏的时刻;t=0.53 s时,由于水流与右侧墙壁碰撞出现“分流”现象,一部分水体沿右侧的水槽底边壁出现回流,在弹性结构右下角出现液体压强的较大区域,此时弹性结构的顶端出现向左回弹;t=0.65 s时,右侧的水体沿弹性结构右侧边壁爬升,与结构左侧流动的水体产生“对流现象”,在弹性结构顶端形成了一个漩涡,弹性结构物开始晃动,解释了图8中本文方法得到的位移曲线在t=0.65 s附近出现轻微震荡的原因,而其他3种方法没有捕捉到该现象,同时可以观察到在弹性结构的顶端形成局部压力较大的区域;t=0.84 s和t=0.94 s时,受弹性结构运动的影响,左侧的水体被不断挤压,右侧坠落的水体和回流的水体开始发生冲撞,并且这种不规则非线性作用不断加强,导致图8中4种方法所得的位移曲线在t=0.65 s以后存在较大的差异。

图10 不同时刻流体的压力和弹性结构的应力分布

图11给出了具有相同布置的弹性结构(图7)和刚性结构在结构左侧中点处压力值的对比。从图11可以看到,在水流第一次冲击结构时,由于弹性结构产生弯曲,相比于刚性结构减小了水流的冲击作用,所以弹性结构左侧中点处所受到的压力更小。图11中弹性结构左侧中点处压力值震荡较大,最终趋向于刚性结构左侧中点处的压力值。