摘要：

针对黄土颗粒流模型细观参数标定过程中试错法计算量大的现状，采用敏感性分析与回归分析相结合的方法，在分析了黄土细观因素与宏观力学响应对应关系和多因素非线性联合反演的基础上，快速确定了黄土的细观参数，并得出如下结论：黄土的抗剪强度τ¯τ¯c随半径乘子λ¯λ¯的增大呈指数函数规律增大，随摩擦系数μ的增大呈线性增大，随法向黏结强度σ¯σ¯c的增大呈幂函数规律增大；内摩擦角φ′随摩擦系数μ的增大呈对数函数规律增大；黏聚力C随半径乘子λ¯λ¯的增大呈线性增大，随法向黏结强度σ¯σ¯c的增大呈对数函数规律增大。数值试验与室内试验结果误差较小，证明了此方法在黄土细观参数标定中的可行性，为其他材料的细观力学参数标定提供了新思路，也为黄土地区滑坡、泥石流和融水降水侵蚀等地质灾害的颗粒流模拟分析提供了参数参考。

关键词：

细观参数标定;正交试验;敏感性分析;颗粒流;

董建鹏(1987—),男，硕士研究生，主要从事岩土工程数值模拟研究。

*李辉(1973—),男，正高级工程师，博士，主要从事岩土防灾减灾研究。

基金：

青海省科技厅项目(2020-ZJ-718);

引用：

董建鹏,李辉.黄土颗粒流宏细观对应关系与参数标定方法研究[J].水利水电技术(中英文),2022,53(4):180-191.DONGJianpeng,ibrationmethodofloessparticleflow[J].WaterResourcesandHydropowerEngineering,2022,53(4):180-191.

0引言

黄土高原地区水土流失、滑坡、泥石流等自然灾害对沿线路桥和人民群众生命财产安全造成巨大威胁，对黄土力学性质研究具有重要的实际意义。颗粒流方法具有成本低、可重复等优点，往往被作为物理试验的补充和延伸。在对黄土的力学性质进行颗粒流数值模拟研究时，一般采用试错法通过不断对比宏观力学与变形响应来标定模型的细观参数，由于宏细观参数之间关系的复杂性，试错的过程枯燥且耗时，国内外学者对此进行了深入的探索和改进。在参数标定方法的探索方面：郭鸿等建立了二维规则排列颗粒的平行黏结模型，通过数学推导得到宏观弹性模量、抗拉强度与颗粒黏结强度、颗粒行列数、颗粒半径等参数的定量关系；郭万里等提出了一种同时确定粗粒料模拟孔隙率和摩擦系数的方法。邓树新等采用试验设计法研究了硬岩的细观参数的标定问题；王晋伟等提出了正交-等值线法对堆石料的细观参数进行了标定；刘佳明等根据细观参数与宏观力学参数的拟合关系提出了适用于类岩材料的细观参数标定方法；陈亚东等在总结了大量数值试验结果的基础上，提出了与砂土孔隙率、内摩擦角和压缩模量相匹配的细观参数的标定方法。在宏细观对应关系研究方面：学者对岩石类材料和粗粒土的研究较多，普遍通过模拟单轴拉伸、压缩试验、巴西劈裂试验、大型三轴试验等分析细观参数与宏观响应的关系，得到宏观弹性模量与细观有效模量呈线性关系、单轴抗压强度与泊松比主要受法向与切向黏结强度和摩擦系数控制等结论。对黏性土的研究相对较少，一般通过分析三轴试验的应力-应变曲线研究峰值强度、残余强度、剪缩剪胀特性、泊松比等指标与细观参数的关系，认为峰值强度和残余强度主要与围压、孔隙率、摩擦系数和黏结强度有关；应变软化性质主要受摩擦系数、黏结强度和围压影响；剪缩、剪胀量随摩擦系数和黏结强度的增大而增加。也有学者研究了颗粒形状和结构性对强度变形特征的影响，认为形状单一会导致数值试验结果偏低,力链分布、裂缝数量、应变软化硬化特征与结构性有密切关系。对黄土的研究方面，同霄等分析了重力条件、颗粒密度、刚度和模型尺寸效应对强度曲线的影响，蒋明镜等提出了非饱和结构性黄土三维接触模型，分析了平面应变条件下有效围压与黄土宏观响应的关系。

