张晓平 ，王思敬，李 黎，王彦兵

（1. 敦煌研究院 古代壁画保护国家文物局重点科研基地，甘肃 敦煌 736200；2. 中国科学院地质与地球物理研究所 中国科学院工程地质力学重点实验室，北京 100029；3. 中国文化遗产研究院，北京 100871）

1 引 言

我国是世界上文明起源最早的国家之一，保存在西北古丝绸之路沿线数千公里的战国、汉代、明代长城及新疆交河故城遗址等土遗址文物都是我国珍贵文化遗产的重要组成部分，具有非常重要的历史、科学和艺术价值，其中交河故城则是西北干旱地区土遗址的代表。现存土建筑20余万平方米，是研究中西文化交流，中国西部和中亚地区历史以及城市建筑史、宗教史、美术史、民族史等多种学科不可多得的珍贵实物资料，具有极大的保存价值。但由于遗址区常年大风、沙暴，集中式强降雨等，造成遗址严重风化风蚀、雨水冲刷侵蚀、墙体开裂、坍塌，急需进行抢修，目前土遗址的保护加固工作仅步入初始阶段[1]。其中加固所需的 PS(potassium silicate)材料[1]及其加固机制研究是最基础的研究内容。

土遗址 PS加固机制研究是土遗址加固的理论部分，在进行大量土遗址抢修加固应用[2－6]的同时，也已经开展了相关的研究探讨[7－8]。通过PS加固遗址土过程的热分析、加固前后物相分析、可溶盐反应、X射线光电子能谱分析、热重分析等得出，交河故城遗址土加固前后物相没有发生变化，PS材料固化的主要物理作用包括：一方面填塞遗址土颗粒间的空隙；另一方面增加土颗粒间的胶结强度,形成较致密的网状结构。此外，PS在遗址土颗粒之间的化学键合作用、PS与遗址土矿物的羟基发生缩合反应，把遗址土进一步结合起来，形成了整体的联结[8]。遗址土风蚀试验结果表明：风蚀主要发生在挟沙风条件下，未加固样的风蚀态势表现为随挟沙风风速增大，风蚀强度增加；加固样品的抗风蚀强度（抗风蚀强度指试样抵抗风蚀的能力。在相同风蚀条件下，抗风蚀强度与风蚀量成反比例关系，即试样风蚀量的减小，表征抗风蚀强度的提高）明显高于未加固样品，其抗风蚀强度是原样的 2.9～38倍[9]。

以上研究结果初步揭示了 PS加固遗址土的加固机制和加固前后的风蚀特性，给颗粒元模拟提出了建模概念。而接下来的研究方向应该是针对大量不同赋存条件、不同特性遗址土分别展开专项研究，验证并深入研究 PS加固机制及加固后各耐环境因素。将土力学研究领域中的颗粒元模拟应用到土遗址保护加固中，为 PS加固机制分析乃至现场加固工程提供理论依据。颗粒元的模拟应用将增加 PS加固机制研究的途径，增加试验的针对性和施工的预见性，减小试验取样和施工反复对遗址的不利影响。但颗粒元在遗址土中要得到有效运用，首要工作是进行颗粒元与遗址土室内试验结果的拟合与分析，这是本文拟进行的首要研究内容，本文选取具有代表性的交河故城遗址土作为研究对象，进行颗粒元模拟研究工作，最终对加固前后遗址土进行风蚀模拟，研究加固前后的风蚀特征。

颗粒元方法最早由Cundall提出[10]，采用颗粒集合模拟岩土体在近 10 年受到广泛关注。该模型不需要输入岩土体的本构模型，只需要指定颗粒间的接触刚度和接触强度参数，颗粒间的计算只是采用简单的牛顿运动定律，便能模拟出岩土体的非线性力学特征和微观破坏机制，与连续网格的有限元计算相比具有独特的优势。通过在颗粒表面增加压力，能模拟风压、液压的作用。张晓平等[11－12]采用颗粒元模拟含了软弱夹层土样的变形破坏机制，并尝试模拟滑坡过程[13－15]及抗滑桩加固[16]；吴顺川[17]利用颗粒元模拟了压力注浆；刘洋[18]利用颗粒元流-固耦合分析模型模拟了饱和砂土的液化现象；张刚[19]利用颗粒元模拟了管涌现象等。遗址土多为人工夯土、版筑土等，不同于自然环境的土体，颗粒元在土遗址 PS加固模拟中的运用，必将拓宽其应用范围，也为土遗址加固机制研究及各耐风蚀环境影响分析提供一种新的思路和方法。

