湖南省大学生研究性学习和创新性实验计划
项 目 申 报 表
如果您无法在线浏览此 PDF 文件，则可以 下载免费小巧的 福昕(Foxit) PDF 阅读器,安装后即可在线浏览 或 下载免费的 Adobe Reader PDF 阅读器,安装后即可在线浏览 或 下载此 PDF 文件
|
学校名称
|
欧博官网
|
学生姓名
|
学 号
|
专 业
|
性 别
|
入 学 年 份
|
王创业
|
201508021028
|
城市轨道
|
男
|
2015
|
周伟
|
201508021023
|
岩土工程
|
男
|
2015
|
王子建
|
201508020827
|
城市轨道
|
男
|
2015
|
罗跃龙
|
201508021016
|
城市轨道
|
男
|
2015
|
郭剑雄
|
201508020922
|
城市轨道
|
男
|
2015
|
指导教师
|
胡敏
|
职称
|
讲师
|
项目所属
一级学科
|
土木工程
|
项目科类(理科/文科)
|
理科
|
学生曾经参与科研的情况
1、参与了指导老师主持的欧博官网桥梁工程领域研究基地开放基金资助项目“基于分形结构的砂卵石土地层桥梁桩基承载力分析”的实验工作。
2、针对张桑高速有关膨胀岩的化学改良进行了调查研究与初步实验。
3、参与了侧限压缩下级配砂砾石料的强度特性及其颗粒破碎特征的相关研究。
4、部分队员参加中南地区结构设计竞赛并获特等奖。
5、参与国家自然科学面上基金“基于块体模型细胞自动机的节理岩质边坡稳定性和演化机理研究”的项目数值模拟。
6、在校期间曾参加过多次科技立项并获奖。
7、参与“互联网+”、“创青春”全国大学生创业大赛并获奖。
|
指导教师承担科研课题情况
1.欧博官网桥梁工程领域研究基地开放基金资助项目，15KC03，基于分形结构的砂卵石土地层桥梁桩基承载力分析，2015/12-2017/12，3万元，在研，主持
2. 国家自然科学基金面上项目，11072095，对碎石桩和堆石路基中应力链的拱效应及其崩塌破坏模式的研究，2011/01-2013/12，25万元，已结题，参与
3．欧博官网引进博士启动基金，侧限压缩下级配砂砾石料的强度特性及其颗粒破碎特征研究，2017/12-2019/12，6万元，在研，主持
|
项目研究和实验的目的、内容和要解决的主要问题
项目研究和实验的目的：
膨胀岩是一种似岩非岩、似土非土的特殊土体，又被称为膨胀土。其具有吸水易膨胀崩解的特性，在公路、边坡、房建工程具有极大的危害[1][2]。目前国内对膨胀岩的治理改良大多使用水泥石灰、二灰土、化学改性剂和物理加固的方法，成本较为昂贵，我们提出采用水泥黏土这一新的改良方法，同时做了初步的实验来论证这个想法的可行性。依托张桑高速公路施工开挖遇到的膨胀岩，现场采集样品后判别其膨胀特性，以此为背景从室内实验与数值模拟的角度来探究膨胀岩作为路基路堤填料的化学改良方案，室内实验将调整不同水泥黏土配合比并测定一定天数后的抗压强度，从满足一定改良强度、满足不同配合比下两个角度下来对改良方案进行比较；同时基于两款离散元软件对膨胀土改性前后的无侧限压缩实验进行模拟，尝试模拟颗粒遇水膨胀破碎崩解机理的过程，为膨胀岩的研究提供新的研究方法。
项目研究和实验的内容：
采用水泥和粘土对膨胀岩进行化学改良，调整水泥、粘土、膨胀岩按不同配合比装填物料入试样中，经过 7天、28天洒水养护后评价各个配合比下试件的无侧限抗压强度。
干湿循环实验反映了改良后击实膨胀岩下自然条件下的变化情况及强度特征，而初步的探究中尚未涉及，因此在进一步的实验中将进行多次洒水与自由风干来模拟干湿循环，并测其干湿循环稳定后的抗压强度。
对实验数据借由origin等专业绘图软件进行二次处理，得出不同配合比下
的改性无侧限抗压强度，确定一定抗压强度下的最优配合比，结合工程实际确定相应改良方案。
借助PFC与EDEM软件进行相应配合比下的实验过程的数值模拟，以此得出更为普遍的规律。
要解决的主要问题：
在实验方面：
1.寻找合适的改良材料组合，该组合可以使改良后的膨胀岩强度达到要求且作为路基使用保证不易受水的影响，同时需满足原料经济、有工程应用价值、施工便利。
2.研究改良材料之间的配比方案，找到最佳的配比方案，使得改良后的膨胀岩在满足强度要求的情况下经济节约。
3.试样添加水分的多少对试样强度有一定影响，因此需要定量探究含水率对形成试样后的强度影响并确定最佳含水率。
4.干湿循环强度是膨胀岩实际应用中的一项重要指标，而初步的探究中尚未涉及，因此我们还需探究改良后的膨胀岩的无侧限抗压强度与干湿循环稳定后的抗压强度。
在数值模拟方面：
1.如何利用PFC的clump功能来模拟颗粒的聚合，精确模拟复杂外形颗粒，设其接触模型为粘性接触模型。
2.模拟颗粒聚合团遇水破碎、软化、崩解、膨胀的现象。
3从颗粒三轴实验模拟程序出发，尝试对改良膨胀岩的聚合物进行无侧向抗压并测其模拟强度。目前我们学习开发了一些相应的程序代码来支撑我们的数值模拟研究，限于篇幅，我们在附件中给出了一部分代码。利用EDEM和PFC软件做改性后模型的无侧向抗压强度数值模拟，其结论做改良试样室内实验的对照，
为研究膨胀岩的软化崩解提供新的思路。
在数据处理方面：
初研的结果显示同一组试样的强度也有偏差，在实验允许的情况下探究含水率和养护条件对强度的影响。同时对数据曲线是否考虑因自然因素引起的强度波动，我们将结合进一步的实验确定级配-强度的曲线。
|
国内外研究现状和发展动态
一、膨胀岩的背景介绍
膨胀岩属于软岩，不同程度上含有亲水矿物，其含水率产生变化时体积发生剧烈改变，其应变、强度等力学特性与含水率、温度等有复杂的关系。
膨胀岩的治理是工程地质领域的一大难题，目前我们对膨胀岩的研究还较为初步。早于1956年Holtz 等[3]便开始了涉及膨胀土的工程研究，而在国内20世纪90年代起，由于其对基建的危害，引发了国内学者的广泛关注。而在南昆铁路工程，南水北调中线工程中均分布有膨胀岩，在一些西部地区如云南、四川等,膨胀岩分布分布面积较大且地质条件复杂,由于膨胀岩的研究不足不断出现相应的施工维护等问题，造成了人员伤亡和重大财产损失。
我们在张桑高速采集了膨胀岩，针对其做了定量的研究。
图1
图2
图3
如上图1图2图3是张桑路段膨胀岩对道路边坡的破坏现场。由此可见膨胀岩对公路边坡等工程施工的危害。
二、实验部分的研究进展
1．膨胀岩的分类判别研究现状与发展动态
膨胀岩易发生胀缩，胀缩是由于岩石吸水、脱水而引起的，这与膨胀岩的含水量密切相关。但不同地区的膨胀岩, 由于其粘土矿物成分与含量的差异, 其破坏机理与治理原则也有所不同。易发生吸水膨胀现象的岩石种类极为复杂, 如泥岩、页岩、粘土岩等。目前根据膨胀岩的膨胀机理将其分为化学膨胀和粘土质膨胀。化学膨胀是由于吸水或结晶体积变化等, 如硬石膏等，另一类是含有粘土矿物, 如泥岩等，其粘土类膨胀岩矿物含量越高, 膨胀率就越高。由于其分类与膨胀机理的复杂性，不少人对此做了一定的研究。
何山等[4]综合考虑影响膨胀的因素，采用数理概率统计方法对膨胀岩的膨胀判定给出了定量的表达方法，对不同影响参数进行了加权处理， 避免了单因素分类对膨胀岩进行分类的不准确性，但由于因素众多而不适合在工程中推广。
朱训国等[5]提出了一种以膨胀岩成分含量、吸水率等为主要控制指标的新标准，并将其应用于工程实际。
2.岩石浸水膨胀的应力-应变本构关系
由于膨胀岩组成成分与破坏机制的复杂性，不少学者尝试从其本构关系来研究其物理力学性能。1970 年Huder 和Amberg 提出Huder-Amberg 模型，他们认为其膨胀应变与应力呈半对数线性关系,而刘晓丽等[5]基于实验结果, 提出了修正的Huder-Amberg 膨胀本构模型，考虑到膨胀岩应变随时间变化的特性，给出了侧限膨胀应变随时间的变化关系，刘等还将岩石一维膨胀本构关系扩展到三维。
罗鸿禧[7]采用高真空冷冻升华干燥法, 对膨胀岩进行了因水份变化而引起的孔径分布变化规律的研究, 揭示胀岩膨胀过程的实质是结构单元内部的细孔和微孔不断吸水膨胀, 产生一定的膨胀力, 引起总体结构重新调整的过程。
朱珍德[8-10]借助南京水利枢纽工程所提供的岩样研究膨胀变形及力学特性，实验结果表明红砂岩膨胀应变随吸水率增加而呈对数增长；借助岩伺服机和应变仪进行实验研究，在Zhong Wanxie等[11]的基础上提出考虑膨胀岩膨胀及软化特性的弹塑性本构模型，得出不同荷载下膨胀特性随吸水率的变化关系 ，并进一步提出用图像处理技术定量描述红砂岩细观结构状态及裂纹扩展变化过程。
非膨胀岩水理性质较稳定，可用于路段填料，如于春江等[12]对实际工程中大量的软岩、泥岩的工程特性进行分析进行研究，发现相应路段的泥岩等非膨胀岩可作为路基填料。付宏渊等[13]通过对广西六寨至宜州高速公路炭质页岩作为路基填料路用性能研究表明：并不是所有的炭质页岩都能用于路堤填料，它取决于炭质页岩的物理力学性能，因此膨胀岩常常需要进行改良改性处理才能作为填料。
流变特性是岩石的重要力学性质之一，郑西贵等[14]在西原体流变模型基础上，考虑湿度效应，建立了膨胀岩在应力和湿度耦合作用下的流变本构关系。
参考Genuchten模型[15]模型用土水特征曲线基础上赵二平等[16]结合工程实际开展了循环加卸荷实验并尝试用指数关系来刻画应变增量随循环次数的变化规律。
张巍等[17]则对崩解物粒径和岩样膨胀性关系展开研究，他提出膨胀岩膨胀性与其崩解物的最大含量粒组颗粒粒径、有效粒径呈负相关，同时针对南水北调中线河南安阳段的膨胀岩设计了饱和与非饱和渗透特性实验．发现膨胀岩的饱和渗透系数、试样下渗速度与含水量的关系。
何建新等[18]对新疆某地膨胀岩重塑土在不同含水率干密度，以及干湿循环状态下的强度指标通过直剪实验进行研究，得出膨胀岩的强度指标随着含水率的变化趋势。
益圣强[19]结合实践对击实条件下膨胀岩含水率与膨胀率进行了探究，得出了膨胀率与起始含水率关系曲线，膨胀率与掺入石灰率、击实功关系曲线；
李国富[20]通过实验分析预测巷道开外后膨胀岩的力学行为并基于实例建立膨胀岩的实验模型，得出了巷道深度与膨胀率等参数的关系。同时发现粘土矿物比表面积的变化是粘土矿物强吸水特征的主要原因，膨胀岩粘土矿物颗粒与水相互作用形成的扩散层厚度及塑性变形是岩石膨胀变形的主要原因.
Lin[21]对泥岩不同含水量下的晶体、形态等做了研究，基于此杨成祥[22]借助CT 技术观察泥岩动态软化并提出了水－岩作用过程中泥岩软化演化的三个阶段。
Heggheim[22]研究了不同介质溶液下水岩的相互作用，杨春和等[24]通过电子扫描电镜等电子设备研究了板岩遇水软化的过程与机理。
车平[25]从矿物含量等方面，对粉砂质泥岩软化实验研究，实验表明水、可溶盐胶结、裂隙综合作用导致岩石软化。
胡博聆[26]依托滇中路段场地探究不同风化程度下泥质粉砂岩的崩解实验。实验结果表明: 岩石风化程度与崩解程度呈正相关，泥质粉砂岩的崩解是由浸水过程、胶结物溶解、干湿循环共同作用引起的。
谭罗荣[27]从力学强度特性的角度讨论了粘土岩泥岩的泥化崩解问题
综上所述，粘土岩类泥化崩解，最初是宏观的结构破坏，然后是岩体吸水矿物水化，这使得岩体不均匀膨胀又加剧岩体的拉裂，粘土矿物进一步水化，逐步使受损软岩产生泥化物，最终形成形成高含水量的泥化物。在整个泥化过程中，实际上是岩体结构由宏观破坏引起微观破坏并相互影响的过程，最终粘土矿物水花泥化崩解强度急剧降低。
杨建林[28]利用粉晶X 射线衍射等设备对泥岩进行测试，研究表明：泥岩的崩解主要由吸水产生的次生孔隙及表面拉应力控制。
黄宏伟等[29]通过扫描电镜等对煤系地层泥岩进行分析,研究了该泥岩软化时微观结构随时间变化的特征。
张建智等[30]认为渗水膨胀岩隧洞变形机制是一个受渗流、膨胀等多因素影响的复杂时变问题，探讨了其变形机制，建立了渗水膨胀围岩受力平衡方程。
卢义玉等[31]针对煤矿穿膨胀岩的缩径问题，在考虑湿度影响的基础上结合ANSYS 和模型实验研究分析了湿度与地应力耦合对钻孔缩径变形的影响。
杨军平等[32]通过控制水温观察膨胀岩在水中崩解模式，探究了温度对膨胀岩崩解软化的影响规律。他发现: 随着温度的升高，粘土矿物膨胀作用增强，崩解速率也渐渐增加。
吴珺华[33]采用压缩仪和收缩仪，对重塑膨胀土反复进行膨胀和收缩，探究了在干湿循环作用下膨胀土的胀缩变形特征。
Terzaghi[34]，Gamble[34]从水分与孔隙的角度探究了红砂岩的崩解机制，而刘晓明等[34]则从能量变化的角度建立了基于能量耗散原理建立红砂岩崩解的能量耗散模型。研究结果表明：红砂岩崩解能量利用率随崩解呈指数衰减。
3.处理方法
为治理膨胀岩这种特殊土体，国内外不少学者做了较多的研究，
罗勇[37]在对南水北调南阳实验段中，水泥改性土的应用做了探究。他对于膨胀土的研究不仅使得实际工程中避免了特殊土体的运走，也节省了路基填料的费用。
目前对于膨胀岩的处理，运用得最广泛的是石灰与水泥，粉煤灰这三种材料。这三种材料廉价易得且性质稳定，适合施工处理。不仅实现了资源的苒利用’ 还减少了对资源的开采，对环境起到了很大的保护作用。
杨和平[38]提出采用微生物技术对西部膨胀土地区公路进行土质改良，并从多个方面考虑了利用微生物技术对膨胀土进行改良,
王小军[39]尝试用生石灰改良膨胀岩，通过室内实验 他的研究得出掺入生石灰改良后膨胀岩改善了其胀缩特性。膨胀岩改良最佳生石灰掺入比在8 %～ 12 %范围内。
何爱华[40]采用高分子量有机阳离子聚丙烯酰胺（CPAM）对矿膨胀性凝灰岩进行改性，模拟工程实践自制模拟注浆仪，发现CPAM改性液抗崩解性能优良，且浓度越高，对于细小颗粒的吸附、孔隙的封堵作用越强。
邹开学[41]依托南昆铁路膨胀岩开展实验,采用日本粉体材料和二灰土、水泥( 加粉煤灰)改良膨胀岩, 为增强韧性添加了高分子材料。实验结果表明, 二灰土、水泥( 加粉煤灰)强度满足要求, 但需增加韧性;
徐晗等[42]过物理性质实验、击实实验、力学性质实验以及膨胀特性实验探究了新乡潞王坟段岩层的压实度，初始含水率与其工程特性之间的关系,
王建军[43]针对路堑边坡工程中常遇的膨胀岩问题，探讨了其破坏机理，总结了设计施工的基本原则，提出膨胀岩边坡防治技术措施。
傅毅静[44]针对南钦铁路沿线风化软岩填料进行水泥改良土实验，结果表明: 强风化粉砂岩、粉砂质泥岩填料掺4 %水泥改良后满足铁路路堤填料要求，比远运合格填料方案更为合理。
于春江[12]对实际工程中大量的软岩、泥岩的工程特性进行分析进行研究，发现相应路段的泥岩可作为路基填料。
陈花林等[45]结合武广客运线对泥质粉砂岩化学改良进行实验。研究结论:泥质粉砂岩化学改良后力学强度及水理性得到较大提高,能满足铁路路基填料的要求。
毛忠良[46]介绍了 CMA 膨胀土生态改性剂在衡茶吉铁路工程中改性膨胀土的施工技术，且改性后的土体基本上消除了膨胀土的胀缩特性，土体的强度大幅度提高，水稳性好，土体稳定。
王浩[47]提出了一种用DAH改良剂与石灰混合溶液渗透改良的新方法。研究表明,改良膨胀土胀缩性降低,同时强度和水稳定性大幅提高。
杨国录等[47]提出了用HPZT 膨胀土改性剂固结膨胀土进行膨胀土渠道施工。该方法抑制了膨胀土的膨胀性能。
颜春[49]探究了树根桩与CMA 混合溶液共同作用下膨胀土的加固防护效果。
4.施工方向的具体措施研究
针对膨胀岩施工的过程，相关学者也做了一定的研究：
青兆麟[50]探究了泥岩用于石路堤填料施工技术
成海坡[51]的研究表明泥质软岩可以作为应用在填泥质粉灰岩的施工过程中， 但要注意压实次数、铺设厚度和压实度。
韩强等[52]提出采用一种化学材料（GCL）等做为渠基隔水材料，设置隔水层，以控制渠基膨胀岩自身稳定和不发生变形
物理改良：王福和等[53]依托贺家坪连接线K2+710～+890 段老路拓宽路堤的工程实例，从满足整体稳定性和均匀沉降的要求出发, 应用Geo- slope软件对土工格栅加固填石路堤进行了分析,结果表明: 土工格栅对填石材料提供了侧向压力，限制了路堤的侧向变形, 因此路堤的整体稳定性得到了较大提高;。
5.室内实验进展
国际岩石力学学会关于泥质膨胀岩室内实验的建议由四部分组成，对膨胀岩的取样和实验处理提供了指导。
惠会清等[54]通过石灰、粉煤灰混合料改良膨胀土实验, 得到了石灰、粉煤灰混合料在改良膨胀土中的最佳添加量;发现膨胀土的自由膨胀率随石灰量的增加而减小, 无侧限抗压强度随石灰量的增加而增大。
陈涛等[55]针对广西已建高速公路部分路段路面开裂等路基病害，从基本物理性质和水稳定性角度对风化泥岩进行了室内实验研究。实验结果表明: 该路段风化泥岩水稳定性较差，经过一定的干湿循环次数后，其强度已不能满足规范要求。
施建振等[55]对南京绕越高速公路风化泥岩作为路基填料的问题进行研究。发现该路段泥岩开挖料水稳定性较差，若利用其作为填筑材料，需要在外面用其他土料包裹以增强路基稳定性。
由于软岩填筑路基的广泛应用，有学者对软岩填筑路基可行性判定与应用做了定量的研究，郑明新[57]通过X-衍射分析矿物，室内实验结合实践提出了风化软岩软岩填筑路基可行性的初步判定方法。
刘新喜等[58]在探究风化软岩用于高等级公路路基填筑时得出压实特性变形特性是评价填料性能的重要指标，对填料在不同压实度下进行承载比实验，实验结果表明，承载比值随填料压实度的增大而增大。
周小顺[59]针对某高速公路施工中所遇到的大面积膨胀土的情况，借由以石灰为核心的改性土室内实验，发现在膨胀土中掺加掺入比为3%的石灰改性后膨胀土的物理力学性质有明显改善。
动态：目前对膨胀土的改良大多属于化学改良
化学改良原理：由土力学经典理论晶格扩张膨胀理论认为凡是具有膨胀晶格构造的粘土矿物都具有一定的膨胀性，引发膨胀岩胀缩变形的主要矿物成分蒙脱石具有两层硅氧四面体夹一层硅氧八面体组成的三层相叠的晶格构造，层间以氧原子形成公共弱键结合。这种结构具有膨胀晶格构造和晶格内离子交换的特性，容易吸附大量的极性水分子，从而使晶层间距增大产生体积膨胀，当有侧限约束时则会产生膨胀力。基于胶体化学原理的扩散双电层理论则认为粘土颗粒表面带负电荷，在静电场作用下必然会吸附带正电的可交换性阳离子来平衡电场力。这种交换性阳离子在水溶液中以水化离子的形式存在，与粘土颗粒表面的负电荷形成的双电层可以吸附极性水分子。蒙脱石颗粒表面双电层中的扩散层厚度较大，结合水膜变厚则会使水膜“模开”粘土颗粒的间距也增大。能提供大量有机阳离化合物的水溶液可减小塑性土样的塑性、收缩及膨胀性,
从上述土力学理论背景可知化学改良膨胀岩是行之有效的。
6.工程实例，对膨胀岩的新认识
龚壁卫等[59][60]指出：人们常常把干缩裂隙与原生裂隙混为一谈而忽视了原生裂隙，发现原生裂隙对滑坡也有重大影响。
范秋雁[61]采用原状膨胀岩进行室内边坡模型实验，揭示在连续降雨和湿-干循环模式下膨胀岩边坡的变形破坏模式。
程永辉等[62]为研究降雨条件下膨胀土边坡失稳的机理，研制相似室内模型，模拟了降雨条件下土边坡失稳破坏过程。研究表明膨胀土边坡失稳具有渐进性且属于边膨胀边失稳。
王军等[63]发现随着干湿气候不断循环, 裂隙逐渐发育,因此利用Matlab数字图像处理方法提出图像灰度化及二值化的数值处理技术,来定量描述膨胀岩裂隙度。
童学军等[64]依托武康二线工程开展膨胀土填料的室内石灰改良实验研究，研究结果表明：经石灰改良处理后膨胀土的胀缩性得到了有效抑制且浸水稳定性和强度都有了明显改善，满足铁路路基填料要求。同时发现石灰改良土强度随养护龄期的增加和压实系数的提高而增加，随着掺灰比的增加，石灰改良土的强度相应增大。
辛文栋等[64]结合兰新线膨胀性泥岩路段，开展了路堑边坡受膨胀性及风蚀的影响研究，提出了相应的防护措施。
二．数值模拟方面的进展
膨胀岩的本构关系的复杂性导致其模型建立较为困难，在工程实践中往往通过现场取样室内实验的方法来开展处理与进一步的研究，地应力和岩层成分的不同导致其特征物性参数的标定难题显得尤为突出，目前在膨胀岩的治理应用中没有较为成熟的数值模拟案例，但仍有一些学者从不同方向做了相关的研究来验证对照室内实验。
王凯等[66][66]为探究膨胀岩湿度应力场对围岩变形与次生应力状态的影响，结合工程实例，基于FLAC3D 软件进行本构模型的二次开发并建立了膨胀岩湿度应力场的本构模型，他发现岩体中水分扩散属于动态扩散且膨胀岩的体积膨胀具有流变性。
刘新喜等[58]利用饱和与非饱和渗流原理，将有限元软件用于边坡的极限平衡分析，探讨降雨强度、持时和压实度等外界参数对边坡稳定性的影响。数值模拟分析表明：强风化软岩用于路堤压实度不能低于90 %；当暴雨强度高于200 mm/d 时，边坡容易失稳。
离散单元法系cundall[67]于1971年最早提出后，广泛应用于岩土等离散介质总，在借鉴Thpompson[68]，calvetti[69]应用离散单元法于岩体的的计算分析基础上，郑立宁等[70]运用数值模拟分析膨胀土路基在干湿胀缩循环下的破坏特征及过程。合理再现膨胀土的室内膨胀率，膨胀里等测试。 通过对膨胀土路基多次胀缩模拟，从细观角度揭示了膨胀土路基在胀缩循环下的浅表层塌滑现象。
参考文献
[1] GATTERMANN, WITTKE J, ERICHSEN W. Modelling water uptake in highly compacted bentonite in environmental sealing barriers[J].Clay Minerals,2001,36(3):435-446.
[2] CHIJIMATSU M, FUJITA T, KOBAYASHI A, et al. Experiment and validation of numerical simulation of coupled thermal, hydraulic and mechanical behaviour in the engineered buffer materials[J].International Journal for Numerical & Analytical Methods in Geomechanics,2000,24(4):403-424.
[3] HOLTZ W G, GIBBS H J. Engineering Properties of Expansive Clays[J].Theoretical Biology & Medical Modelling,1956,8(1):269-276.
[4] 何山，朱珍德，王思敬. 膨胀岩的判别与分类方法探讨[J]. 水利水电科技进展, 2006, 26(4): 62-64, 86.
[5] 朱训国，杨庆. 膨胀岩的判别与分类标准[J]. 岩土力学, 2009(S2): 174-177.
[6] 刘晓丽，王思敬，王恩志，等. 含时间效应的膨胀岩膨胀本构关系[J]. 水利学报, 2006(02): 195-199.
[7] 罗鸿禧，康哲良，徐昌伟. 膨胀岩的胀缩与孔径分布的试验研究[J]. 水文地质工程地质, 1991(02): 27-30.
[8] 朱珍德，邢福东，刘汉龙，等. 南京红山窑第三系红砂岩膨胀变形性质试验研究[J]. 岩土力学, 2004(07): 1041-1044.
[9] 朱珍德，张爱军，邢福东. 红山窑膨胀岩的膨胀和软化特性及模型研究[J]. 岩石力学与工程学报, 2005(03): 389-393.
[10] 朱珍德，杨永杰，蒋志坚，等. 用数字图像处理技术进行膨胀红砂岩细观结构动态劣变特征研究[J]. 岩石力学与工程学报, 2007(10): 2007-2013.
[11] ZHONG W, ZHANG R. THE PARAMETRIC VARIATIONAL PRINCIPLE FOR ELASTOPLASTICITY[J].Acta Mechanica Sinica,1988,4(2):134-137.
[12] 于春江，陈发根，谢添. 泥岩路基填料膨胀工程特性研究[J]. 公路交通科技(应用技术版), 2013(12): 16-19.
[13] 付宏渊，王意明，刘新喜. 炭质页岩路堤变形特性研究[J]. 中外公路, 2012(01): 19-23.
[14] 郑西贵，季明，张农. 考虑湿度效应的膨胀岩流变模型[J]. 煤炭学报, 2012(03): 396-401.
[15] GENUCHTEN M T V. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.[J].Soil Science Society of America Journal,1980,44(44):892-898.
[16] 赵二平，李建林，王瑞红. 不同应力路径下膨胀岩力学特性试验研究[J]. 人民长江, 2014(03): 83-86.
[17] 张巍，尚彦军，曲永新，等. 泥质膨胀岩崩解物粒径分布与膨胀性关系试验研究[J]. 岩土力学, 2013(01): 66-72.
[18] 何建新，刘亮，热依拉，等. 膨胀岩抗剪强度特性研究[J]. 水资源与水工程学报, 2013(06): 92-94.
[19] 益圣强. 膨胀土地基的综合改良法处理试验研究[J]. 山西建筑, 2014(21): 83-85.
[20] 李国富，李珠，戴铁丁. 膨胀岩力学性质试验与巷道支护参数的预测研究[J]. 工程力学, 2010(02): 96-101.
[21] LIN T T, SHEU C, CHANG J E, et al. Slaking mechanisms of mudstone liner immersed in water[J].Journal of Hazardous Materials,1998,58(1–3):261-273.
[22] 杨成祥，宋磊博，王刚，等. CT实时观察下泥岩遇水软化过程的机理[J]. 东北大学学报(自然科学版), 2015(10): 1461-1465.
[23] HEGGHEIM T, MADLAND M V, RISNES R, et al. A chemical induced enhanced weakening of chalk by seawater[J].Journal of Petroleum Science & Engineering,2005,46(3):171-184.
[24] 杨春和，冒海军，王学潮，等. 板岩遇水软化的微观结构及力学特性研究[J]. 岩土力学, 2006(12): 2090-2098.
[25] 车平，宋翔东，虞翔，等. 巢湖地区坟头组泥岩遇水软化特性与机理试验[J]. 同济大学学报(自然科学版), 2012(03): 396-401.
[26] 胡博聆，王继华，赵春宏. 滇中泥质粉砂岩崩解特性试验研究[J]. 工程勘察, 2010(07): 13-17.
[27] 谭罗荣. 关于粘土岩崩解、泥化机理的讨论[J]. 岩土力学, 2001(01): 1-5.
[28] 杨建林，王来贵，李喜林，等. 水作用下泥岩崩解过程与微观机理研究[J]. 水资源与水工程学报, 2014(01): 67-70.
[29] 黄宏伟，车平. 泥岩遇水软化微观机理研究[J]. 同济大学学报(自然科学版), 2007(07): 866-870.
[30] 张建智，俞缙，蔡燕燕，等. 渗水膨胀岩隧洞黏弹塑性蠕变解及变形特性分析[J]. 岩土工程学报, 2014(12): 2195-2202.
[31] 卢义玉，侯吉峰，尤祎，等. 湿度应力场作用下煤矿穿膨胀岩钻孔缩径规律研究[J]. 采矿与安全工程学报, 2014(03): 469-475.
[32] 杨军平，樊永华. 温湿效应对膨胀岩吸水软化作用的影响[J]. 科学技术与工程, 2016(20): 259-263.
[33] 吴珺华，袁俊平，杨松，等. 干湿循环下膨胀土胀缩性能试验[J]. 水利水电科技进展, 2013(01): 62-65.
[34] HANDY R L.:Soil Mechanics in Engineering Practice[J].Journal of Geology,1967,76(5):149-150.
[35] GAMBLE J C. Durability - Plasticity Classification of Shales and Other Argillaceous Rocks[J].1971.
[36] 刘晓明，熊力，刘建华，等. 基于能量耗散原理的红砂岩崩解机制研究[J]. 中南大学学报(自然科学版), 2011(10): 3143-3149.
[37] 罗勇. 浅谈水泥改性土在南水北调南阳膨胀土试验段工程中的运用[J]. 河南水利与南水北调, 2013(13): 46-47.
[38] 杨和平，贺迎喜，江唯伟. 微生物影响岩土工程的现状及用微生物技术改良膨胀土的思考[J]. 中外公路, 2007(04): 228-231.
[39] 王小军，赵中秀. 用生石灰改良膨胀岩基床的室内试验研究[J]. 岩土工程学报, 1998(01): 39-43.
[40] 何爱华，蒋志明，夏鑫，等. CPAM黏土改性液性能试验研究[J]. 公路与汽运, 2014(06): 98-102.
[41] 邹开学，李玉锋，韩会增. 南昆铁路膨胀岩(土)改良研究[J]. 路基工程, 2000(4): 1-4.
[42] Han Xu，Bin Huang，Xiaomin He(徐晗，黄斌，何晓民)(膨胀岩工程特性试验研究)[Z]. 南京: 2007716-722.
[43] 王建军，何小林，王涛. 膨胀岩路堑边坡病害机理与防护技术探讨[J]. 科技创新导报, 2014(26): 105-108.
[44] 傅毅静，张强. 风化软岩路堤填筑与填料改良试验研究[J]. 路基工程, 2013(01): 88-91.
[45] 陈花林. 石灰改良泥质粉砂岩路基填料的试验研究[J]. 铁道工程学报, 2008(03): 10-14.
[46] 毛忠良. CMA改性剂及其在衡茶吉铁路膨胀土路基工程中的应用[J]. 路基工程, 2012(02): 25-28.
[47] Hao Wang(王浩)(DAH石灰混合溶液渗透方法改良膨胀土技术研究)长安大学, 2003.
[48] 杨国录，叶建民，陈士强，等. 南水北调中线工程膨胀土改性施工技术[J]. 节水灌溉, 2008(01): 54-57.
[49] 颜春，刘彦，黄中文. 树根桩+CMA混合溶液在膨胀土路堑边坡处治中的综合应用[J]. 重庆交通学院学报, 2006(03): 69-73.
[50] 青兆麟. 泥岩填石路堤施工技术应用[J]. 内蒙古公路与运输, 2012(02): 12-14.
[51] 成海坡. 泥质软岩填筑路基施工技术探讨[J]. 工程建设与设计, 2015(03): 115-117.
[52] 韩强. 膨胀岩渠段施工中防水毯的应用[J]. 中国水运(下半月), 2012(09): 249-250.
[53] 王福和. 土工格栅增强泥岩填料路基稳定性分析[J]. 公路交通科技(应用技术版), 2008(03): 22-24.
[54] 惠会清，胡同康，王新东. 石灰、粉煤灰改良膨胀土性质机理[J]. 长安大学学报(自然科学版), 2006(02): 34-37.
[55] 陈涛，张文慧，蓝日彦，等. 风化泥岩路基填料水稳定性室内试验研究[J]. 科学技术与工程, 2012(35): 9580-9584.
[56] 施建振，陈勇，王保田. 风化泥岩作为路基填料的试验研究[J]. 现代交通技术, 2010, 7(1): 29-30.
[57] 郑明新，方焘，刁心宏，等. 风化软岩填筑路基可行性室内试验研究[J]. 岩土力学, 2005(S1): 53-56.
[58] 刘新喜，夏元友，刘祖德，等. 复杂应力下强风化软岩湿化变形试验研究[J]. 岩石力学与工程学报, 2006(05): 925-930.
[59] 周小顺，王岗，覃朗. 石灰改性膨胀土机理及施工质量关键点研究[J]. 科学之友, 2012(02): 69-71.
[60] 龚壁卫，程展林，胡波，等. 膨胀土裂隙的工程特性研究[J]. 岩土力学, 2014(07): 1825-1830.
[61] 范秋雁，刘金泉，杨典森，等. 不同降雨模式下膨胀岩边坡模型试验研究[J]. 岩土力学, 2016(12): 3401-3409.
[62] 程永辉，程展林，张元斌. 降雨条件下膨胀土边坡失稳机理的离心模型试验研究[J]. 岩土工程学报, 2011(S1): 416-421.
[63] 王军，龚壁卫，张家俊，等. 膨胀岩裂隙发育的现场观测及描述方法研究[J]. 长江科学院院报, 2010(09): 74-78.
[64] 童学军，张文学，杨有海. 铁路路基膨胀土填料改良试验研究[J]. 交通运输研究, 2011(20): 125-129.
[65] 辛文栋. 风蚀环境下膨胀性泥岩路堑边坡防护措施应用研究[J]. 铁道标准设计, 2012(12): 22-26.
[66] 王凯，刁心宏. 膨胀岩湿度应力场本构模型二次开发研究[J]. 岩石力学与工程学报, 2015(S2): 3781-3792.
[67] CUNDALL P. A computer model for simulating progressive large scale movement in block rock systems[C]. 1971.
[68] THOMPSON N, BENNETT M R, PETFORD N. Development of characteristic volcanic debris avalanche deposit structures: New insight from distinct element simulations[J].Journal of Volcanology & Geothermal Research,2010,192(3):191-200.
[69] CALVETTI F, CROSTA G, TATARELLA M. Numerical simulation of dry granular flows: From the reproduction of small-scale experiments to the prediction of rock avalanches[J].Rivista Italiana Di Geotecnica,2000,2(2):21-38.
[70] 郑立宁，谢强，柳堰龙，等. 基于颗粒流的膨胀土路基破坏特征分析[J]. 铁道工程学报, 2011(03): 32-36.
|
本项目学生有关的研究积累和已取得的成绩
为避免实验开展的盲目性，相关人员收集并阅读了大量有关膨胀岩化学改良的文献，分类整理后主要归结为如下十类：膨胀岩的分类；本构关系；室内实验；数值模拟；处理方法；工程实例；边坡失稳研究；岩石风化程度；对膨胀岩的新认识。如图4，对不同分类下的文献
我们也做了较为深入的学习与总结。如下图5是一部分本构关系的整理分类，按照时间、作者、研究内容的顺序进行了阅读与整理。
 图4 文献分类
