模拟流体冲击致结构破坏问题的SPH-PD耦合方法

时浩天，郭 力

(东南大学 土木工程学院，南京 210000)

流固耦合问题具有多场、多物理的特性，广泛存在于工程应用与生产生活中。在滑坡、泥石流冲击下游建筑，海啸冲击堤坝、海岸结构等自然灾害场景，海上平台受强海浪冲击作用等工程问题中，往往涉及到自由表面流体与结构的剧烈相互作用。由于流体包含自由表面破碎、浪花飞溅和融合等高度非线性流动过程，同时结构在强冲击作用下会发生大变形乃至破坏失效等复杂力学行为，使得这类问题的求解极具挑战性。

在流固耦合问题中，基本方程包括流体和结构各自区域内的控制方程及其边界条件，以及流固界面上力的平衡条件和几何相容条件。将这些方程组纳入到统一的框架中进行求解非常困难，已有的一些统一求解策略仅适用于简单的流固耦合问题。相应地，分区求解方法越来越受到人们的关注。分区求解法将流体和结构在各自区域分别求解，通过在流固界面上传递数据(压力、速度、位移等)来实现耦合计算。分区求解法的优势在于可以最大程度地利用已有的流体和结构求解方法，综合现有的流体和结构求解模块，通过求解器之间的数据传递实现耦合求解，从而显著降低了开发成本。

构建流固界面上数据的有效传递策略是分区求解法中的关键问题。流固界面往往是时变的，为了精确捕捉流固界面的位置，已经发展了任意拉格朗日-欧拉方法[1](arbitrary Lagrangian-Euler，ALE)和欧拉-拉格朗日耦合方法[2](coupled Euler-Lagrangian,CEL)。这些方法主要采用动网格技术，捕捉流固界面的位形，同时保证每个时间步内流体和结构在各自网格区域内正确运算。动网格技术能够有效处理较小变形下的流固耦合问题，但是在处理任意复杂的流固界面时，追踪流固界面需要网格进行大量的处理与重构，过程非常繁琐。浸入边界法[3](immersed boundary,IB)采用固定网格，避免了动网格技术的缺点，但是在流固界面上的分辨率和精度较低。IB在高雷诺数条件下的应用较为有限，多适用于流体粘性较大的生物医学问题。

对于分区求解方法，通过组合现有计算流体力学和固体力学中的求解框架，发展了众多的流固耦合求解方法。Wick[4]将ALE方法和相场法相耦合，用于求解流固耦合问题，同时可以考虑结构的脆性断裂行为；Zheng 等[5]将有限元法与边界元法耦合，研究了浸没在无限可压缩流体中弹性结构的自由振动问题。上述基于网格的方法可以有效地处理流固耦合问题，但是当流体域包含自由表面、结构域涉及大变形时，对网格的处理会导致高昂的计算成本。此外，当结构域涉及到裂纹的萌生及扩展和材料失效时，会产生几何上的不连续和奇异性问题，基于网格的方法将存在固有的困难。针对这类问题的求解，无网格法具有独特的优势。因此，发展一套具有鲁棒性和高精度的无网格类方法，来求解包含自由流体表面、结构大变形和材料断裂失效的复杂流固耦合问题，无疑具有重要的意义。

光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)是一种具有完全拉格朗日属性的无网格粒子类方法，广泛应用于自由表面流、多相流、非牛顿流体、湍流、考虑表面张力的流体流动的模拟以及化学扩散和沉淀分析等[6]。此外，SPH也可以用来模拟固体，比如弹性体、柔性体[7]。沈雁鸣等[8]采用SPH统一算法模拟了自由流体冲击弹性体问题，证明了SPH具有独立处理流固耦合问题的能力；施书文等[9]将SPH和光滑点插值法结合，模拟了弹性体入水问题。在泥石流、滑坡等灾难仿真以及近海结构受强海浪冲击作用等场景下，考虑结构的破坏和失效是十分有必要的，但是考虑材料和结构断裂失效过程的SPH流固耦合模型还很少。

