¹成都理工大学环境与土木工程学院，四川 成都

²西南交通大学土木工程学院，四川 成都

收稿日期：2021年4月28日；录用日期：2021年9月23日；发布日期：2021年9月30日

摘要

西部丘陵区采用开挖泥质岩土石料填筑高路堤、高填方机场跑道、高土石坝和高填方场地等填筑体情况越来越多；泥质岩土石料碾压填筑体的二元结构及其变形特性与填筑体长期变形趋势和稳定性相关，目前对此问题的研究还不够深入系统。本文利用理论分析与数值模拟方法，在采用赫兹接触模型建立的泥质岩二元结构及简化球形接触模型基础上，就分层碾压填筑体变形特征、机制和影响因素进行了系统研究，结果表明：所建理论模型宏观变形总趋势与相似数值分析模型的变形及实际填筑体监测获得的变形总趋势基本一致；二元结构填筑体沉降变形与土石料颗粒平均粒径、单层细观结构和二元结构颗粒粒径组合等因素有关；填筑体平均粒径减小其宏观总变形量增大；二元结构层间的平均粒径差异变大则结构界面附近的局部变形差异明显，具有变形突变特征；各单层颗粒排列方式对填筑体宏观总变形也有影响，其中孔隙率较大的等效球形矩形排列填筑体的变形大于孔隙率较小的等效球形三角形排列填筑体的变形。该文研究结果为进一步研究泥质岩土石料填筑体细、宏观二元结构特征及其变形机制和工程应用参考。

关键词

泥质岩土石料，碾压填筑体变形，二元结构模型，赫兹接触模型理论，颗粒大小，PFC^2D

Study on the Deformation Characteristics and Mechanism of Argillaceous Earth-Rock Compacted Fillings

Zhibin Zhong¹, Ronggui Deng², Yi Sun², Yuanyuan Wang², Ling Dong²

¹College of Environment and Civil Engineering, Chengdu University of Technology, Chengdu Sichuan

²School of Civil Engineering, Southwest Jiaotong University, Chengdu Sichuan

Received: Apr. 28^th, 2021; accepted: Sep. 23^rd, 2021; published: Sep. 30^th, 2021

ABSTRACT

In western hilly region, increasing numbers of excavated argillaceous earth-rock materials were used to fill compacted bodies, such as high embankment, high filled airport runway, high earth-rock dam and high filled construction sites, whose dualistic structure and deformation characteristics were related to their long-term deformation trend and stability. There were hardly any researches on this issue recently. Based on the established argillaceous dualistic structure and simplified spherical contact model, the de-formation characteristics, mechanism and influencing factors of layered compacted fillings were systematically studied by means of Hertz contact theory and compared with discrete element numerical model. The results showed that the macroscopic deformation trend derived from established theoretical model was basically similar to that of numerical model and monitoring of actual compacted filling. The settlement of fillings was related to such factors as the average particle size of earth-rock materials, single-layer micro-structure and particle size combination model of dualistic structure. The total macroscopic deformation increased with the decreasing of average particle size. The local deformation difference was more significant when the average particle size difference between two layers be-came larger. As well, the particle arrangement of each monolayer had influence on the overall macro-scopic deformation of fillings, among which the deformation of rectangular arrangement of equivalent spherical with larger porosity was larger than that of triangular arrangement with lower porosity. The results can provide references for further study on the micro and macro deformation characteristics and mechanism of argillaceous earth-rock fillings with dualistic structure as well as engineering application.

Keywords:Argillaceous Earth-Rock Material, Deformation of Compacted Filling, Dualistic Structure Model, Hertz Contact Theory, Particle Size, PFC^2D

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

红层在中国西南、西北及中部地区广泛分布，其工程地质特性极其复杂 [1]。红层多为以软弱泥质岩夹泥质砂岩的软硬互层结构 [2]。

四川盆地丘陵区红层泥质岩自然开挖土石广泛用于机场跑道及高速公路路堤填筑料，经分层填筑碾压后成为一种土石混合体，其变形特性不同于一般的填土路堤和填石路堤 [3]。泥质岩开挖土石料的变形特性一方面和泥质岩本身的性质有关，另一方面也与填筑的碾压施工工艺有关。其中，由于碾压工艺特性及作用机理，在分层碾压中常呈现出每填筑碾压单层表部密实性大于底部，有时泥质砂或砂质泥岩块体较大，某填筑单层底部还会出现孔隙较大“架空”现象，使某一填筑单层具有如图1所示的上密下疏的“二元结构”特征。

泥质岩开挖土石料的变形特性对于工程建设稳定性至关重要，目前国内外针对红层软岩土石料展开了广泛的研究，主要集中于工程地质特性、物理力学特性、水理特性等方面。周应华 [4] 等利用MTS815Teststar程控伺服岩石力学试验探究了川东地区某红层边坡中的砂岩、粉砂岩和泥岩的弹性模量与围压的关系。章清叙 [5] 等对周期荷载作用下红砂岩的三轴剪切变形特性进行了较系统的试验研究，得到三轴周期荷载的上限应力和幅值对疲劳破坏特性有显著影响的结论。蒋关鲁 [2] 等通过动三轴试验，分析了红层泥岩土石混合体作为路基基床填料，在低围压条件下的累积变形特性。郭永春 [6] [7]、J. Hadizadeh [8]、ME Chenever [9] 等分析了红层软岩水岩作用的物理化学效应。宋磊 [10]、章翠英 [11] 等对软岩的水岩作用下的微细观机理进行了研究。邓华锋、周美玲 [12] [13] [14] 等研究了水岩作用下红层软岩的蠕变特性、劣化规律及机理。