上述研究在参数标定方法，以及颗粒流细观参数与宏观强度和变形响应的关系方面取得了诸多重要成果，然而仍有以下方面亟待更细致地探讨：第一，目前学者较多采用三轴试验的颗粒流模型，直剪试验的颗粒流模拟较少，且直剪试验的颗粒流模型均是针对PFC2D,事实上只有极少数工程问题能简化为二维平面的应力和应变问题；第二，目前的研究多集中在岩石、粗粒土等材料，对黄土的研究较少且偏重于对泊松比、应力-应变曲线等变形特征的拟合，对抗剪强度、内摩擦角、黏聚力等强度特征研究较少，对黄土细观参数的标定仍停留在试错法上，缺少更快速准确的标定方法；第三，数值模拟试验分析仅采用正交设计，由于正交设计因素分散、均匀等特点，无法单独从正交数值试验结果中得到细观参数与宏观响应之间明确的回归关系，需要在正交试验的基础上辅以其他手段；第四，颗粒流细观参数众多、模型建立较复杂，以往学者对模型的简化、待标定细观参数的筛选、细观参数取值范围等方面内容的讨论较少，给初学者造成很大困扰。

基于以上讨论，本文建立了直剪试验的PFC3D颗粒流模型，针对黄土数值模拟研究中的细观参数标定问题，提出了敏感性分析和回归分析相结合的方法，在分析了单细观因素与宏观力学响应对应关系和多因素非线性联合反演的基础上快速确定了黄土的细观参数，为黄土高原地区边坡的颗粒流模拟分析提供了参数基础；同时也详细介绍了模型的简化、细观参数的筛选，讨论了细观参数取值范围与宏观响应的关系，为其他材料的细观力学参数标定提供了新思路。

1直接剪切试验的颗粒流模拟

1.1建立模型

颗粒流细观参数的标定，采用数值试验与室内试验宏观响应相比较的方法，标定流程如图1所示。即：第1步，在室内通过级配实验、密度试验、弹性模量试验等确定颗粒的粒径、密度和弹性模量；第2步，通过假设与理论分析减少待标定细观参数的个数，降低标定难度；第3步，设计正交试验并利用试验结果对细观参数进行宏观响应的敏感性分析，找出起控制作用的细观参数；第4步，固定其他参数不变，依次改变控制性参数带入模型计算，总结出单因素回归公式和反演公式并分析其与宏观响应的关系；第5步，对调自变量与因变量位置，将宏观响应视为自变量、控制性参数视为因变量，在参考单因素反演公式特点的基础上，总结出多因素非线性联合反演公式；第6步，室内试验结果带入联合反演公式得到对应的控制性参数，其余细观参数按参照值带入模型计算得到数值模拟结果；第7步，对数值模拟结果和室内试验结果进行误差分析，参考单因素回归关系微调结果。

图1细观参数标定流程

PB模型是根据Potyondy和Cundall等提出的模型建立的，由非活动阻尼器和参考间隙为零的平行黏结组成(见图2)。与线性接触模型的点接触不同，PB模型在两个相互胶结的接触面(piece)之间可以传递力与力矩，其中胶结宽度有半径乘子λ¯控制，该特点能模拟形状不规则黄土颗粒间相互错动产生的粒间弯矩。由于法向黏结强度σ¯c的存在，PB模型中相互胶结的接触面能承受拉力，使接触面的间隙gs大于零时接触依然存在，该特点能模拟黄土颗粒间因范德华力、分子力和毛细力等形成的胶结作用。PB模型中的线性分量和黏结分量是平行关系，两者相互独立。当拉力超过法向黏结强度σ¯c或剪切力超过切向黏结强度c¯时胶结断裂，PB模型退化为线性接触模型(见图3)。

