一种管涌现象的数值分析方法

摘要

本发明公开了一种管涌现象的数值分析方法，包括步骤：基于流体动力学相似准则结合介质对流体的相间力作用项，获得渗透变形相似准则；基于渗透变形相似准则结合颗粒流方法建立细观颗粒模型；所述模型满足重力水平与原型一致，以及颗粒级配缩放且保持相同孔隙率条件下随机生成骨架及填充颗粒；对模型进行管涌过程模拟，在模型内填充颗粒达到平衡状态或骨架颗粒发生破坏时将模型内参数信息输出；将所述渗透变形相似准则与模型内参数信息结合获得实际数值，反馈至原型以分析管涌现象。本发明符合达西和非达西定律，在同时存在两种流态时模型均能够对原型进行真实描述，克服了常规模型或离心试验不能实现颗粒粒径与模型同比缩放的难题、分析准确。

说明书

技术领域

本发明涉及一种管涌现象的数值分析方法，属于土木工程技术领域。

背景技术

中国是世界上遭受洪水灾害最为严重的国家之一，长江中下游平原区洪灾尤为频 繁、严重。堤防的破坏包括管涌、漫顶、崩岸和整体失稳等各种形式，大量洪灾资料表 明，堤基管涌危害最大，堤坝发生管涌破坏时，渗流场具有强烈的空间特性，目前对管 涌发展过程的渗流场空间性状方面的研究较少，归其原因，一方面，因为当前对管涌发 生机理的认识还不够深入；另一方面，由于工程水力条件的复杂性，对渗流场的分析较 为困难。

当前，工程界和理论界对管涌问题的研究多偏重于对土的抗渗性能方面研究，然而 随着管涌研究工作的开展，对渗流场具有强烈的时空特性已逐渐达成一种共识。如能描 述渗流过程中材料几何特性和水力特性的复杂变化，得到水头、水力梯度、颗粒位移场 和孔隙率等量的动态变化，则将为管涌过程的研究提供强有力的支持。

本发明之前，对管涌的研究分为理论研究、室内试验以及数值试验研究，而渗透变 形试验中一般假定渗流过程符合达西定律，当域内出现非达西渗流时，采用常规模型或 离心模型获得的结论与原型存在一定的偏离。在离心试验中，如采用缩小的颗粒粒径进 行试验会面临一个问题：原型颗粒直径范围在砂土范围内，而粒径与模型长度比尺同比 缩小后，粒径可能已经处于粘土范围，而粘土和砂土中发生渗流现象时，土与水之间的 作用力面临较大的差别，粘土表面可能存在强弱结合水的作用，而砂土中则不会存在这 种情况，也即采用同比缩小粒径的方法虽然机理上实现了相似，但是细观范畴的物理、 化学作用的影响很可能使得结果和预期差别较大；同时采用缩尺关系选择，在保证模型 与原型级配相似的同时并不能消除原型土与模型之间的力学性质差异；特别是在选用满 足特定条件的流体时面临较大的困难。针对管涌的数值分析大多采用有限元模拟渗流 场，连续性方法对管涌发生前相对稳定渗流阶段的模拟可行，当达到临界状态后，无法 考虑颗粒流失引起的土体几何特性和水力特性相互影响的这种复杂水土相互作用，故无 法全面解释管涌发生机理。管涌产生过程中复杂的水土相互作用决定了管涌的发展是一 个非线性的动态过程，目前尚无公认的最合适的研究方法与理论。已有学者开始尝试利 用显微摄像可视化跟踪技术结合数字图像识别分析手段，对砂土管涌的整个过程进行跟 踪记录，从细观层面揭示水土相互作用贯穿于管涌发展的全过程，然而该技术仅局限于 室内试验层面，对于实际管涌渗透破坏过程的分析尚不具备适应性。

发明内容

本发明所要解决的技术问题在于克服现有技术的不足，提供一种准确便捷的、能真 实反映实际管涌渗透破坏过程的一种管涌现象的数值分析方法，基于现有流体动力学相 似准则，考虑达西-非达西状态下介质与流体间的相间力作用，基于完备的多孔介质流 体动力学方程推导渗透破坏相似准则，给出砂土发生渗透变形时遵循的相似准则，使得 当计算区域内出现达西、非达西渗流，或域内同时存在该两种流态时模型均能够对原型 进行真实描述。

