离散元模拟中三维试验孔隙率向二维模拟孔隙率转换方法

技术领域

本发明涉及离散元数值仿真模拟技术领域，尤其涉及一种三维试验孔隙率向二维模拟孔隙率的转换方法。

背景技术

离散元法因其在解决非连续介质问题上的独特优势，被广泛应用于颗粒材料的细观行为研究，并且被认为是对传统模型试验和有限元数值模拟的重要补充。一般认为，对于离散元数值模拟，其模拟结果的可靠性很大程度取决于模型的选择及模型参数的选取，如颗粒与墙体的接触刚度和摩擦系数、模量、阻尼和泊松比等。而实际上，对上述参数影响的讨论只有针对某一已确定孔隙率的颗粒堆积体时，才是有意义的。

相比于三维模型，二维离散元模型在计算效率上有着巨大优势，能够以较小的计算成本完成对研究对象机理性问题的研究。在三维离散元法研究中，颗粒堆积体孔隙率可以直接由实测孔隙率确定，而二维状态下的孔隙率则需要以实测孔隙率为依据，通过适当的方法转换为二维孔隙率。因为基于面积计算的二维孔隙率与基于体积计算的三维孔隙率存在着极大差异。至今为止，还没有一种三维试验孔隙率与二维模拟孔隙率的转换方法能够广泛应用于处理土力学问题。

为了解决上述问题，我们通过对比分析现有的一些转换方法，提出了一种新的三维试验孔隙率向二维模拟孔隙率转化的方法。

发明内容

本发明的目的是为解决背景技术中提到的问题，为此提出了一种新的三维试验孔隙率向二维模拟孔隙率转换的方法。

为解决上述技术问题，本发明采用的技术方案是：首先通过对比分析已有方法，提出了计算三维试验孔隙率所对应二维模拟孔隙率初值的抛物线方程，随后在已经得到的二维模拟孔隙率初值基础上设计了一套流程，以迭代的方式得到合适的二维模拟孔隙率终值，具体内容如下：

1提出抛物线方程，确定二维模拟孔隙率初值

目前，离散元研究中有六种常见的用于联系三维试验孔隙率与二维模拟孔隙率的方法：1.1现有方法

1.1.1最密堆积法

假设颗粒在二维和三维种具有相同的直径，通过匹配二维六边形晶格结构和三维面心立方体结构得出公式(1)，由于两者分别代表了二维和三维状态下的最密堆积形式，因此在这里被命名为最密堆积法。

式中：n

1.1.2最松散堆积法

与最密堆积法的推到过程类似，最松散堆积法通过匹配二维正方形结构和三维简单立方体结构得出公式(2)。

1.1.3区间映射法

区间映射法如式(3)所示，通过比较二维六边形晶格结构和三维六边形填充结构可以对该方法进行验证，两种结构都基于单一颗粒尺寸系统。

1.1.4组合方法

该方法结合了最密堆积法与区间映射法，并引入了相对密度D

式中：

当D

1.1.5线性插值法

土力学中的密度定义也可以用于三维孔隙率向二维孔隙率的转换。线性插值法基于这样的假设：二维状态与三维状态(试验状态)的密度是相同的，通过式(5)得到二维孔隙率。

式中：n

对式(5)进行变换，可得式(6)：

n

式中：

1.1.6准三维孔隙率法

该方法假设二维圆盘的直径和厚度等于球体的直径，通过比较圆盘和球体的体积得到计算准三维最大和最小孔隙率的方程。

/>

随后，使用线性插值法计算准三维孔隙率：

最终，二维孔隙率可由式(10)求得：

该方法可以看作是最松散堆积法和线性插值法的结合。

1.2现有方法的理论评估

通过对比基于不同方法得到的二维模拟孔隙率和三维试验孔隙率之间的关系，对上述方法进行了理论评估。实验室测得的砂土最大孔隙率和最小孔隙率分别为0.4571和0.3253。基于组合法得到的二维孔隙率随着试验孔隙率的增加而降低，这与二维孔隙率和三维试验孔隙率之间的一般规律相反；基于最松散状态法和区间映射法得到的二维孔隙率存在负值，不能用于处理土力学问题；最密堆积法是基于单一尺寸系统的假设提出的，因此无法广泛应用于多尺寸颗粒系统；准三维孔隙率法只适用于二维圆盘的直径和厚度等于球体直径的情况；线性插值法高度依赖于通过试验得到的最大孔隙率和最小孔隙率，然而，理论上线性插值法应当将二维正方形结构和二维六边形晶格结构与三维简单立方结构和三维面心立方结构匹配。

当线性插值法采用上述理论假设，二维孔隙率的计算将与实验测得的最大孔隙率与最小孔隙率无关，对应于试验孔隙率n

1.3抛物线方程