纵观国内外对土石料的细观结构研究，所采用的细观模型多为随机结构模型，几乎没有针对泥质岩开挖土石料分层碾压二元结构填筑体变形特性及其机理进行的研究。简化的细观结构模型对正确认识其力学特性具有重要作用；可从细观结构上认识土石混合料填筑体的变形力学机理。

图1. 二元结构特征

针对泥质岩开挖料分层碾压填筑体的变形特性，本文基于二元结构特点的泥质岩开挖土石料分层填筑体进行二维建模，以成都某在建国际机场跑道泥质岩土石料填筑体为对象，生成二维细观结构试样模型，进行变形机理分析，并用二维颗粒流离散元分析方法对土石料填筑体的细观变形力学破裂机理进行对比研究。

2. 土石料填筑体变形机理分析

2.1. 等效圆形颗粒体结构模型分析

普通土岩填筑体的松密程度和颗粒常为渐变的。为了从理论上深入探讨填土变形特性及其机理，将土岩填筑体颗粒上下叠合的三个颗粒简化成等效圆柱截面平面应变问题的颗粒弹性体模型。此时可以借助Hertz颗粒群体接触解，从细观上考虑两颗粒体之间的接触变形关系(见图2)，坐标y垂直纸面向内。根据Hertz理论的要求对填筑体做如下简化：1) 表面都是连续并且是非协调的，小应变范围， $a_i\ll R_i$ ( $i=1,2,3,\cdots ,n$ )；2) 每个球形颗粒体可看作是一个弹性半空间， $a_i\ll R_{i-1},R_i$ ( $i=1,2,3,\cdots ,n$ )；3) 表面摩擦忽略，即接触法面内 $q_x=q_y=0$ ；4) 剪切变形对竖向变形的影响忽略。

根据Hertz理论 [15]，作用在圆柱径向上的压缩线荷载 $p_i$ 在O₁点引起的Hertz面分布压力为：

${\sigma }_i=\frac{2p_i}{\pi a_i}{\left(1-\frac{x^2}{a_i^2}\right)}^{\frac{1}{2}}$ (1)

式中，a_i为接触半宽度， $i=1,2,\cdots ,n$ (n由填土颗粒大小及总厚度确定)；由式(2)给出：

$a_i^2=\frac{4p_i{R}_i}{\pi E_i^{*}}$ (2)

式中， $E_i^{*}$ 为两接触球体的综合模量(平面应变)，由式(3)计算，

$\frac{1-{\mu }_i^2}{E_i^{*}}=\frac{1-{\mu }_{i-1}^2}{E_{i-1}}+\frac{1-{\mu }_i^2}{E_i}$ (3)

式中其余参数见图2。

图2. 圆柱物体的截面线接触体系

同时，为了满足圆柱的圆形边界应力为零，Timoshenko与Goodier给出了x向基本均匀双轴拉应力。

${\sigma }_{xi}=-\frac{p_i}{\pi R_i}$ (4)

圆形截面任意一点A (见图2)的应力由以下三部分组成：1) 在O₁点引起的接触应力，可由式(1)计算；2) 在O₂点引起的接触应力，因O₂点离A点的距离较大，可以看作是由集中力p_i+1所引起的应力；3) 由公式(4)给出的双轴拉力。三者叠加得到A点的应力为，

${\sigma }_{xi}=-\frac{p_i}{\pi }\left\{\frac{1}{R_i}-\frac{2\left(a_i^2+2a_{i+1}^2\right)}{a_i^2{\left(a_i^2+z^2\right)}^{1/2}}+\frac{4z}{R_i^2}\right\}$ (5)

${\sigma }_{zi}=-\frac{p_i}{\pi }\left\{\frac{1}{R_i}-\frac{2}{2R_i-z}-\frac{2}{{\left(a_i^2+z^2\right)}^{1/2}}\right\}$ (6)

其中，

$p_i=\left({\displaystyle \underset{k=0}{\overset{i}{\sum }}2\gamma R_k}+{\sigma }_{i0}\right)\times 2R_i$ (6-1)

${\sigma }_{i0}=\frac{3p{\left(2R_1\right)}^2{\left({\displaystyle \underset{k=1}{\overset{i}{\sum }}2R_k}\right)}^3}{2\pi {\left({\left(2{\displaystyle \underset{k=1}{\overset{i}{\sum }}{R}_k}\right)}^2+x^2\right)}^{5/2}}$ (7)

式中，R_i为自上而下第i层颗粒的等效半径，p为土岩填筑体表面等效均布荷载， $\gamma$ 为土岩填筑体平均重度。对于土岩“半无限”填筑体，认为基本符合弹性平面应变条件，这样可得到z方向的应变为，

${\epsilon }_{zi}=\frac{1-{\mu }_i^2}{E_i}\times \left\{{\sigma }_{zi}-\frac{{\mu }_i}{1-{\mu }_i}{\sigma }_{xi}\right\}$ (8)

z从0至R_i，对 ${\epsilon }_{zi}$ 积分得到圆柱体上半部O₁C的压缩变形量，因 $a_i\ll R_i$ 忽略其微量项，可以得到z方向的位移量，