本发明具体采用以下技术方案解决上述技术问题：

一种管涌现象的数值分析方法，包括以下步骤：

步骤(1)、基于流体动力学相似准则结合介质对流体的相间力作用项，通过引入 扩展的D-B-F方程考虑渗流的达西-非达西效应，获得砂土发生渗透变形时遵循的相似准 则；

步骤(2)、基于渗透变形相似准则结合颗粒流方法建立与原型相似的细观颗粒模型； 所述细观颗粒模型满足重力水平与原型一致，以及在细观颗粒模型中按照与原型相同的 颗粒级配曲线将原型土样颗粒按渗透变形相似准则缩放且保持相同孔隙率条件下随机 生成骨架及填充颗粒；

步骤(3)、对所述细观颗粒模型进行管涌过程模拟，在模型内填充颗粒达到平衡状 态或骨架颗粒发生破坏时将模型内参数信息输出；

步骤(4)、将所述渗透变形相似准则与模型内参数信息结合获得实际数值，并将 实际数值反馈至原型，以分析管涌现象。

进一步地，作为本发明的一种优选技术方案：所述步骤(1)砂土发生渗透变形时遵 循的相似准则为：$\left(\begin{array}{cc}{\alpha }_n=1& \left(a\right)\\ {\alpha }_{\sigma }/{{\alpha }_{\gamma }\alpha }_l=1& \left(b\right)\\ {\alpha }_v{\alpha }_t/{\alpha }_l=1& \left(c\right)\\ {\alpha }_v^2/{\alpha }_g{a}_l=1& \left(d\right)\\ {\alpha }_{\mathrm{\rho f}}/{\alpha }_{\mathrm{\rho f}}{\alpha }_v^2=1& \left(e\right)\\ {\alpha }_v{\alpha }_l/{\alpha }_{\mu }=1& \left(f\right)\\ {\alpha }_v/{\alpha }_k=1& \left(g\right)\\ {\alpha }_{\bar{k}}={\alpha }_g{\alpha }_l={\alpha }_g{\alpha }_d& \left(h\right)\end{array}\right)$

式中α_n、α_σ、α_γ、α_t、α_v、α_l、α_g、α_ρf、α_μ、α_k、α_d分别表示：孔隙 率、应力、重度、时间、速度、长度、加速度、流体密度、运动粘滞系数、渗透系数相 似常数、等效渗透系数相似常数、粒径缩放系数；(a)式表征模型与原型中孔隙率必须一 致；(b)式表征模型与原型应力状态一致，包含有效应力及总应力对应一致；(c)～(f)式分 别表征流体动力学相似中的谐时准则、重力相似准则、压力相似准则、粘滞力相似准则； (g)、(h)表征孔隙流体处于层流及紊流状态时相间动力交换系数对多孔介质中流体动力 相似的影响。

进一步地，作为本发明的一种优选技术方案：所述步骤(2)中细观颗粒模型按渗透变 形相似准则缩放满足：颗粒粒径与模型同比缩放和流体密度与模型尺寸同比缩放、粘滞 系数比尺为模型比尺的平方。

进一步地，作为本发明的一种优选技术方案：所述步骤(3)中模型内参数变化数据包 括模型内颗粒速度、应力分布、流量、流速、压力数据。

进一步地，作为本发明的一种优选技术方案：所述步骤(4)将所述渗透变形相似 准则与模型内参数变化数据结合获得获得与原型一致的临界水力梯度，并参照通过渗透 变形相似准则获得的与惯性时间一致的时间比尺分析管涌过程。

本发明采用上述技术方案，能产生如下技术效果：

本发明管涌现象的数值分析方法的优点和效果在于：

(1)采用常规模型及离心模型研究渗流问题时均是基于达西定律，在模型中出现 非达西渗流时仍采用基于达西定律的相关结论是不合适的。本发明克服了采用常规模型 试验或离心模型试验研究砂土渗透变形问题时面临的局限性。