近来，Silling等[10]提出的近场动力学(Peridynamics,PD)采用非局部积分方程，不连续位移场可以自然地包含在控制方程中，因而在模拟材料中裂纹萌生、扩展方面具有天然的优势，已经成为研究热点。目前，PD已经成功用于混凝土板侵彻问题[11]、充液裂缝分析[12]和准脆性材料的冲击响应和损伤扩展模拟[13]。此外，Gao等[14]利用近场动力学微分算子离散Navier-Stokes方程(简称N-S方程)，实现了基于近场动力学统一求解框架的流固耦合数值模型，但是其计算效率还有待于进一步研究。近场动力学的理论基础不同于无网格法，但是其数值实施策略与无网格法十分类似。Ganzenmüller等[15]详细讨论了在粒子离散条件下PD与SPH的相似性，结果表明，运用SPH或PD模拟基于变形梯度的经典材料模型，采用节点积分时，两者的数值控制方程是一致的，这意味着在某种情况下PD退化为无网格方法。这为联合SPH和PD来求解流固耦合问题提供了理论上的依据。

本文将SPH和PD相结合，发展了一套能够考虑流体自由表面和结构断裂失效的无网格数值模型。基于动态边界条件思想，构建了流固界面上的数值传递策略，最终达到耦合求解的目的。该模型对于结构断裂过程中产生的新流固界面无需额外的处理，计算过程简洁。通过应用人工状态方程、人工黏性、预测-校正时间步进格式和结构区域的时间冻结技术，改善了计算精度和效率。以溃坝冲击下游弹性挡板问题为例，验证了模型在捕捉自由液面和结构破坏方面的有效性，为模拟一定规模的流固耦合问题提供了基础。

1 物理控制方程

采用分区法求解流固耦合问题时，流体域和固体域采用各自独立的控制方程，下面分别给出流体域的SPH控制方程和结构域的PD控制方程。

1.1 SPH控制方程

SPH最初用于解决天体运动过程之间的相互作用问题，Monaghan将其引入流体仿真领域，并率先用SPH来模拟自由表面流[16]。SPH采用如下形式来近似场函数f(x)及其在x位置的空间导数∇·f(x)

(1)

(2)

式中，Ω为流体区域，W(x-x′,h)为权函数或核函数，h为光滑长度。当定义了核函数的支持域后，支持域半径一般取光滑长度的k倍。核函数W一般选择偶函数，并且具有以下三个性质：

归一化条件

(3)

δ函数性质

(4)

紧支属性

W(x-x′,h)=0, |x-x′|>kh

(5)

这里采用Wendland[17]提出的五次样条核函数：

(6)

式中,R=x/h，对于2D问题αd=7/(4πh2)。

上面的〈f(x)〉称为函数f(x)的一个核近似。通过核近似将原场函数及其空间导数近似为核函数支持域上的积分，由此导出SPH粒子近似：

(7)

(8)

式中,mj为粒子质量，ρj为粒子密度，N为粒子i的邻域粒子数。粒子近似将空间上连续的积分运算离散为域内有限个粒子的求和运算。

拉格朗日描述下的流体控制方程为非守恒形式的N-S方程，采用指标记法为：

连续性方程

(9)

动量方程(无外力作用)

(10)

能量方程

(11)

利用SPH核近似和粒子近似可以得到SPH离散形式的流体控制方程，简称SPH控制方程，其中ρ、m、v、P、e分别表示粒子的密度、质量、速度、压力、内能。式中，下标a表示与流体粒子a相关的物理量，下标b表示粒子a的支持域内的其他流体粒子b的相关物理量，Wab表示粒子a的核函数在粒子b处的值，vab=va-vb，rab=ra-rb：

连续性方程

(12)

动量方程

(13)

式中，g为重力加速度，大小为9.81 m/s2。P为流体压力，此处将不可压缩流体视为弱可压缩流体，通过引入人工状态方程来显式求解压力项，避免了求解压力泊松方程。流体粒子的压力由下式计算

(14)

(15)

Πab为黏性项，此处采用人工黏性，其值由下式确定

(16)

式中，α为第一黏度系数，上划线和下标ab表示粒子a和b相应物理量的平均值。

(17)

能量方程

(18)

在通过SPH控制方程移动粒子时，需要引入XSPH修正[16]来有效避免粒子堆积，使得粒子分布更加均匀，同时防止粒子之间发生非物理穿透。粒子规则移动由下式定义