图2PB模型

图3PB模型的退化模型

PB模型细观参数众多，为减小参数标定的难度，通过假设和理论分析减少待标定的细观参数数量，并对计算模型做出适当简化。PFC计算程序中，计算时步Δt(timestep)的估算方法如下

式中，M为颗粒的总质量；K为颗粒的总刚度；k为颗粒刚度；ρs为颗粒密度；r为颗粒半径；β为调整系数。

计算时步Δt表示颗粒的力与位移状态每更新一次所间隔的机器时间。程序每向前推进一个时间步长Δt,所有颗粒的力与位移状态都要被更新一次，颗粒数量越多，每个时间步长Δt消耗的现实时间越久。由式(1)可知，时间步长Δt随颗粒半径r和颗粒密度ρs的减小而减小，随颗粒刚度k的增大而减小。粒径比指颗粒最大粒径与最小粒径的比值，当粒径比很大时，小颗粒游弋于大颗粒之间的缝隙中形成“悬浮”颗粒(见图4)。“悬浮”颗粒的存在导致模型在预压成样阶段的不平衡力比aratio较难达到稳定的值，且对模型的宏观力学性质影响较小[20]。因此，应从增大颗粒半径r、增大颗粒密度ρs、降低颗粒刚度k、减少颗粒数量、减少“悬浮”颗粒数量等方面入手，在保证结果准确性的前提下提高计算效率。

图4大颗粒间的悬浮颗粒

黄土以粗粉粒为主体构成土体骨架，粒径范围在0.01～0.05mm之间。研究表明，当模型的特征长度比值L/R足够大(L/R40)时，L/R对粗粒土抗剪强度没有影响。经过多次试算，选择粒径放大50倍，在半径为0.0005～0.00125m之间采用均匀概率分布生成颗粒，此时颗粒数目较多、计算流畅，且计算结果与粒径放大40倍、20倍、5倍时相比差别较小。颗粒密度ρs指单位体积颗粒的质量，不随颗粒尺寸改变而改变，因此，颗粒密度ρs采用土颗粒密度经验值2.7×103kg/m3。孔隙率n与室内试验结果一致。

刚度比K*指颗粒法向刚度kn与切向刚度ks的比值，即K*=kn/ks。大量学者的研究成果[20,21,22,23]表明，刚度比与泊松比γ之间存在线性关系，即

式中，a、b为待定系数。

在土体的宏观力学性质中，弹性模量E、剪切模量G与泊松比γ存在如下关系

由于泊松比γ是无量纲的量，刚度比K*与泊松比γ之间可建立类似的关系[24],即

黄土的泊松比γ取经验值0.3,带入式(4)得到刚度比K*=2.6α。为简化起见，本文取α=1。参考以往学者的研究成果,令K∗=K¯¯¯∗,E∗=E¯¯¯∗Κ*=Κ¯*,E*=E¯*。至此，待标定的细观参数减少至6个(见表1)。

直剪试验的数值模拟如图5所示，为确保准静态剪切，设定加载速率为0.001m/s,以降低颗粒的惯性效应,通过墙体径向伺服保证加载时围压偏差在1%以内。

图5直剪试验的数值模拟

1.2正交试验设计

正交试验设计(orthogonalexperimentaldesign)是研究多因素多水平的一种试验设计方法，它是根据正交性从全面试验中挑选出部分有代表性的点进行试验，这些有代表性的点具备了“均匀分散，齐整可比”的特点，大幅度减少了试验次数，是一种高效率、快速、经济的试验设计方法。本次正交试验涉及6个因素，每个因素设置5个水平，正交试验设计如表2所列，正交矩阵序列及数值模拟结果如表3所列，其中抗剪强度τ¯τ¯c采用法向应力为100kPa时的值。