2 遗址土的物理力学性质

2.1 遗址土的物理特性

交河故城本体土的类型有4种：生土、夯土、垛泥和土坯。主要为粉土和粉质黏土。颗粒分析结果表明，粉粒含量普遍在80%左右，垛泥的粒度成分变化极大，级配不佳，而生土和夯土的级配良好。李最雄等[1]曾统计了故城遗址建筑材料的物理性质指标；2008年和法国等[20]又进行了大量的室内常规试验。比较两作者所给出的物性指标数据基本吻合，现仅给出文献[20]的最新试验资料，如表 1所示。

2.2 PS加固前后遗址土的力学性质

将原状土用5 %浓度PS溶液喷洒3次，每次喷洒时间间隔为24 h，全部喷洒完毕后在室内养护至自然风干，测试其单轴抗压和抗拉强度，选取典型的土样 JH-029、JH-030、JH-011和 JH–012，其测试结果如表2。可以看出，经过PS材料加固后，生土和垛泥的强度得到了明显提高，提高幅度120%～180%，生土的强度要高于垛泥[20]。

3 PS加固前后遗址土的力学性质颗粒元(PFC)模拟

本研究首先采用颗粒元对生土加固前后力学性质进行拟合，然后模拟加固前后生土的风蚀过程。

3.1 颗粒元模型简介

遗址土模拟的基础是颗粒元中的连接颗粒模型(bonded particle model，简称BPM)，因此，首先对BPM模型及其参数进行介绍。

本次模拟采用 Itasca公司开发的颗粒元程序PFC2D，该程序自1995年首次发布以来，在许多研究项目中得到了应用。PFC2D通过有限厚度的圆盘集合或单层（垂直于平面方向）球来描述土体或岩石，颗粒间通过黏结和摩擦连接，连接为一定刚度和颗粒间的有限区域。当指定的连接强度超过施加在连接上的应力（如拉力、剪力或由于颗粒旋转而形成的弯矩），连接发生破裂形成破裂面，如裂纹。在这些过程中，PFC不需要任何本构关系来描述屈服后的响应和断裂韧度来控制破裂行为，而只需要颗粒运动定律、简单的颗粒变形公式和颗粒连接断裂的准则，有关PFC的详细情况参见文献[21－23]。平行连接模型是PFC中连接模型(BPM)的一种，近似为两个颗粒在接触点处胶结作用，能承受拉力、剪力和弯矩，是模拟岩土材料的理想模型。确定平行连接模型的主要参数如下：

描述矿物颗粒参数有 { Rmin, Rmax/Rmin, ρ,Ec，(kn/ks), μ}；描述胶结作用的平行连接参数有

其中：Rmin和 Rmax/ Rmin描述颗粒半径，在模型中满足以Rmax和Rmin为上下限的均匀分布；ρ为颗粒的密度；Ec和分别为颗粒和平行连接的弹性模量；和分别为颗粒和平行连接的法向刚度与切向刚度的比率；为设置平行连接半径R的半径系数，为颗粒半径；μ为颗粒间的摩擦系数；和分别为平行连接的抗拉强度和抗剪强度。在给定弹性模量和法向刚度与切向刚度的比率后，通过式（1）、（2）来计算颗粒或平行连接的法向和切向刚度。

对法向刚度来讲，可以想象成如图1所示的弹性梁。

图1 颗粒间接触和平行连接的等效连续Fig.1 Equivalent continuum material of parallel bond particle system

平行连接颗粒间相对运动会形成轴向力 T、切向力V和弯矩M。连接材料能承担的最大法向应力和切向应力由式（3）、（4）表达。

