流固耦合问题的PD-SPH 建模与分析

姚学昊，黄 丹

(河海大学工程力学系，江苏，南京 211100)

流固耦合(Fluid-structure interaction, FSI)问题广泛存在于科学和工程领域，常涉及结构变形破坏以及流体破碎等复杂流动现象[1]，往往难以获得解析解。随着计算力学研究发展，任意拉格朗日-欧拉(Arbitrary lagrangian-eulerian, ALE)方法[2-3]、浸入边界法(Immersed boundary method, IBM)[4]等各种数值方法被用于此类问题分析，然而在处理结构大变形、准确快速地追踪自由表面和运动边界等方面依然存在困难。

与网格类方法不同，光滑粒子动力学(Smoothed particle hydrodynamics, SPH)方法[1,5]是一种拉格朗日无网格法，它适合于追踪自由表面和运动边界，且在处理大变形问题时不会产生网格畸变。ANTOCI等[6]采用SPH 方法模拟流体与弹性结构的相互作用问题，获得了较理想的结果。KHAYYER 等[7]和ZHANG 等[8]分别提出了一种不可压缩SPH 方法以及多分辨率SPH 方法来提高流固耦合问题的求解精度与效率。但SPH 方法在模拟固体时仍需要一些修正技术以克服拉伸不稳定问题[1,5]等固有缺陷。近年来，研究人员开发了诸多SPH 与其他方法耦合的算法，如SPH 与有限元法(Finite element method, FEM)[9-10]、光滑点插值法(Smoothed point interpolation method, SPIM)[11]的耦合等。然而，上述算法主要针对含结构变形的流固耦合问题，在涉及固体损伤与破坏的相关问题中具有局限性。近场动力学(Peridynamics, PD)[12-14]是一种新兴的非局部方法，采用空间积分代替传统微分方程，在分析不连续力学问题时展现出独特优势，可应用于多种大变形和结构破坏[15-17]、流固耦合结构破坏[18-19]等模拟。采用PD-SPH 耦合方法来分析流固耦合问题，能同时发挥PD 方法模拟固体破坏和SPH 方法追踪流体运动界面的优势。

近来，部分学者开展了PD 与SPH 方法结合的初步研究。ZHOU 等[20]结合两种方法提出SPD(Smoothed peridynamics)方法并应用于大变形与断裂问题。FAN 等[21-22]提出一种PD-SPH 耦合模型模拟爆炸引起的土壤破碎问题；SUN 等[23]最近提出一种基于移动虚粒子的PD-SPH 耦合方法，并应用于流固耦合问题。

本文建立一种能求解流固耦合及结构破坏的新型PD-SPH 耦合模型。针对流-固界面处理，提出一种基于虚粒子和排斥力的耦合算法，既能防止流体粒子穿透交界面，又能对交界面处流体粒子进行边界缺陷修正，提高计算精度。通过模拟流体作用下的弹性结构变形问题，验证了耦合方法的可行性和有效性，并进一步开展了涉及结构断裂破坏的复杂流固耦合问题全过程仿真分析。

1 数值模型

1.1 固体模型

1.1.1 常规态型PD 基本方程

如图1 所示，将空间域R离散为一系列包含物性信息且具有相互作用f的PD 粒子，在任意时刻t，粒子Xi的运动方程可写为：

图1 PD 模型示意图Fig. 1 Sketch for PD model

PD 理论分为键型(bond-based)PD 理论、常规态型(ordinary state-based)PD 理论以及非常规态型(non-ordinary state-based)PD 理论。由于键型与非常规态型PD 理论通常分别存在泊松比限制及数值不稳定问题，本文基于常规态型PD 理论进行固体域求解。此时，点对相互作用f可表示为：

1.2 流体模型

1.2.1 控制方程

1.3 流-固交界面处理方案

PD 和SPH 方法均通过粒子离散计算域，故可通过基于粒子-粒子接触的耦合方案处理流-固交界面。图3 为本文提出的基于虚粒子和排斥力的PDSPH 耦合方案示意图。在该方案中，PD 粒子将作为两种类型的耦合虚粒子(边界虚粒子和内部虚粒子)参与SPH 计算。同时，反作用力将作用于PD粒子从而实现SPH 粒子对PD 粒子的影响。

该耦合方案具体计算流程如图4 所示。

图4 基于虚粒子和排斥力的PD-SPH 耦合方案流程Fig. 4 Flow chart for PD-SPH coupling scheme based on virtual particle and repulsive force

2 算例分析

2.1 静水压力作用下的铝板变形