2宏细观参数对应关系分析

2.1敏感性分析

使用SPSS(statisticalpackageforthesocialsciences)统计分析软件对各因素的主效应进行多因素方差分析，得到F统计量和相伴概率Sig.值。F统计量是指在零假设成立的情况下，符合F分布的统计量，F统计量等于参数组间差异与组内差异的比值，当比值超过某一限定值，认为参数组间存在显著差异；相伴概率Sig.值指原假设发生的概率，F统计量越大，则相伴概率Sig.值越小。根据F统计量大小及相伴概率Sig.值，可以判断各细观参数对宏观参数的影响程度，若,则产生了非常显著的影响；若0.01Sig.≤0.05,则产生了显著影响；若,则影响不显著，分析结果如图6所示。

图6敏感性分析

由图6(a)可知，半径乘子λ¯、切向黏结强度c¯、法向黏结强度σ¯c、有效模量E*的概率相伴值Sig.均小于0.05,对抗剪强度τ¯c产生显著影响。从F统计量看，摩擦系数μ与摩擦角ϕ¯对抗剪强度τ¯c的影响程度较弱。按照对抗剪强度τ¯c的影响程度排序：半径乘子λ¯法向黏结强度σ¯c切向黏结强度c¯有效模量E*摩擦系数μ摩擦角ϕ¯。

由图6(b)可知，有效模量E*、法向黏结强度σ¯c、半径乘子λ¯、切向黏结强度c¯、摩擦系数μ、摩擦角ϕ¯的概率相伴值Sig.均小于0.05,对内摩擦角φ′均产生显著影响；摩擦系数μ、摩擦角ϕ¯、切向黏结强度c¯、半径乘子λ¯的概率相伴值Sig.均小于0.01,对内摩擦角φ′的影响非常显著。按照对内摩擦角φ′的影响程度排序：摩擦系数μ摩擦角ϕ¯切向黏结强度c¯半径乘子λ¯有效模量E*法向黏结强度σ¯c。

由图6(c)可知，半径乘子λ¯、切向黏结强度c¯、法向黏结强度σ¯c、有效模量E*的概率相伴值Sig.均小于0.05,对黏聚力C产生显著影响，从F统计量看，摩擦系数μ与摩擦角ϕ¯的影响程度较弱，摩擦系数μ几乎不对黏聚力C产生影响。按照对黏聚力C的影响显著程度排序：半径乘子λ¯法向黏结强度σ¯c切向黏结强度c¯有效模量E*摩擦角ϕ¯摩擦系数μ。

2.2单因素回归分析

为进一步获取宏观参数与其控制因素之间的定量关系，以编号17的正交试验结果为参照组，仅选用Sig.≤0.01所对应的细观参数，对宏观响应进行单因素回归分析。

以半径乘子λ¯为例，保持其余细观参数不变，使半径乘子λ¯从0.5等步长增至1.5并分别代入模型计算，为了得到更精确的结果，将半径乘子λ¯的水平数增加至9个，得到半径乘子λ¯单因素回归分析试验结果如表4所列，半径乘子λ¯与宏观响应的关系如图7所示。

图7半径乘子λ¯与宏观响应的对应关系

由图7(a)可以看出，抗剪强度τ¯c随半径乘子λ¯的增大大致呈指数函数规律增大，按指数函数对其进行拟合，得到拟合关系式为τ¯c=50.837λ¯,R2=0.9801τ¯c=50.837λ¯,R2=0.9801。由图7(b)可以看出，黏聚力C随半径乘子λ¯的增大大致呈线性增大，拟合关系式为C=27.37-49.91,R2=0.9465。由图7(c)可以看出，半径乘子λ¯与内摩擦角φ′没有明显的拟合关系。