A和I是平行连接截面（如图1(b)）的面积和转动惯量，T为正值表示抗拉。当或超过相应的连接强度连接就会破坏，在两个颗粒间形成拉裂纹和剪裂纹。

PFC的颗粒参数和连接参数可以通过单轴压缩试验和巴西劈裂试验来拟合。通过单轴压缩试验可得到径向应变、轴向应变和峰值强度。抗拉强度通过巴西劈裂计算得到。PFC试样的生成和加载方法在文献[21－23]中有详细论述。由于PFC中引入了随机模型，因而颗粒半径和强度参数也都服从统计的随机模型，虽然同样的命令程序，但没有两次PFC的计算结果会完全相同。因此，在参数研究中，每种参数计算10次，求平均值。这样就会增加很大的计算工作量，但如同实际的岩石力学试验一样，可避免由于个别样品的差异而得到不相干的结论。

3.2 颗粒元模型

颗粒元模拟土体，颗粒粒径不能与土体粒径一一对应，但可以按照级配曲线进行近似粒径比例模拟。实际上交河故城粉粒（粒径0.005～0.075 mm）含量普遍在80%左右，因此，出于减小计算量，使得计算具有实用性，故不考虑颗粒级配，而拟采用粒径0.005～0.075 mm范围内均匀分布颗粒模型来模拟遗址土。实际颗粒元计算中发现，粒径0.005～0.075 mm范围内均匀分布颗粒，计算50 mm ×50 mm的原状土压缩试样或 70.7 mm×70.7 mm的重塑土样，仍然存在颗粒数巨大计算机无法执行的情况。故考虑将土体粒径扩大 20倍进行模拟，即在 PFC模型中，颗粒直径在 0.1～1.5 mm 范围内均匀分布。如图2所示，模型尺寸50 mm×50 mm，颗粒直径0.1～1.5 mm范围均匀分布，模型包含4 177个颗粒。

PFC颗粒模型中，指定颗粒密度是指颗粒本身的密度。由于颗粒之间有孔隙，计算模型的密度比颗粒密度要小。因此，在拟合实际土样密度中，先将实际土样密度乘以计算模型体积（二维模型中为单位厚度），计算出总质量；然后用总质量除以模型颗粒的实际体积，能得到应赋颗粒密度值。这一段通过独立编写扩展Fish语言来实现，即只需输入实际土样的密度，便能计算出PFC颗粒密度值，并自动赋值。

3.3 生土PS加固前后抗压、抗拉强度模拟

生土压缩试样计算模型尺寸为50 mm×50 mm，巴西劈裂试样为φ50 mm的圆盘。生土取天然密度平均值，参照表1生土密度范围为1.51～1.77 g/cm3，故近似取生土PFC模型密度为1.64 g/cm3。采用Fish计算得到对应颗粒密度为1.833 g/cm3。

通过反复试算逐步逼近生土的抗压和抗拉强度参数，最后取表3所示参数作为生土加固前、后对应的颗粒元模型参数。对比表中的颗粒参数可以看出，PS加固后生土的强度参数比加固前要大。也就是说，PS加固土体的作用是通过增加颗粒间平行连接强度来实现的。如3.1节所述，出于统计考虑，每种试样抗压和抗拉试验分别计算10次（见表4），取10次计算结果的平均值列于表 5。模拟得到加固前后抗压和抗拉强度参数与试验结果非常接近。

4 PS加固前后遗址土风蚀颗粒元(PFC)模拟

4.1 风蚀模型

遗址土风蚀试验结果表明，风蚀主要发生在挟沙风条件下。目前对风沙运动的力学机制研究，在沙漠化过程研究中较多，且集中在风沙跃移运动发展过程。大量的研究者从离散的角度对其进行了数值模拟[24－28]，并对风沙的启动机制和颗粒间的相互作用进行了相关研究，这些研究吻合于自然风沙跃移运动过程的基本特征。本文的研究内容主要集中在分析加固前、后遗址土的抗风蚀能力变化，因此，将风沙的运动简化为一部分以一定初速度撞向土体表面的颗粒。撞击的结果使得土体表面部分颗粒间连接断开，从土体中剥落，形成风蚀（见图3）。图中左侧为土体，右侧竖线处为随机产生不同粒径的挟沙风颗粒，以一定的水平初速度（箭头所示）撞向土体。