(19)

ε为修正系数，取ε=0.5。

1.2 PD控制方程

根据牛顿第二定律，近场动力学的控制方程可以表达为如下形式

b(x,t)

(20)

式中，H为空间内物质点x的近场范围，u为物质点的位移矢量场，b为一个预定义的体力密度场，ρ为物质密度，f是PD模型中的点对力函数，它包含了材料常数，类似于传统理论中的本构关系，此处的f反映了物质点x′对物质点x的作用效果。

为了能够模拟裂纹在变形体中自然演化的过程，Silling[10]在2005年提出了经典微观弹脆性(prototype microelastic brittle，PMB)本构模型，适用于各向同性材料。记ξ=x′-x，η=u′-u，则PMB模型中点对力函数f(η,ξ)可表达为如下形式

(21)

式中,c是微弹模量。通过对比PD理论中的应变能密度表达式和经典连续介质力学中的应变能密度表达式可得

(22)

(17)中S是物质点对的伸长率，其定义式为

(23)

键型近场动力学中，材料的损伤可以通过物质点之间“键”的断裂来描述。当两个物质点之间的伸长率S达到临界值Scr时，“键”即断裂，物质点之间的相互作用消失。该过程是不可逆的，临界伸长率的具体数值可以通过材料宏观的抗压强度Sc0、抗拉强度St0来定义

(24)

将上述断键准则引入点对力函数的表达式，可得更新后的点对力函数

(25)

当键的破坏累计成一个面的时候，也就形成了宏观的裂纹。每个点定义了局部损伤φ(x,t)

(26)

(27)

通过统计一点近场范围内的断键数目来计算φ(x,t)，从而模拟裂纹的萌生和扩展。键型近场动力学求解的数值过程不依赖网格，通过离散节点携带的信息，求解域被离散为有限的物质点，边界条件可以直接施加到材料的边界物质点层上，其施加边界条件的具体形式如下：

u(x,t)=U0,x∈Rc

b(x,t)=p(x,t),x∈Rc

(28)

2 基于动态边界条件的耦合策略

2.1 动态边界条件

根据动态边界条件的思想，图1(a)中的刚体粒子与流体粒子满足同样的SPH控制方程。对于靠近边界的流体粒子a1，当其支持域内出现刚体粒子时，刚体粒子对其产生的影响与支持域内其他流体粒子对其产生的影响按照同一套控制方程来计算，即按照式(12)更新密度，依据式(13)更新加速度。这时刚体粒子被看作具备流体粒子的特性，这类粒子称为动态边界粒子。与流体粒子不同，刚体粒子的位置不随时间改变，速度始终保持为零，从而实现固定边界的效果。由于刚体边界粒子被当做流体粒子参与计算，不需要做额外的判别，动态边界条件可以在程序中高效统一实现。

(a) 流体粒子受动态边界粒子作用 (b) 弹性体粒子受流体粒子的反作用图1 SPH-PD耦合策略Fig.1 SPH-PD coupling strategy

2.2 流固界面耦合策略

考虑到动态边界条件的便捷性，在SPH-PD模型中同时将其应用到流固界面的处理上。在图1(a)中，对于靠近流体-弹性体界面的流体粒子a2，将弹性体粒子也作为动态边界粒子处理，此时弹性体粒子对流体粒子的影响与刚体粒子对流体粒子的影响完全一样。与刚体粒子不同的是，弹性体粒子的位置和速度随时间变化，其变化量由弹性体区域内的PD控制方程所确定。在弹性体粒子作为动态边界粒子对流体粒子产生作用力时，流体粒子对弹性体粒子有大小相等、方向相反的反作用力，如图1(b)所示。为简明起见，这里只给出了一个流体粒子a2对弹性体的反作用力，实际弹性体受到的反作用力应该是全部受其影响的流体粒子所产生反作用力的总和。计算该作用力，并在一个时间步内作为载荷施加到弹性体上，参与PD控制方程的运算，从而实现流固界面上的数据交换。