图5 本构关系相关文献
对膨胀岩的化学改良原理有了一定了解。同时利用暑假时间做了膨胀岩的化学改良实验，获得了不少实验数据，下面是一些实验的结果。如图6
图6 不同水泥粘土比改良后试样抗压强度
在实验方面：我们已经找到了膨胀岩路用的化学改良材料，分别是粘土和水泥。确定水泥和粘土占总质量的25%左右，试件在无侧限压强实验下有最大强度，并且膨胀岩、粘土、水泥、水四种物质配合在一起无侧限压强普遍在400KPa以上，能达到路用要求。另外还总结了配合比的变化对改良后试件强度的影响，水泥粘土占总质量的比例不变，粘土比水泥为5：1，4：1，3：1，试件强度增加，如果固定水泥和粘土的比例不变，大致呈现水泥和粘土占总质量的增加，试件强度先增加后减小的规律。
首先进行纵向对比，反应的是水泥和黏土比值一定的条件下，膨胀岩含量不同的情况下，7天无侧限抗压强度的对比。从图中可以看出，在水泥与黏土含量1：3和1：5的比值下，不添加任何膨胀岩时强度最高，膨胀岩含量从75%-60%的过程中，强度呈现从高到底再从低到高的趋势，在水泥与黏土含量1:4和1:6的比值下，强度由低到高又由高到低最后由低再到高。由此可见，改良后试件的强度在60%-75%范围内是先增加后降低再增加，而最佳的膨胀岩配比是70%-75%之间。
而在数值模拟方面：能较为熟练地使用颗粒流程序PFC2D,PFC3D模拟一些实验过程中的力学性质与颗粒膨胀，挤压的变化，如图7、图8；同时也能较好地使用EDEM软件对颗粒的简化关系进行相应的模型建立与模拟，下一步即引入clump用多个颗粒组成模板模拟实际颗粒探究相应配合比下的力学特性，同时在附件中给出了相应的模拟程序例子。
图7 PFC模拟颗粒的堆积碰撞
图8 EDEM对球形颗粒的填充
由于膨胀岩本身的特殊性，我们初步的学习是将实际模型简化成球体来处理，但进一步的实验与模拟需要考虑颗粒的破碎和复杂颗粒模型的替换来取代球体颗粒的简化。
经学习，我们掌握了origin等专业绘图软件对颗粒数据进行二次处理如图9：
图9 origin绘图
|
项目的创新点和特色
1. 膨胀岩的化学改良有不少方案，但大都针对具体案例，较少有系统的研究，我们尝试提出水泥黏土这种新的膨胀岩改良方案并从一定强度下最佳配合比、一定配合比下的对应强度两个角度对弱膨胀岩、中膨胀岩展开研究，且初步的研究表明在弱膨胀岩中这种改良方法是满足路用要求的。
2.有不少学者使用生石灰对膨胀岩进行改良，生石灰遇水反应生成碳酸钙，而碳酸钙与水、空气三相作用下生成可溶性碳酸氢钙，不仅破坏改良土体结构，还导致水分渗入土体，使用水泥黏土则能很好的避免这种情况，由于黏土压实后渗透性极差能避免水分的渗入导致改良土的破坏，能更稳定地提高改良膨胀岩的强度。
3.积极研究膨胀岩的化学改良方案可以给工程带来极大的经济效益，并且减少工程中膨胀岩的浪费。目前大多学者使用石灰与水泥或者二灰土进行改良，我们尝试着用水泥与粘土来对膨胀岩进行改良，石灰相对于粘土较为昂贵，这有利于节约经济，与其他改良膨胀岩的方法相比, 水泥粘土改良是成本低、原料易得且绿色环保。
4.实验已有初步的结论，同时也验证了新改良方法的可行性。初步的实验采用的设计方案：水泥比粘土为1：3，1：4，1：5,1:6四组（工程中常用的用量），及粘土比水泥为3：1，4：1，5：1，1：6四组不加膨胀岩的实验组下7天，28天改良膨胀岩的无侧限抗压强度，如下图，在不考虑压实系数因素影响情况下，部分改良后膨胀土的28天无侧限抗压强度达到0.7MPa，完全符合高速公路路基填料的要求。
图10 水泥黏土比与28天抗压强度关系
5.由于从新的角度提出了膨胀岩的改良方法，相应改良工艺与原理可申请关专利。
6.取样自张桑高速，初研实验成果较好，进一步的研究成果可用于相应地段的处理，因此相应研究具有较强实用性。
7.有数值模拟作为参照实验，目前较多学者使用ansys等对膨胀岩进行模拟，但是结果不够理想，出于岩体属于离散元的考虑，我们采用两种离散元程序来对颗粒进行三轴无侧限抗压模拟，EDEM能给出初步而又规律的结论，而PFC程序由于软件自身模型功能的可操作性较大适合精确模拟。
8.在已有的PFC、EDEM数值模拟基础上基于离散元原理结合工程实际建立膨胀岩考虑湿度效应下的软化崩解模型，编写满足工程实际与室内实验对照需求的复杂颗粒替换代码，从新的角度模拟泥岩崩解宏微观变化的破碎崩解情况。
|
项目的技术路线及预期成果:
项目的技术路线:
我们采用目前较为普遍的方法是朱训国提出的一种以膨胀岩中亲水矿物成分的含量为基础控制指标，以岩块的干燥饱和吸水率、极限膨胀量、极限膨胀力为主要控制指标新的标准。我们先对原样进行初步的判断（如表），再进行相应配合比下的实验。经过初步判断，取自张桑高速（如图11）的膨胀岩是属于弱膨胀岩再对膨胀岩进行初研时进行水泥粘土的改良。
表一 膨胀岩判断依据指标
图11张桑高速
在工程中应用较为广泛的改良物料是水泥，石灰和其他化学改性剂，生石灰吸水后变成碳酸钙，碳酸钙暴露在空气与土壤中会产生化学反应生成碳酸氢钙，影响改性效果，同时易诱发灾害，而化学改性剂成本较高，综合经济与施工的考虑，我们采用水泥与粘土做原材料来研究对膨胀岩的改良。但针对膨胀岩的研究大多停留在针对具体的工程而采取的措施上，缺乏系统的探究。因此我们提出一种新的改良膨胀岩方法：水泥粘土的化学改良，试探究不同配合比下对应的石灰粘土改良膨胀岩无侧限抗压强度。
碎石类土和掺水泥的级配碎石直径和高均为150mm，因此我们将试件制成如下模型（图12），即直径与高均为150mm的PVC管；
图12 装填后的PVC管试样
张巍等[17]提出膨胀岩崩解物粒径分布的差异性对其膨胀性具有一定的指示意义。由此在实验中我们：使膨胀岩破碎之后根据规范最大粒径不能超过37.5mm。因此在破碎膨胀岩之后应将膨胀岩过37.5mm筛孔。采用的振动筛，如图13、图14
图13 标准振动筛
图14振动电机
振动筛为ZBS-92A型，摇动次数221次/min,摆动行程25mm,摇摆次数147次/min,电机功率0.37KW.
对待不同的膨胀岩应该换不同的配比，应控制不同的水泥粘土比，在预研中我们结合王小军等的研究成果考虑引入水泥的因素，弱膨胀岩的水泥掺入量可以定为1%，2%，3%，4%四组方案，中膨胀岩的水泥掺入量可以设为2%，3%，4%，5%四组，强膨胀岩的水泥掺入量可以设为2%，3%，4%，5%四组。
由于确定我们采集的张桑高速膨胀岩试样为弱膨胀岩，我们调整水泥和粘土占总质量的25%，30%，35%，40%，水泥比粘土为1：3，1：4，1：5,1:6四组，除此之外还增加了粘土比水泥为3：1，4：1，5：1，1：6四组不加膨胀岩的实验组，
膨胀岩破碎后应该风干，测试过风干后膨胀岩的含水率在进行实验。风干后膨胀岩的含水率用来计算膨胀岩的真实使用量。加入水的剂量为膨胀岩的最优含水率，将粘土加入拌合均匀浸润12-24小时。浸润后将水泥加入拌合均匀，并在一小时做成试件。按照膨胀岩干密度和试件体积算出试样需要的拌合料质量，据此给试件模子加料。模子内壁先涂油，分2到3次加料，每次加料后用捣实棒捣实，之后用千斤顶压实试件。压实过程中要有两个盖子盖在试件上下两侧，离试件2cm左右，将填料压实至设计尺寸。每一个配比做成4个试样，试件做完后称量保证试件质量在30g的范围内波动。
试样养护条件为每天中午与下午各洒一次水，用量以使试样完全湿润但不能漫出为宜。洒水7天养护测一次无侧限压强，洒水28天养护测一次无侧限压强。通过数据分析得到最佳膨胀岩、粘土、水泥配合比，同时评价各个配合比下试件的强度。养护7天后可用总体一半拿出来做无侧限抗压强度实验，试件放在实验机上，应使试件变形等速增加，保持速率约为1mm/min。
下图15-17为膨胀岩65%黏土,30%,水泥5%,水6.3%的改良膨胀岩与粘土和水泥做成的试件在三轴仪上破坏后的图片。
粘土和水泥做成的试件破坏截面与试件底面成45度角，主要是剪切破坏，而膨胀岩、粘土、水泥做成的试件在无侧限压缩下先向侧向膨胀最后沿着膨胀岩和粘土接触界面破坏，膨胀岩自身完好。
图15水泥粘土试样破坏形式
图16水泥粘土试样破坏形式
图17改良膨胀岩试样破坏形式
三轴实验仪器如图18：
图18三轴实验仪
仪器参数为TSZ-1型：试样直径：39.1,61.8mm,轴向压力：0-10KN;围压反压：0-2MPa;孔隙压力：0-2MPa;体积变化：0-50ml;剪切速率：0.001-2.4mm/min.
进一步的研究中，而在进一步的实验中，在结合邹开学等实验的基础上，设定水泥四种配比，粘土七种配比，混合后有28种组合，每种组合做12个试件，每个级配试做6个模型，其中三个为7天，三个为28天。一半在水中养护，一半在低温下养护，有两种养护方式，每种要有三个用来做七天强度实验，三个用来做28天强度实验，总计336个试样。
混合料采取5次加入，每一次填料后用击实锤砸实25次，试件做完后称量保证试件质量在30g的范围内波动。洒水7天养护测一次无侧限压强，洒水28天养护测一次无侧限压强。通过数据分析得到最佳膨胀岩、粘土、水泥配合比，同时评价各个配合比下试件的强度。
在数据处理时无侧限抗压强度小于2mpa,计算至0.01mpa,大于2mpa,计算至0.1mpa.另一半试件在龄期28天后在进行无侧限抗压强度实验。目前的实验只是实验性的对特定级配下的膨胀岩做了初步的探究，由于经费和实验仪器的限制，
尚待更进一步开展对一定围压下的实验及膨胀岩的风化，模拟软件分析，应力本构关系等进一步的研究。
图19 应力应变曲线
在实验数据方面，由于实验数据较多，我们选取部分数据来补充支撑我们的申请书，膨胀岩的改良基本步骤已在申请书中得到体现，而改良效果的鉴定则分为如下两个方面：测试改良后试样的7天无侧限压缩抗压强度与28天无侧限压缩抗压强度。调整膨胀岩水泥黏土三者的含量比后装填试样进行实验，统计不同的实验数据，同时我们也给出了一部分经Excel、origin数据处理后的配合比-强度图，如图20-22。
图20 配合比与抗压强度关系综合图
图21 水泥粘土配合比与28天抗压强度关系
图22水泥粘土配合比与7天抗压强度关系
当我们的实验有较为完整的结论后我们计划将相应水泥粘土配合比下的改性应用于实际工程中。
在数值模拟方面，结合对离散元软件的使用，从三个角度来探究水泥粘土改良膨胀岩的特性：1.利用PFC多个颗粒组合形成的clump来模拟颗粒的聚合，设置其接触模型为粘性接触模型。2.模拟颗粒聚合稳定后遇水膨胀的现象。3.从颗粒三轴实验模拟程序出发，尝试对改良膨胀岩的聚合物进行无侧向抗压测得相应强度。目前我们学习开发了一些相应的程序代码来支撑我们的数值模拟研究，限于篇幅，我们在附件中给出了一部分代码。
数据分析：
将膨胀岩作为一个变量，粘土和水泥作为一个变量，同时也可控制水的含量，用控制变量法的原理来研究膨胀岩在哪种配合比下有最大强度和膨胀岩的试件强度随着配合比的变化是如何变化。
预期成果：
1.撰写项目“膨胀岩的化学改良及其路用特性研究”的结题报告；
2.本项目相关研究成果拟发表核心及以上期刊学术论文2~3篇 。
3.借助PFC、EDEM基于离散元结合工程实际建立膨胀岩考虑湿度效应下的软化崩解模型，编写满足工程实际与室内实验对照需求的复杂颗粒替换代码，结合泥岩崩解宏微观变化撰写颗粒破碎崩解的运行代码，
4.确定综合满足张桑高速公路路基与边坡填料工程与节约耗资要求下的最佳水泥黏土配合比，给出不同强度要求下最佳配合比。
5.刻画不同水泥黏土比改良后膨胀岩7天、28天无侧限抗压强度关系，由此预测不同配合比的改良膨胀岩所对应的抗压强度。
6.依托水泥粘土这一新型膨胀岩改良方法，对实验成果的应用措施与关键技术申请相应专利，针对构建的岩石崩解等破坏形式撰写代码并申请专利。
|
年度目标和工作内容
第一年度：
工作目标：以实验探索和完善膨胀岩的化学改良方案，使膨胀岩的水泥粘土化学改良方案确实可行。
工作内容：
2018年3月-4月
收集整理相关文献、资料，设计实验方案；
整理包括实验相关的文献的本构关系；室内实验；数值模拟；处理方法；工程实例；边坡失稳研究；岩石风化程度；对膨胀岩的新认识相关文献
2018年5月
采购实验装置及原料；在膨胀岩路段等购回一些相应的粘土并运回相应膨胀岩，三轴实验仪可向学校申请，购买振动筛以筛分不同粒径的膨胀岩碎石。在网上商城或附近建材市场，购买PVC管，在附近的批发市场购置一定量的水泥，
2018年6月-9月
购回膨胀岩后开展实验，实验分为两个方面：
1.确定一定比例的膨胀岩的含量情况下针对不同配合比的水泥粘土比来开展膨胀岩的改良研究
2.确定一定配合比的水泥粘土比调整膨胀岩的含量，经过水分养护后测量7天与28天的无侧限抗压强度
2018年10月
分析实验结果；当我们的实验有较为完整的结论后我们计划将相应水泥粘土配合比下的改性应用于实际工程中。
3.撰写一定强度下一定质量比膨胀岩下最优配合比的成果论文；撰写初步探究数值模拟应用于颗粒替换与三轴实验模拟的成果论文。
2018年11-12月
1.对实验结果进行分析整理并进行对比，得出相关结论；将膨胀岩碎块，粘土，水泥，水按照一定比例配合制成试样，每一个配比制作多个个试件，洒水7天养护测一次无侧限压强，洒水28天养护测一次无侧限压强。通过数据分析得到最佳膨胀岩、粘土、水泥配合比，同时评价各个配合比下试件的强度。
2.学习颗粒流数值模拟软件，主要表现在对PFC与EDEM的复杂实际颗粒模型的导入与颗粒的破碎，崩解泥化的模拟机制与代码学习上。
依托水泥粘土这一新型膨胀岩改良方法，对实验成果的应用措施与关键技术申请相应专利，针对构建的岩石崩解等破坏形式撰写代码并申请专利。
第二年度：
工作目标：采用控制变量法来对比实验数据，得出膨胀岩在不同的配比下强度之间的关系，给出张桑高速公路路基填料与路堤填料对应的最佳配比方案。
工作内容：
2019年1月
借助数值模拟软件设计确定一定比例的膨胀岩的含量情况下针对不同配合比的水泥粘土比来开展膨胀岩的改良研究确定一定配合比的水泥粘土比调整膨胀岩的含量，经过水分养护后测量7天与28天的无侧限抗压强度。
撰写满足一定强度下不同质量比膨胀岩下最优配合比的论文；撰写构建膨胀土的数值模拟分析的成果撰写论文；
2019年2月
撰写膨胀岩的化学改良及其路用特性研究结题报告并提交相应论文专利的成果；给出膨胀岩如何用clump来模拟实际颗粒变形受力的程序；给出膨胀岩如何三轴实验下颗粒受力破坏的程序与构建膨胀岩三轴实验下应力应变的数值模拟分析程序的代码。
2019年3月-4月
结合项目申请的专利、论文及相关数值模拟结果，申报结题。
|
指导教师意见
膨胀岩遇水膨胀失水收缩的特性常常诱发各种工程事故，甚至重大安全事故。本项目以张桑高速某膨胀岩边坡为背景，研究膨胀岩碎石体的路用性能改良。本小组同学在前期查阅了大量文献，设计出了可行的研究方案并取得了部分研究成果，研究方案可行，建议优先资助
签字： 胡敏 日期：2018.4.25
|
附件一：基于PFC利用clump生成复杂颗粒模拟实际颗粒力学特性
;新执行
new
SET random ;reset random-number generator
set echo=off
;***************** 定 义 *****************
;*****************************************************
def setup
;**定义空间大小-------------------------------------
height= 0.200 ;
width= 0.100 ;
vel_rad_z = 0.14; 0.1-0.5; ratio of velocity of rad and z
;==========需要修改的参数==========
gradeN = 2 ; d>=16
ballN = 360; 颗粒个数
vol_expand = 30; 体积扩大系数
;多面体类型
polyType= 'ELLIP'
;最小填充球的直径
filledballd= '3'
;多面体最小包容盒三个边长，按大到小编号为L，B，T
Ratiostr= '[1-0.5-0.25]'
Ratio12= 2/1 ; = L/B
Ratio23= 2/1 ; = B/T
;压缩参数
_Wzvel = -5.0 ; 5#墙体速度
Target_press = -1e5 ; 压缩停止条件：围压大小
b_acc = 2000 ; 球体加速度
;==================================
s_stiff=1e8 ;initial stiffnesses
n_stiff=1e8
w_stiff=1e100
R_critical = 0.0095/2 ;需要替换的球体的最小半径(临界半径)
;径向与轴向的扩大系数
yy= 4*vol_expand*vel_rad_z*height/width
nw= exp(ln(abs(yy))/3) ;=yy^(1/3)
nh= vol_expand/nw^2
height_cy= height*nh
width_cy = width *nw
cmpstage =1 ;1-第一阶段，2-第二阶段
;-----------------------------------------------------
;**定义级配信息
rMult = 1.1
;alph = sqrt(Ratio12*Ratio12+1+(1/Ratio23/Ratio23)); for polyhedron
alph = Ratio12 * rMult ; for ellipsoid
array p(10) R(10) v(10) Rs(7)
Rs(1) = '19'
Rs(2) = '16'
Rs(3) = '13.2'
Rs(4) = '9.5'
Rs(5) = '4.75'
Rs(6) = '2.36'
Rs(7) = '1.18'
R(1) = 0.019/2
R(2) = 0.016/2
R(3) = 0.0132/2
R(4) = 0.0095/2
R(5) = 0.00475/2
R(6) = 0.00236/2
R(7) = 0.00118/2
loop k (1,7)
if R(k)>= R_critical
R(k)= alph * R(k) ;粒径扩大
end_if
end_loop
p(1) = 0 ;19
p(2) = 5 ;16
p(3) = 11 ;13.2
p(4) = 14 ;9.5
p(5) = 22 ;4.75
p(6) = 14 ;2.36
p(7) = 9.5 ;1.18
v(1) = 0 ; v=1/(r^3) r=(R1+R2)/2
v(2) = 0.0015
v(3) = 0.0026
v(4) = 0.0055
v(5) = 0.0221
v(6) = 0.1781
v(7) = 1.4427
pvsum=0 ;p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4)+p(5)*v(5)+p(6)*v(6)+p(7)*v(7)
loop k (1, gradeN)
pvsum = pvsum+ p(k)*v(k)
end_loop
;-----------------------------------------------------
; **孔隙率参数
vpoly = 0 ;total volume of polyhedrons or ellipsoids
;-----------------------------------------------------
; **定义存取参数
a_size = ballN+1
IO_READ = 0
IO_WRITE= 1
IO_FISH = 0
IO_ASCII= 1
;读取数据参数
a_size = 200
array ab(200) x_(200) y_(200) z_(200) r_(200) ;
array asz(5000) asr(5000)
nsz = 5000 ; 应力数据平滑长度 建议1000~5000
avsz= 0
avsr= 0
idx = 0
;保存文件名头
filehead= 'R16D100H200'
saveFile= filehead+polyType+Ratiostr+'.SAV'
;-----------------------------------------------------
end
;----------------------------------------------------
;生成墙面
def make_walls ; create walls: a cylinder and two plates
extend=1.0 ;0.2
rad_cy=0.5* width_cy
;圆柱面
_z0=-extend*height
_z1=height_cy*(1.0+extend)
command
wall type cylinder id=1 kn=w_stiff end1 0.0 0.0 _z0 end2 0.0 0.0 _z1 &
rad rad_cy rad_cy
end_command
;底面
_x0=-rad_cy*(1.0+extend)
_y0=-rad_cy*(1.0+extend)
_z0= 0.0
_x1= rad_cy*(1.0+extend)
_y1=-rad_cy*(1.0+extend)
_z1= 0.0
_x2= rad_cy*(1.0+extend)
_y2= rad_cy*(1.0+extend)
_z2= 0.0
_x3=-rad_cy*(1.0+extend)
_y3= rad_cy*(1.0+extend)
_z3= 0.0
command
wall id=5 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
;顶面
_x0=-rad_cy*(1.0+extend)
_y0=-rad_cy*(1.0+extend)
_z0= height_cy
_x1=-rad_cy*(1.0+extend)
_y1= rad_cy*(1.0+extend)
_z1= height_cy
_x2= rad_cy*(1.0+extend)
_y2= rad_cy*(1.0+extend)
_z2= height_cy
_x3= rad_cy*(1.0+extend)
_y3=-rad_cy*(1.0+extend)
_z3= height_cy
command
wall id=6 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
end
;----------------------------------------------------
;生成颗粒
def create_ball
n1 = 0
n2 = 0
w1 = -rad_cy
w2 = rad_cy
h1 = 0.0
h2 = height_cy
loop k (2, gradeN)
n1 = n2+1
n2 = n2+round(ballN*p(k)*v(k)/pvsum)
if n2>n1 then
r1 = R(k)
r2 = R(k-1)
clr=k-2
command
gen id=n1,n2 rad=r1,r2 x=w1,w2 y=w1,w2 z=h1,h2 filter ff_cylinder
prop dens=2700 ks=s_stiff kn=n_stiff color=clr range id=n1,n2
end_command
end_if
end_loop
ballmax = n2
ii=out(string(ballmax)+'particles were created')
end
;----------------------------------------------------
;修改墙体刚度 change lateral wall stiffnesses
def cws
command
wall type cylinder id 1 kn=w_stiff
end_command
end
;----------------------------------------------------
;过滤函数
def ff_cylinder
ff_cylinder=0
_brad= fc_arg(0)
_bx= fc_arg(1)
_by= fc_arg(2)
_bz= fc_arg(3)
_rad= sqrt(_bx^2+_by^2)
if _rad+_brad > rad_cy then
ff_cylinder= 1
end_if
end
;----------------------------------------------------
def read_balls
status = open(txtfn, IO_READ, IO_ASCII)
status = read(ab, a_size)
status = close
nball= parse(ab(1), 1)
loop n (1, nball)
x_(n)= parse(ab(n+1), 1)
y_(n)= parse(ab(n+1), 2)
z_(n)= parse(ab(n+1), 3)
r_(n)= parse(ab(n+1), 4)
end_loop
vpoly = vpoly+ parse(ab(nball+4), 4)
;ii=out(' nball='+string(nball)+' vpoly='+string(vpoly))
end
;-----------------------------------------------------
def clump_make
bp= ball_head
clmp_id= 10000
n3a= ballmax+1
n3b= ballmax+1
loop while bp # null
if b_id(bp)<=ballmax ;get rid of all new balls
r_base = b_rad(bp)
x_base = b_x(bp)
y_base = b_y(bp)
z_base = b_z(bp)
section
loop k (1,7)
if r_base>= R(k)
if r_base>= R_critical
txtfn= polyType+Rs(k)+'M'+filledballd+'R'+Ratiostr+'T1.txt' ;+string(int(1+20*urand))+'.txt'
ii=out(' Read from... '+ txtfn)
read_balls ;读取数据（x，y，z，r）
ii=b_delete(bp)
; gen all the balls
n3a= n3b+1
scl= r_base/R(k)
cal_transpara
loop i (1,nball)
x1= x_base + scl*(axx*x_(i)+axy*y_(i)+axz*z_(i))
y1= y_base + scl*(ayx*x_(i)+ayy*y_(i)+ayz*z_(i))
z1= z_base + scl*(azx*x_(i)+azy*y_(i)+azz*z_(i))
r1= scl*r_(i)
n3b=n3b+1
command
;pause
ball id=n3b rad=r1 x=x1 y=y1 z=z1
end_command
end_loop
; make into a clump
clmp_id= clmp_id+1
command
clump id=clmp_id perm range id=n3a,n3b ;not broken
end_command
else
vpoly = vpoly+(4.0/3.0)*pi*b_rad(bp)^3
endif
exit section
endif
end_loop
end_section
endif
bp=b_next(bp)
end_loop
command
property density=2650 ks=s_stiff kn=n_stiff fric=0
end_command
cmpstage = 2
end
;----------------------------------------------------
def cal_transpara
;绕x、y、z轴旋转角度
ax= urand *2.0*pi
ay= urand *2.0*pi
az= urand *2.0*pi
axx= cos(ay)*cos(az)-sin(ax)*sin(ay)*sin(az)
axy=-cos(ax)*sin(az)
axz= sin(ay)*cos(az)+sin(ax)*cos(ay)*sin(az)
ayx= cos(ay)*sin(az)+sin(ax)*sin(ay)*cos(az)
ayy= cos(ax)*cos(az)
ayz= sin(ay)*sin(az)-sin(ax)*cos(ay)*cos(az)
azx=-cos(ax)*sin(ay)
azy= sin(ax)
azz= cos(ax)*cos(ay)
end
;----------------------------------------------------
; 墙体索引
def wall_addr
wadd1=find_wall(1)
wadd5=find_wall(5)
wadd6=find_wall(6)
end
;----------------------------------------------------
; 应力数组清零
def aszero
idx = 0
avsz= 0
avsr= 0
loop i (1,nsz)
asz(i)= 0
asr(i)= 0
end_loop
end
;----------------------------------------------------
; 试件压缩
def compress
Wrvel= vel_rad_z*_Wzvel
Wzvel= _Wzvel
loopk= 0
cmpbefore = 0 ; rad=1, z=-1, unknown=0
lpkn= 8
if cmpstage =1
lpkn= 4
end_if
loop k (1,lpkn)
;消除原有应力
;aszero
zeroall
command
wall id 6 zvel= 0
cycle 5000 ; to release the old stress
end_command
;aszero
cycn = 0
;新压缩循环
section ;+++++++++
loop while 1 # 0
;调整速度
if new_rad > 0.5*width
w_radvel(wadd1)= Wrvel
Wzvel0 =0
if cmpbefore<= 0
bar_apply ;only rad-acc
end_if
cmpbefore= 1
else
w_radvel(wadd1)= 0
Wzvel0 = Wzvel
if cmpbefore>= 0
baz_apply ;only z-acc
end_if
cmpbefore=-1
end_if
;循环计算
command
wall id 6 zvel= Wzvel0
cycle 50
end_command
cycn = cycn+1
;显示结果
statereq
;退出
if abs(avsr) >= abs(0.9*Target_press)
if cycn <= 5
loopk= loopk+1
if loopk >= 2 ;change the velocity
if new_rad > 0.5*width
Wrvel = Wrvel*0.5
else
Wzvel = Wzvel*0.5
end_if
end_if
if loopk >= 4 ;change the velocity
w_radvel(wadd1)= 0
ba_release
exit
end_if
else
loopk= 0
end_if
exit section
end_if
end_loop
end_section ;+++++++++
w_radvel(wadd1)= 0
end_loop
;清除加速度
ba_release
end
;----------------------------------------------------
; 球体加载卸载
def bar_apply
bp = ball_head
loop while bp # null
bforce = b_acc * 11100.29404 * b_rad(bp)^3 ;(4.0/3.0) * pi * b_rad(bp)^3 * 2650
brad = sqrt( b_x(bp)^2 + b_y(bp)^2 )
;if brad >= width/3 ;width/2 * (2/3)
b_xfap(bp) = -1.0 * bforce * b_x(bp) / brad
b_yfap(bp) = -1.0 * bforce * b_y(bp) / brad
b_zfap(bp) = 0
;else
; b_xfap(bp) = 0
; b_yfap(bp) = 0
; b_zfap(bp) = 0
;end_if
bp = b_next(bp)
end_loop
end
def baz_apply
bp = ball_head
loop while bp # null
bforce = b_acc * 11100.29404 * b_rad(bp)^3 ;(4.0/3.0) * pi * b_rad(bp)^3 * 2650
b_xfap(bp) = 0
b_yfap(bp) = 0
b_zfap(bp) = -1.0 * bforce
bp = b_next(bp)
end_loop
end
def ba_release
bp = ball_head
loop while bp # null
b_xfap(bp) = 0
b_yfap(bp) = 0
b_zfap(bp) = 0
bp = b_next(bp)
end_loop
end
;----------------------------------------------------
; 求取平均应力和应变
def get_ss
get_stress
get_strain
end
def get_ass
;计算瞬间应力应变
get_stress
get_strain
;平均化处理
idx = idx+1
idx0= idx+1
if idx > nsz
idx = idx -nsz
end_if
if idx0 > nsz
idx0 = idx0 -nsz
end_if
asz(idx)= wszz
asr(idx)= wsrr
avsz = avsz + (asz(idx)-asz(idx0))/nsz
avsr = avsr + (asr(idx)-asr(idx0))/nsz
;ii=out(' avsz='+string(avsz)+' avsr='+string(avsr))
end
def get_stress
new_rad=w_radend1(wadd1)
zdif=w_z(wadd6)-w_z(wadd5) ;注意：w_z( )函数返回值不是wall的z坐标。
new_height= height_cy+zdif ; 为墙体的*累积位移*，包括以前所有cycle的位移量
wsrr=-w_radfob(wadd1)/(new_height*2.0*pi*new_rad) ;径向应力
wszz=0.5*(w_zfob(wadd5)-w_zfob(wadd6))/(pi*new_rad^2.0) ;轴向应力
;ii=out(' wradfob='+string(w_radfob(wadd1)))
end
def get_strain
rdif=new_rad-rad_real
werr=2.0*rdif/(rad_real+new_rad) ;径向应变
wezz=(new_height-height_real)/height_real ;2.0*zdif/(height+new_height) ;轴向应变
wevol=wezz+2.0*werr ;体积应变 = x应变+y应变+z应变
end
;----------------------------------------------------
; 将height_real重置为当前高度
def reset_real
height_real= abs(height_cy+w_z(wadd6)-w_z(wadd5))
rad_real = w_radend1(wadd1)
end
;----------------------------------------------------
;统计空隙率 (height_real may be modified)
def poroscal
reset_real ;reset the height and rad
tot_vol= height_real*pi*rad_real^2.0
sum=0.0 ;get actual porosity
bp=ball_head
loop while bp # null
sum= sum+4.0/3.0*pi*b_rad(bp)^3
bp = b_next(bp)
end_loop
pballs= (1.0-sum/tot_vol)*100
ppolys= (1.0-vpoly/tot_vol)*100
ii=out(' The last height is '+string(height_real)+',radius is '+string(rad_real))
ii=out(' ball-porosity is '+string(pballs)+'%'+',poly-porosity is '+string(ppolys)+'%')
end
;----------------------------------------------------
; 查询当前高度/直径
def statereq
height_cur=abs(height_cy+w_z(wadd6)-w_z(wadd5))
diam_cur= w_radend1(wadd1)*2
ii=out(' Nows height is '+string(height_cur)+',Diameter is '+string(diam_cur))
ii=out(' avsz ='+string(avsz)+', avsr ='+string(avsr))
ii=out(' Wzvel='+string(Wzvel)+', Wrvel='+string(Wrvel))
end
;----------------------------------------------------
; 清零
def zeroall
;zero all balls
bp = ball_head
loop while bp # null
loop k (1,3)
b_vvel(bp,k) =0
b_vrvel(bp,k)=0
end_loop
bp = b_next(bp)
end_loop
ii= out(' ok to zero all balls!')
;zero all clumps
clp = clump_head
loop while clp # null
loop k (1,3)
cl_vvel(clp,k) =0
cl_vrvel(clp,k)=0
end_loop
clp = cl_next(clp)
end_loop
ii= out(' ok to zero all clumps!')
end
;----------------------------------------------------
;***************** 运 行 *****************
;*****************************************************
plot create model ; view 1
plot set cap size 25
plot set mag 1.5
plot set rot 0 0 0 ;30 0 40
plot add ball blue green cyan orange red lblue lgreen lred
;plot add wall lgray
;plot set plane normal 1 0 0 origin 0 0 0
;plot close
setup
;生成墙体
make_walls
wall_addr
;生成颗粒
create_ball
;设置height/rad_real为初始值
reset_real ;为方便get_ss中wezz变量计算
history id=1 avsz ;wszz
history id=2 avsr ;wsrr
plot create StressStation ; view 2
plot hist 1 2
;history id=3 Wzvel
;history id=4 Wrvel
;plot create VelocityStation ; view 3
;plot hist -3 -4
;开启应力跟踪
set fishcall 0 get_ass
;压缩球体颗粒
compress
wall id 6 zvel=0
cyc 5000 ;20000
zeroall
;生成多面体颗粒
clump_make
;压缩多面体颗粒
aszero
compress
wall id 6 zvel=0
cyc 20000
zeroall
;关闭应力跟踪
set fishcall 0 remove get_ass
SET w_stiff=1e7 ;make lateral wall stiffness=1/10 of ball stiffness
cws
poroscal
save saveFile
return
;*****************************************************
附件二：基于PFC利用clump生成复杂颗粒模拟实际颗粒受力变形
;新执行
new
SET random ;reset random-number generator
set echo=off
;***************** 定 义 *****************
;*****************************************************
def setup
;**定义空间大小-------------------------------------
height= 0.100 ;
width= 0.0750 ;
vel_rad_z = 0.4 ; 0.1-0.5; ratio of velocity of rad and z
;==========需要修改的参数==========
gradeN = 6 ; d>=1.18
MaterialV = 0.0002 ; 颗粒总体积 m3
vol_expand = 10 ; 体积扩大系数
;压缩参数
_Wzvel = -5.0 ; 5#墙体速度
Target_press = -1e6 ; 压缩停止条件：围压大小
b_acc = 2000 ; 球体加速度
Name = 'ZhouChanghong'
;==================================
s_stiff=1e8 ;initial stiffnesses
n_stiff=1e8
w_stiff=1e10
;径向与轴向的扩大系数
yy= 4*vol_expand*vel_rad_z*height/width
nw= exp(ln(abs(yy))/3) ;=yy^(1/3)
nh= vol_expand/nw^2
height_cy= height*nh
width_cy = width *nw
;ii=out('height_cy='+string(height_cy)+'height='+string(height)+'nh'+string(nh))
;command
; pause;
;end_command
cmpstage = 1 ;1-第一阶段（压缩基球），2-第二阶段（压缩多面体）
;-----------------------------------------------------
;**定义级配信息
rMult = 1.0
array p(10) R(10)
R(1) = 0.016/2
R(2) = 0.0132/2
R(3) = 0.0095/2
R(4) = 0.00475/2
R(5) = 0.00236/2
R(6) = 0.00118/2
R(7) = 0.0006/2
;各粒径分计筛余：
p(1) = 0 ;16
p(2) = 5 ;13.2
p(3) = 22 ;9.5
p(4) = 55 ;4.75
p(5) = 0 ;2.36
p(6) = 8 ;1.18
p(7) = 10 ;0.6
psum=0.0
loop k (1, gradeN)
psum = psum+ p(k)/0.1*0.1
end_loop
;-----------------------------------------------------
; **孔隙率参数
vpoly = 0 ;total volume of polyhedrons or ellipsoids
;-----------------------------------------------------
; **定义存取参数
IO_READ = 0
IO_WRITE= 1
IO_FISH = 0
IO_ASCII= 1
;读取数据参数
array asz(5000) asr(5000)
nsz = 5000 ; 应力数据平滑长度 建议1000~5000
avsz= 0
avsr= 0
idx = 0
;统计数据
array angles(1000)
;保存文件名头
filehead= 'VMA'
saveFile= filehead+'-'+Name+'.SAV'
;-----------------------------------------------------
end
;----------------------------------------------------
;生成墙面
def make_walls ; create walls: a cylinder and two plates
extend= 0.5 ;0.2
rad_cy=0.5* width_cy
;圆柱面
_z0=-extend*height
_z1=height_cy*(1.0+extend)
command
wall type cylinder id=1 kn=w_stiff end1 0.0 0.0 _z0 end2 0.0 0.0 _z1 &
rad rad_cy rad_cy
end_command
;底面
_x0=-rad_cy*(1.0+extend)
_y0=-rad_cy*(1.0+extend)
_z0= 0.0
_x1= rad_cy*(1.0+extend)
_y1=-rad_cy*(1.0+extend)
_z1= 0.0
_x2= rad_cy*(1.0+extend)
_y2= rad_cy*(1.0+extend)
_z2= 0.0
_x3=-rad_cy*(1.0+extend)
_y3= rad_cy*(1.0+extend)
_z3= 0.0
command
wall id=5 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
;顶面
_x0=-rad_cy*(1.0+extend)
_y0=-rad_cy*(1.0+extend)
_z0= height_cy
_x1=-rad_cy*(1.0+extend)
_y1= rad_cy*(1.0+extend)
_z1= height_cy
_x2= rad_cy*(1.0+extend)
_y2= rad_cy*(1.0+extend)
_z2= height_cy
_x3= rad_cy*(1.0+extend)
_y3=-rad_cy*(1.0+extend)
_z3= height_cy
command
wall id=6 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
end
;----------------------------------------------------
;生成颗粒
def create_ball
n1 = 0
n2 = 0
w1 = -rad_cy
w2 = rad_cy
h1 = 0.0
h2 = height_cy
loop k (2, gradeN)
v= 4.0/3.0*pi*(0.5*(R(k-1)+R(k)))^3
n1 = n2+1
n2 = n2+round(MaterialV*p(k)/psum/v)
if n2>n1 then
r1 = R(k)
r2 = R(k-1)
clr=k-2
command
gen id=n1,n2 rad=r1,r2 x=w1,w2 y=w1,w2 z=h1,h2 filter ff_cylinder
prop dens=2700 ks=s_stiff kn=n_stiff color=clr range id=n1,n2
end_command
end_if
end_loop
ballmax = n2
ii=out(string(ballmax)+' particles were created')
command
pause 3
end_command
end
;----------------------------------------------------
;修改墙体刚度 change lateral wall stiffnesses
def cws
command
wall type cylinder id 1 kn=w_stiff
end_command
end
;----------------------------------------------------
;过滤函数
def ff_cylinder
ff_cylinder=0
_brad= fc_arg(0)
_bx= fc_arg(1)
_by= fc_arg(2)
_bz= fc_arg(3)
_rad= sqrt(_bx^2+_by^2)
if _rad+_brad > rad_cy then
ff_cylinder= 1
end_if
end
;----------------------------------------------------
; 墙体索引
def wall_addr
wadd1=find_wall(1)
wadd5=find_wall(5)
wadd6=find_wall(6)
end
;----------------------------------------------------
; 应力数组清零
def aszero
idx = 0
avsz= 0
avsr= 0
loop i (1,nsz)
asz(i)= 0
asr(i)= 0
end_loop
end
;----------------------------------------------------
; 试件压缩
def compress
Wrvel= vel_rad_z*_Wzvel
Wzvel= _Wzvel
nlpk= 8
if cmpstage=1
nlpk= 4
end_if
loopk= 0
cmpbefore = 0 ; rad=1, z=-1, unknown=0
loop mlpk (1,nlpk)
;消除原有应力
;aszero
zeroall
command
wall id 6 zvel= 0
wall id 5 zvel= 0
cycle 5000 ; to release the old stress
end_command
;aszero
cycn = 0
;新压缩循环
section ;+++++++++
loop while 1 # 0
;调整速度
if new_rad > 0.5*width
Wzvel0 = 0
w_radvel(wadd1)= 0.5*Wrvel
bar_apply ;only rad-acc
else
w_radvel(wadd1)= 0
Wzvel0 = Wzvel
baz_apply ;only z-acc
end_if
;循环计算
Wzvel0_= Wzvel0/2.0
_Wzvel0= -Wzvel0/2.0
command
wall id 6 zvel= Wzvel0_
wall id 5 zvel= _Wzvel0
cycle 50
end_command
cycn = cycn+1
;显示结果
statereq
;退出
if abs(avsr) >= abs(0.9*Target_press)
if cycn <= 5
loopk= loopk+1
if loopk >= 2 ;change the velocity
if new_rad > 0.5*width
Wrvel = Wrvel*0.5
else
Wzvel = Wzvel*0.5
end_if
end_if
if loopk >= 4 ;change the velocity
w_radvel(wadd1)= 0
;ba_release
exit
end_if
else
loopk= 0
end_if
exit section
end_if
end_loop
end_section ;+++++++++
w_radvel(wadd1)= 0
ii=out('The '+string(mlpk)+'-th loop is over!')
statereq
command
pause 2
end_command
end_loop
;清除加速度
ba_release
end
;----------------------------------------------------
; 球体加载卸载
def bvr_apply
bp = ball_head
loop while bp # null
brad = sqrt( b_x(bp)^2 + b_y(bp)^2 )
if brad >= width/2
b_xvel(bp) = Wrvel * b_x(bp) / brad
b_yvel(bp) = Wrvel * b_y(bp) / brad
b_zvel(bp) = 0
else
b_xvel(bp) = 0
b_yvel(bp) = 0
b_zvel(bp) = 0
end_if
bp = b_next(bp)
end_loop
end
def bar_apply
bp = ball_head
loop while bp # null
bforce = b_acc * 11100.29404 * b_rad(bp)^3 ;(4.0/3.0) * pi * b_rad(bp)^3 * 2650
brad = sqrt( b_x(bp)^2 + b_y(bp)^2 )
if brad >= width/3 ;width/2 * (2/3)
b_xfap(bp) = -1.0 * bforce * b_x(bp) / brad
b_yfap(bp) = -1.0 * bforce * b_y(bp) / brad
b_zfap(bp) = 0
else
b_xfap(bp) = 0
b_yfap(bp) = 0
b_zfap(bp) = 0
end_if
bp = b_next(bp)
end_loop
end
def baz_apply
bp = ball_head
loop while bp # null
bforce = b_acc * 11100.29404 * b_rad(bp)^3 ;(4.0/3.0) * pi * b_rad(bp)^3 * 2650
b_xfap(bp) = 0
b_yfap(bp) = 0
if b_z(bp)> height_cy/2.0
b_zfap(bp) = -1.0 * bforce
else
b_zfap(bp) = 1.0 * bforce
end_if
bp = b_next(bp)
end_loop
end
def ba_release
bp = ball_head
loop while bp # null
b_xfap(bp) = 0
b_yfap(bp) = 0
b_zfap(bp) = 0
bp = b_next(bp)
end_loop
end
;----------------------------------------------------
; 求取平均应力和应变
def get_ss
get_stress
get_strain
end
def get_ass
;计算瞬间应力应变
get_stress
get_strain
;平均化处理
idx = idx+1
idx0= idx+1
if idx > nsz
idx = idx -nsz
end_if
if idx0 > nsz
idx0 = idx0 -nsz
end_if
asz(idx)= wszz
asr(idx)= wsrr
avsz = avsz + (asz(idx)-asz(idx0))/nsz
avsr = avsr + (asr(idx)-asr(idx0))/nsz
;ii=out(' avsz='+string(avsz)+' avsr='+string(avsr))
end
def get_stress
new_rad=w_radend1(wadd1)
zdif=w_z(wadd6)-w_z(wadd5) ;注意：w_z( )函数返回值不是wall的z坐标。
new_height= height_cy+zdif ; 为墙体的*累积位移*，包括以前所有cycle的位移量
wsrr=-w_radfob(wadd1)/(new_height*2.0*pi*new_rad) ;径向应力
wszz=0.5*(w_zfob(wadd5)-w_zfob(wadd6))/(pi*new_rad^2.0) ;轴向应力
;ii=out(' wradfob='+string(w_radfob(wadd1)))
end
def get_strain
rdif=new_rad-rad_real
werr=2.0*rdif/(rad_real+new_rad) ;径向应变
wezz=(new_height-height_real)/height_real ;2.0*zdif/(height+new_height) ;轴向应变
wevol=wezz+2.0*werr ;体积应变 = x应变+y应变+z应变
end
;----------------------------------------------------
; 将height_real重置为当前高度
def reset_real
height_real= abs(height_cy+w_z(wadd6)-w_z(wadd5))
rad_real = w_radend1(wadd1)
end
;----------------------------------------------------
;统计空隙率 (height_real may be modified)
def poroscal
reset_real ;reset the height and rad
tot_vol= height_real*pi*rad_real^2.0
sum=0.0 ;get actual porosity
bp=ball_head
loop while bp # null
sum= sum+4.0/3.0*pi*b_rad(bp)^3
bp = b_next(bp)
end_loop
pballs= (1.0-sum/tot_vol)*100
ii=out(' The last height is '+string(height_real)+',radius is '+string(rad_real))
ii=out(' porosity is '+string(pballs)+'%')
end
;----------------------------------------------------
; 查询当前高度/直径
def statereq
height_cur=abs(height_cy+w_z(wadd6)-w_z(wadd5))
diam_cur= w_radend1(wadd1)*2
ii=out(' Nows height is '+string(height_cur)+',Diameter is '+string(diam_cur))
ii=out(' avsz ='+string(avsz)+', avsr ='+string(avsr))
ii=out(' Wzvel='+string(Wzvel)+', Wrvel='+string(Wrvel))
end
;----------------------------------------------------
; 清零
def zeroall
;zero all balls
bp = ball_head
loop while bp # null
loop k (1,3)
b_vvel(bp,k) =0
b_vrvel(bp,k)=0
end_loop
bp = b_next(bp)
end_loop
ii= out(' ok to zero all balls!')
;zero all clumps
clp = clump_head
loop while clp # null
loop k (1,3)
cl_vvel(clp,k) =0
cl_vrvel(clp,k)=0
end_loop
clp = cl_next(clp)
end_loop
ii= out(' ok to zero all clumps!')
ba_release
end
;----------------------------------------------------
;***************** 运 行 *****************
;*****************************************************
plot create model ; view 1
plot set cap size 25
plot set mag 1.5
plot set rot 0 0 0 ;30 0 40
plot add ball blue green cyan orange red lblue lgreen lred
;plot add wall lgray
;plot set plane normal 1 0 0 origin 0 0 0
;plot close
setup
;生成墙体
make_walls
wall_addr
;生成颗粒
create_ball
;设置height/rad_real为初始值
reset_real ;为方便get_ss中wezz变量计算
history id=1 avsz ;wszz
history id=2 avsr ;wsrr
plot create StressStation ; view 2
plot hist 1 2
;history id=3 Wzvel
;history id=4 Wrvel
;plot create VelocityStation ; view 3
;plot hist -3 -4
;开启应力跟踪
set fishcall 0 get_ass
;压缩球体颗粒
compress
wall id 6 zvel=0
wall id 5 zvel=0
cyc 20000
statereq
zeroall
;关闭应力跟踪
set fishcall 0 remove get_ass
SET w_stiff=1e7 ;make lateral wall stiffness=1/10 of ball stiffness
cws
poroscal
save saveFile
return
;*****************************************************
附件三：PFC三轴实验程序
new
SET random ; 随机生成
plo wall
plo add ball
caLL 消除浮点.txt ;
def make_walls ;生成墙
wextend =0.1
hextend =0.1
w_stiff= 1.6e10 ;墙的刚度，取球的1.1倍
_width=-width
_x0 = _width*(1.0 + wextend) ;width是长方体的宽度的一半
_y0 = _width*(1.0 + wextend)
_z0 = 0
_x1 = width*(1.0 + wextend)
_y1 = _width*(1.0 + wextend)
_z1 = 0
_x2 = width*(1.0 + wextend)
_y2 = width*(1.0 +wextend)
_z2 = 0
_x3 = _width*(1.0 + wextend)
_y3 = width*(1.0 + wextend)
_z3 = 0
command
wall id=1 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
_x0 = _width*(1.0 + wextend)
_y0 = _width*(1.0 + wextend)
_z0 = height ;height是长方体的高度
_x1 = _width*(1.0 + wextend)
_y1 = width*(1.0 + wextend)
_z1 = height
_x2 = width*(1.0 + wextend)
_y2 = width*(1.0 + wextend)
_z2 = height
_x3 = width*(1.0 + wextend)
_y3 = _width*(1.0 + wextend)
_z3 = height
command
wall id=2 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
_x0 = width
_y0 = _width*(1.0 + wextend)
_z0 = height*(1.0 + hextend)
_x1 = width
_y1 = width*(1.0 + wextend)
_z1 = height*(1.0 + hextend)
_x2 = width
_y2 = width*(1.0 + wextend)
_z2 = -hextend*height
_x3 = width
_y3 = _width*(1.0 + wextend)
_z3 = -hextend*height
command
wall id=3 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
_x0 = _width
_y0 = _width*(1.0 + wextend)
_z0 = height*(1.0 + hextend)
_x1 = _width
_y1 = _width*(1.0 + wextend)
_z1 = -hextend*height
_x2 = _width
_y2 = width*(1.0 + wextend)
_z2 = -hextend*height
_x3 = _width
_y3 = width*(1.0 + wextend)
_z3 = height*(1.0 + hextend)
command
wall id=4 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
_x0 = _width*(1.0 + wextend)
_y0 = width
_z0 = -hextend*height
_x1 = width*(1.0 + wextend)
_y1 = width
_z1 = -hextend*height
_x2 = width*(1.0 + wextend)
_y2 = width
_z2 = height*(1.0 + hextend)
_x3 = _width*(1.0 + wextend)
_y3 = width
_z3 = height*(1.0 + hextend)
command
wall id=5 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
_x0 = _width*(1.0 + wextend)
_y0 = _width
_z0 = height*(1.0 + hextend)
_x1 = width*(1.0 + wextend)
_y1 = _width
_z1 = height*(1.0 + hextend)
_x2 = width*(1.0 + wextend)
_y2 = _width
_z2 = -hextend*height
_x3 = _width*(1.0 + wextend)
_y3 = _width
_z3 = -hextend*height
command
wall id=6 kn=w_stiff face (_x0,_y0,_z0) (_x1,_y1,_z1) (_x2,_y2,_z2) &
(_x3,_y3,_z3)
end_command
end
; ----------------------------------------------------
def assemble ; 颗粒生成
s_stiff=0.0 ; 初始切向刚度
n_stiff=1.4e10 ;初始法向刚度
tot_vol = 4 * height * width^2.0 ;长方体的体积
rbar = 0.5 * (rlo + rhi) ;rlo最小球半径、rhi最大球半径
num = int((1.0 - poros) * tot_vol / (4.0 / 3.0 * pi * rbar^3))
mult = 1.6
rlo_0 = rlo / mult
rhi_0 = rhi / mult
command
gen id=1,num rad=rlo_0,rhi_0 x=_width,width y=_width,width z=0,height tries 10000000;方形的话就不用过滤器了
;上面球的半径是缩小了mult倍的
prop dens=2700 ks=s_stiff kn=n_stiff ;dens即是岩石颗粒的密度
end_command
ii = out(string(num)+' particles were created')
sum = 0.0 ;获得有效孔隙度
bp = ball_head
loop while bp # null
sum = sum + 4.0 / 3.0 * pi * b_rad(bp)^3
bp = b_next(bp)
end_loop
pmeas = 1.0 - sum / tot_vol
mult = ((1.0 - poros) / (1.0 - pmeas))^(1.0/3.0)
command
ini rad mul mult ;将球的半径放大mult倍
cycle 30000
prop ks=1.4e10 fric 0.5
cycle 3100
end_commandend
; ----------------------------------------------------
; ----------------------------------------------------
macro zero 'ini xvel 0 yvel 0 zvel 0 xspin 0 yspin 0 zspin 0' ;宏，现在定义了后面就可以直接引用了xspin可能是x方向的角速度
SET height=0.1 width=0.025 rlo=0.001 rhi=0.002 poros=0.40
make_walls
assemble
SET flt_def=3 flt_remain=0.0 ;上面把球生成完了以后就消除浮点，其中 flt_def=3表示球的连接小于3的就是浮点，需要删除的
flt_eliminate ;浮点的命令之前有的
zero
save tt_ass.SAV
;fname: triax_2.DAT Servo-control and initial stress state - triax sample
res tt_ass.SAV ; restore compacted assembly
; ----------------------------------------------------
def get_ss
zdif = w_z(wadd2) - w_z(wadd1)
ydif = w_y(wadd5) - w_y(wadd6)
xdif = w_x(wadd3) - w_x(wadd4)
new_width1 = 2*width+xdif
new_width2 = 2*width+ydif
new_height = height + zdif
wsxx = 0.5*(w_xfob(wadd4) - w_xfob(wadd3))/(new_width2*new_height) ;
wsyy = 0.5*(w_yfob(wadd6) - w_yfob(wadd5))/(new_width1*new_height) ;y
wszz = 0.5*(w_zfob(wadd1) - w_zfob(wadd2))/(new_width1*new_width2) ;z
wexx = 2.0 * xdif / (2*width + new_width1) ;x向应变参看手册186页3.16
weyy = 2.0 * ydif / (2*width + new_width2) ;y向应变
wezz = 2.0 * zdif / (height + new_height) ;z向应变
wevol = wexx + weyy + wezz ;体积应变
end
; ----------------------------------------------------
def get_gain ; 选取轴向和径向获得的伺服参数
alpha = 0.5 ; 松弛参数
count = 0
avg_stiff = 0
cp = contact_head ; 找到上、下墙上接触的平均数量
loop while cp # null
if c_gobj2(cp) = wadd1
count = count + 1
avg_stiff = avg_stiff + c_kn(cp)
end_if
if c_gobj2(cp) = wadd2
count = count + 1
avg_stiff = avg_stiff + c_kn(cp)
end_if
cp = c_next(cp)
end_loop
ncount = count / 2.0
avg_stiff = avg_stiff / count
gz = 4 * alpha * width^2.0/ (avg_stiff * ncount * tdel) ;可以参看187页公式3.22，是增益参数
count = 0
avg_stiff = 0
cp = contact_head ; 找到左、右墙上接触的平均数量
loop while cp # null
if c_gobj2(cp) = wadd3
count = count + 1
avg_stiff = avg_stiff + c_kn(cp)
end_if
if c_gobj2(cp) = wadd4
count = count + 1
avg_stiff = avg_stiff + c_kn(cp)
end_if
cp = c_next(cp)
end_loop
ncount = count / 2.0
avg_stiff = avg_stiff / count
gx = 2 * alpha * width * height/ (avg_stiff * ncount * tdel) ;可以参看187页公式3.22，是增益参数
count = 0
avg_stiff = 0
cp = contact_head ; 找到前、后墙上接触的平均数量
loop while cp # null
if c_gobj2(cp) = wadd5
count = count + 1
avg_stiff = avg_stiff + c_kn(cp)
end_if
if c_gobj2(cp) = wadd6
count = count + 1
avg_stiff = avg_stiff + c_kn(cp)
end_if
cp = c_next(cp)
end_loop
ncount = count / 2.0
avg_stiff = avg_stiff / count
gy = 2 * alpha * width * height/ (avg_stiff * ncount * tdel) ;可以参看187页公式3.22，是增益参数
end
; ----------------------------------------------------
def servo ;伺服控制
while_stepping
get_ss ;计算应力和应变
udx = gx * (wsxx - sxxreq) ;参看187页3.17udx是墙的速度wsxx是x向实际压力，sxxreq是我们需要的x向压力
udy = gy * (wsyy - syyreq)
w_xvel(wadd4) = udx
w_xvel(wadd3) = -udx
w_yvel(wadd6) = udy
w_yvel(wadd5) = -udy
if z_servo = 1 ;选择轴压伺服的开关，是1的时候就开，否则就关
udz = gz * (wszz - szzreq)
w_zvel(wadd1) = udz
w_zvel(wadd2) = -udz
end_if
end
; ----------------------------------------------------
def iterate
loop while 1 # 0
get_gain
if abs((wsxx - sxxreq)/sxxreq) < sig_tol then ;实际的x围压和需要的围压达到一致就不加载了
if abs((wsyy - syyreq)/syyreq) < sig_tol then ;实际的y围压
if abs((wszz - szzreq)/szzreq) < sig_tol then ;实际的轴压和需要的轴压达到一致就不加载了
exit
end_if
end_if
end_if
command
cycle 100
end_command
end_loop
end
; ----------------------------------------------------
def wall_addr
wadd1 = find_wall(1)
wadd2 = find_wall(2)
wadd3 = find_wall(3)
wadd4 = find_wall(4)
wadd5 = find_wall(5)
wadd6 = find_wall(6)
end
wall_addr
zero
prop fric=1.3 ;球表面的摩擦系数，不是摩擦角
prop pb_kn=1.4e10 pb_ks=1.4e10 pb_rad=1 ;法向刚度系数、切向刚度系数、半径乘数
prop pb_nstren=7e10 pb_sstren=7e10 ;法向平行粘结力、切向平行粘结力
SET sxxreq=-1e3 syyreq=-1e3 szzreq=-1e3 sig_tol=1 z_servo=1 ;单轴压缩，围压数量级很低
iterate ; 使所有方向的压力都达到静水压力状态，上面达到的是所有的方向都是达到1MPa
sav 第一静水压力.sav
res 第一静水压力.sav ; restore initial stressed assembly
; ----------------------------------------------------
def set_ini ; 设置初始应变
wezz_0 = wezz ;轴向应变，前面计算了的
wevol_0 = wevol ;轴向速度
end
; ----------------------------------------------------
def conf ; 记录的变量
devi = wszz - wsxx ; 偏应力，最大主应力减最小主应力，也就是说这里x方向取的是最小值
deax = wezz - wezz_0 ; 轴向应变
devol = wevol - wevol_0 ; 体积应变
confx = wsxx ; x压应力
confy = wsyy ; y压应力
confz = wszz ; z压应力
end
; ----------------------------------------------------
def accel_platens
; -----在一定的时步内室加载板达到相应的速度 Accelerates the platens to achieve vel of _vfinal in _nsteps,
; using _nchunks
_niter = _nsteps / _nchunks
loop _chnk (1,_nchunks)
if _close = 1 then
_vel = _chnk*(_vfinal/_nchunks)
else
_vel = -_chnk*(_vfinal/_nchunks)
end_if
_mvel = -_vel
command
wall id 1 zvel= _vel
wall id 2 zvel= _mvel
cycle _niter
end_command
end_loop
end
; ----------------------------------------------------
set_ini
history id=5 conf ;监测
history id=6 devi ;监测偏应力
history id=7 deax ;监测轴向应变
history id=8 devol ;监测体积应变
history id=9 wexx ;监测x方向应变
history id=10 deax ;监测z方向应变
history id=11 confz ;z压应力
history id=12 confx ;z压应力
history id=13 confy ;z压应力
SET hist_rep=50 ;每50步监测一次
SET z_servo=0 ;围压都还保持伺服状态，取消轴向加压的伺服控制
zero ;设定x、y、z方向所有的速度都等于0
sav tt_init.SAV ; ready for modulus and failure tests
;fname: triax_9.DAT (friction and parallel bonds)
res tt_init.sav
set _vfinal= 0.1 _nsteps= 4000 _nchunks= 80 ;_vfinal最终的加载速率，
set _close= 1 ; 加载 _close=1时候是加载
accel_platens
plo his -11 vs -10
cyc 25000000