表3 模拟加固前生土颗粒参数Table 3 Particle parameters of raw soil before reinforcement

表4 生土加固前、后抗压强度和抗拉强度模拟结果统计表Table 4 Numerical simulation of compressive strength and tensile strength before and after reinforcement of raw soil

表5 生土试验结果和模拟结果对比表Table 5 Comparison between experimental and numerical outputs

图3 挟沙风模拟示意图Fig.3 Illustration of sand-driving wind flow simulation

4.2 风蚀模拟计算及分析

文献[9]中采用颗粒速度为 6、10、15、20 m/s的挟沙风进行了风洞吹蚀试验。发现20 m/s挟沙风条件下，加固前后土体的抗风蚀能力差别显著。为了便于比较，本颗粒元模拟中，也选取20 m/s的颗粒速度模拟挟沙风。挟沙风颗粒从图3所示右侧竖线位置随机产生，每个颗粒的初始位置都在该竖线上，但高低位置随机分布。挟沙风颗粒直径服从0.1～1.5 mm之间的均匀分布，与土体颗粒分布保持一致。挟沙风颗粒密度及弹性参数与生土相同：ρ=1.833 g/cm3，Ec=0.08 GPa，kn/ks=1，μ=0.5。计算过程设置每循环5 000个时步，产生一个新的挟沙风颗粒。因此，挟沙风颗粒数与循环步数成正比例的关系，可以用挟沙风颗粒数来代表吹蚀时间的长短。

图4为生土加固前、后不同吹蚀时间（用不同挟沙风颗粒数表示）条件下土样颗粒元模型的吹蚀图片。加固前，在挟沙风颗粒为300和600个时，试样右侧面已经有较明显的风蚀；挟沙风颗粒为900个时，风蚀变得非常显著；挟沙风颗粒继续增加到1 200个和1 500个时，风蚀面继续加大、风蚀掏蚀深度加深。

加固后，在挟沙风颗粒数为300、600、900个时，试样右侧表面风蚀还非常轻微；到挟沙风颗粒达到1 200个以后，表面才出现较显著的风蚀现象。

图4 生土加固前、后不同吹蚀时间（即不同挟沙风颗粒数）条件下土样颗粒元模型的吹蚀Fig.4 Snapshots of erosion under different times (i. e. different numbers of sand-driving wind particles)before and after reinforcement of raw soil

为了定量比较加固前后的风蚀量，即风蚀掉的颗粒面积，图5列出生土加固前、后风蚀掉的颗粒面积与挟沙风颗粒数的对比曲线。可以看出，随着挟沙风颗粒数的增加，即风蚀时间的增长，加固前土样的风蚀量一直呈较快的增长趋势，比加固后的显著偏大。加固后的土样随挟沙风颗粒数的增加，增幅缓慢而且趋于稳定。可见加固后生土土样抗风蚀能力得到了显著的增强。与风洞试验结果统计的规律一致[9]。即加固样品的抗风蚀能力明显高于未加固样品，其抗风蚀强度是原样的2.9～38倍[9]。

图5 生土加固前、后风蚀量与挟沙风颗粒数对比曲线Fig.5 Erosion mass of soil vs. erosion time before and after reinforcement of raw soil

5 结 论

（1）颗粒元模拟通过增加颗粒间平行连接强度，实现对 PS加固的模拟。模拟得到土样加固前后抗压和抗拉强度参数与试验结果接近。

（2）随机生成挟沙风颗粒，以一定的速度撞向土体，模拟挟沙风的吹蚀作用。挟沙风颗粒数与循环步数成正比例，可以用挟沙风颗粒数来代表吹蚀时间的长短。而挟沙风颗粒的速度则代表挟沙风风速。