由图6可知，除了半径乘子λ¯,Sig.≤0.01λ¯,Sig.≤0.01的细观参数还有摩擦系数μ、摩擦角ϕ¯、切向黏结强度c¯和法向黏结强度σ¯c。同样的方法，保持其余细观参数不变，分别等步长改变摩擦系数μ、摩擦角ϕ¯、切向黏结强度c¯和法向黏结强度σ¯c的值带入模型计算，得到其与宏观参数之间的变化曲线如图8～图11所示。

图8摩擦系数μ与宏观响应的对应关系

图9摩擦角ϕ¯与宏观响应的对应关系

图10切向黏结强度c¯与宏观响应的对应关系

图11法向黏结强度σ¯c与宏观响应的对应关系

由图8(a)可以看出：内摩擦角φ′随摩擦系数μ的增大而增大，增大速度逐渐减小并趋于稳定。选择对数函数对其进行拟合，得到拟合关系式为φ′=9.8722lnμ+12.85,R2=0.9962。由图8(b)可以看出，抗剪强度τ¯cτ¯c随摩擦系数μ的增大呈线性增大，其拟合关系式为τ¯c=5.0042μ+36.675,R2=0.9898τ¯c=5.0042μ+36.675,R2=0.9898。由图8(c)可以看出，黏聚力C与摩擦系数μ没有明显的相关性。由图9、10可以看出：抗剪强度τ¯c、内摩擦角φ′和黏聚力C与摩擦角ϕ¯和切向黏结强度c¯均无明显相关性。由图11(a)可以看出，内摩擦角φ′与法向黏结强度σ¯c无明显相关性；由图11(b)可以看出，抗剪强度τ¯c随法向黏结强度σ¯c增大而增大，其拟合关系式为τ¯c=252.48σ¯,R2=0.9897τ¯c=252.48σ¯,R2=0.9897;由图11(c)可以看出，黏聚力C随法向黏结强度σ¯σ¯c的增大而增大，增大速度逐渐减小并趋于稳定，按对数函数对其进行拟合，C=130.59lnσ¯c+218.89,R2=0.9935C=130.59lnσ¯c+218.89,R2=0.9935。

3数值试验与室内试验结果对比

3.1多因素非线性联合反演

对前文单因素拟合公式做逆运算，得到细观参数的单因素反演公式如表5所列。

由前文的敏感性分析可知，一个宏观参数同时受多个细观参数影响，为了能直接通过宏观参数计算得到细观参数，参考细观参数单因素反演公式的规律，将宏观参数作为自变量、细观参数作为因变量，利用SPSS软件仅对Sig.≤0.01的细观参数进行多因素非线性回归分析，得到细观参数的多因素非线性联合反演公式如表6所列。

3.2试验结果对比与讨论

为验证上述参数标定方法的正确性，对张建华等所撰写论文中1#试坑1.1m深度处室内直剪试验结果进行标定。将室内直剪试验测得黄土宏观力学参数代入表6的反演公式，得到半径乘子λ¯、摩擦系数μ和法向黏结强度σ¯c,其余细观参数取参照值，将细观参数代入模型得到宏观响应的数值模拟结果。将数值模拟结果的强度特征和变形特征与室内试验结果对比，结果如图12和图13所示。

图12剪切强度的结果对比

图13剪切应力-剪切应变的结果对比

由图12可以看出，剪切强度的数值试验与室内试验结果吻合较好。对数值试验的强度结果用下式进行误差分析

由表7所列误差分析结果可知，抗剪强度τ¯c、内摩擦角φ′和黏聚力C与室内试验结果的误差在10%以内，证明了敏感性分析与回归分析相结合的方法在黄土细观参数标定中的可行性。若要进一步降低误差，可参考单因素回归分析中揭露的宏细观参数关系调整计算。例如：内摩擦角φ′和抗剪强度τ¯c的数值模拟结果相比室内试验结果偏小、黏聚力C结果偏大，可通过减小半径乘子λ¯的方式减小黏聚力C,但随着半径乘子λ¯的减小抗剪强度τ¯c也在减小。由于内摩擦角φ′和抗剪强度τ¯c均随摩擦系数μ的增大而增大，此时应优先选择微调摩擦系数μ,但两者随摩擦系数μ的增大速度不同，应适当控制摩擦系数μ的增大量。另外，在细观参数微调过程中，应参考其对宏观参数的影响力序列，按F值从大到小顺序依次调整。

