流固耦合破坏分析的多分辨率PD-SPH 方法1)

姚学昊 陈丁 武立伟 黄丹

(河海大学力学与材料学院,南京 211100)

引言

流-固耦合(fluid-structure interaction,FSI)是一种流体与固体相互作用和影响的力学问题,其广泛存在于诸多工程领域,如水上降落、液舱晃荡以及海浪砰击舰船与海洋平台等.由于流体流动、结构变形破坏以及二相间相互作用问题的复杂性,流-固耦合相关问题,特别是流-固耦合导致结构破坏问题的数值模拟一直是研究难点.已有诸多数值方法用于求解FSI 问题,其中,基于网格的数值方法因计算精度和效率高而被广泛应用,如有限单元法(finite element method,FEM)[1-2]、有限体积法(finite volume method,FVM)[3]、有限差分法(finite difference method,FDM)[4]及其混合方法[5-6]等.然而,基于网格的数值方法在求解FSI 问题时往往需要网格重构、界面追踪等特殊数值技术,增加了算法复杂性,并且在处理三维裂纹萌生、动态多裂纹扩展等强不连续问题时依然面临挑战.无网格方法是另一类被广泛使用于FSI 问题的数值方法,其不受网格约束,能在拉格朗日框架下求解流体域或固体域,在复杂流动问题模拟以及固体结构的大变形与断裂分析等方面展现出独特的优势,已成为研究FSI 问题的重要手段.

近场动力学(peridynamics,PD)[7-9]与光滑粒子流体动力学(smoothed particle hydrodynamics,SPH)[10]是两种基于非局部假定的无网格方法,通常将计算域离散为一系列自由移动的粒子,通过求解控制每个粒子运动的积分方程或偏微分方程获得整个系统的力学行为.PD 方法无需预先确定裂纹位置,可自然地描述裂纹萌生和扩展,近年来成功应用于固体断裂[11]、冲击破坏[12-13]等不连续力学问题模拟.针对FSI 问题,部分学者建立了不同的PD 模型,如Zhang等[14]和Zhou等[15]分别建立了用于水力压裂模拟的PD 模型,Ren等[16]建立了模拟海冰在海浪作用下力学行为的PD 模型.然而,当前的PD 模型通常不适用于含自由表面流体的FSI 问题,涉及复杂流动的PD 相关算法仍有待进一步研究[17-18].SPH 方法适用于流体运动模拟,且能处理结构动态响应问题.近年来,Oger等[19]应用弱可压SPH 方法对刚性楔形体入水问题进行研究,获得了较为理想的数值结果,Sun等[20]结合δ+-SPH 与完全拉格朗日格式的SPH 方法研究了多种流体与弹性结构相互作用问题,Khayyer等[21]建立了不可压缩SPH 与传统SPH 混合方法,将其应用于波浪冲击可变形结构问题.尽管SPH 可以单独处理FSI 问题,但目前很少有研究考虑强耦合作用导致的固体破坏.

在PD-SPH 混合方法研究方面,Fan等[22-23]针对土壤的爆炸破碎问题,采用非常规态型PD 模型建立了一种固体粒子与流体粒子互为虚粒子的混合方法,Sun等[24]、姚学昊等[25]针对FSI 问题分别采用基于接触力的混合方法.上述方法均使用相同的粒子分辨率离散固体域和流体域,在求解结构空间尺度远小于流体的问题时需要较高的计算成本.Rahimi等[26]通过一种基于拉格朗日插值的PD-SPH 混合方法,实现了空间多分辨率离散,但未考虑多时间步长问题.

基于上述研究现状,本文提出一种用于流-固耦合结构破坏模拟的多时间和空间分辨率PD-SPH 混合方法.该方法采用不同的粒子分辨率和时间步长对固体域和流体域进行计算以降低计算成本.通过典型的自由表面流体与可变形结构相互作用问题算例验证了方法的精度和效率,并应用于流体冲击混凝土板破坏过程分析.

1 数值模拟方法

1.1 PD 方法

传统连续介质力学理论中物体内部的相互作用由应力张量的散度 ∇ ·σ 表示[9],而PD 方法中固体域被离散为一系列包含质量m和密度ρ等基本信息的无拓扑关系的粒子Xi,每个粒子与其近场范围HX(半径为δ的球形区域)内的其他粒子Xj间存在相互作用f