当计算弹性体粒子对流体粒子的影响时，依据动态边界条件的思想，处于边界的弹性体粒子被视为流体粒子。因此，在此过程中弹性体粒子的质量、密度、体积均取SPH粒子的值，但其位移和速度不按照SPH控制方程计算；当任一时间步内流体粒子的运动计算完成后，再开始计算弹性体粒子在流体粒子反作用下的运动。这时，弹性体粒子的参数应取实际的固相材料参数，通过求解PD控制方程，即可计算得到结构在流体反作用力下的变形和位移。

3 数值方法

3.1 流体域预测-校正时间步进格式

流体域由SPH离散形式的N-S方程组控制。在第n+1个时间步内，流体在第n步的密度、速度、能量的变化率由式(12)、(13)、(18)计算得到：

(29)

根据第n步的物理量来预测第n+1/2步的物理量

(30)

由预测得到的第n+1/2步的物理量来校正第n+1/2步的物理量

(31)

最终计算得到n+1步的物理量

(32)

3.2 流体域基于CFL条件的变时间步长法

为了防止粒子间的非物理穿透，流体域内的时间步长不能任意选取，必须满足时间步长CFL条件(courant-friedrichs-levy)。进行时间积分时，采用满足CFL条件的最大时间步长。每一步的时间步长Δt由下式确定：

Δt=0.3·min(Δtf,Δtcv)

(33)

3.3 结构域时间冻结技术

结构域由PD积分方程控制，在流体域和结构域选择近似粒子间距的情况下，依据CFL条件确定的结构域最大时间步长ΔtPD往往远小于流体域的最大时间步长ΔtSPH。进行时间积分时，应该取较小的时间步长ΔtPD作为耦合算法整体的时间步长。由于流体域的计算量占总体计算量的较大部分，而且耗时比结构域更长，因此对于流体域，选择时间步长ΔtPD进行运算，会增加很多不必要的迭代次数，降低了算法整体的效率。为了提高计算效率，这里在流体域和结构域分别采用不同的时间步长进行计算。具体的时间步长选择和计算策略为：

在第n个时间步内，若流体域和结构域根据CFL条件计算得到的最大时间步长分别为ΔtSPH和ΔtPD，则流体域按照ΔtSPH迭代一次，结构域按照ΔtPD迭代M次。

(34)

这样，相对于结构域，流体域内的时间就仿佛被冻结了一样，因此称为时间冻结技术。在流体域时间冻结期间，流体粒子对弹性体粒子的反作用保持不变。结构在该载荷作用下，通过求解PD控制方程，迭代M步来计算变形和位移。值得注意的是，在这M次迭代过程中，若任一弹性体粒子的位移累计超过光滑长度h，则结构域的迭代终止，流体域解除时间冻结，进行下一步迭代，重新计算对结构的反作用。

4 数值算例

4.1 模型参数

计算参数的选取对计算结果往往会产生一定影响，这里的算例采用统一模型参数，详见表1。

表1 SPH-PD模型统一参数设定Tab.1 Common Parameters of SPH-PD model

计算中流体域和结构域采用相同粒子间距，流体域时间步长由式(33)确定；结构域时间步长在初始时刻通过CFL条件确定，迭代中不动态改变，通过迭代次数M的变化来间接调整结构域迭代的总次数。

4.2 弹性挡板在水流作用下的变形模拟

为了验证所提出的SPH-PD模型的有效性，首先模拟流体流动和结构变形的相互作用。Antoci等[18]试验研究了弹性挡板控制下的排水问题，并用SPH-SPH统一法进行了数值模拟。采用所提出的SPH-PD模型对该问题进行模拟，并与实验结果数据和SPH-SPH统一法模拟结果进行对比。模型的初始条件设置如图2所示，高为H1的水柱被限制在宽为W1的水槽内，水槽右侧和下方为刚性边界，左侧有一高为L1的弹性挡板，弹性挡板的上端固定，下端自由。模型的几何参数，水的密度ρF，弹性板的密度ρS、弹性模量E、泊松比μ见表2。

图2 弹性挡板下排水问题初始条件设置Fig.2 Initial condition settings of drainage problem with elastic plate

表2 弹性挡板下排水问题参数设定Tab.2 Parameters of drainage problem with elastic plate

