基于B样条物质点法的溃坝流模拟研究

徐云卿, 周晓敏,2, 赵世一, 徐聖飞, 孙 政,2

(1. 江西理工大学 土木与测绘工程学院, 江西 赣州 341000;2. 江西理工大学 江西省环境岩土与工程灾害控制重点实验室, 江西 赣州 341000)

0 引 言

水库大坝是土木工程中常见的工程结构,在防洪、能源生产、蓄水等方面都发挥着重要作用．但大坝在具有社会经济效益的同时,也存在着溃坝的潜在风险,一旦发生溃坝,将会对人民生命财产安全造成极大的威胁．因此,开展水坝的溃决过程及其演进规律研究,对于溃坝洪水灾害的预测和防治具有重要意义．

溃坝流是一种典型的自由表面流动问题,研究内容主要包括波前到达时间、波前流速、波前位置和沿程各处的水位等[1]．目前,研究方法主要有理论分析[2-4]、物理试验[5-6]、数值模拟[7-8]．随着数学理论和计算机技术的高速发展,数值模拟广泛应用于溃坝水流的研究．数值模拟研究相比较于物理试验具有更高的灵活性,不局限于场地设置,计算周期相对较短,能够大幅度降低试验成本等优点．Han等[8]采用Preissmann法的一维模型耦合采用有限差分法的溃坝二维模型求解溃坝扩散波动情况;胡四一和谭维炎[9]以TVD格式预测了高坝瞬溃的溃坝波演进过程;王大国等[10]采用CBOS有限元法建立了水波模型,精确模拟和分析了下游河床无水时溃坝模型的自由水面运动特征;张建伟等[11]采用光滑粒子流体动力学方法(SPH)建立了包含下游水位的溃坝模型,对于溃坝流冲击过程中波前到达时间、波前流速、波前位置能够做到精确模拟．

但溃坝流涉及到水波演进大变形问题,有网格类的有限差分法(FDM)[8]、有限体积法(FVM)[7]、有限元法(FEM)[10]等方法在模拟时会出现网格畸变并引起误差．尽管无网格化的SPH[11]避免了网格畸变,但由于追踪的是网格边界上质量、动量和能量的通量流动,需求解非线性对流项,增加了求解难度,且不易于追踪各质点的运动时间历程,使得计算效率低．因此,需要开发一种可以减少误差和提高计算效率的数值模拟方法．

物质点法(MPM)[12]作为一种相对新型的粒子型算法,其具有Lagrange 粒子和 Euler背景网格双重描述,既包含了Lagrange粒子型无网格算法的优势,又可避免网格畸变和求解非线性对流项,可方便地追踪自由表面和粒子信息的时间历程; 同时利用Euler型背景网格求解动量方程和梯度,具有计算效率和处理基本边界条件[13-14]的优势．物质点法已广泛应用各个工程领域问题[15-19]中．然而,传统的物质点法采用线性插值形函数来实现物质点和背景网格节点之间的信息映射．由于网格边界处的形函数梯度不连续,当物质点穿越背景网格边界时,将产生网格穿越误差,引起数值振荡．Gan等[20]采用高阶B样条基函数[21]替换线性插值基函数,节点自由度空间代替背景网格节点,减少了网格穿越误差,提高了物质点法的计算精度和收敛性．目前B样条物质点法(BSMPM)已经用于模拟自由表面流体问题,如具有自由表面[22-23]的复杂流动问题,流体-结构相互作用问题[24]等．

1 计 算 方 法

1.1 弱可压缩B样条物质点法

1.1.1 B样条基函数

B样条基函数的定义有很多种,尽管所得的表达式不同,但其本质是一致的．本文选用计算稳定、易于理解的Cox-de Boor递归公式[25]来计算B样条基函数．在Cox-de Boor递归公式中,零阶基函数Si,0(q=0)如下所示:

(1)

对于q≥1,基函数由Cox-de Boor递归公式定义为

(2)

式中,ξi是第i个节点,i=1,2,…,n+q+1,n是基函数的总数目,q是基函数的阶数．由式(1)可以看出,如果基函数的阶数为零(q=0),则基函数呈阶梯分布,即在第i节点区间[ξi,ξi+1)上基函数Si,0(ξ)等于1,而其他节点区间等于0．

采用Si,0(ξ)和Si+1,0(ξ),通过式(2)可以得到Si,1(ξ);当所有Si,1(ξ)计算完成,采取同样的方法可以得到Si,2(ξ),持续这个计算过程直到所有需要的Si,q(ξ)计算完毕,其递归过程如图1所示．

图1 Cox-de Boor递归示意图Fig. 1 Schematic diagram of the Cox-de Boor recursion

(a) 参数网格空间 (b) 节点自由度空间(a) The parametric grid space(b) The node degree of freedom space图2 网格空间示意图Fig. 2 Illustration of grid spaces

假设0/0=0,Si,q(ξ)的导数也可以用递归公式定义,即