（3）风蚀模拟结果表明，在20 m/s的挟沙风吹蚀作用下，风蚀程度随吹蚀时间的增加而增大，未加固土样的风蚀程度增幅度远大于加固土样；同样吹蚀时间条件下，加固土样的抗风蚀强度明显高于未加固土样。PS加固后生土的抗风蚀能力得到显著增强，与风洞试验结果的统计规律一致。

（4）本研究拟合的颗粒元模型，可以应用于遗址土加固的进一步研究，如加固机制研究及耐风蚀、雨蚀、冻融等诸环境影响分析研究。

[1]李最雄. 丝绸之路古遗址保护[M]. 北京: 科学出版社,2003.

[2]李最雄, 王旭东, 张志军, 等. 秦佣坑土遗址的加固试验[J]. 敦煌研究, 1998, (4): 151－158.LI Zui-xiong, WANG Xu-dong, ZHANG Zhi-jun, et al.Consolidation tests of Qinyong site[J]. Dunhuang Research, 1998, (4): 151－158.

[3]李最雄, 王旭东, 郝利民. 室内土建筑遗址的加固试验—半坡土建筑遗址的加固试验[J]. 敦煌研究, 1998, (4):144－150.LI Zui-xiong, WANG Xu-dong, HAO Li-min. Research on the conservation of indoor earth structure sites—Experimentation of chemical consolidation on ancient earth—Structure Site of Ban Po and Qin Yong Keng[J]. Dunhuang Research, 1998, (4): 144－150.

[4]孙满丽. 吐鲁番交河故城保护加固研究[D]. 兰州: 兰州大学, 2006.

[5]王旭东, 张鲁, 李最雄, 等. 银川西夏 3号陵的现状及保护加固研究[J]. 敦煌研究, 2002, (4): 64－72.WANG Xu-dong, ZHANG Lu, LI Zui-xiong, et al. The study of existing condition and consolidation project of NO.3 tomb of the Western Xia Mausoleums[J].Dunhuang Research, 2002, (4): 64－72.

[6]苏伯民, 李最雄, 胡之德. PS与土遗址作用机理的初步探讨[J]. 敦煌研究, 2000, 63(1): 3－35.SU Bo-min, LI Zui-xiong, HU Zhi-de. Study on consolidating mechanism between PS and earth relics[J].Dunhuang Research, 2000, 63(1): 3－31.

[7]王旭东等. 国家十一五科技支撑计划课题: 大遗址保护关键技术研究与开发—土遗址保护关键技术研究报告[R]. 敦煌: 敦煌研究院, 2009.

[8]LI Li, WANG Si-jing, SHAO Ming-shen, et al. Impact of environment factors on consolidating earthen architecture sites with PS[C]//Proceedings of International Symposium Conservation of Ancient Sites &ISRM-Sponsored Regional Symposium. Beijing: Science Press, 2010.

[9]赵海英, 李最雄, 汪稔, 等. PS 材料加固土遗址风蚀试验研究[J]. 岩土力学, 2008, 29(2): 392－396.ZHAO Hai-ying, LI Zui-xiong, WANG Ren, et al. Wind erosion experiment of ancient earthen site consolidated by PS material[J]. Rock and Soil Mechanics, 2008, 29(2):392－396.

[10]CUNDALL P A. A computer model for simulating progressive large scale movements in blocky rock system[C]//Proceedings of the Symposium of the International Society of Rock Mechanics. Nancy, France:[s. n.], 1971.

[11]张晓平, 吴顺川, 张志增, 等. 含软弱夹层土样变形破坏过程细观数值模拟及分析[J]. 岩土力学, 2008, 29(5):1200－1204.ZHANG Xiao-ping, WU Shun-chuan, ZHANG Zhi-zeng,et al. Numerical simulation and analysis of failure process of soil with weak intercalated layer[J]. Rock and Soil Mechanics, 2008, 29(5): 1200－1204.