(2)基于现有多孔介质流体动力学相似准则，考虑介质对流体的相间力作用项， 给出了砂土发生渗透变形时遵循的相似准则，其核心思想是确保颗粒粒径与模型同比缩 放，并保证粘滞系数比尺与长度比尺一致，以此为前提建立的渗流模型与原型相似。当 计算区域内出现达西、非达西渗流，或域内同时存在该两种流态时模型均能够对原型进 行真实描述，因而当原型中流态由达西向非达西转化时，基于渗透变形相似准则建立的 模型对于该过程也能准确地描述。

(3)颗粒流方法以颗粒为基本计算单元，模拟颗粒的缩放具有独到的优势，结合 其研究渗透变形的优势，克服了常规模型试验或离心试验不能实现颗粒粒径与模型同比 缩放的难题。

(4)所述研究在分析平台进行，约束条件极少，应用范围广。避免了室内式样制 备、仪器设备操作等繁琐过程，提高了试验效率，节约了研究成本。

本发明工序简明，在颗粒流理论的基础上，基于多孔介质流体动力学相似准则，考 虑介质对流体的作用项，给出了砂土发生渗透变形时遵循的相似准则，建立模型操作方 便、计算简单、经济高效、分析更加准确。

附图说明

图1为本发明管涌现象的数值分析方法中建立细观颗粒模型示意图。

图2为本发明细观颗粒模型边界条件示意图。

图3为本发明管涌现象的数值分析方法的流程图。

图4A为细观颗粒模型中孔口A发生管涌时的临界水力梯度与室内试验的对比计算 结果图。

图4B为细观颗粒模型中孔口B发生管涌时的临界水力梯度与室内试验的对比计算 结果图。

图4C为细观颗粒模型中孔口C发生管涌时的临界水力梯度与室内试验的对比计算 结果图。

图5A为本发明试样I-1在C号孔口开启时发生渗透破坏初始时刻的示意图。

图5B为本发明试样I-1在C号孔口开启时发生渗透破坏细颗粒流失时刻的示意图。

图5C为本发明试样I-1在C号孔口开启时发生渗透破坏粗颗粒流失时刻的示意图。

图5D为本发明试样I-1在C号孔口开启时模型破坏时刻的示意图。

具体实施方式

下面结合说明书附图对本发明的实施方式进行描述。

本发明设计了一种管涌现象的数值分析方法，实施的具体步骤如图3所示，包括：

步骤(1)、基于流体动力学相似准则结合介质对流体的作用项，推导获得砂土发 生渗透变形时遵循的相似准则。具体如下：

(1.1)基于现有多孔介质流体动力学相似准则，考虑介质对流体的相间动力作用 项，引入扩展的D-B-F方程；

(1.2)考虑渗流的达西-非达西效应，根据几何相似，运动相似，动力相似，推导 砂土发生渗透变形时遵循的相似准则。具体推导如下：

在流体动力学中，假设不考虑温度场的变化，对于不可压缩流体，求解由散粒体组 成的多孔介质中流场和压力场时必须遵守质量守恒方程：

$\frac{\partial n}{\partial t}=-(▿·n\bar{v})---\left(1\right)$

公式(1)中，n为多孔介质的孔隙率，为流体速度，

平均动量守恒方程为：

$\left(\begin{array}{c}\frac{\partial \left(n\bar{v}\right)}{\partial t}+(▿·n\bar{v}\bar{v})=\mathrm{ng}-\frac{n}{{\rho }_f}▿p\\ -\frac{n}{{\rho }_f}▿·\tau -\frac{\beta (\bar{v}-\bar{u})}{{\rho }_f}\end{array}\right)---\left(2\right)$

区别于N-S方程，本发明考虑了介质对流体的作用项，详见方程右端最后一项。

本发明将式(2)称为扩展的Darcy-Brinkman-Forchheimer方程，左端第一项表示局部惯性 力，第二项表示变位惯性力，两项之和表征惯性力，也即加速度项；方程右端各项分别 表征体积力、压力、粘滞力及固液相间作用力对加速度的贡献。