利用体积膨胀标量可将拉伸标量状态分解为球量部分es=θx/3和偏量部分ed=e-es.此时,借鉴经典理论中弹性应变能密度表达式,可得到适用于简单弹性材料的常规态型PD 变形能密度函数

式中,k和α是与体积模量K和剪切模量G相关的PD 材料参数,可根据PD 理论中变形能密度与传统理论的应变能密度等效获得[27].计算变形能密度函数关于拉伸标量状态e的弗雷谢(Fréchet)导数则可得到力标量状态t

针对断裂问题,假定两个粒子Xi和Xj之间的相对位置满足一定关系时,粒子对内粒子间的相互作用消失.可引入间断函数表示

1.2 SPH 方法

在SPH 方法中,同样将流体域离散为一系列粒子xi,通过核函数插值和粒子近似技术可得到离散形式的不可压缩流体控制方程

式中,D /Dt表示物质导数;xji=xj-xi为粒子在当前构形中的相对位置矢量,vij=vi-vj为相对速度矢量;h为光滑长度;p为压力;fg,fΠ和fυ分别为体力、人工黏性力和物理黏性力;c0为参考声速,其值应足够大以满足弱可压缩假定[28]

式中,|v|max,pmax和ρ0分别为流体的最大流速、最大压力以及参考密度.∇i表示对粒子xi坐标的空间导数;Wij=W(xij,h)为核函数,本文选用B-样条核函数[10]

式中,qs=|xij|/h;空间维度dim=1,2,3 分别对应αw=1,1 5/7π,3 /2π.控制方程组(15)中的连续性方程右侧第二项为人工耗散项[28],用以抑制因数值噪声产生的压力振荡,其中 δs为耗散强度控制常数,常取值为0.1; 〈 ∇ρ〉L为修正密度梯度

对于不可压缩流体,粒子xi受到的物理黏性力可简化为

式中,µ为流体动力黏度;φ=0.1h用于防止奇点的出现.为了提高数值稳定性,Monaghan[29]在动量方程中加入了人工黏性项fΠ以抑制数值振荡,其表达式为

2 多分辨率PD-SPH 耦合方案

本文采用具有不同初始间距的固体粒子和流体粒子离散计算域,并通过流体粒子-虚粒子接触耦合方案处理流-固交界面.为了更显著地提高计算效率,利用不同的时间步长对固体和流体控制方程进行时间积分.

2.1 多分辨率空间离散

如图1 所示,当SPH 流体粒子a靠近流-固界面时,PD 固体粒子b和c均位于SPH 粒子a的支持域Sa内,此时粒子b和c将以虚粒子b′和c′的形式加入粒子a的邻近搜索列表,从而保证粒子a的支持域不被流-固界面截断,实现SPH 边界缺陷修正.考虑到PD 粒子与SPH 粒子不同的初始间距会导致虚粒子与流体粒子光滑长度不同,本文中令虚粒子的光滑长度等于SPH 粒子的光滑长度h,这也使得流体粒子a的支持域Sa的大小与虚粒子b′,c′的支持域Sb′,Sc′ 完全相同.此时,SPH 粒子a的控制方程变为

图1 多分辨率PD-SPH 耦合方案Fig.1 Multi-resolution PD-SPH coupling scheme

式中,Kr和nr为用户自定义参数[10].

需要注意的是,在该耦合方法中,虚粒子 κ′根据PD 粒子在初始构形中的位置分为表面虚粒子(如b′)和内部虚粒子(如c′),仅被动地被流体粒子搜索,自身不进行SPH 数值积分,其位置信息来自对应的PD 粒子,密度等场变量信息则通过插值计算获得[25].如图1 所示,对于表面虚粒子b′,其场变量信息来自其支持域Sb′ 内的流体粒子;内部虚粒子c′的场变量信息则来自其支持域Sc′ 内的流体粒子和表面虚粒子.虚粒子速度(无滑移边界条件)和密度表达式为

对于PD 部分,界面处PD 粒子受到SPH 粒子的作用力,相当于在耦合界面处对PD 粒子施加力边界条件.为了在流-固界面处满足动力学条件,保证系统动量守恒,令SPH 粒子a对PD 粒子 κ 的作用力Fatoκ