[12]张晓平, 吴顺川, 张兵, 等. 软弱夹层几何参数对试样力学行为影响颗粒元模拟研究[J]. 工程地质学报, 2008,16(4): 539－545.ZHANG Xiao-ping, WU Shun-chuan, ZHANG Bing,et al. Particle flow code simulation of specimen mechanical behavior with different geometric parameters of weak seams[J]. Journal of Engineering Geology,2008, 16(4): 539－545.

[13]张晓平, 吴顺川, 王思敬. 类土质路堑边坡动态监测及数值模拟分析[J]. 岩石力学与工程学报, 2008, 27(增刊2): 3431－3440.ZHANG Xiao-ping, WU Shun-chuan, WANG Si-jing.Dynamic monitoring and numerical analysis of soil-like cut slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(Supp.2): 3431－3440.

[14]张晓平. 类土质边坡破坏机理分析及工程应用研究[D].北京: 北京科技大学, 2007.

[15]吴顺川, 张晓平, 刘洋. 基于颗粒元模拟的含软弱夹层类土质边坡变形破坏过程分析[J]. 岩土力学, 2008,29(11): 16－21.WU Shun-chuan, ZHANG Xiao-ping, LIU Yang.Analysis of failure process of similar soil slope with weak intercalated layer based on particle flow simulation[J].Rock and Soil Mechanics, 2008, 29(11): 16－21.

[16]张晓平, 王思敬, 王幼明, 等. 二维离散元模拟抗滑桩的折算方法研究[J]. 岩土工程学报, 2010, 32(2): 271－278.ZHANG Xiao-ping, WANG Si-jing, WANG You-ming,et al. Conversion of anti-sliding piles into 2-dimensional discrete element simulation[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(2): 271－278.

[17]吴顺川. 压力注浆复合锚固桩地基处治理论研究及工程应用[D]. 北京: 北京科技大学, 2004.

[18]刘洋. 砂土液化破坏的细观力学机制与数值模拟[D].上海: 同济大学, 2006.

[19]张刚. 管涌现象细观机理的模型试验与颗粒流数值模拟研究[D]. 上海: 同济大学, 2007.

[20]和法国, 谌文武, 张景科, 等. PS 材料加固交河故城土体试验研究[J]. 敦煌研究, 2007, (5): 38－41.HE Fa-guo, CHEN Wen-wu, ZHANG Jing-ke, et al.Experimental research on PS reinforcing Jiaohe cultural heritage soil[J]. Dunhuang Research, 2007, (5): 38－41.

[21]POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329－1364.

[22]CHO N, MARTIN C D, SEGO D C, et al. A clumped particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(7): 997－1010.

[23]Itasca Consulting Group, Inc. PFC2D (Particle Flow Code in 2 Dimension)Version 3.1[M]. Minneapolis,Minnesota: Itasca Consulting Group, Inc., 2004.

[24]李万清. 风沙跃移运动发展过程的离散动力学模拟[J].中国沙漠, 2006, 26(1): 47－53.LI Wan-qing. Discrete dynamics simulation on developing process of aeolian sand saltation[J]. Journal of Desert Research, 2006, 26(1): 47－53.

[25]HAFF P K, ANDERSON R S. Grain scale simulations of loose sedimentary beds: The example of grain-bed impacts in aeolian saltation[J]. Sedimentology, 1993,40(2): 175－198.

[26]亢力强, 郭烈锦. 风沙运动的 DPM 数值模拟[J]. 工程热物理学报, 2006, 27(3): 441－444.KANG Li-qiang, GUO Lie-jin. Simulation of windblown sand movement by discrete particle model[J]. Journal of Engineering Thermophysics, 2006, 27(3): 441－444.

[27]亢力强, 郭烈锦. 风沙跃移中颗粒与多粒径床面碰撞的数值模拟[J]. 工程热物理学报, 2006, 27(1): 82－84.KANG Li-qiang, GUO Lie-jin. Numerical simulation of particles impacting multiple-grain size bed in saltation[J].Journal of Engineering Thermophysics, 2006, 27(1):82－84.

[28]亢力强, 郭烈锦. 风沙跃移中颗粒冲击起动的数值模拟[J]. 自然科学进展, 2005, 15(2): 252－256.