流固耦合(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，粒子X_i的运动方程可写为：

图  1  PD模型示意图

Figure  1.  Sketch for PD model

下载: 全尺寸图片 幻灯片

式中：$\rho$为材料密度；u为位移矢量；$\tilde {\boldsymbol{b}}$为作用在粒子X_i上的外载荷密度；L_u为t时刻粒子X_i所受合内力密度，与点对相互作用f相关，其表达式为：

式中：H_X为半径为δ的近场范围；${V_j}$为粒子X_j的体积。

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

式中：${\boldsymbol{\xi }} = {{\boldsymbol{X}}_j} - {{\boldsymbol{X}}_i}$为初始构形中PD粒子间的相对位置矢量；$\underset{\raisebox{0.5em}{--}}{{\boldsymbol{T}}}$为力矢量状态：

式中：$\underset{\raisebox{0.5em}{--}}{{\boldsymbol{M}}} = ({\boldsymbol{\xi}} + {\boldsymbol{\eta}} )/| {{\boldsymbol{\xi}} + {\boldsymbol{\eta}} } |$为单位方向矢量状态，${\boldsymbol{\eta}} = {{\boldsymbol{u}}_j} - {{\boldsymbol{u}}_i}$为相对位移矢量；$\underset{\raisebox{0.5em}{--}}{t}$为力标量状态。

“状态”是态型PD理论引入的用于描述粒子物理行为的概念，分别定义参考位置矢量状态$\underset{\raisebox{0.5em}{--}}{{\boldsymbol{X}}} \langle {\boldsymbol{\xi}} \rangle = {\boldsymbol{\xi}}$、形变矢量状态$\underset{\raisebox{0.5em}{--}}{{\boldsymbol{Y}}} \langle {\boldsymbol{\xi}} \rangle = {\boldsymbol{\xi}} + {\boldsymbol{\eta}}$以及拉伸标量状态$\underset{\raisebox{0.5em}{--}}{e} \langle {\boldsymbol{\xi}} \rangle = | {\underset{\raisebox{0.5em}{--}}{{\boldsymbol{Y}}} \langle {\boldsymbol{\xi}} \rangle } | - | {\underset{\raisebox{0.5em}{--}}{{\boldsymbol{X}}} \langle {\boldsymbol{\xi}} \rangle } |$，并根据力状态与能量密度的关系^[12]，可得到线弹性常规态型PD模型的力标量状态表达式：

式中：k'和α为PD材料参数^[24]；$\underset{\raisebox{0.5em}{--}}{x} = | {\underset{\raisebox{0.5em}{--}}{{\boldsymbol{X}}} \langle {\boldsymbol{\xi}} \rangle } |$；$\underset{\raisebox{0.5em}{--}}{\omega }$为影响函数^[15]，其值只与粒子对长度$\left| {\boldsymbol{\xi}} \right|$有关；${q_{\rm{v}}}$和$\theta$分别为体积权重和体积膨胀标量，定义为：

${\underset{\raisebox{0.5em}{--}}{e} ^d}$为拉伸标量状态$\underset{\raisebox{0.5em}{--}}{e} \langle {\boldsymbol{\xi}} \rangle$的偏量部分，即${\underset{\raisebox{0.5em}{--}}{e} ^d} = \underset{\raisebox{0.5em}{--}}{e} - \theta \underset{\raisebox{0.5em}{--}}{x} /3$。

1.1.2   损伤与断裂描述

基于粒子X_i与X_j之间的伸长率${s_{ij}}$判断点对相互作用f的有无，即判断粒子间“键”的断裂与否，可引入间断函数$\mu$描述：

式中：${s_{ij}} = \left( {\left| {{\boldsymbol{\xi}} + {\boldsymbol{\eta}} } \right| - \left| {\boldsymbol{\xi}} \right|} \right)/\left| {\boldsymbol{\xi}} \right|$；${s_0}$为临界伸长率，可根据断裂参数定义为${s_0} = \sqrt {5{G_0}/9K\delta }$，其中，${G_0}$为断裂能密度，K为体积模量。

根据每个粒子近场范围内断键情况，可定义局部损伤D($0 \leqslant D \leqslant 1$)为^{[15, 25]}：

1.2   流体模型

1.2.1   控制方程

如图2所示，SPH方法同样采用一系列粒子x_i离散流体域Ω。此时，所有粒子遵照控制方程的规律运动即可描述流场的变化。在SPH方法中，基于核函数W，通过核近似和粒子近似可将场函数或其导数的表达式转化为对支持域S(半径为κh，$\kappa$为常数，h为光滑长度)内所有粒子加权求和的形式。由此可得不可压缩流体控制方程的离散形式：

图  2  SPH模型示意图

Figure  2.  Sketch for SPH model

下载: 全尺寸图片 幻灯片

式中：v为流体速度矢量；${{\boldsymbol{v}}_{ij}} = {{\boldsymbol{v}}_i} - {{\boldsymbol{v}}_j}$为粒子间相对速度；m为粒子质量；p为流体压力；${c_0}$为参考声速^[5]；g为重力加速度；${\nabla _i}$为梯度算子；核函数${W_{ij}} = W({|{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}|,h)$，本文选用B-样条型核函数^[1]，此时$\kappa = 2$；${\delta _{\rm{s}}}$为控制耗散强度的可调参数，常取值为0.2；${{\boldsymbol{\psi}} _{ij}}$的表达式为：

式中，$\left\langle {\nabla \rho } \right\rangle _i^L$为修正密度梯度^[5]。${{\boldsymbol{\varPi}} _{ij}}$为人工粘度^[11]，用于抑制数值振荡：

式中：${\alpha _\pi }$为粘度控制参数；${{\boldsymbol{x}}_{ij}} = {{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}$为现时构形中粒子间相对位置；${\bar c_{ij}} = ({c_i} + {c_j})/2$和${\bar \rho _{ij}} = ({\rho _i} + {\rho _j})/2$分别表示平均声速和平均密度；$\varphi = 0.1h$用于避免因粒子过近产生的数值发散。物理粘性项${\boldsymbol{f}}_i^\upsilon$可表示为：

式中，${\upsilon _{\rm{o}}}$为流体运动粘度系数。

1.2.2   粒子位移修正

粒子分布不规则会引起数值不稳定，通过引入粒子位移技术(Particle shifting technique, PST)重新分布粒子，能够有效提高SPH方法的精度和鲁棒性^[26]。在PST方法中，粒子x_i的位置需要额外修正至${{\boldsymbol{x}}_i} + \delta {{\boldsymbol{r}}_i}$，其中位置修正矢量$\delta {{\boldsymbol{r}}_i}$记为：

式中：${v_{\max }}$为最大流速；$\Delta t$为计算时间步长；$\Delta x$为初始粒子间距。在求得$\delta {\boldsymbol{r}}$后，可根据泰勒公式对流体粒子的场变量进行修正^[26]：

式中，${\phi _i}$和${\phi _{i'}}$分别为更新前后的粒子场变量。

1.3   流-固交界面处理方案

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

图  3  PD-SPH耦合方案示意图

Figure  3.  Sketch for coupled scheme of PD-SPH

下载: 全尺寸图片 幻灯片

在耦合处理过程中，首先根据固体的初始构形确定每个PD粒子对应的耦合虚粒子类型。如图3(a)所示，若PD粒子位于结构表面，则将其伪装为边界虚粒子(粗虚线圆)，其它PD粒子则伪装为内部虚粒子(细虚线圆)。然后，在每个时间步的初始时刻，将当前PD粒子${\boldsymbol{X}}_n^{{\text{PD}}}$的位置信息传递给耦合虚粒子${\boldsymbol{x}}_{n'}^{{\text{SPH}}}$，即${\boldsymbol{x}}_{n'}^{{\text{SPH}}} = {\boldsymbol{X}}_n^{{\text{PD}}} + {\boldsymbol{u}}_n^{{\text{PD}}}$(${{\boldsymbol{x}}_{n'}}$包含两种类型的耦合虚粒子)。同时，耦合虚粒子参与SPH计算所需的场变量信息，如密度、速度，需要根据其当前支持域(半径为两倍光滑长度2h)内邻近粒子的对应场变量信息插值得到^[1]。如图3(a)所示，边界虚粒子(如粒子$b'$)的场变量信息来自其支持域内的流体粒子，而内部虚粒子(如粒子$c'$)的场变量信息则来自其支持域内流体粒子和边界虚粒子。耦合虚粒子的密度计算式为：

对于耦合虚粒子速度，则需根据边界条件获得^[8]。

如图3(b)所示，当SPH流体粒子a与耦合虚粒子(如粒子$b'$和$d'$)的间距小于两倍光滑长度(2h)，粒子产生接触，此时认为粒子a、$b'$和$d'$均位于流固耦合区域。耦合虚粒子$n'$(包含粒子$b'$和$d'$)将参与粒子$a$的SPH控制方程计算，从而实现SPH边界缺陷的修正以及速度边界的施加。此外，为防止流体粒子穿透交界面，耦合虚粒子$n'$同时会对流体粒子$a$施加排斥力${{\boldsymbol{f}}_{\rm{c}}}$，其表达式为：

式中，排斥力${{\boldsymbol{f}}_{{c}}}$根据接触势函数$\phi ({{\boldsymbol{x}}_a})$^[27-28]的梯度求得，势函数$\phi ({{\boldsymbol{x}}_a})$的表达式为：

式中：$\Delta {h_{{\rm{avg}}}}$为平均光滑长度；${K_c}$和${n_c}$为用户自定义参数。因此，PD粒子对SPH粒子$a$的合作用力为${{\boldsymbol{F}}_{{\rm{Pt}}{S_a}}} = \displaystyle\sum\limits_{n' \in \,\,{S_a}} {\left[ {{{\boldsymbol{F}}_{n'a}} - {{\boldsymbol{f}}_c}\left( {{{\boldsymbol{x}}_a}} \right)} \right]}$，其中${{\boldsymbol{F}}_{n'a}}$由式(11)获得。根据牛顿第三定律，将PD粒子$n$受到的流体粒子的反作用力(即耦合虚粒子$n'$受到的作用力)${{\boldsymbol{F}}_{{{{\rm{St{P}}}}_n}}} = - \displaystyle\sum\limits_a {\left[ {{{\boldsymbol{F}}_{n'a}} - {{\boldsymbol{f}}_c}\left( {{{\boldsymbol{x}}_a}} \right)} \right]}$代入式(1)可得考虑流体作用的PD运动方程：

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

图  4  基于虚粒子和排斥力的PD-SPH耦合方案流程

Figure  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\;{\text{kg/}}{{\text{m}}^3}$；水箱底部为厚度d = 0.05 m的铝板，其密度是$2700\;{\text{kg/}}{{\text{m}}^3}$，弹性模量67.5 GPa，泊松比0.34。

图  5  静水压力下的铝板示意图

Figure  5.  Sketch for aluminium plate under hydrostatic pressure

下载: 全尺寸图片 幻灯片

初始时刻，铝板突然受到水柱的静水压力载荷而发生变形。经过短时间的初始振荡，系统最终达到具有特定静态变形的平衡状态。根据理论解^[9]，2.0 m高水柱产生的静水压力载荷作用下铝板跨中挠度大小为$- 6.85 \times {10^{ - 5}}\;{\text{m}}$。在模拟计算中：人工粘性参数${\alpha _\pi }{\text{ = }}\;{\text{0}}{\text{.03}}$；参考声速${c_0} = 190\;{\text{m/s}}$；光滑长度$h = 1.{\text{33}}\,\Delta x$；近场范围半径$\delta = 3\Delta x$；时间步长$\Delta t = 1 \times {10^{ - 6}}\;{\text{s}}$；总模拟时间为1 s。流体粒子与固体粒子的初始间距一致，粒子空间分辨率分别取$L/\Delta x = 100$、$L/\Delta x = 80$和$L/\Delta 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时不同空间分辨率下流体压力云图和铝板变形

Figure  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/\Delta x = 60$时，由于粒子数较少，PD和SPH方法本身的计算精度较低；随着空间分辨率的增加，PD-SPH耦合方法的跨中挠度计算结果逐渐向静态理论解析解收敛，并最终在$L/\Delta x = 100$时，与解析解基本一致，相对误差仅为3.07%。由此表明：本文提出的PD-SPH耦合算法能有效、准确模拟准静态流固耦合问题。

图  7  铝板跨中挠度

Figure  7.  Mid-span deflection of aluminium plate

下载: 全尺寸图片 幻灯片

表  1  不同空间分辨率下铝板跨中挠度计算误差

Table  1.  Computational error of mid-span deflection of aluminium plate with different spatial resolutions

分辨率	跨中挠度/m	相对误差/(%)
L/Δx=60	–8.22×10–5	20.0
L/Δx=80	–7.40×10–5	8.03
L/Δx=100	–7.06×10–5	3.07

下载: 导出CSV 

| 显示表格

2.2   溃坝水流冲击弹性板

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

图  8  溃坝水流冲击弹性板的初始条件示意图

Figure  8.  Sketch for initial conditions in dam-break flows impacting the elastic plate

下载: 全尺寸图片 幻灯片

为了对PD-SPH结果进行定量验证，在PD-SPH模拟中计算了弹性板顶端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点水平方向位移时程

Figure  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结果的比较

Figure  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  流体作用下结构破坏的初始条件示意图

Figure  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  不同时刻结构变形破坏和流体飞溅的模拟结果

Figure  12.  Simulation results of structural deformation and failure and fluid splashing at different times

下载: 全尺寸图片 幻灯片

图  13  结构质心A的坐标时间历程

Figure  13.  Time history of the coordinate of centroid A

下载: 全尺寸图片 幻灯片

3   结论

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

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

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

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