2.2 多分辨率时间离散

考虑到PD 与SPH 的完全显式计算,分别定义结构和流体的计算时间步ΔtS和ΔtF,其中,固体计算的时间步长应满足[31]

通常,固体计算的时间步长ΔtS会小于流体计算的时间步长ΔtF,若采用单一时间步长会出现流体计算耗时过多的问题.本文将采用多时间步长方案,并通过执行一次流体时间积分同时进行 ℓ=ΔtF/ΔtS次固体时间积分的方式保证流体与固体计算的时间一致性.该方案示意图如图2 所示.

图2 ℓ=6 时多时间步长方案Fig.2 Multi-timestep scheme(ℓ=6)

多分辨率PD-SPH 耦合计算流程如图3 所示.

图3 多分辨率PD-SPH 耦合计算流程Fig.3 Flow chart for multi-resolution PD-SPH coupling calculation

3 算例验证

为验证多分辨率PD-SPH 混合方法的精度和计算效率,首先对经典的FSI 问题进行数值分析.随后对流固相互作用导致的结构断裂破坏问题进行模拟,进一步验证耦合方案对于FSI 问题的适用性.本文算例均在Intel i7-7700CPU(主频 3.6GHz)个人计算机上完成,所取光滑长度和近场范围半径分别为h=1.3dxF和δ=3.0dxS,其中dxF和dxS分别表示SPH粒子与PD 粒子的初始间距.

3.1 液柱静压力下的铝板变形

首先对液柱静压力作用下的铝板变形问题进行数值模拟.如图4 所示,长L=1.0m、厚d=0.05 m,密度 ρS=2700kg/m3,弹性模量E=67.5 GPa,泊松比ν=0.34 的铝板突然受到高H=2 m,密度ρF=1000kg/m3的水柱作用.根据文献[33],铝板经过初始振荡后达到平衡状态,其跨中挠度解析解为 - 6.85×10-5m.模拟中,参考声速c0=60m/s,总模拟时长为 5 s.

图4 液柱静压力下铝板变形的初始条件(单位:m)Fig.4 Initial condition for the aluminum plate under hydrostatic pressure(unit:m)

针对当前问题,考虑到过大的流-固分辨率比会对计算精度和效率会产生一定的影响,设置了固定结构空间分辨率条件下的3 种多分辨率工况,具体参数设置如表1 所示.图5 给出了3 种工况下铝板跨中挠度时程,由图可见,经过初始振荡,3 种工况下铝板跨中挠度均收敛至约为 - 6.70×10-5m,与解析解相比误差仅为2.19%.图6 给出了t=5 s 时工况III 的流体粒子分布与压力场以及铝板挠度场,可以看出,此时流体压力场与铝板挠度场均是光滑的.上述结果表明,本文多分辨率PD-SPH 混合方法计算精度较高.

表1 液柱静压力下铝板变形问题工况布置Table 1 The setup of different cases for aluminum plate deformation with hydrostatic pressure

图5 铝板跨中挠度时间历程Fig.5 Time histories for mid-span deflection of aluminum plate

图6 t=5 s 时工况III 的流体压力场与板的挠度场Fig.6 Fluid pressure and plate deflection field in case-III at t=5 s

图7为3 种工况下FSI 系统归一化总能量[21](分别为当前和初始时刻系统总能量)时间历程.由图可见,归一化总能量同样存在振动变化过程,随着FSI 系统进入稳定状态,系统总能量最终保持稳定.t=5 s 时,3 种工况的系统总能量损失分别为0.23%,0.34%以及0.49%,表明多分辨率混合方法具有较好的能量守恒特性.综合图5与图7 可以发现,铝板位移初始振荡的持续时间与FSI 系统总能量损失有关,系统能量损伤越少,其位移振荡时间越长.该发现与Zhang等[34]所得结论一致,进一步验证了本文的多分辨率混合方法.

图7 FSI 系统归一化总能量时间历程Fig.7 Time histories for normalized total energy of FSI system

图8 给出了3 种工况下全仿真过程计算时长,可见,随着流体粒子间距增大,流体计算所需时间步长也增大,总模拟时长显著减少.工况III 的计算速度约为工况I 的11.3 倍.由此表明,多分辨率PD-SPH混合方法能够在保证计算精度的同时有效提高计算效率.