${\delta }_i=p_i\frac{1-{\mu }_i^2}{\pi E_i}\left\{\left[2\ln \left(\frac{4R_i}{a_i}\right)-1\right]+\frac{{\mu }_i}{1-{\mu }_i}\left[3+\frac{2\left(a_i^2+2a_{i+1}^2\right)}{a{}_i{}^2}\ln \left(\frac{a_i}{2R_i}\right)\right]\right\}$ (9)

对于圆形颗粒体下半部的压缩量，可通过类似分析得到，

${\delta }_{i+1}=p_{i+1}\frac{1-{\mu }_i^2}{\pi E_i}\left\{\left[2\ln \left(\frac{4R_i}{a_{i+1}}\right)-1\right]+\frac{{\mu }_i}{1-{\mu }_i}\left[3+\frac{2\left(a_{i+1}^2+2a_i^2\right)}{a{}_{i+1}{}^2}\ln \left(\frac{a_{i+1}}{2R_i}\right)\right]\right\}$ (10)

上、下球体接触区域产生的压缩变形为，

${\delta }_{ii+1}=2R_i-\left({\left(R_i^2-a_i^2\right)}^{1/2}+{\left(R_i^2-a_{i+1}^2\right)}^{1/2}\right)$ (11)

圆形接触体直径上的总压缩量为，

$\delta ={\displaystyle \underset{1}{\overset{n}{\sum }}\left({\delta }_i+{\delta }_{i+1}+{\delta }_{ii+1}\right)}$ (12)

等径简化情况， $R_1=R_i=R_n=R$ ， $n=\mathrm{int}\left(\frac{h}{2R}+0.5\right)$ ，h为填筑体厚度。

例如，成都某在建国际机场场道某段土岩料填筑体总厚度平均h = 5.0 m，颗粒等效半径R = 0.025 m，n为100。填筑体综合弹性模量E = 45 MPa，泊松比μ = 0.30，填筑体重度γ = 0.02 MN/m³。土面区地面均布荷载为0；飞行区道床层 + 水稳层 + 结构层总厚度1.18 m均布净荷载约为5.9 kPa，最大的A380飞机引起跑道地基面等效均布为38.34 kPa，综合为44.24 kPa。根据上述方法计算得到，土面区和道槽区的工后沉降量分别为34.930 mm和35.111 mm，飞机等附加荷载对填土的沉降变形影响不大，主要是填土自重作用引起的变形沉降。

2.2. 矩形排列等效圆形颗粒体结构模型分析

对于空间集合在一起的土岩料填筑体颗粒，可以简化为等效等径球形颗粒群。在此情况下颗粒之间除上、下球体排列接触外，左、右及前、后也有类似排列。根据几何原理，颗粒矩形(体)排列时孔隙率最大，三角形(或空间四面体)排列时孔隙率最小。对于等效等径颗粒平面问题，其矩形排列情况如图3所示。

由于前、后及左、右均有颗粒排列接触约束，也可以用上述的平面型圆柱Hertz接触理论进行分析。在图3中，颗粒聚合体中取微小单元进行分析，其高度为h，土岩颗粒的弹性模量为E，泊松比为μ，土岩颗粒重度为 $\gamma$ ，颗粒受到竖向外荷载为p_i，分析单元上覆土岩厚度为h_i。

在线荷载p_i作用下(见图4)，土岩填土颗粒群的竖向压缩变形量，仍可由式(4)至式(11)计算得到。

图3. 等效球形矩形排列

图4. 矩形排列模型受力分析

其中， $a_i^2=\frac{4p_i{R}_i}{\pi E_i^{*}}$ ，当 $R_i=R$ 时， $E_i^{*}=E_i/2$ ，则式(5)和式(6)化简为：

${\sigma }_{xi}=-\frac{p_i}{\pi }\left\{\frac{1}{R_i}-\frac{2\left(a_i^2+2a_{i+1}^2\right)}{a_i^2{\left(a_i^2+z^2\right)}^{1/2}}+\frac{4z}{R_{{}^i}^2}\right\}$ (13)

${\sigma }_{zi}=-\frac{p_i}{\pi }\left\{\frac{1}{R_i}-\frac{2}{2R_i-z}-\frac{2}{{\left(a_i^2+z^2\right)}^{1/2}}\right\}$ (14)

$p_i$ 作用下的压缩量为：

${\delta }_i=p_i\frac{1-{\mu }_i^2}{\pi E_i}\left\{2\ln \left(\frac{4R_i}{a_i}\right)-1\right\}$ (15)

同样， $p_{i+1}$ 作用下的压缩量为：

${\delta }_{i+1}=p_{i+1}\frac{1-{\mu }_i^2}{\pi E_i}\left\{2\ln \left(\frac{4R_i}{a_{i+1}}\right)-1\right\}$ (16)

${\delta }_{i2}=p_{i2}\frac{1-{\mu }_i^2}{\pi E_i}\left\{2\ln \left(\frac{4R_i}{a_{i2}}\right)-1\right\}+2R_i-2{\left(R_i^2-a_{i2}^2\right)}^{1/2}$ (17)

其中， $a_{i2}=\frac{4p_{i2}{R}_i}{\pi E_i^{*}}$ ， $p_{i2}=\frac{{\mu }_i}{1-{\mu }_i}\left(p_i+p_{i+1}\right)$ ，故矩形排列等效球形颗粒填土层总变形为，

