Cosserat-particle finite element method for large deformation analysis of rock and soil
-
摘要: 粒子有限元法（PFEM）既继承了有限元法坚实的数学基础，又具有模拟大变形、复杂边界问题的能力，在流固耦合、岩土工程领域有广泛的应用。但另一方面，岩土体在大变形过程中往往具有应变软化特性和应变局部化现象，为保持问题求解的适定性，需要在本构方程中引入正则化机制，采用Cosserat连续体理论是引入正则化机制有效方法之一。将PFEM计算方法与Cosserat连续体理论结合，发展了Cosserat-PFEM方法。与传统PFEM中使用三角形单元不同，提出的新方法将边界识别、网格划分相互独立进行，使得四边形等单元的使用成为可能，以提高数值求解的精度与克服三角形单元对应变局部化问题模拟的倾向性。算例表明，本文发展的Cosserat-PFEM方法及基于ABAQUS软件开发的程序是可靠和有效的，拓展了PFEM的应用范围，具有模拟大变形问题并保持问题适定性的能力，适用于大变形渐进破坏问题的模拟。
-
关键词:
- 粒子有限元 /
- Cosserat连续体 /
- 边界识别 /
- 大变形 /
- 四边形单元
Abstract: The particle finite element method (PFEM) inherits the solid mathematical foundation of the finite element method and possesses the capability of modeling the problems with large deformation and complex boundary, so it has been widely used in the fields of fluid-structure coupling and geotechnical engineering. While in the process of large deformation of rock and soil, strain softening and strain localization often occur. In order to keep the well-posedness for the large deformation problem, it is necessary to introduce the regularization mechanism into the constitutive equation. The Cosserat continuum theory is one of the effective methods to introduce the regularization mechanism. Combining the PFEM with the Cosserat continuum theory, the Cosserat-PFEM method is developed. Besides, unlike the traditional PFEM using triangular elements, in the Cosserat-PFEM method the boundary recognition and the mesh generation are carried out independently, which makes it possible to use quadrilateral elements, so as to improve the numerical accuracy and overcome the tendency of triangular elements to simulate the strain localization. The examples show that the developed Cosserat-PFEM method and the programme explored based on the ABAQUS software are reliable and efficient, and expand the application scope of the PFEM. It is also demonstrated that the Cosserat-PFEM has the capability to simulate large deformation problems and keep the well-posedness of the problems, and is suitable for the simulation of large deformation and progressive failure problems. -
0. 引言
堆石料是土石坝的主要筑坝材料，一般占坝体工程70%~90%[1]。堆石料的力学性质将直接影响坝体的安全运行，随着筑坝技术的提高以及调节性大水库建设的需求，中国已建和在建的高坝数量日益增多，例如双江口、两河口、如美、古水等水利工程，最大坝高可达300 m。坝高的提升首先引起堆石料所受应力的增大，导致长期变形中颗粒破碎现象明显并影响堆石料的强度、渗透、变形等性质[2-5]。
针对堆石料的颗粒破碎现象，国内外学者进行了大量研究[6-8]。颗粒破碎指标的选择是研究堆石料颗粒破碎规律的基础，Marsal[6]将各个粒组试验前后的质量百分数的差值，取绝对值求和定义了破碎指标Bg。Hardin[7]根据初始级配与理想状态下颗粒完全破碎时级配所围成的面积作为破碎潜能，定义了破碎指标Br。Einav[8]认为颗粒破碎存在极限级配，将极限级配与初始级配所围成的面积作为破碎潜能，定义了破碎指标BE。郭万里[9]考虑到颗粒破碎的极限级配不易获取，将初始级配的级配面积作为破碎潜能，定义了破碎指标Bw。
对于堆石料在三轴剪切过程中颗粒破碎演化规律的研究，现存的研究方法主要有两种：①构建考虑颗粒破碎的本构模型。此类模型大部分是用孔隙比e的变化来反映颗粒破碎情况[10-11]，并构建包含孔隙比e的状态参数，将新的状态参数与弹塑性模型中的参数进行关联。此类模型虽然在建模过程中考虑了颗粒破碎的影响，但是对于破碎演化仅仅体现在孔隙比e的变化，不能实现堆石料级配的演化。②建立反映颗粒破碎率随应力和应变变化的数学模型。基于不同的颗粒破碎指标，众多学者研究了颗粒破碎率在不同力学状态下的变化规律[12-14]。蔡正银等[15]研究了不同级配、不同密度、不同围压条件下堆石料的颗粒破碎现象，发现破碎颗粒主要集中粒径在20 mm以上的颗粒。刘汉龙等[16]发现颗粒破碎率随围压的增大而增大，二者近似满足双曲线关系。高玉峰等[17]对不同岩性的堆石料展开大型三轴剪切试验，发现颗粒破碎率与围压呈线性增长关系。Salim等[18]基于三轴试验，建立了破碎指标Bg与剪应变和平均正应力的关系。Jia等[19]认为颗粒破碎不仅发生在三轴剪切过程中，试样振动压实放入三轴仪同样会产生明显的颗粒破碎，并基于塑性功建立了三轴剪切作用下破碎指标Br的变化规律。
以上的研究主要针对不同力学状态下破碎指标的变化，而就堆石料本身而言，颗粒破碎引起级配分布的变化。很少有涉及级配与剪应变、平均正应力的变化规律。朱俊高等[20]、郭万里等[21]基于一个双参数的级配方程，提出了实现“应力应变—破碎指标—级配分布（SBG）”转换的模型。其模型理论清晰，但是在实现“破碎指标—级配分布”转换时，双级配参数造成计算量大不易求解。本文在其基础上改进，将分形分布用于堆石料颗粒破碎的研究，建立一个便于运算的“应力应变—破碎指标—级配分布（SBG）”模型。
1. 破碎指标与级配分布的关系
1.1 级配的定量表示方法
关于堆石料颗粒破碎的研究已经取得众多成果，但大多数研究是针对“应力应变—破碎指标”的关系，通过引入适当的破碎指标来定量表示不同力学状态下颗粒破碎的程度。而对于堆石料整体而言，其最直观的物理变化就是级配的改变。如果可以建立“应力应变—级配分布”的转换模型，就可以从本质上反映颗粒破碎作用下堆石料性能的变化。
模型的建立首先需要解决级配分布的定量表示，大量学者根据级配的分布特性建立了众多的级配方程，其中以Talbot等[22]提出的分形级配方程和朱俊高等[20]提出的双参数级配方程应用最为广泛。相比之下，双参数级配方程由于受两个级配参数共同控制，在进行级配拟合时，其精度略高于分形级配方程。
郭万里通过引入双参数级配方程，构建了一个“应力应变—破碎指标—级配分布”的转换模型，其模型在理论上是可以应用的。然而，该模型在进行“破碎指标—级配分布”转化时，需要构建两个破碎指标从而建立两个方程来求解级配参数。无论是破碎指标Bw还是Bg，与级配参数的转换公式都是十分复杂的，方程组不便求解，而造成求解问题出现的主要原因是级配方程形式复杂且存在两个参数。
考虑到分形级配方程形式简单、参数单一，且在研究土体级配分布中得到广泛应用。本文将基于分形级配方程，并结合已有三轴剪切下堆石料的颗粒破碎成果，构建一个便于应用的“应力应变—破碎指标—级配分布”模型。
基于分形理论，Talbot等提出的分形级配方程[22]：
P=(ddmax)3−D×100%, (1) 式中D为分形维数，P为小于某粒径d的累计含量，dmax为最大粒径。
为了说明用分形级配方程定量表示堆石料级配时，具有充足的准确性，本文对古水、马吉、如美、茨哈峡堆石料分别进行了拟合，如图 1所示。4种堆石料的级配分布与拟合曲线吻合良好，仅仅是茨哈峡堆石料出现略微偏差，拟合参数见表 1。茨哈峡堆石料的相关系数的平方R2为0.98，其他堆石料为0.99。
表 1 不同堆石料的拟合参数Table 1. Fitting parameters of different rockfills土体类型 dmax/mm D R2 古水堆石料 800 2.46 0.99 马吉堆石料 700 2.59 0.99 如美堆石料 600 2.60 0.99 茨哈峡堆石料 300 2.67 0.98 虽然以上的拟合效果良好，但有学者发现分形级配方程仅仅适合拟合双曲线型级配分布，对反S型级配分布不适用[14, 20]。基于此，本文对两种形状的级配分布进行拟合，如图 2所示。分形级配方程对反S型级配分布的拟合效果较差。然而，笔者对实际工程中的堆石料级配进行筛查后，发现并没有反S型级配分布。主要原因是该形状的级配存在明显的粒组缺失，而实际工程中的堆石料为了获得高密实度，设计时不会出现粒组缺失现象。反S型级配分布仅仅出现在反滤料、砂土、黏土的设计中，故本文无需考虑拟合出现偏差的问题。
1.2 破碎指标的选择
颗粒破碎指标的选择是研究堆石料颗粒破碎的基础，意义明确的破碎指标可以很清晰地反映颗粒的破碎程度。本节对常用的破碎指标进行对比总结，并确定新模型所涉及的破碎指标。
常见的破碎指标基本可以分成2类：第1类为特定粒径（d10，d15等）和组合系数（Cu，Cc），该类指标形式简单、便于计算，但所反映的破碎程度过于片面，难免以偏概全。第2类用各个粒组质量分数的总体变化衡量破碎程度。Marsal[6]将各个粒组颗粒破碎前后质量分数的差值（ΔWk）取绝对值求和，定义了破碎指标Bg：
Bg = ∑|ΔWk|。 (2) Hardin[7]根据颗粒破碎后级配的变化，认为所有颗粒均具有破碎成微小颗粒的潜能，即级配曲线与P=100%完全重合，将初始级配与P=100%所围成的面积作为破碎潜能Sp，用试验后级配与初始级配所围成的面积作为颗粒破碎量S1，如图 3所示，将S1与Sp的比值作为破碎指标Br。Einav[8]认为颗粒破碎存在极限级配，将极限级配与初始级配所围成的面积作为破碎潜能，如图 4所示，定义了破碎指标BE：
BE=S1/Sp。 (3) 理论上，破碎指标Br与BE的范围为0~1，然而Br不可能取到1甚至远小于1，其反映的破碎程度仅仅在一个较小的区间内。相比之下，破碎指标BE可以取到0和1，指标值为0时表示颗粒未发生破碎，指标值为1时表示颗粒破碎达到最大程度。对于极限级配的确定，已有研究表明[23-24]，就实际工程中常规堆石料的级配范围而言，极限级配的分形维数趋近于2.7。
1.3 BE与分形维数的关系
根据破碎指标BE的定义，只要可以计算不同状态下级配在P–lgd坐标下的面积，就可以实现直接求解。由于粒径d取不到0 mm，故对P在区间[k, 1]进行积分，然后将k趋近于0得到极限值，如图 5所示。
级配曲线在P–lgd坐标下的面积：
S=limk→0SA=∫1k(lgdmax−lgdk)dP, (4) 式中，dk为P=k时所对应的粒径。
对分形级配方程进行变形，可得
lgdmax−lgdk=−13−DlgP。 (5) 将式（5）代入式（4），可得
S=1ln10⋅13−D。 (6) 将式（6）代入式（3），可得
BE=(Dc−D0)(3−Du)(Du−D0)(3−Dc)。 (7) 式中Dc为当前级配的分形维数，D0为初始级配的分形维数，Du为极限级配的分形维数，其中D0，Du为已知。故，破碎指标BE与当前级配的分形维数Dc之间一一对应的转换关系建立。
2. 破碎指标与应力应变的关系
Jia等[19]通过对古水堆石料进行三轴固结排水剪切试验，研究了堆石料在剪切过程中的颗粒破碎规律，发现颗粒破碎率同时受剪应变εs和平均正应力p控制，但没有提出相应的预测模型。本节基于文献[19]的试验数据，构建一个用于预测不同力学状态下破碎指标的模型。试验数据见表 2（此处的破碎指标BE与文献[14]有所差别，是采用不同级配方程求取级配面积造成的结果，总体而言趋势相同）。
围压/kPa 剪应变εs/% 平均正应力p/kPa D BE/% 0 0 0 2.2742 0 500 1.87 1242 2.2777 0.34 500 4.56 1422 2.3070 3.34 500 8.44 1494 2.3416 7.21 500 12.49 1469 2.3378 6.77 500 16.01 1420 2.3418 7.24 1000 1.75 2046 2.2891 1.48 1000 4.60 2343 2.3294 5.80 1000 7.20 2552 2.3657 10.16 1000 10.66 2534 2.3490 8.09 1000 15.19 2621 2.4026 15.14 1500 1.90 2760 2.3000 2.60 1500 7.34 3504 2.3913 13.55 1500 9.31 3596 2.3929 13.78 1500 12.47 3543 2.4180 17.40 1500 14.35 3602 2.4325 19.66 2000 1.85 3336 2.3009 2.70 2000 4.55 4064 2.3549 8.81 2000 7.26 4484 2.3901 13.39 2000 10.75 4653 2.4235 18.25 2000 13.76 4619 2.4386 20.63 通过分析剪应变和平均正应力与BE的关系，得到一个可以定量表示BE随剪应变和平均正应力变化的公式：
BE=a[1−exp(−bεs)]arctan(cp/pa)。 (8) 式中，3个参数可以直接通过拟合获得。另外，在剪应力–剪应变试验曲线完整的情况下，参数b可以通过临界状态直接确定。当土体受剪达到临界状态时，平均正应力p保持不变而剪应变εs持续增大，颗粒破碎量将趋近于恒定值。根据式（8）可知，要想破碎指标BE趋近于恒定值，[1-exp(-bεs)]需要无限接近于1。此处假设[1-exp(-bεs)]≥0.99时，土体达到临界状态。此时需满足bεs≥4.6，换句话说，用εcs表示临界剪应变：
b=4.6/εsc。 (9) 受试验仪器的限制，三轴剪切下的剪应变常常达不到临界值。根据已有研究[25-26]，可以按下式计算临界状态的应力比MC及相应的临界剪应变：
η=ACε2s+εsε2s+1,Mc=AC。} (10) 式中，A，C为拟合参数。
本文提取了Jia等[19]的三轴固结排水剪切试验的剪应力–剪应变曲线，将应力比η和剪应变ε绘制成散点；并根据式（10）进行拟合，如图 6所示。剪应变为18%时，不同围压下的应力比η趋近于水平线，且临界应力比MC较为接近，说明式（10）的正确性。临界剪应变εcs=18%，根据式（9），参数b=25.6。
将b=25.6代入式（8），并对表 2中的数据进行拟合，如图 7所示。图 7中，a=0.3257，c=0.0149，相关系数的平方R2达到了0.90，表明模型具有较高的预测精度。
以上是对Jia等[19]的数据进行建模和验证，仅仅是对1组数据进行拟合不代表模型广泛的适用性，此处又引用了郭万里[27]的试验数据，见表 3。利用式（8）对三轴剪切试验的数据直接进行拟合，如图 8所示。拟合参数a=0.3284，b=24.6，c=0.0727，相关系数的平方R2为0.95，进一步论证了模型的适用性。
围压/kPa 剪应变εs/% 平均正应力p/kPa BE/% 0 0 0 0 200 4.41 570 10.6 200 9.95 580 14.4 200 13.40 560 15.6 500 3.37 1140 11.6 500 7.53 1300 18.7 500 14.30 1320 23.1 1000 2.60 1800 12.2 1000 8.56 2370 26.0 1000 15.20 2480 32.1 1500 1.63 2160 13.2 1500 5.34 3010 30.2 1500 16.30 3570 41.4 3. SBG模型的验证
以上分别研究了“破碎指标—级配分布”、“应力应变–破碎指标”的关系，本节将对以上转换关系联立，对“应力应变—破碎指标—级配分布”模型进行验证。将表 2中的数据通过式（8）进行拟合，得到不同力学状态下预测的破碎率BE，再根据式（7）计算相应的分形维数。对围压为2000 kPa，剪应变为7.26%，13.76%的级配进行拟合，结果如图 9所示，预测级配与试验后级配拟合良好。为了说明模型对所有数据的总体预测精度良好，将所有20组实测级配的分形维数和模型预测得到的分形维数绘制于图 10，可以看出实测值与拟合值的最大误差仅仅在0.02左右。故，本文建立的“应力应变—破碎指标—级配分布”模型拥有良好的适用性。
4. 结论
本文基于分形级配方程，研究了三轴剪切过程中颗粒破碎率和级配变化的一般规律，建立了相应的“应力应变—破碎指标—级配分布（SBG）”转换模型，模型成功预测了不同力学状态下级配的变化。
（1）用分形级配方程定量表示堆石料级配分布具有较高的准确性。虽然分形级配方程对反S型级配分布不适用，但反S型级配分布存在明显的粒组缺失，实际堆石料设计时不存在该形状的级配分布。
（2）对比了常见的破碎指标，认为Einav提出的破碎指标BE最适合表示颗粒破碎率，BE的取值范围为0~1，指标值取0时表示未发生颗粒破碎，指标值取1时表示颗粒破碎达到最大程度。
（3）建立了一个可以定量表示破碎指标随剪应变和平均正应力变化的模型，模型共有3个参数a，b，c，其中参数b与土体的临界状态有关，由临界剪应变可以直接求得。对2组不同的三轴剪切试验数据进行拟合，相关系数的平方均在0.9以上。
（4）本文建立的SBG模型成功预测了不同力学状态下堆石料级配的变化。其使用方法简洁：首先通过式（8）对试验数据进行拟合，获得“应力应变—破碎指标”的转换；然而按式（7）实现“破碎指标—级配分布”的转变。
-
![]()
图 7 不同单元类型对应承载力曲线图
Figure 7. Corresponding bearing capacity curves of different element types
![]()
图 10 Cosserat-PFEM计算终止时等效塑性应变图
Figure 10. Diagram of equivalent plastic strain at the end of calculation by Cosserat-PFEM
表 1 材料参数表
Table 1 Parameters of materials
弹性模量
E/kPa泊松比
ν黏聚力
c/kPa内摩擦角
φ/(°)密度
ρ/(kg⋅m-3)30000 0.35 0.1 19.0 1490 
下载: 导出CSV
表 3 材料参数表
Table 3 Parameters of soil
材料参数 值 材料参数 值 弹性模量E/kPa 50000 内摩擦角φ/(°) 25 泊松比ν 0.30 剪胀角ψ/(°) 0 初始黏聚力c0/
kPa50 Cosserat剪切模量Gc/
kPa10000 黏聚力软化模量hp/
kPa-30 内部长度参数lc/
m0.06
下载: 导出CSV
-
[1] HIRT C W, AMSDEN A A, COOK J L. An arbitrary Lagrangian-eulerian computing method for all flow speeds[J]. Journal of Computational Physics, 1997, 135(2): 203-216. doi: 10.1006/jcph.1997.5702
[2] HU Y, RANDOLPH M F. A practical numerical approach for large deformation problems in soil[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1998, 22(5): 327-350. doi: 10.1002/(SICI)1096-9853(199805)22:5<327::AID-NAG920>3.0.CO;2-X
[3] YU L, HU Y X, LIU J, et al. Numerical study of spudcan penetration in loose sand overlying clay[J]. Computers and Geotechnics, 2012, 46: 1-12. doi: 10.1016/j.compgeo.2012.05.012
[4] NOH B W F. CEL: A time dependent two space-dimensional, coupled Eulerian Lagrangian code[C]// Methods of Computational Physics, New York, 1964.
[5] PUCKER T, GRABE J. Numerical simulation of the installation process of full displacement piles[J]. Computers and Geotechnics, 2012, 45: 93-106. doi: 10.1016/j.compgeo.2012.05.006
[6] LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. The Astronomical Journal, 1977, 82: 1013-1024. doi: 10.1086/112164
[7] 廉艳平, 张帆, 刘岩, 等. 物质点法的理论和应用[J]. 力学进展, 2013, 43(2): 237-264. https://www.cnki.com.cn/Article/CJFDTOTAL-LXJZ201302003.htm LIAN Yanping, ZHANG Fan, LIU Yan, et al. Material point method and its applications[J]. Advances in Mechanics, 2013, 43(2): 237-264. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXJZ201302003.htm
[8] SULSKY D, CHEN Z, SCHREYER H L. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1/2): 179-196.
[9] 孙玉进, 宋二祥. 大位移滑坡形态的物质点法模拟[J]. 岩土工程学报, 2015, 37(7): 1218-1225. doi: 10.11779/CJGE201507007 SUN Yujin, SONG Erxiang. Simulation of large-displacement landslide by material point method[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(7): 1218-1225. (in Chinese) doi: 10.11779/CJGE201507007
[10] 高宇新, 朱鸿鹄, 张春新, 等. 砂土中锚板上拔三维物质点法模拟研究[J]. 岩土工程学报, 2022, 44(2): 295-304. doi: 10.11779/CJGE202202011 GAO Yuxin, ZHU Honghu, ZHANG Chunxin, et al. Three- dimensional uplift simulation of anchor plates in sand using material point method[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(2): 295-304. (in Chinese) doi: 10.11779/CJGE202202011
[11] OÑATE E, IDELSOHN S R, DEL PIN F, et al. The particle finite element method—an overview[J]. International Journal of Computational Methods, 2004, 1(2): 267-307. doi: 10.1142/S0219876204000204
[12] 张雪, 盛岱超. 一种模拟土体流动的连续体数值方法[J]. 岩土工程学报, 2016, 38(3): 562-569. doi: 10.11779/CJGE201603021 ZHANG Xue, SHENG Daichao. Continuum approach for modelling soil flow in geotechnical engineering[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(3): 562-569. (in Chinese) doi: 10.11779/CJGE201603021
[13] ZHANG X, SLOAN S W, OÑATE E. Dynamic modelling of retrogressive landslides with emphasis on the role of clay sensitivity[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42(15): 1806-1822. doi: 10.1002/nag.2815
[14] ZHANG X, OÑATE E, TORRES S A G, et al. A unified Lagrangian formulation for solid and fluid dynamics and its possibility for modelling submarine landslides and their consequences[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 343: 314-338. doi: 10.1016/j.cma.2018.07.043
[15] LI H C, ZHANG S. Implement the particle finite element method in ABAQUS[J]. Japanese Geotechnical Society Special Publication, 2020, 8(3): 70-75. doi: 10.3208/jgssp.v08.c16
[16] YUAN W H, WANG H C, ZHANG W, et al. Particle finite element method implementation for large deformation analysis using Abaqus[J]. Acta Geotechnica, 2021, 16(8): 2449-2462.
[17] TANG H X, GUAN Y H, ZHANG X, et al. Low-order mixed finite element analysis of progressive failure in pressure-dependent materials within the framework of the Cosserat continuum[J]. Engineering Computations, 2017, 34(2): 251-271.
[18] EDELSBRUNNER H, MÜCKE E P. Three-dimensional alpha shapes[J]. ACM Transactions on Graphics, 1994, 13(1): 43-72.
[19] HU Y, RANDOLPH M F. H-adaptive FE analysis of elasto-plastic non-homogeneous soil with large deformation[J]. Computers and Geotechnics, 1998, 23(1/2): 61-83.
[20] NAZEM M, CARTER J P, AIREY D W. Arbitrary Lagrangian-Eulerian method for dynamic analysis of geotechnical problems[J]. Computers and Geotechnics, 2009, 36(4): 549-557.
[21] TIAN Y H, CASSIDY M J, RANDOLPH M F, et al. A simple implementation of RITSS and its application in large deformation analysis[J]. Computers and Geotechnics, 2014, 56: 160-167.
-
其他相关附件
计量
- 文章访问数: 0
- HTML全文浏览量: 0
- PDF下载量: 0
