免费下载！
[!--downpath--]概括：
针对黄土颗粒流模型细观参数标定过程中试错法估计量较大的现状，采用敏感性分析和回归分析相结合的方法，分析细观诱因与细观诱因之间的对应关系。黄土的宏观热响应和多诱因非线性。 基于线性联合反演，快速确定了黄土的细观参数，并得出以下推论：黄土的剪切硬度τ¯τ¯c随着直径乘数λ¯λ¯的减小呈指数下降，减小随摩擦系数μ的减小呈线性关系，随法向胶粘剂硬度σ¯σ¯c的减小呈幂函数规律减小； 内摩擦角φ′随摩擦系数μ的减小呈对数函数规律减小内聚力C随直径倍数λ¯λ¯的减小呈线性减小，随法向结合硬度σ¯σ¯的减小而减小c 中的对数函数定律。 数值试验与室外试验结果偏差较小，证明了该方法在黄土细观参数标定中的可行性，为其他材料的介温参数标定提供了新思路。 为侵蚀等地质洪水的颗粒流模拟分析提供参数参考。
关键词：
细观参数标定； 正交检验； 敏感性分析; 粒子流；
关于作者：
董建鹏（1987—），男，硕士研究生，主要从事岩土工程数值模拟研究。
*李辉（1973—），男，中级工程师，博士，主要从事岩土防灾减灾研究。
基金：
湖南省科技厅项目（2020-ZJ-718）；
引用：
董建鹏，李辉． 黄土颗粒流宏微观对应及参数标定方法研究[J]. 水利水电技术（中英文），2022,53(4):180-191. ]. 沃特兰, 2022, 53(4): 180-191。
0 前言
黄土高原水土流失、滑坡、泥石流等自然灾害对路桥沿线人民的生命财产安全造成极大威胁，对黄土热学性质的研究具有重要的现实意义。 粒子流法具有成本低、重复性好等优点，常被用作化学实验的补充和延伸。 在对黄土热特性进行粒子流数值模拟研究时，通常采用试错法，通过不断比较宏观热特性和变形响应来标定模型的细观参数。 由于宏观和微观参数之间的关系复杂，试错过程繁琐费时，国内外学者对此进行了深入的探索和改进。 在参数标定方法的探索中：郭红等。 完善了二维规则排列颗粒的平行键合模型，通过物理推理得到了宏观弹性挠度、拉伸硬度和颗粒键合硬度、颗粒行列数和颗粒直径等定量参数。 关系; 郭万里等。 提出了一种同时测定粗粒材料模拟孔隙率和摩擦系数的方法。 邓树新等。 通过实验设计研究了硬岩细观参数的标定； 王金伟等。 提出了正交等值线法标定砂岩桩的细观参数； 热参数的拟合关系提出了适用于类岩石材料的细观参数标定方法； 陈亚东等。 在总结大量数值试验结果的基础上，提出了与泥岩孔隙度、内摩擦角和压缩挠度相匹配的细观参数。 查看参数的校准技术。 在宏观与细观对应关系研究方面：学者们对岩石材料和粗粒土做了大量研究，一般通过模拟单轴拉伸、压缩试验、巴西劈裂试验来分析细观参数与宏观参数之间的关系, 以及大规模的三轴试验。 根据响应关系，宏观弹性挠度与细观有效挠度呈线性关系，单轴压缩硬度和模量主要受法向和切向结合硬度和摩擦系数控制。 对粘性土的研究相对较少。 通常通过分析三轴试验的挠度-应变曲线来研究峰值硬度、残余硬度、剪切膨胀特性、泊松比等指标与细观参数的关系。 硬度主要与围压、孔隙率、摩擦系数和结合硬度有关； 应变软化性能主要受摩擦系数、结合硬度和围压的影响； 剪切收缩和剪胀性随着摩擦系数和结合硬度的降低而降低。 也有学者研究了颗粒形状和结构对硬度变形特性的影响，认为形状单一会导致数值试验结果偏高，力链分布、裂纹数量、特性应变软化和硬化与结构密切相关。 在黄土研究方面，童晓等。 分析了重力条件、颗粒密度、刚度和模型尺寸对硬度曲线的影响。 江明镜等。 提出了非饱和结构性黄土的三维接触模型，并分析了平面应变条件下的有效范围。 压力与黄土宏观响应的关系。
上述研究在参数标定方法、颗粒流细观参数与宏观硬度和变形响应的关系等方面取得了许多重要成果。 试验的粒子流模型，直剪试验的粒子流模拟较少，直剪试验的粒子流模型为PFC2D。 事实上，只有少数工程问题可以简化为二维平面的挠度和应变问题； 二是目前的研究多集中在岩石、粗粒土等材料上，而对黄土的研究较少，主要集中在模量与应力应变曲线、剪切硬度、内摩擦角、粘聚力等硬度特性的研究较少，黄土细观参数的标定仍然是一种试错法，缺乏更快、更准确的标定方法； 第三，数值模拟实验的分析只采用正交设计，因为正交设计的激励是分散的。 、均匀性等特点，仅从正交数值实验的结果很难得到细观参数与宏观响应之间明确的回归关系，需要在正交实验的基础上辅以其他手段； 四、粒子流的细观参数较多 1.模型构建较为复杂。 以往学者对模型的简化、待标定的细观参数的选取、细观参数的取值范围等问题讨论甚少，给初学者造成很大的困惑。
基于以上讨论，本文构建了直剪试验的PFC3D颗粒流模型。 针对黄土数值模拟中细观参数的标定问题，提出了敏感性分析与回归分析相结合的方法。 基于宏观热响应对应和多原因非线性联合反演，快速确定了黄土的微观参数，为黄土高原地区坝体颗粒流模拟分析提供了参数依据； 同时，模型的简化，细观参数的筛选，探讨了细观参数的取值范围与宏观响应的关系，为其他材料的细观热参数标定提供了新的思路。
1 直剪试验粒子流模拟
1.1 建立模型
颗粒流微观参数的标定采用数值试验与室外试验宏观响应对比的方法。 标定过程如图1所示，即：第一步，在室外通过骨料实验、密度试验、弹性挠度试验等确定颗粒的粒径、密度和弹性挠度； 第二步，通过假设和理论分析，减少需要标定的细观参数数量，增加标定难度； 第三步，设计正交试验，利用试验结果对细观参数的宏观响应进行敏感性分析，找出起控制作用的细观参数； step 4，固定其他参数不变，依次改变将控制参数引入模型估计，总结单激励回归公式和反演公式，分析其与宏观响应的关系； step 5，将自变量和因变量的位置颠倒，将宏观响应作为自变量，将控制参数作为因变量物理二力平衡公式，在参考单激励反演公式特点的基础上，多激励非线性联合归纳反演公式； 步骤6，将室外试验结果带入联合反演公式得到相应的控制参数，其余细节观测参数根据参考值带入模型估计数值模拟结果； 第七步，分析数值模拟结果与实验室测试结果的偏差，并参考单激励回归关系对结果进行微调。
图1 细观参数标定过程
PB模型是根据He等人提出的模型构建的，由非主动减震器和参考间隙为零的平行键组成（见图2）。 不同于线接触模型的点接触，PB模型可以在两个相互胶结的接触面（片）之间传递力和扭矩，胶结长度由直径倍数λ¯控制，可以模拟不规则黄土颗粒的晶间相互位移产生的剪切力。 由于法向结合硬度σ¯c的存在，PB模型中相互结合的接触面可以承受拉力，因此当接触面的间隙gs小于零时，接触仍然存在。 该特性可以模拟范德华力、分子力和毛细管力等产生的胶结效应。PB模型中的线性分量和内聚力分量是平行且相互独立的。 当拉力超过法向结合硬度 σ¯c 或剪切力超过切向结合硬度 c¯ 时，结合断裂，PB 模型退化为线性接触模型（见图 3）。
图2PB模型
图3 PB模型的退化模型
PB模型中有很多细观参数。 为了降低参数标定的难度，通过假设和理论分析减少了需要标定的细观参数数量，适当简化了估计模型。 PFC估算程序中，估算时间步长Δt()的计算方法如下
式中，M为粒子总质量； K 是粒子的总偏转； k 是粒子的偏转; ρs 是颗粒的密度； r为颗粒直径； β 为调整系数。
估计的时间步长 Δt 表示每次更新粒子的力和位移状态之间的机器时间。 程序每前进一个时间步Δt，所有粒子的力和位移状态都要更新一次。 粒子越多，每个时间步长Δt实际消耗的时间就越长。 由式（1）可知，时间步长Δt随着粒径r和粒子密度ρs的减小而减小，随着粒子偏转k的减小而减小。 粒径比是指最大粒径与最小粒径的比值。 当粒径比大时，小颗粒在大颗粒之间的缝隙中移动，产生“悬浮”颗粒（见图4）。 “悬浮”颗粒的存在使得模型的不平衡力在预压阶段难以达到稳定值，对模型的宏观热性能影响不大[20]。 因此，应减小粒子直径r、增大粒子密度ρs、减小粒子偏转k、减少粒子数、减少“悬浮”粒子数等，在保证的前提下提高估计效率结果的准确性。
图4 大颗粒之间的漂浮颗粒
黄土主要由构成基岩骨架的粗粉砂颗粒组成，粒径为0.01～0.05mm。 研究表明，当模型的特征宽度比L/R足够大时（L/R>40），L/R对粗粒土的剪切硬度没有影响。 经过多次试算，选择粒子大小放大50倍，生成的粒子在直径0.0005~0.之间概率分布均匀。 此时粒子数量多，计算顺利，估计结果比粒子尺寸放大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以减小颗粒的惯性作用，外层径向伺服wall用于保证围压误差在1%以内。
图5 直剪试验数值模拟
1.2 正交实验设计
正交实验设计()是一种研究多激励、多水平的实验设计方法。 它根据实验的正交性从综合实验中选择一些有代表性的点进行实验。 这个代表点具有“分散均匀、整齐可比”的特点，大大减少了测试次数。 是一种高效、快速、经济的测试设计方法。 本次正交试验涉及6个诱因，每个诱因设置5个水平。 正交试验设计如表2所示，正交矩阵序列和数值模拟结果如表3所示。其中，剪切硬度τ¯τ¯ c采用法向挠度为时的值。
2 宏观参数与微观参数对应关系分析
2.1 敏感性分析
利用SPSS()统计分析软件对各诱因的主效应进行多诱因残差分析，得到F统计量和Sig。 相对随机率的值。 F统计量是指原假设成立时服从F分布的统计量。 F统计量等于链表之间的差异与组内差异的比率。 当比值超过一定限值时，认为链表之间存在显着差异。 ; 相位随机率 Sig。 值是指原假设的概率。 F统计量越大，相位随机率Sig越小。 价值。 根据 F 统计量和 Sig 的大小。 相位随机率的值，可以确定各个细观参数对宏观参数的影响程度。 如果是Sig.0.05，影响不明显。 分析结果如图 6 所示。
图 6 敏感性分析
从图 6（a）可以看出，直径乘数 λ¯、切向结合硬度 c¯、法向结合硬度 σ¯c 和概率相关值 Sig。 有效挠度E*均大于0.05，对剪切硬度τ¯c形成显着影响。 从 F 统计量来看，摩擦系数 μ 和摩擦角 ϕ¯ 对剪切硬度 τ¯c 的影响微弱。 根据对剪切硬度τ¯c的影响程度排序：直径乘数λ¯>法向结合硬度σ¯c>切向结合硬度c¯>有效挠度E*>摩擦系数μ>摩擦角φ¯。
从图 6（b）可以看出，概率相关值 Sig。 若大于0.05，则对内摩擦角φ′有显着影响； 如果摩擦系数μ、摩擦角φ¯、切向结合硬度c¯、半径乘数λ¯都大于0.01，则它们对φ'的影响非常明显。 按照对内摩擦角φ′的影响程度排序：摩擦系数μ>摩擦角φ¯>切向结合硬度c¯>直径倍数λ¯>有效挠度E*>法向结合硬度σ¯c。
从图6(c)可以看出，直径乘数λ¯、切向结合硬度c¯、法向结合硬度σ¯c和概率值Sig。 的有效挠度E*均大于0.05。 从F统计量来看，摩擦系数μ和摩擦角φ¯的影响程度较弱，摩擦系数μ对凝聚力C几乎没有影响。按对凝聚力C的影响程度排序：直径乘数λ¯ > 法向结合硬度 σ¯c> 切向结合硬度 c¯> 有效挠度 E*> 摩擦角 ϕ¯> 摩擦系数 μ。
2.2 单一激励回归分析
为进一步获取宏观参数与其控制激励之间的定量关系，以正交试验No.17的结果作为参照组，仅取Sig对应的细观参数。
以直径乘数λ¯为例，保持其他细观参数不变，将直径乘数λ¯从0.5增加到1.5，分别代入模型估计。 为了得到更准确的结果，直径乘数λ¯的水平数减少为9，直径乘数λ¯单激励回归分析检验结果如表4所示。直径乘数λ¯与宏观响应如图 7 所示。
图7 直径倍数λ¯与宏观响应的对应关系
从图7(a)可以看出，随着直径倍数λ的减小，剪切硬度τ¯c大致呈指数函数减小，并按指数函数进行拟合，得到拟合关系作为 τ¯ c=50.837e0.1639λ¯，R2=0.9801τ¯c=50.837e0.1639λ¯，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.μ+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σ¯c0.7703，R2=0.9897τ¯c= 252.48σ¯c0.7703, 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的细观参数。
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) For the that have a great on the but have no clear with them, in order to their to the , it is to set them as fixed . For , in this paper, the bond c¯ and angle ϕ¯ have a very on the angle φ′ and C, but they have no clear with the . the value. (4) When the test are , the are often in good , but the are quite , which is to the of the in the stage. The of soil is by the of the . The - of soil clay and loose sand show work , and the - of super-soil clay and dense sand show work . , , the curve of the test to be be first. If it shows , it is to set the and to zero or a small value pre- into . value, so that the off the pre- stage to super-soil clay or dense sand; if it shows work , it is to set the or to a value pre- value, in the under the of or , the be fully and , in loose sand or soil clay with large pores. On this basis, the next work is out.
4。结论
(1) , it can be known that:
The and order of of each on the shear τ¯c and C are , the λ¯ has the on the shear τ¯c and C, and the angle ϕ¯ and μ have the on shear τ¯c and C; each has on angle φ′, among which angle ϕ¯ and μ have The angle φ' has the .
(2) - , it can be known that:
The shear τ¯c of loess with the of λ¯, with the of μ, and with power with the of bond σ¯c The angle φ′ as a with the of the μ; the C with the of the λ¯, and with the bond σ¯c The is a law.
(3) the the test and the test, it can be known that:
The and shear -shear of the test are in good with the test , and the of the is 10%. The of the , new ideas for the of of other , and also for the flow of such as in loess areas, flows, and Base.
Water and ( and )
"Water and ( and )" of the of Water is a () of China's water and . This the in the , , , , and of water in my , as well as the , , , and of water and . 技术。 The main of the are: and water , , , , , , , water and , , , , metal , water , Water , flood and , , new , urban water , rural water , soil and water , , water , water , etc.