图8 不同工况的计算时间Fig.8 Computational times for different cases

3.2 溃坝流体冲击弹性闸门

Antoci等[35]对溃坝水流冲破弹性闸门问题进行了实验研究,该问题的初始条件示意图如图9 所示.实验在宽L=0.1 m、内部水深H1=0.14 m 的水箱中进行.箱体左侧刚性壁面下部放置了上端固定下端自由的弹性橡胶闸门,其宽度为s=0.005 m、高度为H=0.079 m,材料参数为:密度 ρS=1100kg/m3,弹性模量E=7.8 MPa,泊松比 ν=0.47.数值模拟中,参考声速c0=15m/s,水体密度ρF=1000kg/m3,总模拟时长0.4 s.针对该问题的多分辨率工况设置如表2 所示.

表2 溃坝流体冲破弹性闸门问题工况布置Table 2 The setup of different cases for dam-break flow through an elastic gate

图9 溃坝流体冲破弹性闸门的初始条件示意图(单位:mm)Fig.9 Schematic sketch for initial conditions in the dam-break flow through an elastic gate(unit:mm)

图10给出了工况III 中不同时刻的弹性闸门在时变水压作用下的变形以及自由液面变化过程,并与实验观察结果[35]进行比较.由图10可以看出,PD-SPH 所得弹性闸门与自由液面形状与实验结果基本吻合.在t=0.08s 时,闸门自由端在水压力作用下向x轴负方向移动,形成水流出口.在t=0.16 s 时,闸门自由端位移与水流出口显著增大.随后,由于水位的降低,水压力减小,水流出口逐渐减小.值得注意的是,由于本文开展的是二维模拟,实验中的流体飞溅现象不会出现.

图10 不同时刻弹性闸门变形和自由液面演变Fig.10 Elastic gate deformation and change of fluid free surface

图11为弹性闸门自由端水平位移时程曲线.通过与文献[34]的SPH 模拟结果以及文献[35]的实验数据对比可以看出,3 种工况下所得位移结果均吻合较好.在P D-S P H 模拟中,弹性闸门于t=0.122 s时刻达到最大变形,此时其自由端水平位移约为0.043 m.图12 给出了3 种工况的计算总时长,比较可知,仅采用多分辨空间离散方法(工况II)即可有效提高计算效率;若采用多时间与空间分辨率,即工况III,其计算速度为单一分辨率(工况I)的8.5 倍.

图11 弹性闸门自由端水平位移时程Fig.11 Time histories of horizontal displacements at the free end of elastic gate

图12 3 种工况的计算总时长Fig.12 Total calculation time in three cases

综上,本文提出的PD-SPH 混合方法在不同分辨率工况下均能有效传递流-固相互作用信息,可用于涉及复杂流体流动和结构大变形的FSI 问题.

3.3 重力坝水力劈裂

为了验证多分辨率PD-SPH 混合方法对于FSI 致结构破坏问题的适用性,选取经典的带裂缝Koyna 重力坝问题,分析其在超载作用下的断裂响应.初始条件如图13 所示,坝高 103 m、坝顶和坝底宽度分别为w=14.8 m和L=70m.水平裂缝位于大坝上游面 66.5 m 高程处,长度为1.93 m.大坝材料参数为:密度 ρS=2450kg/m3,弹性模量E=25 GPa,泊松比 ν=0.2.考虑大坝承受水头超载作用,根据文献[36],取上游面超高水头为 10m,即大坝上游面水深为113 m.模拟中,PD 粒子间距dxS=0.2 m,时间步长ΔtS=5×10-6s ;水体密度 ρF=1000kg/m3,SPH 粒子间距dxF=0.4 m,时间步长ΔtF=5×10-5s,参考声速c0=1500m/s.为与文献工况保持一致,防止水流溢出,水体高于坝体部分右侧添加了固壁边界;计算中未考虑上游水渗入裂缝产生的缝面水压.

图13 Koyna 重力坝模型(单位:m)Fig.13 Koyna gravity dam model(unit:m)

图14 给出了Koyna 坝在水头超载作用下的裂缝扩展路径.可以看出,水力裂缝首先在初始裂缝顶点萌生,并沿水平方向扩展约3.2 m.随后,裂缝开始向右下方扩展,此时裂缝与水平方向夹角为33°.本文所得裂缝扩展路径与文献[36]结果对比如图15所示,可以看出,两种方法所得裂缝路径基本一致,表明本文提出的多分辨率PD-SPH 混合方法能够较好地再现流体作用下的结构断裂过程.