为了检验耦合方案的有效性、稳定性及计算精度，首先研究了静水压力作用下的铝板变形问题[9]。如图5 所示，在宽1.0 m 的水箱中，箱内水深H= 2.0 m，密度为 1000 kg/m3；水箱底部为厚度d= 0.05 m 的铝板，其密度是 2700 kg/m3，弹性模量67.5 GPa，泊松比0.34。

图5 静水压力下的铝板示意图Fig. 5 Sketch for aluminium plate under hydrostatic pressure

初始时刻，铝板突然受到水柱的静水压力载荷而发生变形。经过短时间的初始振荡，系统最终达到具有特定静态变形的平衡状态。根据理论解[9]，2.0 m 高水柱产生的静水压力载荷作用下铝板跨中挠度大小为 -6.85×10-5m。在模拟计算中：人工粘性参数απ= 0.03 ； 参考声速c0=190 m/s；光滑长度h=1.33Δx；近场范围半径 δ=3Δx；时间步长 Δt=1×10-6s；总模拟时间为1 s。流体粒子与固体粒子的初始间距一致，粒子空间分辨率分别取L/Δx=100、L/Δx=80和L/Δx=60，对应的粒子间距为0.01 m、0.0125 m 以及0.0167 m，此时SPH 粒子总数分别为21060、13648、7836，PD 粒子总数分别为525、340、195。本算例在配有主频为2.7 GHz 的Intel CPU 的计算机上完成，计算总耗时分别为8 h 31 min、6 h 4 min、4 h 32 min。

图6 给出了模拟得到的t=1 s 时不同空间分辨率(粒子间距)下流体压力云图以及放大1000 倍后的铝板变形图。可见，基于PD-SPH 的流固耦合求解方法在不同空间分辨率条件下均可得到光滑的流体压力场和铝板挠度场。

图6 t = 1 s 时不同空间分辨率下流体压力云图和铝板变形Fig. 6 The fluid pressure contours and the deflection of the aluminium plate at t = 1 s with different spatial resolutions

图7 为三种不同空间分辨率下铝板跨中挠度的时间历程，表1 则给出了不同分辨率下跨中挠度数值解及其与解析解的相对误差。结合图7 与表1 可以看出，空间分辨率L/Δx=60时，由于粒子数较少，PD 和SPH 方法本身的计算精度较低；随着空间分辨率的增加，PD-SPH 耦合方法的跨中挠度计算结果逐渐向静态理论解析解收敛，并最终在L/Δx=100时，与解析解基本一致，相对误差仅为3.07%。由此表明：本文提出的PD-SPH 耦合算法能有效、准确模拟准静态流固耦合问题。

表1 不同空间分辨率下铝板跨中挠度计算误差Table 1 Computational error of mid-span deflection of aluminium plate with different spatial resolutions

图7 铝板跨中挠度Fig. 7 Mid-span deflection of aluminium plate

2.2 溃坝水流冲击弹性板

如图8 所示，宽为0.584 m 的水箱左侧存在高H= 0.292 m，宽L= 0.146 m 的水柱，水柱在重力作用下倒塌产生溃坝流，并冲击水箱中高he=0.08 m，厚度为w= 0.012 m 的弹性板，随后发生剧烈的流固耦合作用。该算例中，水的密度为1000 kg/m3；弹性板的密度为 2500 kg/m3，弹性模量为1.0 MPa，泊松比为0.0。模拟中，人工粘性参数设为 απ= 0.02 ； 参考声速c0=35 m/s；粒子间距Δx= 0.002 m，对应的SPH 粒子与PD 粒子总数分别为12744 和246；光滑长度h=1.5Δx；近场范围半径 δ=3Δx； 时间步长为 Δt=5×10-6s；总模拟时间为0.75 s。使用与算例2.1 相同的计算机配置，计算总耗时为2 h 20 min。

图8 溃坝水流冲击弹性板的初始条件示意图Fig. 8 Sketch for initial conditions in dam-break flows impacting the elastic plate

为了对PD-SPH 结果进行定量验证，在PDSPH 模拟中计算了弹性板顶端A点(如图8 所示)的水平位移时间历程，与其它已有文献结果[29-32]的对比如图9 所示。可以看出：在水流冲击弹性板前期(0.4 s 前)，PD-SPH 方法的位移计算结果与其它方法所得结果吻合良好，且弹性板的位移最大值以及最大位移的出现时间基本一致，其中，PD-SPH 方法在0.244 s 时获得弹性板水平位移最大值，约0.046 m；在0.4 s 以后，受自由表面融合现象的影响，不同方法得到的位移结果存在明显差异，PD-SPH 模拟结果与粒子有限元法(PFEM)[29]、SPH-再生核粒子法(RKPM)[32]模拟结果更为接近。