式(2)中n为多孔介质的孔隙率，为流体速度，为颗粒平均运动速度，τ为粘 滞力张量，g为重力加速度，ρf为流体密度，p为压力，β为相间动力交换系数,根据 Ergun方程及Wen和Yu(1966)提出公式确定，具体表达如下:

$\beta =(150\frac{{(1-n)}^2}{n{\bar{d}}^2}{\mu }_f+1.75\frac{(1-n)}{\bar{d}}{\rho }_f|\bar{v}-\bar{u}\left|\right)n\le 0.8---\left(3\right)$

公式(3)中μ_f为表征流体粘滞系数，为颗粒等效粒径。

相应的压力梯度和表观流速可以表示为：

$\frac{\mathrm{dp}}{\mathrm{dx}}=-150\frac{{(1-n)}^2}{n^3{\bar{d}}^2}{v}_f{\rho }_f|\bar{v}-\bar{u}|-1.75\frac{(1-n)}{n^3\bar{d}}{\rho }_f{|\bar{v}-\bar{u}|}^2---\left(4\right)$

公式(4)中dp/dx表征压力梯度。

方程第一项即为层流状态时的Kozeny-Carman公式，第二项表示压力梯度与速度的 平方项成比例。当两相均不占主导时，压力梯度与速度的关系介于1～2之间，此即渗透 系数比值为重力比尺的x次方的原因。完全处于紊流状态时如仅考虑速度平方项，忽略 第一项，可以将式(5)表示成：

$v^2=\frac{n^3\bar{d}g}{1.75(1-n)}i=\bar{k}·i---\left(5\right)$

为方便阐述，将式(5)中定义为等效渗透系数，i为水力梯度，g为重力加速度， 结合层流状态时的达西定律表达式v＝k·i，k为渗透系数。令β₁、β₂分别代表式(3)右端 的两项，β₁、β₂分别为完全层流及完全紊流条件下相间力系数，则式(2)可以表示为：

$\beta ={\beta }_1+{\beta }_2=\frac{{\rho }_f{\mathrm{gn}}^2}{k}+\frac{{\rho }_f{\mathrm{gn}}^3}{\bar{k}}|\bar{v}-\bar{u}|---\left(6\right)$

定义时间、速度、长度、加速度、流体密度、运动粘滞系数、相间动力交换相似常 数分别为：α_t、α_v、α_l、α_g、α_ρf、α_μ、α_β，将相似常数代入动量守恒方程(2)得：

$\left(\begin{array}{c}\frac{{\alpha }_n{\alpha }_v}{{\alpha }_t}\frac{\partial \left(n\bar{v}\right)}{\partial t}+\frac{{\alpha }_n{\alpha }_v^2}{{\alpha }_l}(▿·n\bar{v}\bar{v})={\alpha }_n{\alpha }_g\mathrm{ng}\\ -\frac{{\alpha }_n{\alpha }_p}{{\alpha }_{{\rho }_f}{\alpha }_l{\rho }_f}\frac{n}{{\rho }_f}▿p-{\alpha }_n{\alpha }_v\frac{{\alpha }_v}{{\alpha }_d^2{\alpha }_{{\rho }_f}}\frac{n}{{\rho }_f}▿·\tau -\frac{{\alpha }_{\beta }{\alpha }_v(\bar{v}-\bar{u})}{{\alpha }_{{\rho }_f}{\rho }_f}\end{array}\right)---\left(7\right)$

公式(7)中，α_d为粒径缩放系数，τ的定义为流体剪应力。α_n为孔隙率相似常数。

遵循流体动力学相似基本要求，在相应点所作用的同名力比值以单位体积计算应当 相等，结合式(6)有：

$\frac{{\alpha }_v}{{\alpha }_t}=\frac{{\alpha }_v^2}{{\alpha }_l}{\alpha }_g=\frac{{\alpha }_{{\rho }_f}}{{\alpha }_{{\rho }_f}{\alpha }_l}={\alpha }_v\frac{{\alpha }_v}{{\alpha }_l^2}={\alpha }_{{\beta }_1}{\alpha }_v={\alpha }_{{\beta }_2}{\alpha }_v---\left(8\right)$

公式(8)中分别为完全层流及完全紊流条件下相间力系数相似常数。