将相关的模拟结果与实验结果进行对比，选取了6个典型时刻的结果，如图3所示：第一列为实验结果，第二列为基于本文提出的SPH-PD模型模拟的结果，第三列为Antoci等采用SPH-SPH统一方法模拟的结果。起初弹性挡板在水压作用下发生变形，弹性板自由端向左上方偏移，导致下方出现空隙，从而形成流体排出的通道。在0.16 s以前弹性板的变形在水压作用下不断增大，但随着槽内流体的排出，水位逐渐降低，水压渐渐变小，弹性板的变形开始逐渐减小。通过对比可以看出：基于SPH-PD模型的模拟结果很好地再现了实验中水位的变化和弹性板的变形过程。

(a) 0.04 s

弹性板自由端在水平方向(x方向)和垂直方向(y方向)的位移变化曲线见图4，相应的水位变化过程见图5。可以看出模拟结果与实验数据在变化趋势上吻合良好，因此所提出的SPH-PD模型是有效的。

图4 弹性板自由端水平、垂直方向位移Fig.4 Horizontal and vertical displacement of free end of elastic plate

图5 弹性板右侧水位Fig.5 Water level on the right side of elastic plate

4.3 水坝坍塌冲击作用下弹性挡板的破坏模拟

为了验证SPH-PD模型模拟结构损伤断裂的能力，模拟溃坝冲击下游弹性挡板的破坏过程。在二维平面上有一正方形刚性封闭容器如图6所示。容器边长a=4 m，初始状态下，水柱位于容器内左侧1/4处，水柱的宽度和高度分别为，W=1 m和H=3 m。一弹性挡板位于方形容器中部，底端固定，弹性挡板的厚度和高度分别为，l=0.1 m和h=1 m，与水柱的距离为L=1 m。挡板采用线弹性本构模型，考虑损伤，密度ρ=1 878 kg/m3，弹性模量E=20 MPa，泊松比μ=0.333 3，当键的伸长率超过拉压临界伸长率St0和Sc0时，键即断裂，取St0=0.05和Sc0=-0.09。为了增加流动过程的复杂程度，同时在下游设置一个宽W2=1.9 m、高H2=0.7 m流体区域。水柱在重力作用下自由坍塌冲击中部弹性挡板，利用SPH-PD模型计算流体和结构的运动变形过程，计算总时长为3 s。

图6 溃坝冲击下游弹性挡板模型初始条件设置Fig.6 Initial condition settings of elastic plate with water downstream impacted by dam break

选取了溃坝冲击过程中几个典型时刻的粒子位形及压力分布情况，如图7所示。

图7 弹性挡板在溃坝冲击作用下的断裂过程及t=1.2 s时结构区域损伤场分布Fig.7 Fracture process of elastic plate under dam break impact and damage distribution of structure at t=1.2 s

可以看出，左侧水柱被释放之后撞击弹性挡板，同时形成高达3 m的浪尖。弹性挡板在水流冲击作用下先发生大变形，变形过程中由于局部键的伸长率超出临界值而发生断裂。1.2 s时可以清楚地看到结构区域发生断裂，断裂截面上粒子的损伤量为50%。挡板断裂之后，坍塌的水柱与下游水域发生融合，并携带断裂部分的挡板共同运动。SPH-PD模型可以模拟出流体冲击致结构破坏的典型过程和细节，是有效可靠的。

5 结 论

本文将SPH和PD相结合，提出了一套模拟流固耦合问题的无网格数值方法。所提出的SPH-PD模型采用分区求解策略：SPH用于模拟自由表面流，PD用于模拟结构的变形、断裂行为，通过流固界面上基于动态边界条件的思想实现耦合运算。采用预测-校正时间步进格式，提高了时间积分精度。利用时间冻结技术，改善了流体域内的计算效率，同时采用人工黏性和XSPH修正，改善了粒子分布的均匀性。利用提出的SPH-PD模型，模拟了二维溃坝冲击下游弹性挡板问题，可以高效捕捉到流体自由表面破碎后出现的浪花翻滚、飞溅、融合等复杂现象，同时实现了挡板中裂纹的萌生、扩展直到断裂失效的全过程模拟，从而验证了模型的有效性和鲁棒性。

SPH-PD模型具有无网格粒子类方法的典型优势，可以拓展到更多包含自由表面流、结构损伤失效的流固耦合问题，可以为灾难仿真、影视特效等领域提供了一种新的模拟方法。