(3)

本文研究的是二维溃坝流模拟实验,基函数是双变量基函数,假定节点向量N={ξ1,ξ2,…,ξn+q,ξn+q+1}和H={η1,η2,…,ηm+p,ηm+p+1},ξ和η是参数的方向,n和m分别为沿图2(a)中ξ,η方向的基函数数目．利用张量积结构双变量基函数定义为

S(i,j),(q,p)(ξ,η)=Si,q(ξ)·Sj,p(η),

(4)

式中,i=1,2,…,n;j=1,2,…,m．三变量B样条基函数也可以类似于双变量基函数推导得到．图2(a)参数网格中的空心圆和实心圆分别表示节点向量和物质点．

图2(b)给出的是节点自由度空间,图中(i,j)表示节点自由度空间中对应的(i,j)自由度;物质点(ξ,η)在节点自由度(i,j)上对应的B样条基函数记作S(i,j),(q,p)(ξ,η),其中(q,p)分别表示在节点矢量N和H方向上B样条基函数的阶数．

1.1.2 B样条物质点法

B样条物质点法是基于物质点法改进而来的,B样条插值基函数阶数为一阶时与物质点法中的传统线性插值形函数是一样的,即B样条插值基函数阶数为一阶时的B样条物质点法就是物质点法[26],两者求解步骤基本一致．两者的不同之处在于:在物质点法中,物质点信息采用线性插值形函数映射到相应背景网格节点上,在背景网格节点上施加本质边界条件．而在B样条物质点法中,物质点信息采用B样条基函数映射到节点自由度空间的各节点自由度上,在节点自由度空间施加本质边界条件．其具体算法实现过程如下:

1) 划分参数网格,将研究对象离散成一组物质点．

2) 初始化物质点的位置、质量、动量等物质信息．

3) 采用B样条基函数,在第k个计算时间步,将物质点上的质量和动量信息映射到节点自由度空间的节点自由度上:

(5)

(6)

(7)

(8)

式中,∇S(i,j),(q,p)表示B样条插值基函数的梯度;σ为Cauchy应力张量;b和t分别为体力和表面力矢量;V为体积;h表示面力边界厚度．

5) 计算节点自由度上的加速度和速度:

(9)

(10)

式中,a表示加速度矢量．

6) 在节点自由度空间,施加本质边界条件．

7) 将各节点自由度的加速度和速度信息,通过B样条基函数映射回各物质点,得到第k+1个时间步上物质点的速度和位置:

(11)

(12)

8) 将更新后的各物质点速度重新映射回各节点自由度上,并计算物质点第k+1个时间步上的应变增量Δε:

(13)

(14)

9) 更新物质点的密度、应力、应变、位置等信息,开始下一计算步．

1.2 弱可压缩流体本构模型

溃坝流问题中流体的本构模型可表示为[27]

(15)

(16)

2 计算模型和结果分析

2.1 计算模型

模型布置如图3所示,图3(a)中h0是上游水库水位高度,l0是大坝水库区域的长度,l和h分别是水箱的长度和高度,v是阀门上升速度．图3(b)中L(t)是t时刻流体前缘位置,h(t)是t时刻上游水位高度．边界条件取液柱的左右两侧和上下底面都为可滑移边界条件, 初始下游河床无水．在模拟计算中取重力加速度g=9.81 m/s2,水的密度ρ=1 000 kg/m3．阀门上升速度v=2.0 m/s,人工声速c和时间步长Δt分别为100 m/s和1×10-5s．

(a) 初始时刻(b) t时的水流位置(a) The initial setup (b) The flow position at time t图3 溃坝流问题示意图Fig. 3 Illustration of the dam break flow problem

本文基于三维程序模拟二维平面应变问题, 即在z方向仅设置一个背景网格, 单元网格大小为0.02 m,每个背景单元网格内布置2×2×2个物质点粒子．为了方便进行对比,采用无量纲化的L*-T*,H*-T*曲线表示．L*代表无量纲化后流体波前到达位置,H*代表无量纲化后给定位置水位高程,T*代表流体流动时间:

(17)

(18)

(19)

2.2 模拟结果分析

图4 溃坝流体波前位置随时间的变化Fig. 4 Variation of the dam-break fluid wavefront position with time