将式(6)代入式(8)，并结合重力相似，孔隙率相似，最终得到：

$\left(\begin{array}{cc}{\alpha }_n=1& \left(a\right)\\ {\alpha }_{\sigma }/{\alpha }_{\gamma }{\alpha }_l=1& \left(b\right)\\ {\alpha }_v{\alpha }_t/{\alpha }_l=1& \left(c\right)\\ {\alpha }_v^2/{\alpha }_g{a}_l=1& \left(d\right)\\ {\alpha }_{\mathrm{\rho f}}/{\alpha }_{\mathrm{\rho f}}{\alpha }_v^2=1& \left(e\right)\\ {\alpha }_v{\alpha }_l/{\alpha }_{\mu }=1& \left(f\right)\\ {\alpha }_v/{\alpha }_k=1& \left(g\right)\\ {\alpha }_{\bar{k}}={\alpha }_g{\alpha }_l={\alpha }_g{\alpha }_d& \left(h\right)\end{array}\right)---\left(9\right)$

公式中(9)中，α_k为渗透系数相似常数，为等效渗透系数相似常数。

式(9)中(a)式表征模型与原型中孔隙率必须一致；(b)式表征模型与原型应力状态一 致，包含有效应力及总应力对应一致；(c)～(f)式分别表征流体动力学相似中的谐时准则、 重力相似准则、压力相似准则、粘滞力相似准则；(g)、(h)表征孔隙流体处于层流及紊 流状态时相间动力交换系数对多孔介质中流体动力相似的影响。

步骤(2)、基于渗透变形相似准则结合颗粒流方法建立与原型相似的细观颗粒模型。 具体如下：

(2.1)建立室内细观颗粒模型，如图1所示，上游边界与x＝0的平面重合，在上 游边界施加水头，其它表面均设为不透水边界，上表面预留A、B、C号孔口(开始均墙 体覆盖)，大小遵循模型比尺，以模拟不同的出流口，孔口中心位置设为零压力边界， 模拟不同孔口出流时分别删除相应位置的墙体单元，该模型边界条件示意图如图2所 示，模型左侧为水头边界，孔口为零压力边界，其它边界均设为刚性不透水非滑移边界。

具体参数见表1；

(2.2)按照步骤1所得的渗透变形相似准则，结合原型中颗粒级配，进行细观颗 粒模型中粒径范围确定；

(2.3)然后，在模型中按照与原型相同的颗粒级配曲线将原型土样颗粒按相似准 则缩放，缩放包括颗粒粒径与模型同比缩放、流体密度与模型尺寸同比缩放，粘滞系数 比尺为模型比尺的平方。为避免接触冲刷现象的产生，将模型上边界摩擦系数设置为较 大值，生成边界模型后采用PFC3D颗粒随机生成器先产生粗颗粒，后产生细颗粒，生 成过程中保持孔隙率与原型一致，以保持相同孔隙率条件下随机生成骨架及填充颗粒。

(2.4)由细观颗粒模型可完美满足渗透变形相似准则，模型必须满足重力水平与 原型一致，通过提高重力水平以抵消模型缩尺影响。确保细观颗粒模型达到初始应力平 衡状态，模拟无渗流条件下初始状态，模型达到初始稳定状态后首先使模型处于静水压 力中。

步骤(3)、对细观颗粒模型施加水头边界进行渗透变形，模拟管涌过程，记录模型 内参数变化数据。具体为：

(3.1)对细观颗粒模型施加水头边界进行渗透变形分析，分别开启不同的孔口(删 除相应位置的上覆墙体)，并将其设置为自由面，模拟不同孔口发生渗透破坏，通过监 控手段获取模型内相应参数如孔口流出的颗粒体积、孔口流量及单元内孔隙率，颗粒位 置，颗粒运移速度，压力，流速等，分析模型内颗粒速度、应力分布及渗透特性，分析 流体速度、压力及颗粒与流体之间的耦合力，以在模拟过程中实时获取模型内相应参数 信息。