$\delta ={\displaystyle \underset{1}{\overset{n}{\sum }}{\delta }_{hi}={\displaystyle \underset{1}{\overset{n}{\sum }}\left({\delta }_i+{\delta }_{i+1}+\mu {\delta }_{i2}+{\delta }_{ii+1}\right)}}$ (18)

利用2.1节的实例数据进行计算，土面区和飞行区沉降变形量分别为30.318 mm和30.678 mm。显然，比上述模型沉降值稍小，说明侧向约束会使得竖向变形沉降难度加大，竖向变形沉降量变小。

2.3. 三角形排列等效圆形粒体结构模型分析

等径球形颗粒在平面三角形排列情况下，某个圆形颗粒周边可排六个圆心连线夹角为60˚的等角度圆形面接触点，见图5和图6。接触力分别为 $p_{i1}$ 、 $p_{i3}$ 、 $p_{i2}$ 和 $p_i{}_{+\text{11}}$ 及 $p_i{}_{+\text{13}}$ 。因六个接触点作用的对称性和平衡条件有，

$p_{i1}=p_{i3}$ , $p_{i1}=\frac{\sqrt{3}}{2}\frac{p_i}{2}$ (19)

$p_{i+11}=p_{i+13}$ , $p_{i+11}=\frac{\sqrt{3}}{2}\frac{p_{i+1}}{2}$ (20)

$p_{i2}=\frac{{\mu }_i}{1-{\mu }_i}\left(p_i+p_{i+1}\right)/2$ (21)

土岩料填筑体的压缩变形是因第i层球形颗粒( $p_{i1}$ - $p_{i+11}$ 或 $p_{i3}$ - $p_{i+13}$ )方向接触点变形及其之间土岩颗粒压缩变形在竖直方向的累积。类似于前面的分析， $p_{i1}$ 、 $p_{i3}$ 和 $p_{i2}$ 方向的变形分别为，

${\delta }_{i1}={\delta }_{i3}=\frac{1-{\mu }_i^2}{\pi E_i}\left\{\left[2p_{i1}\ln \left(\frac{4R_i}{a_i}\right)-1\right]+2p_{i+11}\ln \left(\frac{4R_i}{a_{i+1}}\right)-1\right\}+2R_i{\left(R_i^2-a_i^2\right)}^{1/2}-{\left(R_i^2-a_{i+1}^2\right)}^{1/2}$ (22)

${\delta }_{i2}=\frac{2\left(1-{\mu }_i^2\right)p_{i2}}{\pi E_i}\left\{2\ln \left(\frac{4R_i}{a_{i2}}\right)-1\right\}+2R_i-2{\left(R_i^2-a_{i2}^2\right)}^{1/2}$ (23)

其中，

$a_i={\left(\frac{4\sqrt{3}{p}_{i+1}{R}_i}{2\pi E_i^{*}}\right)}^{1/2}$ , $a_{i+1}={\left(\frac{4\sqrt{3}{p}_{i+1}{R}_i}{2\sqrt{3}\pi E_i^{*}}\right)}^{1/2}$ , $a_{i2}={\left(\frac{4p_{i2}{R}_i}{\pi E_i^{*}}\right)}^{1/2}$ (23-1)

图5. 等效球形三角形排列

图6. 三角排列模型受力分析

则总的沉降变形为，

$\delta ={\displaystyle \underset{1}{\overset{n}{\sum }}{\delta }_{hi}={\displaystyle \underset{1}{\overset{n}{\sum }}\left(\frac{\sqrt{3}}{2}{\delta }_{i1}+{\mu }_i{\delta }_{i2}\right)}}$ (24)

同样利用2.1节的实例数据进行计算，土面区和飞行区的沉降变形量分别为9.399 mm和9.511 mm。显然，比以上两种模型的沉降变形量都小，说明三角形排列颗粒结构侧向变形约束比矩形排列颗粒结构侧向变形约束强，使得竖向变形沉降难度较大，产生的竖向变形沉降量较小。

3. 土石料二元结构填筑体变形机理分析

上述介绍了土岩填筑体松密度和粒径渐变情况下的模型与机理分析。如前所述，类似机场等大型土石料碾压填筑体，因碾压静、动应力的扩散作用，填筑体往往具有二元结构特征，出现某填筑层上部碾压相对密实而平均颗粒小，下部相对疏松而平均颗粒较大，甚至存在架空现象。为了分析其变形机理，选取填筑单层的上部相对密实层与下部相对疏松(甚至架空)层厚度相当。上部相对密实层颗粒与下部相对疏松层颗粒结构界面，其下侧大颗粒R₁与上侧小颗粒粒径r_m2之比为n₁，其特征见图7。

图7. 填筑体二元结构界面模型

根据绘图统计分析(图7)，边界接触带里的一个较大圆形颗粒之上，接触排列的圆形小颗粒数n_o与粒径比n₁的统计关系为

$n_o=1.385\times n_1$ (个) (25)

可供小颗粒接触排列的弧长与粒径比n₁的统计关系为

$l=7.2496R_1\times n_1^{-0.093}$ (26)

相邻两个小颗粒圆心与大颗粒圆心连线的夹角与粒径比n₁的统计关系为

${\alpha }_0=60\times n_1^{-0.07455}$ (度) (27)

大颗粒与小颗粒接触排列特征见图7。实际土岩填筑体结构层中，可以根据检测开挖剖面初步确定二元结构的界限。分别取样测得界面上、下颗粒的平均粒径r_i和R_i，或者粒径随高度位置不同的变化情况，计算界面上、下填筑体粒径比n₁ = R₁/r_n，假设所有颗粒变形参数都为E和μ。

