Implicit DEM analyses of size and shape effects on crushing strength of rockfill particles
-
摘要: 采用隐式离散元法（DEM），亦称非光滑接触动力学（NSCD）法，模拟了堆石料单颗粒的一维压缩破碎过程。采用多面体颗粒近似真实堆石料颗粒的不规则棱角状形状，并通过Voronoi空间划分将母颗粒离散成若干多面体子颗粒，子颗粒间通过接触面上的内聚力模型（CZM）胶结发生相互作用，胶结断裂则可反映颗粒的破碎过程，可避免传统DEM中常用的碎片替代法和胶结小球法的缺点。首先使用巴西劈裂试验获得CZM参数，其后开展了不同粒径颗粒的一维压缩试验，结果表明颗粒破碎强度服从Weibull分布，且破碎强度大小和变异性均随着颗粒尺寸的增大而减小。还进一步模拟了不同形状颗粒的压缩试验，发现颗粒形状与加载方向均对颗粒破碎强度有显著影响。长轴方向加载破碎强度通常最低，且在相同等效粒径下，球状、椭球状、扁平状颗粒的平均破碎特征强度依次降低。Abstract: An implicit version of discrete element method (DEM) called non-smooth contact dynamics (NSCD) is used to simulate the crushing process of single rockfill particle under one-dimensional (1D) compression. The angular and irregular shapes of rockfill particles are represented using polyhedrons, which are discretized into smaller elementary cells through the Voronoi tessellation. The interaction between neighboring elementary cells is described by the cohesive zone model (CZM), where hybrid tensile and shear failure modes are considered. Consequently, particle crushing can be captured by the breakage of CZM bonds subjected to external loading. The proposed method for modeling the particle crushing avoids drawbacks of the traditional fragment replacement method and bonded sphere method in DEM. The Brazilian splitting tests on granite are performed to calibrate the parameters for CZM. 1D compression tests on single particles with different sizes and shapes further reveal that the crushing strength of particles follows the Weibull distribution, and both the magnitudes and the variances decrease with the increasing particle size. For the particles with different principal axes (e.g., elongated and platy), the crushing strength loaded in the major axis is the smallest among all directions. It is also shown that the elongated and platy particles have smaller average strength than the spherical ones given the same equivalent particle size.
-
0. 引言
地下洞室对弹性波的散射问题在地震工程、地下结构抗爆、地球物理勘探等领域有着广泛的应用。近年来，许多研究人员基于波函数展开法[1-2]对单相介质中衬砌结构的动力响应进行了探讨。然而，土和岩石具有多孔性和多相性，单相弹性连续介质无法准确描述实际情况。因此，波在两相介质中的传播和动力响应问题受到广泛关注。应用Biot饱和多孔介质波动理论，周香莲等[3]求解了全空间圆柱形衬砌洞室对平面波散射的解析解。Jiang等[4]引入大圆弧假设，解决了平面P波入射下饱和多孔弹性半空间中圆形衬砌洞室周围的散射问题。为了更符合实际工程，张海等[5]使用辅助函数配合波函数展开法和傅立叶级数展开方法，进一步研究了SH波入射下含直边界半圆形衬砌洞室的动力响应。Ding等[6]利用波函数展开法，给出了平面P波入射下饱和多孔弹性介质中复合衬砌动应力集中的解析解。Hasheminejad等[7]解析分析P波和SV波作用下饱和多孔弹性介质中变壁厚内衬的圆柱形洞室的动应力集中问题。Liu等[8]基于平面复变理论得到了多孔弹性半空间中衬砌洞室对P1和SV波散射的解析解。徐长节等[9]，Xu等[10]基于非局部弹性理论分析了孔径效应和孔隙动力学效应对饱和土中圆柱形衬砌动力响应的影响。Ding等[11]基于非局部Biot理论得到了平面P波入射下饱和多孔介质中浅埋复合衬砌洞室动力响应的解析解。
综上，饱和多孔介质中衬砌洞室的动力响应已经得到了广泛而充分的研究。可实际上，土是一种三相介质，即使地下水位以下的土，也可能存在不完全饱和的情况。Li等[12]总结了众多非饱和孔隙弹性模型以及非饱和多孔介质波动特性的研究，表明由于气相的存在，非饱和介质的动态响应在许多方面与传统饱和多孔介质的动态反应显著不同。然而，目前有关非饱和多孔介质中衬砌洞室对地震波散射问题的研究还较为匮乏。Li等[13]研究了非饱和多孔介质中圆形洞室在内部荷载作用下的动力响应。Tan等[14]研究了无限非饱和多孔弹性介质中不完全接触的衬砌洞室对平面P波的散射问题。以考虑混合物理论的非饱和多孔介质动力理论为基础，采用波函数展开法和分离变量法，分别推导了平面P波和SV波入射时衬砌洞室周围位移应力的解析解。分析了不同入射频率下饱和度对衬砌洞室周围动应力集中系数的影响。
1. 分析模型
1.1 数学模型和基本方程
为了简化研究模型，把非饱和多孔介质中圆柱形衬砌洞室，视为无限非饱和多孔介质全空间中一个无限长度的衬砌洞室（图 1）。
衬砌洞室的内径为a，外径为b。本研究基于Wei等[15]提出的非饱和多孔介质动力学理论。控制方程表示为
nS0ρS0¨uS=(MSS+G)∇∇⋅uS+G∇⋅∇uS+MSW∇∇⋅uW+MSN∇∇⋅uN+ˆμW(˙uW−˙uS)+ˆμN(˙uN−˙uS)，nW0ρW0¨uW=MSW∇∇⋅uS+MWW∇∇⋅uW+MWN∇∇⋅uN−ˆμW(˙uW−˙uS)，nN0ρN0u ¨N=MSN∇∇⋅uS+MWN∇∇⋅uW+MNN∇∇⋅uN−ˆμN(˙uN−˙uS)。} (1) 式中：上角标S，W，N分别为非饱和介质中的固体骨架部分、液相部分和气相部分；G为非饱和多孔介质的剪切模量；nα0,ρα0,uα分别为各相的初始体积分数、初始密度和初始位移；MSS，MWW，MNN，MSW，MSN，MWN为与体积模量、剪切模量及模量间耦合相关的弹性系数；ˆμf为渗透率相关的系数，f=W,N。Wei等[15]给出了相关参数的表达式。采用Van等[16]提出的土水特征曲线V-G模型。
非饱和多孔介质中的本构关系表示为
σS=(MSS∇⋅uS+MSW∇⋅uW+MSN∇⋅uN)I+2G[∇uS+(∇uS)T]，σw=(MSW∇⋅uS+MWW∇⋅uW+MWN∇⋅uN)I ，σN=(MSN∇⋅uS+MWN∇⋅uW+MNN∇⋅uN)I 。} (2) 式中：σα为各相的应力张量，α=S,W,N；I为单位张量矩阵。
衬砌介质被视为连续、各向同性的单相弹性介质，衬砌介质位移为uL，衬砌介质的拉梅常数为λL，μL，衬砌密度为ρL。
1.2 动力学方程的求解
根据亥姆霍兹(Helmholtz)矢量分解原理，非饱和多孔介质方程中的位移可以表示为两部分：
uα = ∇φα + ∇×Ψα(α=S,W,N)。 (3) 式中：φα,Ψα分别为非饱和多孔介质中α相的标量势函数和矢量势函数。
本文考虑简谐波入射，势函数可表示为φ= φ(x,y,z)e−iωt。将式（3）代入方程（1）可得
(∇2−k2j)φSj=0，φWj=δWjφSj，φNj=δNjφSj(j=1,2,3)，} (4a) (∇2−k24)ΨS=0，ΨW=ζWΨS，ΨN=ζNΨS。} (4b) 从式（4a），（4b）中可以看出非饱和多孔介质中有3种压缩波1种剪切波，ki为3种压缩波波数, i=1,2,3，k4为剪切波的波数，δWj，δNj分别为液相、气相压缩波参与系数，j=1,2,3; ζW，ζN分别为液相、气相剪切波参与系数。详细表达式参见文献[13]。
将式（4a），（4b）中固相波函数转换到圆柱坐标系下并使用分离变量法求解，可得固体骨架的标量势函数和矢量势函数的通解为
φSj(r,θ)=∞∑n=0An,jKn(kjr)cosnθ (j=1,2,3)， (5) ΨS(r,θ)=∞∑n=0BnKn(k4r)sinnθ。 (6) 式中：An,j，Bn为待定系数；Kn为第二类修正Bessel函数，下标n代表阶数。
同理可以得到衬砌介质中标量势函数和矢量势函数的通解：
φL(r,θ)=∞∑n=0(A(1)nIn(kL1r)+A(2)nKn(kL1r))einθ， (7) ΨL(r,θ)=∞∑n=0(B(1)nIn(kL2r)+B(2)nKn(kL2r))einθ。 (8) 式中：φL，ΨL分别表示衬砌介质中的标量势函数和矢量势函数；A(1)n，A(2)n，B(1)n和B(2)n为待定系数，分别表示P和SV波的振幅；In为第一类虚宗量Bessel函数，下标n代表阶数。
2. 波场分析与求解
2.1 波场分析
假设具有圆频率的平面谐波P1或SV波沿x轴正方向传播，用贝塞尔-傅立叶级数表示为
入射P1波：φi=φ0∞∑n=0εninJn(k0r)cos nθe−iωt， (9) 入射SV波：Ψi=Ψ0∞∑n=0εninJn(k0r)cos nθe−iωt。 (10) 式中：k0=ω/cp(入射P1波)或者k0=ω/cs(入射SV波)表示入射波的波数，φ0(Ψ0)表示入射波振幅，Jn(.)表示第一类n阶贝塞尔函数，若n=0则εn=1，若n⩾则 {\varepsilon _n} = 2 且 i = \sqrt { - 1} ，t表示时间， {\text{e}^{ - \text{i}\omega t}} 为时间因子将从后续表达式中省略。
入射平面波进入非饱和多孔介质，由于衬砌洞室的存在，会在非饱和多孔介质中产生散射波，在衬砌介质中产生透射波和内表面的反射波。以P1波入射为例，非饱和多孔介质中的波场可表示为
\left. \begin{array}{l}{\varphi }^{\text{S}}={\varphi }_{0}{\displaystyle \sum _{n=0}^{\infty }{\varepsilon }_{n}}{i}^{n}{\text{J}}_{n}({k}_{\text{0}}r)\text{cos }n\theta \text{+}{\displaystyle \sum _{n=0}^{\infty }{A}_{n,j}{\text{K}}_{n}}({k}_{j}r)\cos n\theta \text{，}\\ {\varphi }^{\rm{W}}={\delta }_{1}^{\rm{W}}{\varphi }_{0}{\displaystyle \sum _{n=0}^{\infty }{\varepsilon }_{n}}{i}^{n}{\text{J}}_{n}({k}_{\text{0}}r)\text{cos }n\theta \text{+}\\ \;\;\;\;\;\;\;\;{\displaystyle \sum _{j=1}^{3}{\delta }_{j}^{\rm{W}}{\displaystyle \sum _{n=0}^{\infty }{A}_{n,j}{\text{K}}_{n}}({k}_{j}r)\cos n\theta }\text{，}\\ {\varphi }^{\rm{N}}={\delta }_{1}^{\rm{N}}{\varphi }_{0}{\displaystyle \sum _{n=0}^{\infty }{\varepsilon }_{n}}{i}^{n}{\text{J}}_{n}({k}_{\text{0}}r)\text{cos }n\theta \text{+}\\ \;\;\;\;\;\;\;\;{\displaystyle \sum _{j=1}^{3}{\delta }_{j}^{\rm{N}}{\displaystyle \sum _{n=0}^{\infty }{A}_{n,j}{\text{K}}_{n}}({k}_{j}r)\cos n\theta \text{ }}\text{，}\end{array} \right\} (11a) \left. \begin{array}{l}{\varPsi }^{\rm{S}}={\displaystyle \sum _{n=0}^{\infty }{B}_{n}{\text{K}}_{n}}({k}_{4}r)\sin n\theta \text{，}\\ {\varPsi }^{\rm{W}}={\zeta }^{\rm{W}}{\varPsi }^{\rm{S}}\text{，}\\ {\varPsi }^{\rm{N}}={\zeta }^{\rm{N}}{\varPsi }^{\rm{S}}。\end{array} \right\} (11b) 衬砌介质中的总波场可表示为
\left. \begin{array}{l}{\varphi }_\text{L}={\displaystyle \sum _{n=0}^{\infty }\left({A}_{n}^{{}^{(1)}}{\text{I}}_{n}({k}_{\text{L}_{1}}r)+{A}_{n}^{{}^{(2)}}{\text{K}}_{n}({k}_{\text{L}_{1}}r)\right)}\cos n\theta \text{，}\\ {\varPsi }_\text{L}={\displaystyle \sum _{n=0}^{\infty }\left({B}_{n}^{{}^{(1)}}{\text{I}}_{n}({k}_{\text{L}_{2}}r)+{B}_{n}^{{}^{(2)}}{\text{K}}_{n}({k}_{\text{L}_{2}}r)\right)}\sin n\theta 。\end{array} \right\} (12) 2.2 边界条件
当 r = a 时，在衬砌介质中：
\left. \begin{array}{l}{\sigma }_{rr}^{\text{L}}=0\text{，}\\ {\sigma }_{r\theta }^{\text{L}}=0。\end{array} \right\} (13) 当 r = b 时，非饱和多孔介质与衬砌介质满足连续变形条件与界面不透水边界条件：
\left. \begin{array}{l}{u}_{r}{}^{\rm L}={u}_{r}{}^{\rm{S}}={u}_{r}{}^{\rm{W}}={u}_{r}{}^{\rm{N}}\text{，}\\ {u}_{\theta }{}^{\rm L}={u}_{\theta }{}^{\rm{S}}\text{，}\\ {\sigma }_{r\theta }{}^{\text{L}}={\sigma }_{r\theta }{}^{\rm{S}}\text{，}\\ {\sigma }_{rr}{}^{\text{L}}={\sigma }_{rr}{}^{\rm{S}}+{\sigma }^{\rm{W}}+{\sigma }^{\rm{N}}。\end{array} \right\} (14) 式中： u_{^r}^\alpha ， u_\theta ^\alpha 分别为非饱和多孔介质各相位移法向和切向分量， \alpha = \text{S,W,N} ； {u_r}^{\rm L} ， {u_\theta }^{\rm L} 分别为衬砌的法向、切向位移分量； {\sigma _{rr}}^{\rm S} ， {\sigma _{r\theta }}^{\rm S} 分别为非饱和多孔介质固相法向、切向应力分量； {\sigma ^{\rm W}} ， {\sigma ^{\rm N}} 分别为非饱和多孔介质的孔隙水压力和孔隙气压力； {\sigma _{rr}}^{\text{L}} ， {\sigma _{r\theta }}^{\text{L}} 分别为衬砌法向和切向应力分量。
2.3 求解
将非饱和多孔介质和衬砌介质中总波场势函数的式（11），（12），代入应力和位移的式（2）中，根据边界条件，可得到波场中未知幅值系数的方程：
\left[\begin{array}{llllllll} G_1 & G_2 & G_3 & G_4 & G_5 & G_6 & G_7 & G_8 \\ H_1 & H_2 & H_3 & H_4 & H_5 & H_6 & H_7 & H_8 \\ T_1 & T_2 & T_3 & T_4 & 0 & 0 & 0 & 0 \\ U_1 & U_2 & U_3 & U_4 & 0 & 0 & 0 & 0 \\ V_1 & V_2 & V_3 & V_5 & -\left.D_1\right|_{r=b} & -\left.D_2\right|_{r=b} & -\left.D_3\right|_{r=b} & -\left.D_4\right|_{r=b} \\ Y_1 & Y_2 & Y_3 & Y_5 & -\left.F_1\right|_{r=b} & -\left.F_2\right|_{r=b} & -\left.F_3\right|_{r=b} & -\left.F_4\right|_{r=b} \\ 0 & 0 & 0 & 0 & \left.D_1\right|_{r=a} & \left.D_2\right|_{r=a} & \left.D_3\right|_{r=a} & \left.D_4\right|_{r=a} \\ 0 & 0 & 0 & 0 & \left.F_1\right|_{r=a} & \left.F_2\right|_{r=a} & \left.F_3\right|_{r=a} & \left.F_4\right|_{r=a} \end{array}\right]\left[\begin{array}{l} A_{n, 1} \\ A_{n, 2} \\ A_{n, 3} \\ B_n \\ A_n^{(1)} \\ A_n^{(2)} \\ B_n^{(1)} \\ B_n^{(2)} \end{array}\right]=\left[\begin{array}{l} G_9 \\ H_9 \\ T_5 \\ U_5 \\ -V_4 \\ -Y_4 \\ 0 \\ 0 \end{array}\right] . (15) 系数 {G}_{i}，{H}_{i}，{D}_{i}，{F}_{i}，{T}_{i}，{U}_{i}，{V}_{i}和{Y}_{i} 的表达式由上述公式计算得出，表达式略。通过求解式（15），得到未知波幅系数，以确定非饱和多孔介质和衬砌中的波场，从而获得其中位移和应力的解析解答。
3. 验证和算例
3.1 退化验证
将衬砌外内半径之比取1.000001，其结果将会退化为非饱和多孔介质中无衬砌洞室的情况。与文献[13]计算结果对比，取参数（n0=0.34，Sr=0.6），图 2给出了P波入射下衬砌洞室与空腔洞室动应力集中因子对比结果。从图 2中可以看出，当衬砌很薄时，衬砌洞室退化后的结果与洞室的计算结果一致。
![]()
图 2 非饱和多孔介质中衬砌外内半径之比为1.000001时衬砌洞室与空腔洞室集中因子对比结果Figure 2. Comparison results of concentration factors between cavity and lining cavity with 1.000001 of ratio of outer radius to inner radius in unsaturated porous media3.2 算例分析
工程抗震研究的重点是衬砌的动应力集中因子(dynamic stress concentration factor，简称DSCF)，本文中考虑介质中环向应力与入射波产生的环向应力的比值，表达式为 {\sigma _{\theta \theta }}^{\rm S}/{\sigma _0} ，其中
\begin{align} {\sigma _{\theta \theta }}^{\rm S} =& \left( {{M_\text{SS}}{\mathit{\nabla }^2}{\varphi ^{\rm S}} + {M_\text{SW}}{\mathit{\nabla }^2}{\varphi ^{\rm W}} + {M_{\rm{SN}}}{\mathit{\nabla }^2}{\varphi ^{\rm N}}} \right) + \\ &\frac{{2G}}{r}\left( {\frac{{\partial {\varphi ^{\rm S}}}}{{\partial r}} + \frac{1}{r}\frac{{{\partial ^2}{\varphi ^{\rm S}}}}{{\partial {\theta ^2}}}} \right) - 2G\left[ {\frac{\partial }{{\partial r}}\left( {\frac{1}{r}\frac{{\partial {\varPsi ^{\rm S}}}}{{\partial \theta }}} \right)} \right] \text{，} \end{align} (16) \begin{array}{c} {\sigma _0} = [2G + {M_\text{SS}} + {M_\text{SW}} + {M_\text{SN}} + \delta _1^{\rm W}({M_\text{WW}} + {M_\text{SW}} + {M_\text{WN}}) + \\ {\delta }_{1}^{\rm{N}}\left({M}_\text{SN}+{M}_\text{WN}+{M}_\text{NN}\right)]{k}_{1}^{2}\text{ (P}波入射\text{)} \text{，} \end{array} (17a) {\sigma _0} = Gk_4^2\;\;\;\;\;\; (\text{SV}波入射) 。 (17b) 从应力表达式中可以看到动剪切模量是一个关键因素，在以往的研究中很少将饱和度作为动剪切模量的考量，实际饱和度的变化对土中动剪切模量影响很大[17]。采用下式计算非饱和土动剪切模量[17]：
G = {G_0} + \frac{{2050}}{\alpha }\ln (\sqrt {{{({S_e})}^{ - 2}} - 1} + {({S_e})^{ - 1}})\tan \phi ' 。 (18) 式中： {G_0} 为土体饱和时动剪切模量； \phi ' 为土体饱和时内摩擦角。
本文分别分析P波和SV波入射时，不同频率下饱和度对衬砌洞室周围动应力集中因子的影响。为避免孔洞尺寸对分析结果的影响采用无量纲频率 \eta = a \times {k_0} ，式中a为衬砌内半径，k0为入射波波数。衬砌介质的参数取 {\lambda _{\rm{L}}} = 3.5×107 kPa， {\mu _{\rm{L}}} = 3.5×107 kPa， {\rho _{\rm{L}}} = 2500{\text{ }}\rm kg/{m^3} 。非饱和多孔介质的材料参数及土水特征曲线参数 \alpha ， d ， {S_{r\rm{W}}}{\text{ }} ， {S_{r\rm{N}}}{\text{ }} ，如表 1所示[13, 17]。
表 1 非饱和多孔介质材料参数Table 1. Material parameters of unsaturated porous medium孔隙率 {n_0} \rho _0^{\rm S}/(kg·m-3) \rho _0^{\rm W}/(kg·m-3) \rho _0^{\rm N}/(kg·m-3) {K_\rm{S}}/\rm kPa 0.36 2650 1000 1.1 3.6×107 {K_\rm{W}}/\rm kPa {K_\rm{N}}/\rm kPa {\nu ^{\rm W}} /(\rm Pa \cdot s) {\nu ^{\rm N}} /(\rm Pa \cdot s) {G_0}/\rm kPa 2×106 110 1.0×10-3 1.8×10-5 9.23×107 \widehat K{\text{/}} kPa \kappa /(\rm m \cdot {\rm s^{ - 1}} ) {\alpha _\text{B}} {S_{r\rm{W}}}{\text{ }} {S_{r\rm{N}}}{\text{ }} 2.0 \times {10^5} 2.5 \times {10^{ - 12}} 1 0.05 1 \upsilon d \alpha \phi ' /(°) 0.3 2 2.0 \times {10^{ - 5}} 20 取衬砌外内半径之比为 \vartheta = 1.1 ，分别计算 \eta 为0.25，1.0，2.0情况下饱和度为Sr为0.2，0.6，0.8，1.0时衬砌外壁动应力集中因子分布情况。图 3，4分别给出了P波和SV波入射时衬砌洞室周边动力集中因子分布曲线。从图中可以看出，P波和SV波入射时衬砌周边的应力均对称分布在入射波两侧，饱和介质与非饱和介质中动应力分布有明显差异。P波入射时，饱和度变化对动应力分布形状影响较大，饱和时动应力分布形状与非饱和状态时有显著不同；SV波入射时，饱和度变化对动应力分布形状影响较小，饱和时动应力大于非饱和状态。P波入射下，当 \eta =0.25时，Sr < 1.0动应力集中因子的最大值随饱和度的增加逐渐增大，在Sr为0.2，0.6时分布形式较一致呈“8”字形，最大值出现在80°，280°；在Sr=0.8时图像分布形式变为3个凸起，最大值出现在70°，180°和290°。Sr=1.0时动应力集中因子呈椭圆形沿洞中均匀分布，最大值角度出现在80°，280°。当 \eta =1.0时，Sr < 1.0动应力集中因子最大值随饱和度的增加逐渐增大，动应力集中因子分布呈“豌豆状”，随饱和度增加逐渐饱满。Sr=1.0时动应力集中因子呈椭圆形分布在坐标系的中上部，最大值出现在0°，360°。当 \eta =2.0时，动应力集中因子最大值随饱和度的增加逐渐增大，曲线集中分布在坐标系的上半部分，随饱和度增加最大值出现位置逐渐向0°坐标轴靠拢。P波入射时，最小值位置出现在入射面和入射背面，基本不受饱和度和入射频率的影响。
![]()
图 3 P波入射时考虑饱和度对剪切模量的影响下衬砌洞室周边动力集中因子分布曲线Figure 3. Curves of dynamic stress concentration factors around lining cavity considering influences of saturation on shear modulus under incidence of P-waves![]()
图 4 SV波入射时考虑饱和度对剪切模量的影响下衬砌洞室周边动力集中因子分布曲线Figure 4. Curves of dynamic stress concentration factors around lining cavity considering influences of saturation on shear modulus under SV-waves由于P波与SV波入射的振动形式不同，衬砌洞室周边动应力集中因子分布不同，在SV波入射时动应力集中因子呈蝶状分布。对比图 3，4可以发现不同相对频率下，饱和时动应力集中因子最大值要大于入射 {P_{}} 波时的动应力集中因子。在图 4中可以发现，SV波入射时，饱和与非饱和情况存在明显差异，相比于Sr=1.0，Sr < 1.0时环向应力数值更小。在Sr < 1.0时饱和度对动应力集中因子影响较小，但随相对频率变化动应力集中因子变化明显，整体上随着入射频率的增加衬砌周边动应力逐渐增加，由于Sr=1.0与Sr < 1.0增幅不一致，在相对频率为2.0时饱和与非饱和动应力集中因子最接近。Sr < 1.0时最大值出现角度不变为40°，320°；Sr=1时最大值出现角度变化为50°，310°；30°，330°；20°，340°。 SV 波入射时最小值出现在轴线上。
4. 结论
推导了平面P波和SV波入射下非饱和多孔介质中内嵌圆柱形衬砌洞室的动力响应解析解。考虑了饱和度对剪切模量的影响，求解并分析了饱和度和频率对圆柱形衬砌洞室动应力集中因子的影响，得到3点结论。
（1）饱和度对动应力集中因子的影响不容忽视，无论P波还是SV波入射，饱和时动应力分布与非饱和时有显著不同。且P波入射时，饱和度对动应力集中因子的分布形式影响较大，随着饱和度增加，动应力因子最大值出现角度逐渐向0°偏移，动应力集中因子最大值逐渐增加；SV波入射下，饱和度小于1时，饱和度对动应力集中因子影响较小。
（2）由于P波与SV波的振动形式不同，衬砌洞室周边动应力集中因子分布不同。P波入射时，动应力集中因子分布更多呈椭圆形或“豌豆”状；SV波入射时，动应力集中因子分布呈蝶状。
（3）入射频率主要影响动应力集中因子的大小。随着入射频率的增加，P波和SV波入射下洞室周边动应力集中因子最大值均逐渐增加。
致谢: 感谢香港科技大学朱凡博士为本研究内容提出的宝贵意见。 -
![]()
图 1 单边接触Signorini-Coulomb条件
Figure 1. Signorini-Coulomb condition for a unilateral contact
![]()
图 3 巴西劈裂试验荷载-位移曲线（内嵌图为岩柱断裂模式）
Figure 3. Load-displacement curves in Brazilian disk split tests (inset shows fracture pattern of disk)
![]()
图 5 单颗粒一维压缩荷载‒位移曲线（内嵌图为颗粒破碎形态）
Figure 5. Load-displacement curve for single-particle compression tests (inset shows crushed particle)
![]()
图 6 部分生成颗粒的Voronoi离散
Figure 6. Examples of particle discretization using Voronoi tessellation
![]()
图 7 等效直径27.5 mm单颗粒破碎强度Weibull分布
Figure 7. Weibull distribution of crushing strength for particles with equivalent diameter of 27.5 mm
![]()
图 8 不同尺寸颗粒的Voronoi离散
Figure 8. Discretization of particles with different sizes using Voronoi tessellation
![]()
图 9 特征强度与颗粒粒径关系
Figure 9. Relationship between characteristic strength and particle size
![]()
图 11 S2椭球状颗粒的破碎强度Weibull分布
Figure 11. Weibull distribution of crushing strength for ellipsoidal particles (S2)
![]()
图 12 S3扁平状颗粒的破碎强度Weibull分布
Figure 12. Weibull distribution of crushing strength for platy particles (S3)
表 1 不同尺寸颗粒破碎强度的Weibull参数
Table 1 Weibull parameters for particles with different sizes
等效直径D/mm 子颗粒数 特征强度σ0/MPa Weibull模量m 27.5 75 10.50 4.20 31.6 114 10.00 4.77 35.7 165 8.71 6.04 39.7 229 8.04 7.31 43.9 307 7.92 8.49 
下载: 导出CSV
表 2 不同形状颗粒破碎强度的Weibull参数
Table 2 Weibull parameters for particles with different shapes
形状 S1 S2 S3 {\sigma _0} /MPa m {\sigma _0} /MPa m {\sigma _0} /MPa m 长轴 10.5 4.2 8.86 2.35 6.53 3.66 中轴 10.5 4.2 9.18 5.15 8.32 5.26 短轴 10.5 4.2 9.18 5.15 10.6 3.74 平均值 10.5 4.2 9.07 4.2 8.48 4.2
下载: 导出CSV
-
[1] 徐琨, 周伟, 马刚. 颗粒破碎对堆石料填充特性缩尺效应的影响研究[J]. 岩土工程学报, 2020, 42(6): 1013-1022. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202006006.htm XU Kun, ZHOU Wei, MA Gang. Influence of particle breakage on scale effect of filling characteristics of rockfill material[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(6): 1013-1022. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202006006.htm
[2] MCDOWELL G R. On the yielding and plastic compression of sand[J]. Soils and Foundations, 2002, 42(1): 139-145. doi: 10.3208/sandf.42.139
[3] HUANG J, XU S, YI H, et al. Size effect on the compression breakage strengths of glass particles[J]. Powder Technology, 2014, 268: 86-94. doi: 10.1016/j.powtec.2014.08.037
[4] DE BONO J P, MCDOWELL G R. An insight into the yielding and normal compression of sand with irregularly-shaped particles using DEM[J]. Powder Technology, 2015, 271: 270-277. doi: 10.1016/j.powtec.2014.11.013
[5] LAUFER I. Grain crushing and high-pressure oedometer tests simulated with the discrete element method[J]. Granular Matter, 2015, 17(3): 389-412. doi: 10.1007/s10035-015-0559-z
[6] CHENG Y P, NAKATA Y, BOLTON M D. Discrete element simulation of crushable soil[J]. Géotechnique, 2003, 53(7): 633-641. doi: 10.1680/geot.2003.53.7.633
[7] JEAN M. The non-smooth contact dynamics method[J]. Computer Methods in Applied Mechanics and Engineering, 1999, 177(3/4): 235-257.
[8] CAMACHO G T, ORTIZ M. Computational modelling of impact damage in brittle materials[J]. International Journal of Solids and Structures, 1996, 33(20/21/22): 2899-2938.
[9] JIANG H X, MENG D G. 3D numerical modelling of rock fracture with a hybrid finite and cohesive element method[J]. Engineering Fracture Mechanics, 2018, 199: 280-293. doi: 10.1016/j.engfracmech.2018.05.037
[10] 喻勇, 尹健民. 三峡花岗岩在不同加载方式下的能耗特征[J]. 岩石力学与工程学报, 2004, 23(2): 205-208. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200402004.htm YU Yong, YIN Jianmin. Energy dissipation properties of Three Gorges granite under different loading modes[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(2): 205-208. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200402004.htm
[11] QUEY R, DAWSON P R, BARBE F. Large-scale 3D random polycrystals for the finite element method: generation, meshing and remeshing[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(17/18/19/20): 1729-1745.
[12] SAADAT M, TAHERI A. A cohesive grain based model to simulate shear behaviour of rock joints with asperity damage in polycrystalline rock[J]. Computers and Geotechnics, 2020, 117: 103254. doi: 10.1016/j.compgeo.2019.103254
[13] CANTOR D, AZÉMA E, SORNAY P, et al. Three-dimensional bonded-cell model for grain fragmentation[J]. Computational Particle Mechanics, 2017, 4(4): 441-450. doi: 10.1007/s40571-016-0129-0
[14] 周博, 黄润秋, 汪华斌, 等. 基于离散元法的砂土破碎演化规律研究[J]. 岩土力学, 2014, 35(9): 2709-2716. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201409039.htm ZHOU Bo, HUANG Runqiu, WANG Huabin, et al. Study of evolution of sand crushability based on discrete elements method[J]. Rock and Soil Mechanics, 2014, 35(9): 2709-2716. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201409039.htm
[15] MCDOWELL G R, AMON A. The application of weibull statistics to the fracture of soil particles[J]. Soils and Foundations, 2000, 40(5): 133-141. doi: 10.3208/sandf.40.5_133
[16] FU R, HU X L, ZHOU B. Discrete element modeling of crushable sands considering realistic particle shape effect[J]. Computers and Geotechnics, 2017, 91: 179-191. doi: 10.1016/j.compgeo.2017.07.016
[17] HUILLCA Y, SILVA M, OVALLE C, et al. Modelling size effect on rock aggregates strength using a DEM bonded-cell model[J]. Acta Geotechnica, 2021, 16(3): 699-709. doi: 10.1007/s11440-020-01054-z
[18] MCDOWELL G R. Statistics of soil particle strength[J]. Géotechnique, 2001, 51(10): 897-900. doi: 10.1680/geot.2001.51.10.897
[19] XIAO Y, MENG M Q, DAOUADJI A, et al. Effects of particle size on crushing and deformation behaviors of rockfill materials[J]. Geoscience Frontiers, 2020, 11(2): 375-388. doi: 10.1016/j.gsf.2018.10.010
[20] GRIFFITH A A. The phenomena of rupture and flow in solids[J]. Philosophical Transactions of the Royal Society A, 1921, 221: 163–198.
[21] 孙壮壮, 马刚, 周伟, 等. 颗粒形状对堆石颗粒破碎强度尺寸效应的影响[J]. 岩土力学, 2021, 42(2): 430-438. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202102015.htm SUN Zhuangzhuang, MA Gang, ZHOU Wei, et al. Influence of particle shape on size effect of crushing strength of rockfill particles[J]. Rock and Soil Mechanics, 2021, 42(2): 430-438. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202102015.htm
[22] ZHU F, ZHAO J D. Interplays between particle shape and particle breakage in confined continuous crushing of granular media[J]. Powder Technology, 2021, 378: 455-467. doi: 10.1016/j.powtec.2020.10.020
[23] WANG Y H, MA G, MEI J Z, et al. Machine learning reveals the influences of grain morphology on grain crushing strength[J]. Acta Geotechnica, 2021, 16(11): 3617-3630.
[24] MA G, ZHOU W, REGUEIRO R A, et al. Modeling the fragmentation of rock grains using computed tomography and combined FDEM[J]. Powder Technology, 2017, 308: 388-397.
[25] ZHU F, ZHAO J D. Multiscale modeling of continuous crushing of granular media: the role of grain microstructure[J]. Computational Particle Mechanics, 2021, 8(5): 1089-1101.
[26] XIAO Y, DESAI C S, LIU H L. Testing and modeling on particle breakage for granular soils[J]. International Journal of Geomechanics, 2021, 21(11): 02021001. http://www.xueshufan.com/publication/3193815105
-
期刊类型引用(0)
其他类型引用(2)
计量
- 文章访问数: 292
- HTML全文浏览量: 52
- PDF下载量: 112
- 被引次数: 2