(3.2)在每一级水力梯度下，孔口流速达到稳定状态或者不能保持稳定状态时该 级水头下渗透分析结束。如未发生破坏，施加下一级水头，启用处于静水压力中模型， 重新开始分析，以消除上一级水头作用下部分颗粒已经发生运移的影响，直至模拟对象 不能保持稳定或孔口出流已经稳定为止，否则继续缓慢施加下一级水头。直至在模型内 填充颗粒达到平衡状态或模型骨架颗粒发生破坏，将所述模型内参数数据信息输出。

步骤(4)、将所述渗透变形相似准则与模型内参数变化数据结合获得实际数值，根 据实际数值反馈至原型，以分析管涌现象。具体为：

运用步骤(1)推导得到的渗透变形相似准则，对步骤(3)所得模型内相应参数数 据，如颗粒速度、应力分布、流量、流速等数据进行分析，基于推导的相似准则，获得 实际数值，根据实际数值分析管涌现象。获得实际数值在分析过程中根据需要获得对应 参数数值，本发明中优选获得与原型一致的临界水力梯度，并参照通过渗透变形相似准 则获得的与惯性时间一致的时间比尺，根据同一时间比尺下管涌发展过程的变化，对实 际管涌过程进行表征。

所述渗透变形相似准则与模型内参数变化数据结合，具体以模型中孔口出流速度为 例，根据式(9)中重力相似准则及粘滞力相似准则，均可得到模型与原型中流速一致， 如此，模型中流速及代表实际流速；此时如继续流量与流速关系，根据：Q＝VAt，结合 模型与原型流速一致结论，由于流过截面积比尺为1/N²，时间比尺为1/N，可以得到模 型与原型流量比尺为1/N³。由此可获得原型流量的实际数值，再加以分析。

由此，本方法基于多孔介质流体动力学相似准则，考虑介质对流体的相间力作用项， 给出了砂土发生渗透变形时遵循的相似准则，其核心思想是确保颗粒粒径与模型同比缩 放，以此为前提建立的渗流模型与原型基本相似。当计算区域内出现达西、非达西渗流， 或域内同时存在该两种流态时模型均能够对原型进行真实描述，因而当原型中流态由达 西向非达西转化时，基于该相似准则建立的模型对于该过程也能准确地描述，由此形成 的数值分析方法能够解决现有管涌分析过程中无法同时满足西、非达西效应的技术问 题，克服了常规模型试验或离心试验不能实现颗粒粒径与模型同比缩放的难题，能准确 的运用于实际管涌分析过程。

为了验证本发明的管涌现象的数值分析方法能够模拟管涌现象，对模型内的渗透变 形情况进行分析，特以如下实验数据进行验证。

首先，按照表1所示参数，分别通过室内试验获得6组土样在不同孔口发生渗透 破坏时的临界水力梯度。然后按照本发明的方法，建立细观颗粒模型分析得到对应水力 梯度。

表观临界水力梯度随孔口A、B、C所在位置的变化关系详见图4A、图4B、图 4C。可以看出，两者基本吻合，符合临界水力梯度随着细料含量以及填充密实度的增 加而增大的趋势，与工程实际相吻合。图5A为试样I-1在C号孔口开启时发生渗透破 坏初始时刻的示意图，初始阶段，模型内颗粒初始平衡；随后试样中仅仅有少量的细颗 粒流失，如图5B所示，为C号孔口开启时发生渗透破坏中细颗粒流失时刻的示意图； 随着细颗粒向孔口区域聚集，细颗粒流失一定的时间后，个别粗颗粒开始被细颗粒携裹 带出孔口，随着时间推移，粗颗粒开始逐渐流失，如图5C所示，为C号孔口开启时发 生渗透破坏粗颗粒流失时刻的示意图；最终孔口区域出现较大的凹坑，如图5D所示， 为模型破坏时刻的示意图，其与工程实际相吻合。同时，整个数值模拟过程全部由计算 机执行，避免了室内试验中制备试样、操作仪器的繁琐过程，提高了试验效率，节约了 研究成本。

表1模型参数

由此实验验证采用颗粒流方法模拟管涌破坏可以对模型内各参数进行动态观测， 可以反映管涌破坏的细观过程。证明本发明可准确、便捷地对管涌现象进行解释和分析。 应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之 后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定 的范围。