图14 Koyna 大坝裂缝扩展过程Fig.14 Crack propagation process in Koyna dam

图15 重力坝裂缝扩展路径Fig.15 Crack paths in the gravity dam

3.4 水流冲击作用下的混凝土板崩塌问题

通过上述验证,本节将提出的方法应用于工程常见的水流冲击作用下混凝土板崩塌这一复杂流-固耦合结构破坏问题模拟.初始条件如图16 所示,水体经由宽为w=0.33 m 的水流入口并以5 m/s的速度流入;混凝土板长L=2 m,高H=0.3 m,其左右两端顶点受到固定约束,中部含一条长0.1 m、宽0.015 m 的竖向裂缝,板的具体材料参数为:密度ρS=2400kg/m3,弹性模量E=35 GPa,泊松比 ν=0.2.模拟中,令PD 粒子与SPH 粒子的初始间距分别为dxS=0.005 m和dxF=0.01 m,两相计算时间步长为ΔtS=1×10-6s,ΔtF=5×10-6s ;水体密度ρF=1000kg/m3,参考声速c0=60m/s.

图16 水流冲击混凝土板模型(单位:m)Fig.16 Model of water impacting on a concrete slab(unit:m)

图17为水流冲击作用下混凝土板崩塌过程,图18则给出了混凝土板上表面中部A点(如图16 所示)的竖向位移时程曲线.在t=0.12 s 时,部分水体自水流出口流出,并在重力作用下加速运动;t=0.212 s 时,水流开始冲击混凝土板上部;t=0.215 s时,初始裂缝开始扩展,此时A点的竖向位移较小,仅为 - 9.18×10-5m.随着时间推移,裂缝逐渐向上扩展,并在t=0.222 s 时,扩展至混凝土板上表面形成贯穿裂缝,A点位移也于此刻产生巨大变化.t=0.31 s时,左右两块混凝土板在自重及水流冲击作用下分别绕固定点旋转,使得裂缝左右表面呈倒V 形,此时A点的竖向位移为 -0.11 m ;水流则因混凝土板的阻挡作用形成横向射流.最终,水流完全冲开了两块独立运动的混凝土板,并自由向下运动;A点在t=0.45 s 时的竖向位移为 -0.57 m.本文所得混凝土板断裂过程与文献[37]一致,表明本文提出的多分辨率PD-SPH 混合方法适用于复杂FSI 问题求解,能够自然实现流体运动与结构变形破坏全过程连续仿真.

图17 水流冲击作用下混凝土板崩塌过程Fig.17 Collapse process of concrete slab subjected to impact of water flow

图18 A 点竖向位移时程曲线Fig.18 Time history of vertical displacement at point A

4 小结

本文针对流-固耦合破坏问题提出一种多时间和空间分辨率的PD-SPH 混合方法,并将其应用于典型问题模拟,得到以下结论.

(1)该混合方法利用具有不同初始间距的粒子离散固体域和流体域,采用不同的时间步长对两相的控制方程进行时间积分.通过将PD 粒子设置为具有与流体粒子相同光滑长度的虚粒子处理流-固界面,可使得界面处满足动力学和运动学条件,且能够保证系统动量守恒.

(2)通过液柱静压力作用下铝板变形和溃坝流体冲破弹性闸门两个经典考题模拟,结果表明:多分辨率PD-SPH 混合方法具有较高的计算精度和较好的能量守恒特性,针对不同问题合理选择流-固分辨率比以及时间步长比可有效提高计算效率.

(3)方法应用于经典Koyna 重力坝水力劈裂和水流冲击作用下混凝土板的崩塌问题,得到了固体结构的运动、变形与断裂过程以及固体结构影响下流体的运动状态.所得结果验证了多分辨率PD-SPH混合方法在分析流-固耦合破坏问题方面的适用性.在此基础上,后续可扩展至三维流-固耦合破坏问题研究.本文工作希望为分析工程中同时涉及移动自由表面、复杂流动、流-固耦合导致结构大位移和损伤破坏问题的高效数值仿真提供参考.