为了便于分析和阐明分析思路，做以下简化：1) 界面两侧的大颗粒和小颗粒以竖向对称分布；2) 除界面附近颗粒层外，其余大、小颗粒层均呈紧密性的三角形排列，小颗粒排列呈对称性；3) 与某大颗粒接触的小颗粒数取整数，见图7和图8。在分析时将等效大颗粒和等效小颗粒都划分为上、下半圆部分和左、右半圆部分，参考第2部分介绍的分析思路和方法，对大颗粒之上小颗粒单数(奇数)和双数(偶数)两种可能的情况进行分析。

3.1. 奇数小颗粒体结构模型分析

当小颗粒为奇数时，其接触排列模型见图8。设大颗粒填筑体的总厚度为h_d，则层数为 $m_d=h_d/\left(2R_i\right)$ ，小颗粒层的总厚度为h_x，则层数 $m_x=m_2=h_x/\left(2r_i\right)$ ；大、小颗粒层中颗粒间的排列结构均设为2.3节所

述的三角形排列。R_i为大颗粒层的颗粒等效半径，r_i为小颗粒层的颗粒等效半径。

1) 填筑体界面大颗粒层变形分析

根据图8所示模型，分为上半圆变形和下半圆变形两个部分。

上半圆部分的变形，由x及z方向的平衡条件和对称性有：

$q_{m2,1}+2q_{m2,2}\cos \left({\alpha }_2\right)+\cdots +2q_{m2,m}{}_1\cos \left({\alpha }_m{}_1\right)=2q_{m2}{R}_1$ (28)

$q_{m2,k}=q_{m2,1}\cos \left({\alpha }_k\right)$ , $k=1,\cdots ,m_1$ (29)

图8. 奇数小颗粒界面分析模型

其中， $m_1=1,\cdots ,\left(n-1\right)/2$ ， $n=\mathrm{int}\left(1.385n_1\right)$ ， $n_1=R_1/r_{m2}$ ，R₁及r_m2分别为界面两侧大颗粒等效圆半径和小颗粒等效圆半径； $i=1,\cdots ,m_2$ ，m₂由小颗粒填筑层厚度h₁、r_k和 ${\displaystyle \underset{k=1}{\overset{m2}{\sum }}2r_k=h_1}$ 确定。

$q_{m2}={\displaystyle \underset{k1=1}{\overset{m2}{\sum }}2{\gamma }_1}{r}_{k1}+{\sigma }_{i0}$ (30)

式中， ${\sigma }_{i0}$ 为填筑体表面超载在填土层中引起的竖向附加应力，由公式(7)计算；γ₁为小颗粒层等效重度。

${\alpha }_k=\frac{83.1n_1^{0.2545}\left(k-1\right)}{1.385n_1-1}$ (31)

见图8， $k=1,\cdots ,m_1$ 。利用前述的分析方法，

$a_{m2,k}={\left(\frac{4q_{m2,k}{R}_1}{\pi E^{*}}\right)}^{1/2}$ (32)

${\delta }_{m2,k}=\frac{\left(1-{\mu }^2\right)q_{m2,k}}{\pi E}\left\{2\ln \left(\frac{4R_1}{a_{m2,k}}\right)-1\right\}+R_1-{\left(R_1^2-a_{m2.k}^2\right)}^{1/2}$ (33)

上半圆的总变形 ${\delta }_{11}$ 为

${\delta }_{11}={\displaystyle \underset{k=1}{\overset{m1}{\sum }}{\delta }_{m2,k}\left(\cos \left({\alpha }_k\right)-\mu \sin \left({\alpha }_k\right)\right)}$ (34)

下半圆的变形分析，参考式(22)，在p₁作用下有

$a_{12}={\left(\frac{4p_1{R}_1}{\pi E^{*}}\right)}^{1/2}$ (35)

${\delta }_{12}=\frac{\left(1-{\mu }^2\right)p_1}{\pi E}\left\{2\ln \left(\frac{4R_1}{a_{12}}\right)-1\right\}+R_1-{\left(R_1^2-a_{12}^2\right)}^{1/2}$ (36)

其中， $p_1=\left(q_{m2}+2{\gamma }_2{R}_1\right)\cdot R_1$ ， $p_{21}=\frac{\mu }{1-\mu }{p}_1$

在两侧p₂₁作用下的变形 ${\delta }_{21}$ 为

$a_{21}={\left(\frac{4p_{21}{R}_1}{\pi E^{*}}\right)}^{1/2}$ (37)

${\delta }_{21}=\frac{\left(1-{\mu }^2\right)p_{21}}{\pi E}\left\{2\ln \left(\frac{4R_1}{a_{21}}\right)-1\right\}+R_1-{\left(R_1^2-a_{21}^2\right)}^{1/2}$ (38)

下半圆总变形 ${\delta }_{22}$ 为

${\delta }_{22}={\delta }_{12}\cos \left(30˚\right)-\mu \left({\delta }_{12}\sin \left(30˚\right)+{\delta }_{21}\right)$ (39)

填筑体界面大颗粒层的变形为式(34)和式(39)计算值之和。

2) 填筑体大颗粒段变形分析