如图13所示为剪切应力-剪切应变的对比结果。由于在预压成样阶段，模型中的颗粒在轴向和径向伺服墙体的共同作用下压密和平衡，成样后的颗粒模型尺寸与室内试验中的尺寸略有不同，因此在对比两者的剪切位移时，采用剪切位移与试样直径的比值作为横坐标，即剪切应变。由图13可以看出，数值试验与室内试验的剪切应力-剪切应变曲线吻合得较好。

根据以上标定过程，对敏感性分析与回归分析相结合的参数标定方法给出经验建议：(1)在正交试验开始之前，应先对各细观参数的取值范围进行试算，以使室内试验所得到的宏观参数值落在正交试验的数值模拟结果范围内。(2)为了保证标定精度，应选择与室内试验结果较接近的一组数据作为参照组。(3)对宏观参数影响较大却又与其没有明确相关性的细观参数，为消除其对参数标定的干扰，宜将其设为固定值。例如本文中，切向黏结强度c¯和摩擦角ϕ¯对内摩擦角φ′和黏聚力C产生了非常显著的影响，但其与宏观响应却没有明确的相关性，在参数标定时采用了参照值。(4)在试验结果对比时往往会出现强度结果吻合得较好，但变形结果差别较大的情况，这与成样阶段颗粒试样的密实程度有关。土的宏观变形主要由颗粒间位置的变化引起，正常固结黏土和松砂的应力-应变曲线一般表现为加工硬化，超固结黏土和密砂的应力-应变曲线一般表现为加工软化。因此，在参数标定之前，应先观察待标定的室内试验的变形曲线特征，若其表现为加工软化，宜在预压成样之前将颗粒间初始摩擦系数和初始黏结设为零或较小的值，使颗粒在预压成样阶段均匀弹开，形成超固结黏土或密砂；若其表现为加工硬化，宜在预压成样之前将颗粒间初始摩擦系数或初始黏结设为较大值，使得颗粒在摩擦或黏结的作用下嵌固在一起不能充分均匀弹开，形成有较大孔隙的松砂或正常固结黏土。在此基础上，再进行下一步的参数标定工作。

4结论

(1)通过敏感性分析可知：

各细观参数对抗剪强度τ¯c和黏聚力C的影响程度和影响力排序较一致，半径乘子λ¯对抗剪强度τ¯c和黏聚力C的影响最强，摩擦角ϕ¯和摩擦系数μ对抗剪强度τ¯c和黏聚力C的影响最弱；各细观参数对内摩擦角φ′均有显著影响，其中摩擦角ϕ¯和摩擦系数μ对内摩擦角φ′的影响最强。

(2)通过单因素回归分析可知：

黄土的抗剪强度τ¯c随半径乘子λ¯的增大呈指数函数规律增大，随摩擦系数μ的增大呈线性增大，随法向黏结强度σ¯c的增大呈幂函数规律增大；内摩擦角φ′随摩擦系数μ的增大呈对数函数规律增大；黏聚力C随半径乘子λ¯的增大呈线性增大，随法向黏结强度σ¯c的增大呈对数函数规律增大。

(3)通过数值试验与室内试验对比可知：

数值试验的强度结果和剪切应力-剪切应变结果与室内试验结果吻合较好，其中强度结果的误差在10%以内，达到了细观参数标定的目的，证明了此方法在黄土细观力学参数标定中的可行性，大大减少了工作量，为其他材料的细观力学参数标定提供了新思路，也为黄土地区滑坡、泥石流和融水降水侵蚀等地质灾害的颗粒流模拟分析提供了参数基础。

水利水电技术（中英文）

水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。