图9 A 点水平方向位移时程Fig. 9 Time history of horizontal displacement of point A

图10 为不同时刻下，溃坝水流冲击弹性板产生的流体飞溅与结构变形。可以看出：0.16 s 时，水流冲击弹性板并沿板的左侧向上爬升，同时，弹性板在冲击力作用下出现弯曲变形；0.26 s 时，水流完全覆盖弹性板左侧，使其产生较大变形；0.42 s 时，水流撞击右侧壁面形成射流，弹性板则出现较大幅度的回弹现象。同时，图10 还给出了文献中PFEM 方法[29]和SPH 方法[30]的模拟结果，比较看出：每个时刻3 种方法得到的流体自由表面形状和弹性板变形都是相近的。此外，PD-SPH方法也计算得到了光滑的流体压力场，并且与PFEM结果相似。由此表明：本文提出的PD-SPH 耦合方法对剧烈流固耦合问题的模拟是可行的。

图10 不同时刻弹性板变形和流体飞溅的PFEM[29]、SPH[30]和PD-SPH 结果的比较Fig. 10 Comparison of PFEM[29], SPH[30] and PD-SPH results of plate deformation and fluid splashing at different times

2.3 流体作用下的结构破坏

为了验证本文PD-SPH 方法在流固耦合破坏问题模拟方面的可行性，本节模拟结构在溃坝流体作用下的运动、破坏过程。如图11 所示，水箱长0.8 m，高0.4 m，其内部水柱和挡板的几何尺寸、材料参数均与算例2.2 相同。为了实现挡板在水流作用下的断裂破坏，假定挡板材料强度极限较低。本算例中，SPH 和PD 粒子总数分别为14267 和258，总模拟时间为0.8 s，计算总耗时为2 h 27 min。

图11 流体作用下结构破坏的初始条件示意图Fig. 11 Sketch for initial conditions for structural failure under fluid action

图12 为不同时刻下结构变形破坏和流体飞溅的模拟结果，图13 则给出了结构质心A(如图11 所示)的坐标随时间的变化。在0.18 s 前，流体运动过程以及结构变形与算例2.2 完全一致；0.18 s 时，结构底部开始断裂(如图12(a)所示)，质心A的y坐标时程曲线出现瞬时震荡(如图13(b)所示)；在水流的持续作用下，裂纹不断扩展，产生的应力波也在水体内部传播，导致流体压力场振荡，0.19 s 时裂纹呈现明显的V 字型，上部挡板仅剩1 个粒子与底部连接(如图12(b)所示)；最终在0.2 s 时，上部挡板与底部完全断开(如图12(c)所示)；随后，挡板在水流作用下向水箱右侧运动，于0.318 s 时撞击水箱底部(如图12(d)和图13(b)所示)，并出现反弹现象(接触力由式(18)计算)；经过几次接触回弹，挡板开始沿水箱底部滑动，0.524 s时挡板撞击水箱右侧(如图12(f)和图13(a)所示)，并开始向左滑动。最终，挡板将在水箱底部保持静止状态。本文模拟所得挡板运动过程与文献[33]基本一致，由此表明：本文提出的PD-SPH 方法可应用于流固耦合破坏问题求解，且能够较好地预测流体运动过程和结构变形破坏以及破坏后部分结构的运动。

图12 不同时刻结构变形破坏和流体飞溅的模拟结果Fig. 12 Simulation results of structural deformation and failure and fluid splashing at different times

图13 结构质心A 的坐标时间历程Fig. 13 Time history of the coordinate of centroid A

3 结论

本文建立了一种新的PD-SPH 耦合模型求解流固耦合作用下的结构破坏问题，主要工作和结论如下：

(1)分别采用PD 方法与SPH 方法离散固体域和流体域，可充分发挥两种方法的优势。通过基于虚粒子和排斥力的耦合算法处理流-固交界面，能够对流固界面处的流体粒子施加速度边界，修正核函数缺陷，提高计算精度。

(2)采用PD-SPH 耦合方法模拟了静水压力作用下的铝板变形问题以及溃坝水流冲击弹性板问题，所得结构变形和流体运动过程与解析解或文献结果吻合较好，验证了本文PD-SPH 耦合方法在模拟复杂流固耦合问题方面的适用性。

(3)对流体冲击作用下的结构破坏过程分析表明：所提出的PD-SPH 耦合方法易于实现流固耦合作用下结构从变形到破坏乃至破坏后部分结构运动的全过程仿真模拟，可为流-固耦合结构破坏问题的研究提供新参考。