1、 应用数学和力学编委会，：基于 样条物质点法的溃坝流模拟研究徐云卿，周晓敏，赵世一，徐聖飞，孙 政，（江西理工大学 土木与测绘工程学院，江西 赣州；江西理工大学 江西省环境岩土与工程灾害控制重点实验室，江西 赣州）摘要：溃坝流是水利工程中常见的一种自由表面流动问题，准确模拟溃坝流问题具有重要的工程意义 样条物质点法（）作为一种物质点法（，）的改进算法，其提高了物质点的计算精度和收敛性，且在自由表面流动问题中有独特的算法优势 基于 样条物质点法，通过引入人工状态方程，发展了一种弱可压缩 样条物质点法（）；开展溃坝流问题的模拟研究，分析 样条插值基函数阶数对模拟结果的影响 结果表明：模拟所得流体波
2、前位置、波前流速及给定位置处的高程变化与已有实验结果吻合较好；同时，随着 样条基函数阶数的增加，模拟结果与试验结果吻合度逐渐提高 随着基函数阶数的增加，计算耗时呈约 倍增长；不同阶次 样条物质点法的计算耗时随背景网格尺寸的增长率基本一致，约呈线性增长 验证了弱可压缩 样条物质点法模拟溃坝流问题的有效性，为模拟溃坝流问题提供了一种新的思路和方法 关 键 词：溃坝流；样条物质点法；自由表面流；人工状态方程中图分类号：文献标志码：，（，；，）：，（），（），（），应用数学和力学 卷 期 年 月 ，收稿日期：；修订日期：基金项目：国家自然科学基金项目（；）；江西省主要学科学术和技术带头人青年项目（）；
3、国家级大学生创新创业训练项目（）；江西省大学生创新创业训练项目（）作者简介：徐云卿（），男，硕士（：）；孙政（），男，副教授，博士，硕士生导师（通讯作者：）引用格式：徐云卿，周晓敏，赵世一，徐聖飞，孙政 基于 样条物质点法的溃坝流模拟研究 应用数学和力学，（）：，：；引 言水库大坝是土木工程中常见的工程结构，在防洪、能源生产、蓄水等方面都发挥着重要作用 但大坝在具有社会经济效益的同时，也存在着溃坝的潜在风险，一旦发生溃坝，将会对人民生命财产安全造成极大的威胁 因此，开展水坝的溃决过程及其演进规律研究，对于溃坝洪水灾害的预测和防治具有重要意义 溃坝流是一种典型的自由表面流动问题，研究内容主要包括
4、波前到达时间、波前流速、波前位置和沿程各处的水位等 目前，研究方法主要有理论分析、物理试验、数值模拟 随着数学理论和计算机技术的高速发展，数值模拟广泛应用于溃坝水流的研究 数值模拟研究相比较于物理试验具有更高的灵活性，不局限于场地设置，计算周期相对较短，能够大幅度降低试验成本等优点 等采用 法的一维模型耦合采用有限差分法的溃坝二维模型求解溃坝扩散波动情况；胡四一和谭维炎以 格式预测了高坝瞬溃的溃坝波演进过程；王大国等采用 有限元法建立了水波模型，精确模拟和分析了下游河床无水时溃坝模型的自由水面运动特征；张建伟等采用光滑粒子流体动力学方法（）建立了包含下游水位的溃坝模型，对于溃坝流冲击过程中波前
5、到达时间、波前流速、波前位置能够做到精确模拟 但溃坝流涉及到水波演进大变形问题，有网格类的有限差分法（）、有限体积法（）、有限元法（）等方法在模拟时会出现网格畸变并引起误差 尽管无网格化的 避免了网格畸变，但由于追踪的是网格边界上质量、动量和能量的通量流动，需求解非线性对流项，增加了求解难度，且不易于追踪各质点的运动时间历程，使得计算效率低 因此，需要开发一种可以减少误差和提高计算效率的数值模拟方法 物质点法（）作为一种相对新型的粒子型算法，其具有 粒子和 背景网格双重描述，既包含了 粒子型无网格算法的优势，又可避免网格畸变和求解非线性对流项，可方便地追踪自由表面和粒子信息的时间历程；同时利用
6、 型背景网格求解动量方程和梯度，具有计算效率和处理基本边界条件的优势 物质点法已广泛应用各个工程领域问题中 然而，传统的物质点法采用线性插值形函数来实现物质点和背景网格节点之间的信息映射 由于网格边界处的形函数梯度不连续，当物质点穿越背景网格边界时，将产生网格穿越误差，引起数值振荡 等采用高阶 样条基函数替换线性插值基函数，节点自由度空间代替背景网格节点，减少了网格穿越误差，提高了物质点法的计算精度和收敛性 目前 样条物质点法（）已经用于模拟自由表面流体问题，如具有自由表面的复杂流动问题，流体结构相互作用问题等 本文主要研究两个方面：一方面，运用弱可压缩 样条物质点法（），对溃坝流问题进行模拟
7、研究，并将模拟的流体波前位置、波前速度及给定位置的高程变化与 等的实验结果进行对比分析，验证运用 样条物质点法模拟溃坝流问题的有效性；另一方面，基于不同阶次 样条插值形函数对溃坝流问题进行模拟，并与已有实验结果进行对比，分析阶数对模拟精度和计算效率的影响 计 算 方 法 弱可压缩 样条物质点法 样条基函数 样条基函数的定义有很多种，尽管所得的表达式不同，但其本质是一致的 本文选用计算稳定、易于理解的 递归公式来计算 样条基函数 在 递归公式中，零阶基函数，（）如下应 用 数 学 和 力 学 年 第 卷所示：，（）对于 ，基函数由 递归公式定义为，（），（），（），（）式中，是第 个节点，是基函
8、数的总数目，是基函数的阶数 由式（）可以看出，如果基函数的阶数为零（），则基函数呈阶梯分布，即在第 节点区间，）上基函数，（）等于，而其他节点区间等于 采用，（）和，（），通过式（）可以得到，（）；当所有，（）计算完成，采取同样的方法可以得到，（），持续这个计算过程直到所有需要的，（）计算完毕，其递归过程如图 所示 图 递归示意图 （）参数网格空间（）节点自由度空间（）（）图 网格空间示意图 假设，（）的导数也可以用递归公式定义，即，（），（），（）（）本文研究的是二维溃坝流模拟实验，基函数是双变量基函数，假定节点向量 ，第 期 徐云卿，等：基于 样条物质点法的溃坝流模拟研究 和 ，和 是参数
9、的方向，和 分别为沿图（）中，方向的基函数数目 利用张量积结构双变量基函数定义为（，），（，）（，），（），（），（）式中，；，三变量 样条基函数也可以类似于双变量基函数推导得到 图（）参数网格中的空心圆和实心圆分别表示节点向量和物质点 图（）给出的是节点自由度空间，图中（，）表示节点自由度空间中对应的（，）自由度；物质点（，）在节点自由度（，）上对应的 样条基函数记作（，），（，）（，），其中（，）分别表示在节点矢量 和方向上 样条基函数的阶数 样条物质点法 样条物质点法是基于物质点法改进而来的，样条插值基函数阶数为一阶时与物质点法中的传统线性插值形函数是一样的，即 样条插值基函数阶数为一阶
10、时的 样条物质点法就是物质点法，两者求解步骤基本一致 两者的不同之处在于：在物质点法中，物质点信息采用线性插值形函数映射到相应背景网格节点上，在背景网格节点上施加本质边界条件 而在 样条物质点法中，物质点信息采用 样条基函数映射到节点自由度空间的各节点自由度上，在节点自由度空间施加本质边界条件 其具体算法实现过程如下：）划分参数网格，将研究对象离散成一组物质点）初始化物质点的位置、质量、动量等物质信息）采用 样条基函数，在第 个计算时间步，将物质点上的质量和动量信息映射到节点自由度空间的节点自由度上：（，）（，），（，）（），（）（，）（，），（，）（）（，），（）式中上标 表示计算时间步；下
11、标 示物质点的变量；表示物质点 的位置矢量；为质量；为速度矢量 ）计算节点自由度上的内力（，）和外力（，）：（，）（，），（，）（），（）（，）（，），（，）（）（，），（，）（），（）式中，（，），（，）表示 样条插值基函数的梯度；为 应力张量；和 分别为体力和表面力矢量；为体积；表示面力边界厚度）计算节点自由度上的加速度和速度：（，）（，）（，）（，），（）（，）（，）（，），（）式中，表示加速度矢量）在节点自由度空间，施加本质边界条件）将各节点自由度的加速度和速度信息，通过 样条基函数映射回各物质点，得到第 个时间步上物质点的速度和位置：（，），（，）（）（，），（）（，），（，）（）（
12、，）（）将更新后的各物质点速度重新映射回各节点自由度上，并计算物质点第 个时间步上的应变增量：应 用 数 学 和 力 学 年 第 卷（，）（，），（，）（）（，），（）（，）（，），（，）（）（，）（，），（，）（）（）更新物质点的密度、应力、应变、位置等信息，开始下一计算步 弱可压缩流体本构模型溃坝流问题中流体的本构模型可表示为 （），（）式中，为流体动力黏度系数，对于 流体动力黏度系数 为常数；为黏度参数，对没有体积黏度的流体 ；为可以用状态方程（）来确定的流体压力，为 函数张量，为应变率张量 为了加快计算速度，本文通过引入 将不可压缩流体看作弱可压缩流体，求解流体压力 的公式表示为 ，（
13、）式中，是流体的初始密度；（）；表示人工声速，本文选取 虽然小于实际值，但取为流体最大速度的 倍左右 通过这种方法，在确保流体的近似不可压缩条件下，可以允许更大的步长，在保证计算精度的同时可以大大提高计算效率 计算模型和结果分析 计算模型模型布置如图 所示，图（）中 是上游水库水位高度，是大坝水库区域的长度，和 分别是水箱的长度和高度，是阀门上升速度 图（）中（）是 时刻流体前缘位置，（）是 时刻上游水位高度 边界条件取液柱的左右两侧和上下底面都为可滑移边界条件，初始下游河床无水 在模拟计算中取重力加速度 ，水的密度 阀门上升速度 ，人工声速 和时间步长 分别为 和 （）初始时刻（）时的水流位
14、置（）（）图 溃坝流问题示意图 本文基于三维程序模拟二维平面应变问题，即在 方向仅设置一个背景网格，单元网格大小为 ，每个背景单元网格内布置 个物质点粒子 为了方便进行对比，采用无量纲化的，曲线表示 代表无量纲化后流体波前到达位置，代表无量纲化后给定位置水位高程，代表流体流动时间：（），（）第 期 徐云卿，等：基于 样条物质点法的溃坝流模拟研究 ，（），（）模拟结果分析图 给出了流体波前到达位置随时间变化的关系 图中，点线是 等的物理模型实验结果，虚线、实线和点划线分别为 样条插值基函数阶数 阶、阶和 阶时的模拟结果 从图中可明显看出：）模拟结果与实验结果基本吻合，曲线单调性相近，对于流体波前
15、位置到达水箱右壁的时刻模拟与实验结果基本一致）模拟结果与实验结果相比，模拟结果早期的流体波前位置变化的速度略低于实验结果，中后期的速度略大于实验结果 模拟和实验流体波前位置变化的速度产生差异的主要原因是：）实验早期由于快速拆除阀门会产生较大的剪切应力，并产生射流，从而使流体波前位置变化的速度比模拟结果速度快）实验中后期由于流体与水箱两侧壁碰撞使得流体能量损失，从而使流体波前位置变化的速度比模拟结果速度慢 图 溃坝流体波前位置随时间的变化 为了说明 样条插值基函数阶数对结果的影响，图（）、（）、（）和（）分别给出了给定位置，和 水位高程随时间的变化 从图中可以明显看出：）模拟得到的结果与实验结果
16、曲线单调性相近，模拟结果与实验结果基本吻合）随着 样条插值基函数阶数的增加模拟结果曲线变得愈发平稳和光滑，在 样条插值基函数阶数为 阶时，模拟实验所得到的给定位置，和 水位高程随时间变化关系曲线最为平稳和光滑）与实验结果相比，不同 样条插值基函数阶数下，模拟给定位置，和水位高程随时间变化关系曲线过程存在较为明显的锯齿变化 样条插值基函数阶数的增加产生影响和模拟水位高程随时间变化关系曲线结果存在明显锯齿变化的主要原因是：）在 样条物质点法中，随着 样条插值基函数阶数的增加，节点矢量内的节点和基函数个数将会相应地增加）在样条物质点法中，随着 样条插值基函数阶数的增加，基函数更为光滑，同时基函数在边
17、界处具有更大的梯度，提高了 样条物质点法的求解精度和收敛性）在数值模拟中，对水位高度的变化，是通过捕捉固定水平位置处一个背景网格单元范围内物质点 坐标的最大值得到的，然而试验数据是基于 组试验的平均值，因此图中试验结果较光滑，而数值模拟结果呈锯齿变化，但两者吻合较好 由图 和图 可知，本文的模拟方法可以准确地模拟溃坝流体的传播过程以及水库区和下游区的水位高程变化规律 为了进一步说明溃坝流的演进过程和 样条插值基函数阶数对结果的影响，图 给出了不同时刻速度剖面下 样条物质点法 阶、阶、阶的模拟结果和实验结果的对比 从图中可以看出：）模拟结果的速度剖面与实验结果有良好的一致性，且随着 样条插值基函
18、数阶数的增加，速度剖面愈发一致；）随着时间的推移，溃坝水流的波前位置的移动速度加快，右侧水库内水位逐渐降低 应 用 数 学 和 力 学 年 第 卷图 初始水位高度 ，（）、（）、（）和（）分别为给定位置，和 水位高程随时间的变化 ，（），（），（）（），图 初始水位高度 ，（）、（）、（）和（）分别为 阶、阶、阶的模拟结果和实验结果在 （），（），（）下的速度云图 ，（），（），（）（），（），（），（），注 为了解释图中的颜色，读者可以参考本文的电子网页版本 计算效率通过改变网格尺寸，对 小节所述的溃坝流问题，阶、阶和 阶基函数下 样条物质点法求解耗时进行分析 表 给出了 阶、阶和 阶基函数
19、下 样条物质点法，在网格尺寸为 ，和 下的单步 计算耗时，图 绘制了其变化关系曲线 本文程序运算所用计算机为 位 系统、第 期 徐云卿，等：基于 样条物质点法的溃坝流模拟研究 、内存；程序基于 编辑器，串行计算，耗时通过 命令得到 表 不同网格尺寸下，阶、阶和 阶基函数下 样条物质点法的求解耗时 ，由表 和图 可以看出：）随着基函数阶数的增加，样条物质点法求解计算耗时呈约 倍增长；）不同阶次 样条物质点法的计算耗时随背景网格尺寸的增长率基本一致，约呈线性增长 图 单步 计算耗时随网格尺寸和基函数阶数的变化 结 论本文基于 样条物质点法，通过引入人工状态方程，开展了弱可压缩 样条物质点法模拟溃坝
20、流问题的研究，分析了 样条插值基函数阶数对模拟结果的影响 从模拟结果可以得出：）样条物质点法所得溃坝流问题模拟结果与物理实验结果基本吻合）样条物质点法可以很好地模拟溃坝流体的流动特性 溃坝水流的波前速度，随着时间的推进越来越快 对于给定位置的 高程，随着时间的推进逐渐下降，而，和 的高程，随着时间的推进逐渐上升，验证了 样条物质点法模拟溃坝流问题的可行性）在改变 样条插值基函数阶数的条件下，通过对比 阶、阶和 阶基函数下的溃坝流体波前位置、波前速度及给定位置的高程变化的模拟结果和物理实验结果，得出随着基函数阶数的增加，模拟结果与实验结果拟合度提高）在改变背景网格尺寸条件下，通过对比 阶、阶和
21、阶基函数 样条物质点法计算耗时，得出相较物质点法，更高阶 样条物质点法的计算耗时约呈 倍增长；随着背景网格数目的增加，各阶次 样条物质点法的计算耗时增长率与传统物质点法基本一致，约呈线性增长 参考文献（）：谢任之 溃坝水力学 济南：山东科学技术出版社，（：，（）应 用 数 学 和 力 学 年 第 卷 ，（）：，（）：郑素佩，李霄，赵青宇，等 求解二维浅水波方程的旋转混合格式 应用数学和力学，（）：（，（）：（），（）：，：鲍远林，陈秀荣，程冰，等 有限体积 方法在二维溃坝中的应用 应用数学，（）：（，（）：（），（）：胡四一，谭维炎 用 格式预测溃坝洪水波的演进 水利学报，（）：（，（）：（）
22、王大国，水庆象，等 基于 有限元溃坝流数值模型 工程力学，（）：（，（）：（）张建伟，杜宇，陈海舟 溃坝水流冲击水垫塘的 模拟 水电能源科学，（）：（，（）：（），（）：，（）：，：，：王宇新，李晓杰，王小红，等 炸药爆轰物质点法三维数值模拟 应用数学和力学，（）：（，（）：（）张雄，刘岩，张帆，等 极端变形问题的物质点法研究进展 计算力学学报，（）：（，（）：（）马上，张雄，邱信明 超高速碰撞问题的三维物质点法 爆炸与冲击，（）：（，（）：（）郑勇刚，顾元宪，陈震 薄膜破坏过程数值模拟的 方法 力学学报，（）：（，（）：（）周晓敏，孙政 非 流体的物质点法模拟研究 应用数学和力学，（）：（第 期 徐云卿，等：基于 样条物质点法的溃坝流模拟研究，（）：（），（）：，：，（）：，：，（）：，孙政，周晓敏 物质点法改进算法及其工程应用 长沙：中南大学出版社，（，：，（）周光坰，严宗毅，许世雄，等 流体力学（上）北京：高等教育出版社，（，（）：，（）应 用 数 学 和 力 学 年 第 卷