为了说明B样条插值基函数阶数对结果的影响,图5(a)、(b)、(c)和(d)分别给出了给定位置h1,h2,h3和h4水位高程随时间的变化．从图中可以明显看出: 1) 模拟得到的结果与实验结果曲线单调性相近,模拟结果与实验结果基本吻合．2) 随着B样条插值基函数阶数的增加模拟结果曲线变得愈发平稳和光滑,在B样条插值基函数阶数为3阶时,模拟实验所得到的给定位置h1,h2,h3和h4水位高程随时间变化关系曲线最为平稳和光滑．3) 与实验结果相比,不同B样条插值基函数阶数下,模拟给定位置h1,h2,h3和h4水位高程随时间变化关系曲线过程存在较为明显的锯齿变化．B样条插值基函数阶数的增加产生影响和模拟水位高程随时间变化关系曲线结果存在明显锯齿变化的主要原因是: 1) 在B样条物质点法中,随着B样条插值基函数阶数的增加,节点矢量内的节点和基函数个数将会相应地增加．2) 在样条物质点法中,随着B样条插值基函数阶数的增加,基函数更为光滑,同时基函数在边界处具有更大的梯度,提高了B样条物质点法的求解精度和收敛性．3) 在数值模拟中,对水位高度的变化,是通过捕捉固定水平位置处一个背景网格单元范围内物质点y坐标的最大值得到的,然而试验数据是基于100组试验的平均值,因此图中试验结果较光滑,而数值模拟结果呈锯齿变化,但两者吻合较好．由图4和图5可知,本文的模拟方法可以准确地模拟溃坝流体的传播过程以及水库区和下游区的水位高程变化规律．

图5 初始水位高度h0=0.3 m,(a)、(b)、(c)和(d)分别为给定位置h1,h2,h3和h4水位高程随时间的变化Fig. 5 Initial water level height h0=0.3 m, (a),(b),(c) and (d) denoting the changes of the water level elevation at given positions h1,h2,h3 and h4 with time, respectively

为了进一步说明溃坝流的演进过程和B样条插值基函数阶数对结果的影响,图6给出了不同时刻速度剖面下B样条物质点法1阶、2阶、3阶的模拟结果和实验结果的对比．从图中可以看出: 1) 模拟结果的速度剖面与实验结果有良好的一致性,且随着B样条插值基函数阶数的增加,速度剖面愈发一致; 2) 随着时间的推移,溃坝水流的波前位置的移动速度加快,右侧水库内水位逐渐降低．

图6 初始水位高度h0=0.3 m,(a)、(b)、(c)和(d)分别为BSMPM 1阶、2阶、3阶的模拟结果和实验结果在t=0.32 s(T*=1.83),t=0.41 s(T*=2.34),t=0.46 s(T*=2.63)下的速度云图Fig. 6 Initial water level height h0=0.3 m, (a),(b),(c) and (d) denoting the velocity profiles of the BSMPM 1st-order, 2nd-order and 3rd-order simulation results and experimental results at t=0.32 s(T*=1.83),t=0.41 s(T*=2.34),t=0.46 s(T*=2.63), respectively

2.3 计算效率

通过改变网格尺寸,对2.1小节所述的溃坝流问题,1阶、2阶和3阶基函数下B样条物质点法求解耗时进行分析．表1给出了1阶、2阶和3阶基函数下B样条物质点法,在网格尺寸为0.01 m,0.02 m和0.04 m下的单步CPU计算耗时,图7绘制了其变化关系曲线．本文程序运算所用计算机为64位CentOS Linux 7系统、Inter Xeon Gold 6226R @ 2.90 GHz×64 CPU、128 G内存;程序基于InterFortran90编辑器,串行计算,CPU耗时通过CPU—TIME命令得到．

表1 不同网格尺寸下,1阶、2阶和3阶基函数下B样条物质点法的求解耗时

由表1和图7可以看出: 1) 随着基函数阶数的增加,B样条物质点法求解计算耗时呈约1.5倍增长; 2) 不同阶次B样条物质点法的计算耗时随背景网格尺寸的增长率基本一致,约呈线性增长．

图7 单步CPU计算耗时随网格尺寸和基函数阶数的变化Fig. 7 Variation of the single-step CPU computation time with the grid size and the order of basis functions

3 结 论

本文基于B样条物质点法,通过引入人工状态方程,开展了弱可压缩B样条物质点法模拟溃坝流问题的研究,分析了B样条插值基函数阶数对模拟结果的影响．从模拟结果可以得出:

1) B样条物质点法所得溃坝流问题模拟结果与物理实验结果基本吻合．

2) B样条物质点法可以很好地模拟溃坝流体的流动特性．溃坝水流的波前速度,随着时间的推进越来越快．对于给定位置的h1高程,随着时间的推进逐渐下降,而h2,h3和h4的高程,随着时间的推进逐渐上升,验证了B样条物质点法模拟溃坝流问题的可行性．

3) 在改变B样条插值基函数阶数的条件下,通过对比1阶、2阶和3阶基函数下的溃坝流体波前位置、波前速度及给定位置的高程变化的模拟结果和物理实验结果,得出随着基函数阶数的增加,模拟结果与实验结果拟合度提高．

4) 在改变背景网格尺寸条件下,通过对比1阶、2阶和3阶基函数B样条物质点法计算耗时,得出相较物质点法,更高阶B样条物质点法的计算耗时约呈1.5倍增长;随着背景网格数目的增加,各阶次B样条物质点法的计算耗时增长率与传统物质点法基本一致,约呈线性增长．