根据简化说明，大颗粒填筑体除上述分析的界面层变形外，其下部的变形计算方法与2.3节基本相同，由式(24)计算，其中 $n=n_d-\text{1}$ 。再与式(34)和式(39)计算值相加，得到大颗粒填筑层的总变形量。

3) 填筑体小颗粒界面层变形分析

界面上侧小颗粒层的受力与大颗粒层受力是作用力与反作用力的关系，受大颗粒的“局部屏蔽作用”，忽略界面层小颗粒横向的受力，其变形分析可参考式(32)至式(34)，接触面半宽 ${a^{\prime }}_{m2,k}$ 为，

${a^{\prime }}_{m2,k}={\left(\frac{4q_{m2,k}{r}_{m2}}{\pi E^{*}}\right)}^{1/2}$ (40)

${{\delta }^{\prime }}_{m2,k}=\frac{1-{\mu }^2}{\pi E}\left\{\begin{array}{l}2q_{m2,k}\ln \left(\frac{4r_{m2}}{{a^{\prime }}_{m2-1,k}}\right)-1\\ 2q_{m2,k+1}\ln \left(\frac{4r_{m2}}{{a^{\prime }}_{m2,k+1}}\right)-1\end{array}\right\}+2r_1-{\left(r_1^2-{a^{\prime }}_{m2.k}^2\right)}^{1/2}-{\left(r_1^2-{a^{\prime }}_{m2,k+1}^2\right)}^{1/2}$ (41)

界面小颗粒层的变形 ${{\delta }^{\prime }}_{11}$ 为

${{\delta }^{\prime }}_{11}={\displaystyle \underset{k=1}{\overset{m1}{\sum }}{{\delta }^{\prime }}_{m2,k}\left(\cos \left({\alpha }_k\right)-\mu \sin \left({\alpha }_k\right)\right)/m}$ (42)

4) 填筑体小颗粒段变形分析

根据简化条件，小颗粒填筑体段除上述分析的界面层变形外，其上部的变形计算方法与2.3节基本相同，由式(24)计算，其中 $n=n_x-\text{1}$ 。再与式(34)和式(39)计算值相加就可得到小颗粒填筑段的总变形量。

3.2. 偶数小颗粒体结构模型分析

1) 填筑体界面大颗粒层变形分析

当小颗粒为偶数时，其模型见图9。

图9. 偶数小颗粒界面分析模型

所用参数与3.1节相同。此类情况下，大颗粒上半圆的变形计算公式(29)至公式(34)变为

$2{q^{\prime }}_{m2,1}\cos \left({\alpha }_1/2\right)+\cdots +2{q^{\prime }}_{m2,m}\cos \left({\alpha }_m-{\alpha }_1/2\right)=2q_{m2}{R}_1$ (43)

${q^{\prime }}_{m2,k+1}\cos \left({\alpha }_k\right)={q^{\prime }}_{m2,k}\cos \left({\alpha }_{k+1}\right)$ , $k=1,\cdots ,m$ (44)

其中， $m=1,\cdots ,n_1/2$ ，其余参数如上。

$b_{m2,k}={\left(\frac{4{q^{\prime }}_{m2,k}{R}_1}{\pi E}\right)}^{1/2}$ (45)

$u_{m2,k}=\frac{\left(1-{\mu }^2\right){q^{\prime }}_{m2,k}}{\pi E}\left\{2\ln \left(\frac{4R_1}{b_{m2,k}}\right)-1\right\}+R_1-{\left(R_1^2-b_{m2.k}^2\right)}^{1/2}$ (46)

上半圆的总变形 $u_{11}$ 为

$u_{11}={\displaystyle \underset{k=1}{\overset{m}{\sum }}{u}_{m2,k}\left(\cos \left(\left(k-0.5\right){\alpha }_1\right)-\mu \sin \left(\left(k-0.5\right){\alpha }_1\right)\right)}$ (47)

大颗粒界面层的变形应为公式(47)与公式(39)之和。

2) 填筑体大颗粒段变形分析

大颗粒填筑体段总的变形分析与3.1节第(2)段的分析方法相同，这里不再陈述。

3) 填筑体小颗粒界面层变形分析

界面侧小颗粒层的受力与大颗粒层受力是作用力与反作用力的关系，受大颗粒的“局部屏蔽作用”，忽略界面小颗粒层横向的受力，其变形分析可参考式(32)至式(34)，颗粒接触面半宽 ${b^{\prime }}_{m2,k}$ 为，

${b^{\prime }}_{m2,k}={\left(\frac{4{q^{\prime }}_{m2,k}{r}_{m2}}{\pi E}\right)}^{1/2}$ (48)

${u^{\prime }}_{m2,k}=\frac{2\left(1-{\mu }^2\right){q^{\prime }}_{m2,k}}{\pi E}\left\{2\ln \left(\frac{4r_{m2}}{{b^{\prime }}_{m2,k}}\right)-1\right\}+2r_1-2{\left(r_1^2-{b^{\prime }}_{m2.k}^2\right)}^{1/2}$ (49)

小颗粒界面层的变形 ${u^{\prime }}_{11}$ 为

${u^{\prime }}_{11}={\displaystyle \underset{k=1}{\overset{m}{\sum }}{u^{\prime }}_{m2,k}\left(\cos \left(\left(k-0.5\right){\alpha }_1\right)-\mu \sin \left(\left(k-0.5\right){\alpha }_1\right)\right)/m}$ (50)

4) 填筑体小颗粒段变形分析