尽管采用1.2中理论假设的线性插值法也是基于单一尺寸系统，但该方法对应曲线的两端点仍可最为确定抛物线方程的粗略准则，此外，多尺寸颗粒系统的三维孔隙率为0时，对应的二维孔隙率也为0。据此，可以得到初始的抛物线方程，见式(11)。需要注意的是，抛物线方程与实验室得到的最大孔隙率和最小孔隙率无关，它只是将试验孔隙率转化为二维孔隙率的一个初始准则。

n

2确定二维模拟孔隙率终值

步骤1：将试验孔隙率代入抛物线方程，得到二维模拟孔隙率初值；

步骤2：建立简单的二维离散元数值模型，生成四个刚性墙体使其围成一个矩形区域，随后在其内生成孔隙率与所选二维模拟孔隙率一致的颗粒系统；

步骤3：运行至模型稳定后，测量顶板上墙体所受竖向力，如果该竖向力大于0且小于试样重量的1％，认为使用该孔隙率可以使颗粒充满试样空间。同时，观察重力作用下试样内的力链分布是否符合实际情况。如果满足以上条件，认为该孔隙率是合理的，跳出程序，使用该孔隙率作为二维模拟孔隙率终值。如果不满足，则下降或上升顶板，直至顶板上所受竖向力大于0且小于试样重量的1％，计算得到此时试样孔隙率。

步骤4：计算下降或上升顶板后试样的孔隙率与原孔隙率之间差值的绝对值是否小于0.005。如果绝对值小于0.005，认为此时调整前后孔隙率的差异对模拟结果影响较小，选择原孔隙率作为最终二维孔隙率；如果绝对值大于0.005，则需要对孔隙率进行调整。当步骤2.2中对顶板进行了下降操作时，将孔隙率较原孔隙率调小0.01，当步骤2.2中对顶板进行了上升操作时，将孔隙率较原孔隙率调大0.01。

步骤5：以步骤3中调整后的孔隙率为基础，重复步骤2与步骤3，直至得到二维模拟孔隙率终值。

与现有技术相比，本发明提出的三维试验孔隙率向二维模拟孔隙率转换方法的优势在于：

本发明公开了一种将三维试验孔隙率转换为二维模拟孔隙率的方法，主要内容包括：提出了计算二维模拟孔隙率初值的抛物线方程；设计了一套基于二维模拟孔隙率初值得到二维模拟孔隙率终值的程序。

附图说明

图1为现有方法之间三维试验孔隙率与二维模拟孔隙率关系的对比图；

图2为本发明提出的抛物线方程法与线性插值法在三维试验孔隙率与二维模拟孔隙率关系上的对比图；

图3为本发明设计的基于二维模拟孔隙率初值得到二维模拟孔隙率终值的程序流程图；

图4为实施例中所涉及颗粒系统的颗粒级配图；

图5为实施例中各阶段试样的力链分布图；

上述附图中，各标号说明如下：

线性插值法

F

具体实施方式

以下实施例仅出于说明性目的，而不是想要限制本发明的范围。

实施例

参照发明内容中的第二部分内容，具体步骤如下：

步骤1：将三维试验孔隙率代入抛物线方程，得到二维模拟孔隙率初值；

本实施例中三维试验孔隙率为0.3434，代入抛物线方程式(11)得到二维模拟孔隙率初值为0.14。

步骤2：建立简单的二维离散元数值模型，生成四个刚性墙体使其围成一个矩形区域，随后在其内生成孔隙率与所选二维孔隙率一致的颗粒系统；

本实施例中二维离散元数值模型所涉及参数如表1所示：

表1离散元模拟参数

需要注意的是，本实施例为了提高计算效率，对颗粒粒径进行了调整和放大，具体情况见图4。

步骤3：二维模拟孔隙率为0.14时，模型稳定后顶板和底板上所受竖向力非常大，分别为3854kN/m和3855kN/m，远大于试样自重的846N/m，此时的试样力链分布如图5(a)所示。下一步需要对顶板进行上升操作，直至其所受竖向力大于0且小于试样重量的1％，完成上升操作后的试样力链分布如图5(b)所示。可以明显看出此时的力链分布与实际情况更加吻合，计算此时的试样孔隙率。

步骤4：计算移动顶板后试样的孔隙率与原孔隙率之间差值的绝对值为0.008，大于0.005，需要使用调整后的二维孔隙率进行重新模拟。调整后的孔隙率为0.15。

步骤5：重复步骤3和步骤4。在进行步骤3时发现此时的试样力链分布如图5(c)所示，顶板所受竖向力为2.7N/m，符合标准。底板上测得的竖向力为809N/m，与此时的试样自重836N/m较为接近。认为此时的孔隙率是合理的，可以基于此完成进一步的离散元模拟研究。

以上所述均为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。