根据简化条件，偶数小颗粒填筑体段的变形，参考2.3节的方法进行计算。先由式(24)计算，其中 $n=n_x-\text{1}$ 。再与式(34)和式(39)计算值相加就可得到小颗粒填筑段的总变形沉降量。

参考上述成都某新建国际机场场道，某段土岩料填筑体总厚度平均5.0 m的工程实际，调整一下参数。假设上、下厚度各为2.5 m，上部小颗粒层重度0.021 MN/m³，弹性模量50 MPa，泊松比0.28，颗粒等效平均半径r = 0.0125 m；下部大颗粒等重度0.019 MN/m³，弹性模量40 MPa，泊松比0.32，颗粒等效平均半径R = 0.05 m。此时，颗粒粒径比n₁ = 4。一个大颗粒上可排列6个小颗粒，即n₁ = 6为偶数，m = 3，其厚度加权平均参数值与上述均匀颗粒模型相当。利用二元结构分析方法计算结果见表1。显然，二元颗粒结构填筑体的沉降变形量增加了。

3.3. 土岩二元结构填筑体变形因素分析

为了解单层碾压填筑厚度对填筑体沉降变形的影响，考虑三种不同填筑厚度a、b和c，每种层厚均包含上述颗粒矩形排列、三角形排列及二元结构三种排列结构情况，分别用A、B、C表示；其中A、B的土颗粒平均粒径均为25 mm，C的上部密实层平均粒径为10 mm，下部层平均粒径为40 mm，厚度加权平均颗粒半径为25 mm。采取不同的组合方式构成的分析模型情况见表2；根据填筑料特性，确定的分析参数见表3。

表1. 土石填筑体模型及其变形(mm)

表2. 模型特征参数

表3. 模型物理力学参数

根据前文所推导的理论公式，计算得到上述9种组合方式的填筑层沉降量见图10。

图10. 九种组合方式填筑层的变形–深度曲线

从图10可以看出，所建理论模型宏观变形总趋势与普通填筑体宏观固结变形规律吻合。填筑体单层计算厚度增大，其总宏观变形量增加，而二元结构模型由于上部密实层与下部架空层的平均粒径差异较大，结构界面附近的局部变形差异明显，具有变形突变特征，表现为密实层与架空层在层厚相同时，上部密实层总宏观变形量显著大于下部架空层总宏观变形量，这是由于密实层比架空层的颗粒的平均粒径更小，层内的孔隙率更大，可压缩变形的空间更多，最终的沉降量也更大。

图10中由不同模型细观结构得到的结果表明，在模型土颗粒粒径一致且单层厚度相同时，颗粒呈三角形排列时得到的总宏观变形量比矩形排列时的总宏观变形量更小，反映颗粒的细观结构对宏观变形量有重要影响，这里实质上体现的是细观结构中的孔隙率对沉降量的显著影响，即模型颗粒粒径一致时，呈矩形排列时的孔隙率大，因而最终的变形量大，呈三角形排列时的孔隙率小，因而最终的变形量小。二元结构模型因其较为特殊的排列，上部密实层粒径取为下部架空层粒径的1/2且下部架空层颗粒呈三角形排列，此时得到的变形量位于前两者变形量数值区间内，需要指出的是，本文算例中只选取了二元结构模型的一种情况计算，可见二元结构模型的层间组合、层间粒径比、粒径大小等组合变化形式多样，其沉降变形量的计算也随之更加复杂。

4. 等效圆柱形粒填筑体变形数值模拟分析

4.1. 数值模型构建

鉴于细观结构的非连续性，采用PFC^2D颗粒离散元程序进行模拟，以二维圆球(ball)为基本颗粒单元，在设置边界(wall)范围内规则生成呈矩形、三角形排列的模型，此外再生成二元结构的堆积模型，见图11。

图11. A至C模型的数值模型

文中采用以下方法建立A至C模型：首先根据构建的模型层厚选取相应高度同时高宽(径)比为2:1的模型尺寸，分别建立呈矩形、三角形规则排列及二元结构排列的圆盘颗粒模型；其次设定颗粒的物理参数和接触本构参数，物理参数通过模型基本特性确定，接触本构参数通常采用参数反演试算法对细观参数进行确定，即不断调整细观参数直至PFC模型与室内模型试验表现出相同的宏观响应，由于暂缺相应的室内模型试验，同时本文研究的关注点在模型的压缩变形量而不是其强度特征，为了模拟理论计算中的理想小刚球模型的变形量，在模拟过程中颗粒之间采用线性接触刚度模型，表4为本文所选取的PFC模拟的微观参数。

表4. PFC^2D模拟的微观参数

同时，用四个墙单元(wall)模拟上、下边界及两侧边界，通过伺服程序调节侧边界的位置以保持稳定的围压，由于PFC软件中无法通过对墙体单元施加荷载的方式对模型施加荷载，因此本文是通过在模型顶部预设厚度为0.5 mm的颗粒层作为加载区域，监测加载区域中心颗粒的变形量作为受载后的试样的变形量。

4.2. 数值结果与理论结果比较分析

图12为等效球形粒空隙体结构及二元结构模型最终变形量的理论计算与数值模拟结果对比曲线，同时最终变形量随着颗粒划分的粒径不同而有显著差异。

从对比分析结果来看，当细观模型颗粒为规则排列，且颗粒粒径一定时，数值模拟结果仍出现呈矩形排列的球形粒空隙体比呈三角形排列时的最终变形量明显偏大，表明颗粒的细观结构对变形量有重要影响。当细观模型颗粒呈二元结构排列时，上部颗粒粒径一致时，其最终变形量略小于呈矩形排列模型而大于呈三角形排列模型，这是因为虽然试算的二元结构模型上半部分的颗粒是由于重力沉积法生成而导致孔隙率较矩形排列模型偏大，但其下半部分的颗粒粒径不仅是上半部分颗粒粒径的2倍且同时呈三角形排列，因此二元结构模型的最终变形量是模型颗粒粒径及孔隙率对最终变形量的综合影响结果。

同时，当细观模型颗粒的排列模式一定时，随着颗粒粒径减小，理论计算与数值模拟结果均增大，且当颗粒半径小于1.0 mm后，最终变形量随着颗粒粒径的进一步减小会显著增大。进一步说明了模型颗粒粒径及孔隙率对最终变形量的重要影响。

图12. 填筑体的最终变形量

5. 结论

在土石填料碾压填筑体二元结构模型概化模型下，利用理论分析及数值模拟方法进行亚久，得出的结论如下：

1) 等径圆柱体接触问题弹性理论解，可以用于模拟计算土石料分层填筑碾压二元结构体宏观变形的分析，能够较好的描述实际工程问题。

2) 二元结构模型是较为特殊的排列模式，其上部为孔隙率较小的密实层，下部为孔隙率较大的相对松散层，在本文所设置的参数下，其沉降变形与土石料颗粒平均粒径、单层细观结构和二元结构颗粒粒径组合等因素有关；二元结构界面附近局部变形差异明显，具有变形突变特征。

3) 土石料填筑体细观结构对其变形特征有显著影响；表现为填筑体单层计算厚度增大其宏观总变形增大，孔隙率较大的等效球形矩形排列填筑体的变形大于孔隙率较小的等效球形三角形排列填筑体的变形。

4) 数值模拟结果表明：当颗粒半径小于1.0 mm，变形量随细观结构颗粒粒径的继续减小而增大。需要指出的是，由于数值模拟过程中的细观模型颗粒是可以通过平/转动对细观结构进行调整的，因此数值模拟得到的沉降变形量较理论计算值更大。

基金项目

本文受国家自然科学基金(41272321，51808458)，中国博士后科学基金(2018M640934)资助。

文章引用

钟志彬,邓荣贵,孙怡,王园园,董玲. 泥质岩土石料碾压填筑体变形特性及机理分析
Study on the Deformation Characteristics and Mechanism of Argillaceous Earth-Rock Compacted Fillings[J]. 土木工程, 2021, 10(09): 961-976. https://doi.org/10.12677/HJCE.2021.109107

参考文献

[1] 陈从新, 卢海峰, 袁从华, 等. 红层软岩变形特性试验研究[J]. 岩石力学与工程学报, 2010(2): 261-270.

[2] 蒋关鲁, 房立凤, 王智猛, 等. 红层泥岩路基填料动强度和累积变形特性试验研究[J]. 岩土工程学报, 2010(1): 124-129.

[3] 程强, 周永江, 吉随旺, 等. 红层软岩基本特性和公路建设中的主要技术问题[C]//中国科协第五届青年学术年会. 北京: 中国科学技术出版社, 2004: 468-469.

[4] 周应华, 周德培, 封志军. 三种红层岩石常规三轴压缩下的强度与变形特性研究[J]. 工程地质学报, 2005(4): 477-480.

[5] 章清叙, 葛修润, 黄铭, 等. 周期荷载作用下红砂岩三轴疲劳变形特性试验研究[J]. 岩石力学与工程学报, 2006(3): 473-478.

[6] 郭永春, 谢强, 文江泉. 红层岩土水理性质工程判别准则试验研究[J]. 水文地质工程地质, 2008(4): 71-74.

[7] 郭永春. 红层岩土中水的物理化学效应及其工程应用研究[D]: [博士学位论文]. 成都: 西南交通大学, 2007.

[8] Ballard, T.J., et al. (1994) Fundamentals of Shale Stabilization: Water Transport through Shales. SPE Formation Evaluation, 9, 129-134. https://doi.org/10.2118/24974-PA

[9] Chenevert, M.E. (1970) Shale Control with Balanced-Activity Oil-Continuous Muds. Journal of Petroleum Technology, 22, 1309-1316. https://doi.org/10.2118/2559-PA

[10] 宋磊. 红层软岩遇水软化的微细观机理研究[D]: [硕士学位论文]. 成都: 西南交通大学, 2011.

[11] 周翠英, 李伟科, 向中明, 等. 水-应力作用下软岩细观结构摩擦接触分析[J]. 岩土力学, 2015, 36(9): 2458-2466.

[12] 邓华锋, 周美玲, 李建林, 等. 水-岩作用下红层软岩力学特性劣化规律研究[J]. 岩石力学与工程学报, 2016(S2): 3481-3491.

[13] 周美玲. 考虑水-岩作用的三峡库区红层软岩蠕变特性研究[D]: [硕士学位论文]. 宜昌: 三峡大学, 2016.

[14] 邓华锋, 方景成, 李建林, 等. 含水状态对红层软岩力学特性影响机理[J]. 煤炭学报, 2017, 42(8): 1994-2002.

[15] Hertz, H. (1881) On the Contact of Elastic Solids. Journal für die reine und angewandte Mathematik, 92, 156-171.